[Numpy-discussion] Scipy dot

[Frédéric Bastien]

On Fri, Nov 9, 2012 at 3:32 AM, Nathaniel Smith <njs at pobox.com> wrote:

> On Fri, Nov 9, 2012 at 6:18 AM, Nicolas SCHEFFER
> <scheffer.nicolas at gmail.com> wrote:
> > Fred,
> >
> > Thanks for the advice.
> > The code will only affect the part in _dotblas.c where gemm is called.
> > There's tons of check before that make sure both matrices are of ndim 2.
> > We should check though if we can do these tricks in other parts of the
> function.
> >
> > Otherwise:
> > - I've built against ATLAS 3.10
> > - I'm happy to add a couple more test for C and F-contiguous. I'm not
> > sure how to get the third type (strided), would you have an example?
>
> def with_memory_order(a, order):
>     assert order in ("C", "F", "discontig")
>     assert a.ndim == 2
>     if order in ("C", "F"):
>         return np.asarray(a, order=order)
>     else:
>         buf = np.empty((a.shape[0] * 2, a.shape[1] * 2), dtype=a.dtype)
>         buf[::2, ::2] = a
>         # This returns a view onto every other element of 'buf':
>         result = buf[::2, ::2]
>         assert not result.flags.c_contiguous and not
> result.flags.f_contiguous
>         return result
>
> > The following test for instance checks integrity against
> > multiarray.dot, which I believe is default when not compiled with
> > BLAS.
> > Dot is a hard function to test imho, so if anybody has ideas on what
> > kind of test they'd like to see, please let me know.
> >
> > If that's ok I might now be able to:
> > - Check for more bugs, I need to dig a bit more in the gemm call, make
> > sure everything is ok.
> > - Create an issue on github and link to this discussion
> > - Make a commit in a seperate branch
> > - Move forward like that.
> >
> > ==
> > import numpy as np
> > from time import time
> > from numpy.testing import assert_almost_equal
> >
> > def test_dot_regression():
> >     """ Test numpy dot by comparing with multiarray dot
> >     """
> >     np.random.seed(7)
> >     a = np.random.randn(3, 3)
> >     b = np.random.randn(3, 2)
> >     c = np.random.randn(2, 3)
> >
> >     _dot = np.core.multiarray.dot
> >
> >     assert_almost_equal(np.dot(a, a), _dot(a, a))
> >     assert_almost_equal(np.dot(b, c), _dot(b, c))
> >     assert_almost_equal(np.dot(b.T, c.T), _dot(b.T, c.T))
> >
> >     assert_almost_equal(np.dot(a.T, a), _dot(a.T, a))
> >     assert_almost_equal(np.dot(a, a.T), _dot(a, a.T))
> >     assert_almost_equal(np.dot(a.T, a.T), _dot(a.T, a.T))
>
> You should check that the result is C-contiguous in all cases too.
>
> for a_order in ("C", "F", "discontig"):
>     for b_order in ("C", "F", "discontig"):
>         this_a = with_memory_order(a, a_order)
>         this_b = with_memory_order(b, b_order)
>         result = np.dot(this_a, this_b)
>         assert_almost_equal(result, expected)
>         assert result.flags.c_contiguous
>
> You could also wrap the above in yet another loop to try a few
> different combinations of a and b matrices (perhaps after sticking the
> code into a utility function, like run_dot_tests(a, b, expected), so
> the indentation doesn't get out of hand ;-)). Then you can easily test
> some of the edge cases, like Nx1 matrices.
>

I agree that tests are needed the for Nx1 and variant cases. I saw blas
error being raised with some blas version.

You also need to test with the output provided, so there is 3 loops:

for a_order in ("C", "F", "discontig", "neg"):
    for b_order in ("C", "F", "discontig", "neg"):
        for c_order in ("C", "F", "discontig", "neg"):

I also added the stride type "neg", I'm not sure if it is needed, but that
is other corner cases. neg =>  result = buf[::-1, ::-1]

I just looked again at our code and there is another constrain: that the
strides are multiple of the elemsize. Theano do not support not aligned
array, but numpy does, so there is a need for test for this. You can make
an unaligned array like this:

dtype = "b1,f4"
a = numpy.empty(1e4, dtype=dtype)['f1']

I just saw the strides problems that affect only some. I think that the
best explaination is our code:

        /* create appropriate strides for malformed matrices that are
row or column
         * vectors, or empty matrices.
         * In that case, the value of the stride does not really matter, but
         * some versions of BLAS insist that:
         *  - they are not smaller than the number of elements in the array,
         *  - they are not 0.
         */
        sx_0 = (Nx[0] > 1) ? Sx[0]/type_size : (Nx[1] + 1);
        sx_1 = (Nx[1] > 1) ? Sx[1]/type_size : (Nx[0] + 1);
        sy_0 = (Ny[0] > 1) ? Sy[0]/type_size : (Ny[1] + 1);
        sy_1 = (Ny[1] > 1) ? Sy[1]/type_size : (Ny[0] + 1);
        sz_0 = (Nz[0] > 1) ? Sz[0]/type_size : (Nz[1] + 1);
        sz_1 = (Nz[1] > 1) ? Sz[1]/type_size : (Nz[0] + 1);


So this ask for test with empty matrices too.

HTH

Fred
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20121109/94d5e6b2/attachment.html>