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.