On Fri, Nov 25, 2011 at 2:43 PM, vweber <[email protected]> wrote:
> Hi Chris
>
> That looks excellent! I will give it a try.
> I noticed that there might also be a problem with the number of
> continuation lines (for f95 standard it should be 40 lines (20 for
> fixed form),
> it was extended to 256 for the f2003 std).
>
> BTW do you know if there is a way to optimize an equation (e.g. vs
> operation counts)?
> For example:
>
> from f=(x+y)**2+(x+y)**4 (5 operations)
>
> to t1=(x+y); t2=t1**2; f=t1+t2**2 (4 operations)
>
> and at the end I would like to call fcode() for the optimized
> equation.
Here's an equation that could be optimized
>>> exp(x+3)+y*exp(x+3)/(1+exp(x+3))
y*exp(x + 3)/(exp(x + 3) + 1) + exp(x + 3)
Let's create a pair of equations, containing some similar
symbols:
>>> eq = _, _*z
Run them *at the same time* through common substring extraction (cse):
>>> reps, eqs = cse(eq)
Put the fcode lines in a list:
First get numbered symbols to number the equations (if there were not already
lhs, as in this case):
>>> numbered_e = symbols('e:%i'%len(eqs))
>>> flines = []
>>> for ei in reps:
... flines.append(fcode(Eq(*ei)))
...
>>> for i in range(len(eqs)):
... flines.append(fcode(Eq(numbered_e[i], eqs[i])))
...
Let's look at what we've got:
>>> print '\n'.join(flines)
x0 == exp(x + 3)
x1 == x0*y/(x0 + 1) + x0
e0 == x1
e1 == x1*z
If the equations already have preferred lhs symbols, that's fine, too. Let's
run this again defining two equations as 'a' and 'b':
>>> eqs[0] = Eq(a,y*exp(x + 3)/(exp(x + 3) + 1) + exp(x + 3))
>>> eqs[1] = Eq(b,z*(y*exp(x + 3)/(exp(x + 3) + 1) + exp(x + 3)))
>>> for ei in flatten(cse(eqs),levels=1):
... if not isinstance(ei, Expr):
... print fcode(Eq(*ei))
... else:
... print fcode(ei)
...
x0 == exp(x + 3)
x1 == x0*y/(x0 + 1) + x0
a == x1
b == x1*z
HTH,
Chris
--
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.
