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 [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.
