Why can't you just use the SymPy Sum function? Aaron Meurer
On Nov 10, 2010, at 2:32 PM, Øyvind Jensen wrote: > on., 10.11.2010 kl. 10.56 -0800, skrev Omar Awile: >> Hi all, > > Hi > >> >> I've recently started playing around with sympy. I am mainly >> interested in the fcode printer module as I would like to be able to >> write symbolic expressions and generate on the fly fortran code that >> is inserted into a fortran subroutine that is part of a larger >> parallel code. >> >> After trying out a couple simple examples that compiled nicely into >> fortran code I tried an example that is slightly more complicated, but >> I can't seem to figure out how to write this properly in sympy. So >> here is what I would like to express >> > >> $dw_p = \frac{V_p}{e^2} \sum_{q=1}^{N}{(w_p - w_q) \mu (x_p - x_q)} >> \forall p=1,...,N$ >> >> In fortran this would be something like >> >> do p=1,n >> do q=1,n >> dw(p) = dw(p) + (w(p) - w(q))*mu*(x(p)-x(q)) >> enddo >> dw(p) = dw(p)*(V(p)/(e**2)) >> enddo > > This loop is currently too complex for the fcode printer, at least > without some hacking. Maybe we can figure something out, and this may > give us an idea on how to make the printers more sophisticated. > > The exception is raised in a routine that tries to figure out which > indices (that is Idx-instances) imply summations, and which are fixed. > Repeated indices imply a summation, and this is all figured out by a > recursion into the subexpressions. A decision had to be made about how > to treat an expression like > > (y[i] + x[j])*z[i] > > I see two alternatives: > > 1) (sum_i y[i]*z[i]) + x[j]*z[i] > > 2) Undefined > > The reason I did not choose option 1) in the current implementation, is > because this leads to the summation of objects with 0 indices and 2 > indices. Since this may not be expected by the user, i preferred to err > on the cautious side, raising the exception whenever a sum contains > terms with different sets of indices. However, a long term goal is to > implement all broadcasting rules of numpy so that we can directly > produce and verify valid numpy expressions. This will give us a well > defined behaviour for option 1, so at some point we can probably disable > the exception safely. > > This is the reason you see the error. Here is a relly ugly workaround > that works only because functions and arrays in Fortran look the same: > In order to avoid the index conformance requirements, you can use > functions and symbols instead of IndexedBase and Idx. Accordingly, I > modified your script to: > > from sympy import IndexedBase, Idx, Symbol, Function > import sympy > w = IndexedBase('w') > dw = IndexedBase('dw') > x = IndexedBase('x') > N = Symbol('N',integer=True) > q = Idx('q',N) > p = Idx('p',N) > V = IndexedBase('V') > mu = Symbol('mu') > nu = Symbol('nu') > eps = Symbol('eps') > # expr = (w[q]-w[p])*mu*(x[q]-x[p]) > wfunc = Function('w') > xfunc = Function('x') > psymb = Symbol('p') > expr = (w[q]-wfunc(psymb))*mu*(x[q]-xfunc(psymb)) > sympy.printing.print_fcode(expr, assign_to=dw[p]) > > > This gives me the loop: > > do p = 1, N > do q = 1, N > dw(p) = mu*(-w(p) + w(q))*(-x(p) + x(q)) + dw(p) > end do > end do > > In order to put in the V(p)/e**2, you can use the same trick, but for > optimization, you should probably do it in a separate call to fcode. > Now, the workaround is ugly, and it may be difficult to automate. So > maybe what we really need in order to deal with complex indexed > expressions, is a better way to represent summation. > > Øyvind > >> >> I tried to accomplish this with following python/sympy code: >> http://codepad.org/O7L70QJV >> >> But when I try to run this I get an exception: >> IndexConformanceException: Indices are not consistent: -w[p] + w[q] >> >> I must be doing something wrong, but I don't quite get what.. >> (actually the code I tried doesn't even contain the $\frac{V_p}{e^2}$, >> but I don't know how to express this properly, apart from using >> sympy.sum, but the fcode printer doesn't like sum.) >> >> I was looking a bit through the documentation but wasn't sure were to >> look actually and I even bugged the devs in the IRC channel who >> adviced me to post this to the mailinglist. So I hope I'm not asking a >> too obvious question... >> >> Thanks in advance! >> >> omar >> > > > -- > 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. > To unsubscribe from this group, send email to > sympy+unsubscr...@googlegroups.com. > 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. To unsubscribe from this group, send email to sympy+unsubscr...@googlegroups.com. For more options, visit this group at http://groups.google.com/group/sympy?hl=en.