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
*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.
*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.
*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.
*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.
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.
For more options, visit https://groups.google.com/d/optout.
