On Thursday, 16 January 2014 14:24:41 UTC+1, Björn Dahlgren wrote: > > On Saturday, 28 July 2012 21:45:20 UTC+2, Aaron Meurer wrote: >> >> On Fri, Jul 27, 2012 at 7:52 PM, Matthew Turk <[email protected]> >> wrote: >> > Hi Aaron, >> > >> > Thanks very much for the quick reply. I'm fairly new to sympy, and >> > this is my first time looking at its internals. >> > >> > On Fri, Jul 27, 2012 at 5:01 PM, Aaron Meurer <[email protected]> >> wrote: >> >> I did some digging in the code, and I didn't find a way to do it. It >> >> seems that it is pretty hardcoded that Indexed is not just an indexed >> >> symbol, but a tensor. What we really should do is create a subclass >> >> IndexedTensor and only check for that in get_contraction_structure in >> >> sympy/tensor/index_methods.py, and require the use of that if you want >> >> automatic tensor contraction. >> > >> > On your suggestion I dug into the code, and it looks like by >> > implementing IndexedTensor, one could get the behavior I want by >> > performing a check like: >> > >> > if not any(isinstance(a, IndexedTensor) for a in expr.args): >> > return (set([]), {}, ()) >> > >> > inside _get_indices_Mul. It wasn't clear to me that I could get the >> > correct behavior from get_contraction_structure with what I knew of >> > how SymPy worked. This solution seems somewhat unsatisfying, however. >> > >> > Is there a way to determine whether or not the indices should be >> > iterated over just by the state of the variable "inds" inside that >> > routine? By checking if the length of inds was zero and the length of >> > dummies greater than zero (and returning empty values), I was able to >> > get the behavior I wanted, but it causes two tests to fail, so that >> > won't do. >> > >> >> >> >> Somewhat related is >> >> http://code.google.com/p/sympy/issues/detail?id=2659. >> >> Indexed/IndexedBase should just be generic and additional, tensor-like >> >> assumptions should be in subclasses. >> >> >> >> Feel free to open an issue for this, or even give it a shot >> implementing it. >> > >> > I've opened an issue -- I'd be happy to try to contribute, but I am >> > hesitant because of the scope of the changes to replace IndexedBase >> > with IndexedTensor everywhere expansion is necessary. Do you think >> > there's a way to deduce the correct behavior inside _get_indices_Mul? >> >> I suppose another way would be to create an option for it in the code >> printers, and pass the option through to codegen. >> >> Aaron Meurer >> >> > >> > Thanks again, >> > >> > Matt >> > >> >> >> >> Aaron Meurer >> >> >> >> On Fri, Jul 27, 2012 at 1:58 PM, Matthew Turk <[email protected]> >> wrote: >> >>> Hi all, >> >>> >> >>> I'm trying to get some code generation working, and I was hoping to >> >>> get some clarification about how IndexedBase objects work inside >> >>> codegen. For instance, with this code: >> >>> >> >>> import sympy >> >>> from sympy.utilities.codegen import codegen >> >>> A, B, C, D = map(sympy.IndexedBase, "ABCD") >> >>> m = sympy.Symbol("m", integer=True) >> >>> i = sympy.Idx("i", m) >> >>> expr = sympy.Equality(D[i], A[i]*B[i]+C[i]) >> >>> rv = codegen( ("f1", expr), "C", "temp", header=False) >> >>> print rv[0][1] >> >>> >> >>> I would naively expect each element in D to be assigned >> >>> A[i]*B[i]+C[i]. However, what comes out of codegen iterates over the >> >>> index twice: >> >>> >> >>> [snip] >> >>> for (int i=0; i<m; i++){ >> >>> for (int i=0; i<m; i++){ >> >>> D[i] = A[i]*B[i] + D[i]; >> >>> } >> >>> } >> >>> [/snip] >> >>> >> >>> Everything else looks right, including the assignment to the passed >> >>> pointer. If I change the expression to A[i]+C[i], so there's no >> >>> multiplication, the loop looks like what I'd expect: >> >>> >> >>> [snip] >> >>> for (int i=0; i<m; i++){ >> >>> D[i] = A[i] + C[i]; >> >>> } >> >>> [/snip] >> >>> >> >>> It looks as though the first expression is potentially trying to >> >>> evaluate the components as tensors. Is there a way to modify what >> I'm >> >>> doing in the multiplicative expression to indicate I'm actually >> >>> multiplying flat arrays, and get a single loop rather than the double >> >>> nested one? >> >>> >> >>> Thanks very much, >> >>> >> >>> Matt >> >>> >> >>> -- >> >>> You received this message because you are subscribed to the Google >> Groups "sympy" group. >> >>> To post to this group, send email to [email protected]. >> >>> To unsubscribe from this group, send email to >> [email protected]. >> >>> For more options, visit this group at >> http://groups.google.com/group/sympy?hl=en. >> >>> >> >> >> >> -- >> >> You received this message because you are subscribed to the Google >> Groups "sympy" group. >> >> To post to this group, send email to [email protected]. >> >> To unsubscribe from this group, send email to >> [email protected]. >> >> For more options, visit this group at >> http://groups.google.com/group/sympy?hl=en. >> >> >> > >> > -- >> > You received this message because you are subscribed to the Google >> Groups "sympy" group. >> > To post to this group, send email to [email protected]. >> > To unsubscribe from this group, send email to >> [email protected]. >> > For more options, visit this group at >> http://groups.google.com/group/sympy?hl=en. >> > >> > > So I have been pondering about this for a while, > > IndexedBase and Indexed are quite general but the code printers really do > assume them to > be tensors. I can think of tons of cases where this is not the case. For > finite differences e.g. > > Say that I want to express the simple forward difference stencil for the > first derivative: > > len_y=5 > > y = IndexedBase('y', shape=(len_y,)) > > x = IndexedBase('x', shape=(len_y,)) > > Dy = IndexedBase('Dy', shape=(len_y-1,)) > > i = Idx('i', len_A-1) > > > > # Forward difference > > e=Eq(Dy[i], (y[i+1]-y[i])/(x[i+1]-x[i])) > > e > > > Out[16]: > Dy[i]=(y[i+1]−y[i])(x[i+1]−x[i])−1 > > ccode(e.rhs, assign_to=e.lhs) > > IndexConformanceException: Indices are not consistent: y[i + 1] - y[i] > > > So if I want the above to work I would pretty much need to introduce a > SimpleIndexed and SimpleIndexBase classes > and then patch the codeprinters for this case but it is quite backwards if > you ask me. > Have anyone else made any progress in this direction? > > Writing a separate code printer not assuming tensor behaviour would also be > an alternative (or letting tensor behaviour be a flag). > > What would be the preferred approach here? > >
Ok so I made an attempt to allow a user to bypass the current default behaviour: https://github.com/sympy/sympy/pull/2780 Please let me know what you think! Best regards /Björn -- 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. For more options, visit https://groups.google.com/groups/opt_out.
