Thank you for working on this. I think that this is important work.
I would express my views by they happen to be completely in line with Tim's
response which mostly says that your intuition on this problem seems
sensible.
Thanks again,
-Matt
On Mon, Aug 4, 2014 at 4:13 PM, Jason Moore <moorepa...@gmail.com> wrote:
>
>
>
> Jason
> moorepants.info
> +01 530-601-9791
>
>
> On Mon, Aug 4, 2014 at 3:27 PM, James Crist <crist...@umn.edu> 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 sympy+unsubscr...@googlegroups.com.
>> To post to this group, send email to sympy@googlegroups.com.
>> 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 sympy+unsubscr...@googlegroups.com.
> To post to this group, send email to sympy@googlegroups.com.
> 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
> <https://groups.google.com/d/msgid/sympy/CAP7f1AgqKNa5XG1rJbyMZYktJ3_Q4-P%3Di03HTiLSqbBpepQ8xw%40mail.gmail.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 sympy+unsubscr...@googlegroups.com.
To post to this group, send email to sympy@googlegroups.com.
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/CAJ8oX-GZvp8RAxvOUA5ojzcDfTU_iXH73Ju860xX24JnC2_cEA%40mail.gmail.com.
For more options, visit https://groups.google.com/d/optout.
