On Saturday, 28 July 2012 21:45:20 UTC+2, Aaron Meurer wrote:
>
> On Fri, Jul 27, 2012 at 7:52 PM, Matthew Turk
> <matth...@gmail.com<javascript:>>
> 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
> > <asme...@gmail.com<javascript:>>
> 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
> >> <matth...@gmail.com<javascript:>>
> 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 sy...@googlegroups.com<javascript:>.
>
> >>> To unsubscribe from this group, send email to
> sympy+un...@googlegroups.com <javascript:>.
> >>> 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 sy...@googlegroups.com<javascript:>.
>
> >> To unsubscribe from this group, send email to
> sympy+un...@googlegroups.com <javascript:>.
> >> 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 sy...@googlegroups.com<javascript:>.
>
> > To unsubscribe from this group, send email to
> sympy+un...@googlegroups.com <javascript:>.
> > 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?
--
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.
For more options, visit https://groups.google.com/groups/opt_out.
