Jason moorepants.info +01 530-601-9791
On Mon, Aug 4, 2014 at 3:27 PM, James Crist <[email protected]> wrote: > I'm working on adding support ofr codegeneration with Matrix objects. > Currently an `indexed` type is supported that results in low-level > contiguous arrays. These are always converted into loops though (and I > don't really understand what they're for). In contrast, the intent here is > to provide support for generating matrices that are comprised of sympy > expressions (as seen frequently in `mechanics`). For example, the following > should *work* after this update: > > >>> mat = Matrix([sin(a), cos(b), tan(c)]) > >>> func = autowrap(mat) > >>> func(1, 2, 3) > array([ 0.84147098, -0.41614684, -0.14254654]) > > I have some stuff working already, but have some question before I > progress further. > > *1. Sympy Matrices are always 2 dimensional, should this be true of the > generated code as well?* > > In numpy, the default is the minimum number of dimensions required. For > example, `array([1, 2, 3])` has only one dimension. In contrast, with sympy > `Matrix([1, 2, 3])` has 2 dimensions always (no way around). For > expressions that are inherently column/row vectors should a single > dimension array be created? > > *- Pros:* Less nested data, many scipy routines require single dimension > arrays only (odeint) > *- Cons:* Multiplication and indexing will be different. For this reason > I'm kind of against switching, but I could be swayed. For example: > > >>> np_a = array([1, 2, 3]) > >>> sym_a = Matrix([1, 2, 3]) > >>> np_a.dot(np_a.T) # Perform Multiplication as done in > numpy > 14 > >>> sym_a * sym_a.T # Perform Multiplication as done in > sympy > Matrix([ > [1, 2, 3], > [2, 4, 6], > [3, 6, 9]]) > >>> np_a_res = np_a.reshape((3,1)) # Reshape the array to be 2 > dimensional > >>> np_a_res.dot(np_a_res.T) # Multiply the 2 dimensional > numpy arrays. This results in the same as sympy > array([ > [1, 2, 3], > [2, 4, 6], > [3, 6, 9]]) > >>> np_a[1, 0] > IndexError Traceback (most recent call last > ) > <ipython-input-265-9a4f4bee3353> in <module>() > ----> 1 np_a[1,0] > > IndexError: too many indices > >>> np_a_res[1,0] # Indexing works the same in > 2d as it would in sympy normally > 2 > > > SymPy only has 2D arrays (matrices) and column vectors are still represented as 2D arrays (like Matlab). So I think we should stick to that paradigm. You could have a flag that says "make vectors 1D arrays" if you want to support that. But the C code in the backend, for example, is all based on 1D arrays. The indexing in C to nD arrays is just syntactic sugar. I suspect that is the same in fortran too. Many numpy functions will broadcast so it may be the case that 2D column vectors will work in most cases, but otherwise the user may need to use .squeeze() to remove the unneeded dimension. It makes sense to me to retain the same indexing in the numpy array as the sympy matrix. > > *2. Should the default be numpy arrays?* > > I think yes. They are pretty much used everywhere in scientific python, > and should be supported. numpy.matrix is on its way out, and should not be > used in my opinion. Cython offers some support for generality, so that > anything that offers the buffer protocol can be used. f2py may do the same, > but I have little experience with it. Currently autowrap supports lists as > well. I think this is silly, and results in some inefficiencies. As they're > transformed into numpy.array internally in the call, the user must have > numpy installed, and should be able to do this themselves. > Yes, we should default to outputting numpy arrays and not support the numpy matrix type. The numpy matrix type is going to eventually be deprecated in numpy. No reason to use it in anything new. > *3. For inputting matrices to functions, `MatrixSymbol`, or > `DeferredVector`?* > > Sometimes you may have a lot of input variables, have those variables > expressed as a vector of *known* length. Sympy offers two options for this: > `MatrixSymbol` and `DeferredVector`. They both seem to offer the same > intention, although `MatrixSymbol` seems to be used *way* more/have more > functionality. Currently I'm only supporting `MatrixSymbol`. The idea is > that this should be possible: > > >>> x = MatrixSymbol('x', 3, 1) #Acts as a vector > > > > >>> expr = x[1,1] + sin(x[2,1]) + cos(x[3,1]) > > > > >>> func = autowrap(expr) > > > > >>> inp = np.array([1, 2, 3]) > > > > >>> func(inp) > > 0.9193049302252362 > > > Note that this is the inverse of the issue described in #1. This time it's > the function input that will have dimesion mismatch between what is being > input (a single dimension numpy vector, compared to a 2 dimension sympy > MatrixSymbol/Matrix. Again, the numpy array can be reshaped either > externally or internally with the magic of cython (this won't result in a > copy operation either, just a new memoryview). Still, I'd like an opinion > on this. > I don't think this is really a mismatch. This operation seems more like you are treating the things in inp as a list of values. The output should be defined by expr, which in this case is a scalar. This looks fine. > > *4. What is an Indexed type for?* > > Currently I'm ignoring these, and allowing them to stay as they are while > writing functionality around them for Matrix types. Still, the code > supporting them is messy (maybe it has to be?) and I can't refactor it into > something cleaner until I understand what they do. I've read the docs on > this functionality and see that you can make a Matrix Mutliplication > function. Great. But I plan on supporting that later on as well with > MatrixExpr types and calls to blas routines. There has to be a purpose for > these, I just don't understand it yet. > I think Björn Dahlgren has been working on this stuff. Maybe speaking with him could help clarify this. I don't really get it either. > > *5. Would a ctypes CodeWrapper be wanted?* (Not relevant immediately, but > I'm curious) > > Currently only Cython and f2py are supported for wrapper objects. these > make very robust solutions, but neither is in the standard library. > `ctypes` is a standard library module that provides some functionality of > wrappers around shared libraries (.so or .dll). Not nearly as robust, and > no built in type conversion. Pros: everyone has it, Cons: it's not as good. > I think I'll add this in later. > This could certainly be useful, especially if you are using SymPy without the optional dependencies. I think the autowrap module idea is to have a flexible framework that can use a variety of wrappers and language combinations. For cPython there are are several that are popular: cython, swig, ctypes, boost, etc. But the main goal would to have it working for the two currently supported combinations: C/Cython, Fortran/F2py. > > If you have thoughts on any of these, please let me know! I want to make > this useful for everyone, not just myself :) > > -- > You received this message because you are subscribed to the Google Groups > "sympy" group. > To unsubscribe from this group and stop receiving emails from it, send an > email to [email protected]. > To post to this group, send email to [email protected]. > Visit this group at http://groups.google.com/group/sympy. > To view this discussion on the web visit > https://groups.google.com/d/msgid/sympy/b362375c-2caa-41ee-a176-841d14147f12%40googlegroups.com > <https://groups.google.com/d/msgid/sympy/b362375c-2caa-41ee-a176-841d14147f12%40googlegroups.com?utm_medium=email&utm_source=footer> > . > For more options, visit https://groups.google.com/d/optout. > -- You received this message because you are subscribed to the Google Groups "sympy" group. To unsubscribe from this group and stop receiving emails from it, send an email to [email protected]. To post to this group, send email to [email protected]. Visit this group at http://groups.google.com/group/sympy. To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAP7f1AgqKNa5XG1rJbyMZYktJ3_Q4-P%3Di03HTiLSqbBpepQ8xw%40mail.gmail.com. For more options, visit https://groups.google.com/d/optout.
