Re: [sympy] Comparing a cse'd expr with the original

[Jason Moore]

David,

I'm not sure I understand what algebraic manipulations you are referring
to. In the example, the algebraic manipulation is done with sympy
variables. Only the final expressions are evaluated with floating point
numbers substituted in. The substitutions of floating point numbers do have
to occur in the order defined by the underlying tree (or directed acyclic
graph in the case of the cse'd expressions). Those evaluations will have
floating point error accumulation but the magnitude of that error should be
small.

Jason
moorepants.info
+01 530-601-9791


On Fri, Oct 23, 2020 at 10:37 PM David Bailey  wrote:

> On 23/10/2020 11:59, Oscar Benjamin wrote:
>
> Hi Jason,
>
> I'm approaching this from the perspective that this is a bug in
> cse/lambdify and that you want to identify which part of a large
> expression tree triggers the bug. Since evaluating the whole
> expression is slow I would start by testing smaller subexpressions and
> work up from there. You can get something like that with:
>
> from sympy.utilities.iterables import postorder_traversal
>
> def cse_compare(expr, syms, vals):
> result = lambdify(syms, expr)(*vals)
>
> syms_all, vals_all = list(syms), list(vals)
> reps, [reduced_expr] = cse(expr)
> for csesym, cseexpr in reps:
> val = lambdify(syms_all, cseexpr)(*vals_all)
> syms_all.append(csesym)
> vals_all.append(val)
> result_cse = lambdify(syms_all, reduced_expr)(*vals_all)
>
> return abs(result - result_cse) / abs(result)
>
> def cse_compare_bottom_up(expr):
> syms = list(expr.free_symbols)
> vals = list(np.random.random(len(syms)))
>
> def subexprs():
> seen = set()
> for e in postorder_traversal(expr):
> if e not in seen:
> yield e
> seen.add(e)
>
> # possibly better:
> #exprs = sorted(set(postorder_traversal(expr)), key=lambda e: 
> e.count_ops())
> exprs = subexprs()
>
> for subexpr in exprs:
> print()
> print(subexpr)
> print(cse_compare(expr, syms, vals))
>
> test_expr = (sin(x+1) + 2*cos(x-3)).expand(trig=True)
> cse_compare_bottom_up(test_expr)
>
>
> Oscar
>
> On Fri, 23 Oct 2020 at 07:59, Jason Moore  
>  wrote:
>
> Oscar,
>
> Yes I can try evaluating in different orders. That's an interesting idea.
>
> I also may be breaking lambdify with too many arguments. I can store the 
> numerical results of each subexpression in a dictionary, then identify the 
> minimal set of variables in each subsequent subexpression, and only lambdify 
> with that very much smaller set of arguments. Maybe that will fix it.
>
> I'll try both.
>
> Jasonmoorepants.info
> +01 530-601-9791
>
>
> On Fri, Oct 23, 2020 at 8:56 AM Oscar Gustafsson  
>  wrote:
>
> Yes, you are right that the email is about a different topic. However, it can 
> be worth noting that evaluating floating-point expressions in different 
> order, which is basically what you do using cse:s, can give different results 
> due to the floating-point quantization/rounding happening at different 
> places. Hence, it is not always obvious to tell if the results are 
> "identical". If you are unlucky, the error propagation can be quite 
> disastrous at times, but you are probably correct in that there is something 
> else that is the primary issue here.
>
> BR Oscar
>
> Den tors 22 okt. 2020 kl 18:57 skrev Jason Moore  
> :
>
> As far as I understand the results should match to machine precision.
>
> We are able to match these kinds of models executed on different platforms 
> using different formulation techniques to machine precision (e-14 or more). 
> In fact, we use this to be able to verify that two independent modelers have 
> created the same model.
>
> The fact that my model can't get the same results with two numerical 
> evaluations of sympy code is concerning.
>
> My question in this email though, is how to improve the comparison code so 
> that it doesn't crash for large expressions.
>
> Jasonmoorepants.info
> +01 530-601-9791
>
>
> On Thu, Oct 22, 2020 at 6:49 PM Oscar Gustafsson  
>  wrote:
>
> How much do they differ? I checked the issue and the latest result wasn't 
> that bad, although probably too much to blame floating-point errors. However, 
> if you push the numerical range into e.g. denormalized numbers I guess it can 
> simply be that. Not very likely though. Is it still that type of numerical 
> errors?
>
> BR Oscar
>
> Den tors 22 okt. 2020 09:37Jason Moore  
>  skrev:
>
> Howdy,
>
> I'm trying to track down a bug in this PR:
> https://github.com/pydy/pydy/pull/122
>
> I have quite large SymPy expressions which I evaluate with floating point 
> numbers. I also cse those expressions and evaluate those with the same 
> floating point numbers. I'm getting discrepancies in the results. I'd like to 
> eliminate the possibility that the cse'd expressions are different from the 
> original expression. I wrote this naive function to make the comparison:
>
> def 

Re: [sympy] Comparing a cse'd expr with the original

[David Bailey]

On 23/10/2020 11:59, Oscar Benjamin wrote:

Hi Jason,

I'm approaching this from the perspective that this is a bug in
cse/lambdify and that you want to identify which part of a large
expression tree triggers the bug. Since evaluating the whole
expression is slow I would start by testing smaller subexpressions and
work up from there. You can get something like that with:

from sympy.utilities.iterables import postorder_traversal

def cse_compare(expr, syms, vals):
 result = lambdify(syms, expr)(*vals)

 syms_all, vals_all = list(syms), list(vals)
 reps, [reduced_expr] = cse(expr)
 for csesym, cseexpr in reps:
 val = lambdify(syms_all, cseexpr)(*vals_all)
 syms_all.append(csesym)
 vals_all.append(val)
 result_cse = lambdify(syms_all, reduced_expr)(*vals_all)

 return abs(result - result_cse) / abs(result)

def cse_compare_bottom_up(expr):
 syms = list(expr.free_symbols)
 vals = list(np.random.random(len(syms)))

 def subexprs():
 seen = set()
 for e in postorder_traversal(expr):
 if e not in seen:
 yield e
 seen.add(e)

 # possibly better:
 #exprs = sorted(set(postorder_traversal(expr)), key=lambda e: 
e.count_ops())
 exprs = subexprs()

 for subexpr in exprs:
 print()
 print(subexpr)
 print(cse_compare(expr, syms, vals))

test_expr = (sin(x+1) + 2*cos(x-3)).expand(trig=True)
cse_compare_bottom_up(test_expr)


Oscar

On Fri, 23 Oct 2020 at 07:59, Jason Moore  wrote:

Oscar,

Yes I can try evaluating in different orders. That's an interesting idea.

I also may be breaking lambdify with too many arguments. I can store the 
numerical results of each subexpression in a dictionary, then identify the 
minimal set of variables in each subsequent subexpression, and only lambdify 
with that very much smaller set of arguments. Maybe that will fix it.

I'll try both.

Jason
moorepants.info
+01 530-601-9791


On Fri, Oct 23, 2020 at 8:56 AM Oscar Gustafsson  
wrote:

Yes, you are right that the email is about a different topic. However, it can be worth 
noting that evaluating floating-point expressions in different order, which is basically 
what you do using cse:s, can give different results due to the floating-point 
quantization/rounding happening at different places. Hence, it is not always obvious to 
tell if the results are "identical". If you are unlucky, the error propagation 
can be quite disastrous at times, but you are probably correct in that there is something 
else that is the primary issue here.

BR Oscar

Den tors 22 okt. 2020 kl 18:57 skrev Jason Moore :

As far as I understand the results should match to machine precision.

We are able to match these kinds of models executed on different platforms 
using different formulation techniques to machine precision (e-14 or more). In 
fact, we use this to be able to verify that two independent modelers have 
created the same model.

The fact that my model can't get the same results with two numerical 
evaluations of sympy code is concerning.

My question in this email though, is how to improve the comparison code so that 
it doesn't crash for large expressions.

Jason
moorepants.info
+01 530-601-9791


On Thu, Oct 22, 2020 at 6:49 PM Oscar Gustafsson  
wrote:

How much do they differ? I checked the issue and the latest result wasn't that 
bad, although probably too much to blame floating-point errors. However, if you 
push the numerical range into e.g. denormalized numbers I guess it can simply 
be that. Not very likely though. Is it still that type of numerical errors?

BR Oscar

Den tors 22 okt. 2020 09:37Jason Moore  skrev:

Howdy,

I'm trying to track down a bug in this PR:

https://github.com/pydy/pydy/pull/122

I have quite large SymPy expressions which I evaluate with floating point 
numbers. I also cse those expressions and evaluate those with the same floating 
point numbers. I'm getting discrepancies in the results. I'd like to eliminate 
the possibility that the cse'd expressions are different from the original 
expression. I wrote this naive function to make the comparison:

def compare_cse(expr):

 args = list(expr.free_symbols)
 args.remove(TIME)
 args += me.find_dynamicsymbols(expr)
 args = list(sm.ordered(args))
 vals = list(np.random.random(len(args)))

 eval_expr = sm.lambdify(args, expr)
 res = eval_expr(*vals)

 replacements, reduced_expr = sm.cse(expr)
 for repl in replacements:
 var = repl[0]
 sub_expr = repl[1]
 eval_sub_expr = sm.lambdify(args, sub_expr)
 vals.append(eval_sub_expr(*vals))
 args.append(var)
 eval_reduced_expr = sm.lambdify(args, reduced_expr[0])
 cse_res = eval_reduced_expr(*vals)

 np.testing.assert_allclose(res, cse_res)

This seems to work for small expressions, but for the large expressions run 
this actually kills my Python interpreter after some minutes.

Is there a 

Re: [sympy] Comparing a cse'd expr with the original

[Jason Moore]

Also, if anyone is curious what prompted this question, I do now have a
nice example running in the pydy docs for this problem:

https://pydy.readthedocs.io/en/latest/examples/carvallo-whipple.html

There are still some things that are incorrect in the derivation (the
warning boxes) but the simulation and animation run. The sympy expressions
for this problem are quite large and it is tricky to get correct.

Jason
moorepants.info
+01 530-601-9791


On Fri, Oct 23, 2020 at 1:29 PM Jason Moore  wrote:

> Oscar,
>
> Yes, this is what I want to check. Thanks for taking the time to write
> this code. I'll give it a try.
>
> In the meantime I did get my code to run. I think Python was actually
> being "Killed" by trying to lambdify an expression with thousands of
> operations. I then switched to expr.evalf(subs=) and it
> wouldn't complete either. Finally, expr.xreplace()
> completes. I also reduced the number of symbolic arguments for the lambdify
> call by ensuring only the ones that were needed were used.
>
> def compare_cse(expr, args=None):
>
> def find_args(expr):
> args = list(expr.free_symbols)  # get constants and time
> try:  # time not explicit, so remove if there
> args.remove(TIME)
> except ValueError:  # time not present
> pass
> args += me.find_dynamicsymbols(expr)  # get the functions of time
> return list(sm.ordered(args))
>
> if args is None:
> args = find_args(expr)
> vals = list(np.random.random(len(args)))
> value_dict = {s: v for s, v in zip(args, vals)}
>
> # This is where it is getting killed! Can't lambdify M it seems.
> #print('lambdifying the whole expression')
> #eval_expr = sm.lambdify(args, expr)
> #print('Evaluate the whole expression')
> #res = eval_expr(*vals)
> # this also doesn't complete (haven't waited for it to complete)
> #res = sm.matrix2numpy(expr.evalf(subs=value_dict), dtype=float)
>
> print('xreplace whole expr')
> res = sm.matrix2numpy(expr.xreplace(value_dict), dtype=float)
>
> print('Find common sub expressions')
> replacements, reduced_expr = sm.cse(expr)
>
> print('Evaluate each subexpression')
> for (var, sub_expr) in replacements:
> print('Sub expression {}'.format(var))
> sub_expr_args = find_args(sub_expr)
> sub_expr_vals = [value_dict[s] for s in sub_expr_args]
> eval_sub_expr = sm.lambdify(sub_expr_args, sub_expr)
> value_dict[var] = eval_sub_expr(*sub_expr_vals)
>
> red_expr_args = find_args(reduced_expr[0])
> red_expr_vals = [value_dict[s] for s in red_expr_args]
> eval_reduced_expr = sm.lambdify(red_expr_args, reduced_expr[0])
> cse_res = eval_reduced_expr(*red_expr_vals)
>
> np.testing.assert_allclose(res, cse_res)
>
> The good thing is that I am getting matching results from an xreplaced
> large expression and the cse'd version of that expression evaluated using
> lambdify. I'm think the bug is now in PyDy's c-code generation (that use's
> cse).
>
> Maybe I can adapt your code to evaluating the c code in some way.
>
> Jason
> moorepants.info
> +01 530-601-9791
>
>
> On Fri, Oct 23, 2020 at 1:00 PM Oscar Benjamin 
> wrote:
>
>> Hi Jason,
>>
>> I'm approaching this from the perspective that this is a bug in
>> cse/lambdify and that you want to identify which part of a large
>> expression tree triggers the bug. Since evaluating the whole
>> expression is slow I would start by testing smaller subexpressions and
>> work up from there. You can get something like that with:
>>
>> from sympy.utilities.iterables import postorder_traversal
>>
>> def cse_compare(expr, syms, vals):
>> result = lambdify(syms, expr)(*vals)
>>
>> syms_all, vals_all = list(syms), list(vals)
>> reps, [reduced_expr] = cse(expr)
>> for csesym, cseexpr in reps:
>> val = lambdify(syms_all, cseexpr)(*vals_all)
>> syms_all.append(csesym)
>> vals_all.append(val)
>> result_cse = lambdify(syms_all, reduced_expr)(*vals_all)
>>
>> return abs(result - result_cse) / abs(result)
>>
>> def cse_compare_bottom_up(expr):
>> syms = list(expr.free_symbols)
>> vals = list(np.random.random(len(syms)))
>>
>> def subexprs():
>> seen = set()
>> for e in postorder_traversal(expr):
>> if e not in seen:
>> yield e
>> seen.add(e)
>>
>> # possibly better:
>> #exprs = sorted(set(postorder_traversal(expr)), key=lambda e:
>> e.count_ops())
>> exprs = subexprs()
>>
>> for subexpr in exprs:
>> print()
>> print(subexpr)
>> print(cse_compare(expr, syms, vals))
>>
>> test_expr = (sin(x+1) + 2*cos(x-3)).expand(trig=True)
>> cse_compare_bottom_up(test_expr)
>>
>>
>> Oscar
>>
>> On Fri, 23 Oct 2020 at 07:59, Jason Moore  wrote:
>> >
>> > Oscar,
>> >
>> > Yes I can try evaluating in different orders. That's an interesting
>> idea.
>> >
>> > I also may be breaking lambdify 

Re: [sympy] Comparing a cse'd expr with the original

[Jason Moore]

Oscar,

Yes, this is what I want to check. Thanks for taking the time to write this
code. I'll give it a try.

In the meantime I did get my code to run. I think Python was actually being
"Killed" by trying to lambdify an expression with thousands of operations.
I then switched to expr.evalf(subs=) and it wouldn't
complete either. Finally, expr.xreplace() completes. I
also reduced the number of symbolic arguments for the lambdify call by
ensuring only the ones that were needed were used.

def compare_cse(expr, args=None):

def find_args(expr):
args = list(expr.free_symbols)  # get constants and time
try:  # time not explicit, so remove if there
args.remove(TIME)
except ValueError:  # time not present
pass
args += me.find_dynamicsymbols(expr)  # get the functions of time
return list(sm.ordered(args))

if args is None:
args = find_args(expr)
vals = list(np.random.random(len(args)))
value_dict = {s: v for s, v in zip(args, vals)}

# This is where it is getting killed! Can't lambdify M it seems.
#print('lambdifying the whole expression')
#eval_expr = sm.lambdify(args, expr)
#print('Evaluate the whole expression')
#res = eval_expr(*vals)
# this also doesn't complete (haven't waited for it to complete)
#res = sm.matrix2numpy(expr.evalf(subs=value_dict), dtype=float)

print('xreplace whole expr')
res = sm.matrix2numpy(expr.xreplace(value_dict), dtype=float)

print('Find common sub expressions')
replacements, reduced_expr = sm.cse(expr)

print('Evaluate each subexpression')
for (var, sub_expr) in replacements:
print('Sub expression {}'.format(var))
sub_expr_args = find_args(sub_expr)
sub_expr_vals = [value_dict[s] for s in sub_expr_args]
eval_sub_expr = sm.lambdify(sub_expr_args, sub_expr)
value_dict[var] = eval_sub_expr(*sub_expr_vals)

red_expr_args = find_args(reduced_expr[0])
red_expr_vals = [value_dict[s] for s in red_expr_args]
eval_reduced_expr = sm.lambdify(red_expr_args, reduced_expr[0])
cse_res = eval_reduced_expr(*red_expr_vals)

np.testing.assert_allclose(res, cse_res)

The good thing is that I am getting matching results from an xreplaced
large expression and the cse'd version of that expression evaluated using
lambdify. I'm think the bug is now in PyDy's c-code generation (that use's
cse).

Maybe I can adapt your code to evaluating the c code in some way.

Jason
moorepants.info
+01 530-601-9791


On Fri, Oct 23, 2020 at 1:00 PM Oscar Benjamin 
wrote:

> Hi Jason,
>
> I'm approaching this from the perspective that this is a bug in
> cse/lambdify and that you want to identify which part of a large
> expression tree triggers the bug. Since evaluating the whole
> expression is slow I would start by testing smaller subexpressions and
> work up from there. You can get something like that with:
>
> from sympy.utilities.iterables import postorder_traversal
>
> def cse_compare(expr, syms, vals):
> result = lambdify(syms, expr)(*vals)
>
> syms_all, vals_all = list(syms), list(vals)
> reps, [reduced_expr] = cse(expr)
> for csesym, cseexpr in reps:
> val = lambdify(syms_all, cseexpr)(*vals_all)
> syms_all.append(csesym)
> vals_all.append(val)
> result_cse = lambdify(syms_all, reduced_expr)(*vals_all)
>
> return abs(result - result_cse) / abs(result)
>
> def cse_compare_bottom_up(expr):
> syms = list(expr.free_symbols)
> vals = list(np.random.random(len(syms)))
>
> def subexprs():
> seen = set()
> for e in postorder_traversal(expr):
> if e not in seen:
> yield e
> seen.add(e)
>
> # possibly better:
> #exprs = sorted(set(postorder_traversal(expr)), key=lambda e:
> e.count_ops())
> exprs = subexprs()
>
> for subexpr in exprs:
> print()
> print(subexpr)
> print(cse_compare(expr, syms, vals))
>
> test_expr = (sin(x+1) + 2*cos(x-3)).expand(trig=True)
> cse_compare_bottom_up(test_expr)
>
>
> Oscar
>
> On Fri, 23 Oct 2020 at 07:59, Jason Moore  wrote:
> >
> > Oscar,
> >
> > Yes I can try evaluating in different orders. That's an interesting idea.
> >
> > I also may be breaking lambdify with too many arguments. I can store the
> numerical results of each subexpression in a dictionary, then identify the
> minimal set of variables in each subsequent subexpression, and only
> lambdify with that very much smaller set of arguments. Maybe that will fix
> it.
> >
> > I'll try both.
> >
> > Jason
> > moorepants.info
> > +01 530-601-9791
> >
> >
> > On Fri, Oct 23, 2020 at 8:56 AM Oscar Gustafsson <
> oscar.gustafs...@gmail.com> wrote:
> >>
> >> Yes, you are right that the email is about a different topic. However,
> it can be worth noting that evaluating floating-point expressions in
> different order, which is basically what you do using cse:s, can give
> different results 

Re: [sympy] Comparing a cse'd expr with the original

[Oscar Benjamin]

Hi Jason,

I'm approaching this from the perspective that this is a bug in
cse/lambdify and that you want to identify which part of a large
expression tree triggers the bug. Since evaluating the whole
expression is slow I would start by testing smaller subexpressions and
work up from there. You can get something like that with:

from sympy.utilities.iterables import postorder_traversal

def cse_compare(expr, syms, vals):
result = lambdify(syms, expr)(*vals)

syms_all, vals_all = list(syms), list(vals)
reps, [reduced_expr] = cse(expr)
for csesym, cseexpr in reps:
val = lambdify(syms_all, cseexpr)(*vals_all)
syms_all.append(csesym)
vals_all.append(val)
result_cse = lambdify(syms_all, reduced_expr)(*vals_all)

return abs(result - result_cse) / abs(result)

def cse_compare_bottom_up(expr):
syms = list(expr.free_symbols)
vals = list(np.random.random(len(syms)))

def subexprs():
seen = set()
for e in postorder_traversal(expr):
if e not in seen:
yield e
seen.add(e)

# possibly better:
#exprs = sorted(set(postorder_traversal(expr)), key=lambda e: e.count_ops())
exprs = subexprs()

for subexpr in exprs:
print()
print(subexpr)
print(cse_compare(expr, syms, vals))

test_expr = (sin(x+1) + 2*cos(x-3)).expand(trig=True)
cse_compare_bottom_up(test_expr)


Oscar

On Fri, 23 Oct 2020 at 07:59, Jason Moore  wrote:
>
> Oscar,
>
> Yes I can try evaluating in different orders. That's an interesting idea.
>
> I also may be breaking lambdify with too many arguments. I can store the 
> numerical results of each subexpression in a dictionary, then identify the 
> minimal set of variables in each subsequent subexpression, and only lambdify 
> with that very much smaller set of arguments. Maybe that will fix it.
>
> I'll try both.
>
> Jason
> moorepants.info
> +01 530-601-9791
>
>
> On Fri, Oct 23, 2020 at 8:56 AM Oscar Gustafsson  
> wrote:
>>
>> Yes, you are right that the email is about a different topic. However, it 
>> can be worth noting that evaluating floating-point expressions in different 
>> order, which is basically what you do using cse:s, can give different 
>> results due to the floating-point quantization/rounding happening at 
>> different places. Hence, it is not always obvious to tell if the results are 
>> "identical". If you are unlucky, the error propagation can be quite 
>> disastrous at times, but you are probably correct in that there is something 
>> else that is the primary issue here.
>>
>> BR Oscar
>>
>> Den tors 22 okt. 2020 kl 18:57 skrev Jason Moore :
>>>
>>> As far as I understand the results should match to machine precision.
>>>
>>> We are able to match these kinds of models executed on different platforms 
>>> using different formulation techniques to machine precision (e-14 or more). 
>>> In fact, we use this to be able to verify that two independent modelers 
>>> have created the same model.
>>>
>>> The fact that my model can't get the same results with two numerical 
>>> evaluations of sympy code is concerning.
>>>
>>> My question in this email though, is how to improve the comparison code so 
>>> that it doesn't crash for large expressions.
>>>
>>> Jason
>>> moorepants.info
>>> +01 530-601-9791
>>>
>>>
>>> On Thu, Oct 22, 2020 at 6:49 PM Oscar Gustafsson 
>>>  wrote:

 How much do they differ? I checked the issue and the latest result wasn't 
 that bad, although probably too much to blame floating-point errors. 
 However, if you push the numerical range into e.g. denormalized numbers I 
 guess it can simply be that. Not very likely though. Is it still that type 
 of numerical errors?

 BR Oscar

 Den tors 22 okt. 2020 09:37Jason Moore  skrev:
>
> Howdy,
>
> I'm trying to track down a bug in this PR:
>
> https://github.com/pydy/pydy/pull/122
>
> I have quite large SymPy expressions which I evaluate with floating point 
> numbers. I also cse those expressions and evaluate those with the same 
> floating point numbers. I'm getting discrepancies in the results. I'd 
> like to eliminate the possibility that the cse'd expressions are 
> different from the original expression. I wrote this naive function to 
> make the comparison:
>
> def compare_cse(expr):
>
> args = list(expr.free_symbols)
> args.remove(TIME)
> args += me.find_dynamicsymbols(expr)
> args = list(sm.ordered(args))
> vals = list(np.random.random(len(args)))
>
> eval_expr = sm.lambdify(args, expr)
> res = eval_expr(*vals)
>
> replacements, reduced_expr = sm.cse(expr)
> for repl in replacements:
> var = repl[0]
> sub_expr = repl[1]
> eval_sub_expr = sm.lambdify(args, sub_expr)
> vals.append(eval_sub_expr(*vals))
> args.append(var)
>   

Re: [sympy] Comparing a cse'd expr with the original

[Jason Moore]

Oscar,

Yes I can try evaluating in different orders. That's an interesting idea.

I also may be breaking lambdify with too many arguments. I can store the
numerical results of each subexpression in a dictionary, then identify the
minimal set of variables in each subsequent subexpression, and only
lambdify with that very much smaller set of arguments. Maybe that will fix
it.

I'll try both.

Jason
moorepants.info
+01 530-601-9791


On Fri, Oct 23, 2020 at 8:56 AM Oscar Gustafsson 
wrote:

> Yes, you are right that the email is about a different topic. However, it
> can be worth noting that evaluating floating-point expressions in different
> order, which is basically what you do using cse:s, can give different
> results due to the floating-point quantization/rounding happening at
> different places. Hence, it is not always obvious to tell if the results
> are "identical". If you are unlucky, the error propagation can be quite
> disastrous at times, but you are probably correct in that there is
> something else that is the primary issue here.
>
> BR Oscar
>
> Den tors 22 okt. 2020 kl 18:57 skrev Jason Moore :
>
>> As far as I understand the results should match to machine precision.
>>
>> We are able to match these kinds of models executed on different
>> platforms using different formulation techniques to machine precision (e-14
>> or more). In fact, we use this to be able to verify that two independent
>> modelers have created the same model.
>>
>> The fact that my model can't get the same results with two numerical
>> evaluations of sympy code is concerning.
>>
>> My question in this email though, is how to improve the comparison code
>> so that it doesn't crash for large expressions.
>>
>> Jason
>> moorepants.info
>> +01 530-601-9791
>>
>>
>> On Thu, Oct 22, 2020 at 6:49 PM Oscar Gustafsson <
>> oscar.gustafs...@gmail.com> wrote:
>>
>>> How much do they differ? I checked the issue and the latest result
>>> wasn't that bad, although probably too much to blame floating-point errors.
>>> However, if you push the numerical range into e.g. denormalized numbers I
>>> guess it can simply be that. Not very likely though. Is it still that type
>>> of numerical errors?
>>>
>>> BR Oscar
>>>
>>> Den tors 22 okt. 2020 09:37Jason Moore  skrev:
>>>
 Howdy,

 I'm trying to track down a bug in this PR:

 https://github.com/pydy/pydy/pull/122

 I have quite large SymPy expressions which I evaluate with floating
 point numbers. I also cse those expressions and evaluate those with the
 same floating point numbers. I'm getting discrepancies in the results. I'd
 like to eliminate the possibility that the cse'd expressions are different
 from the original expression. I wrote this naive function to make the
 comparison:

 def compare_cse(expr):

 args = list(expr.free_symbols)
 args.remove(TIME)
 args += me.find_dynamicsymbols(expr)
 args = list(sm.ordered(args))
 vals = list(np.random.random(len(args)))

 eval_expr = sm.lambdify(args, expr)
 res = eval_expr(*vals)

 replacements, reduced_expr = sm.cse(expr)
 for repl in replacements:
 var = repl[0]
 sub_expr = repl[1]
 eval_sub_expr = sm.lambdify(args, sub_expr)
 vals.append(eval_sub_expr(*vals))
 args.append(var)
 eval_reduced_expr = sm.lambdify(args, reduced_expr[0])
 cse_res = eval_reduced_expr(*vals)

 np.testing.assert_allclose(res, cse_res)

 This seems to work for small expressions, but for the large expressions
 run this actually kills my Python interpreter after some minutes.

 Is there a better way to write this? I need some way to verify that
 numerical evaluation of a cse'd expression gives the same result as
 numerical evaluation of the expression that can actually complete in a
 reasonable amount of time.

 Jason
 moorepants.info
 +01 530-601-9791

 --
 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 view this discussion on the web visit
 https://groups.google.com/d/msgid/sympy/CAP7f1Aivs53AcwyaoyO5%2Bzkn%2BPrX3TfoJJ5X85qQNZ6JRbo5JQ%40mail.gmail.com
 
 .

>>> --
>>> 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 view this discussion on the web visit
>>> https://groups.google.com/d/msgid/sympy/CAFjzj-LQNBFoF7Q4Rb%2BZG3J35-Et9b__SED-DMJvwUYqyKYADQ%40mail.gmail.com

Re: [sympy] Comparing a cse'd expr with the original

[Oscar Gustafsson]

Yes, you are right that the email is about a different topic. However, it
can be worth noting that evaluating floating-point expressions in different
order, which is basically what you do using cse:s, can give different
results due to the floating-point quantization/rounding happening at
different places. Hence, it is not always obvious to tell if the results
are "identical". If you are unlucky, the error propagation can be quite
disastrous at times, but you are probably correct in that there is
something else that is the primary issue here.

BR Oscar

Den tors 22 okt. 2020 kl 18:57 skrev Jason Moore :

> As far as I understand the results should match to machine precision.
>
> We are able to match these kinds of models executed on different platforms
> using different formulation techniques to machine precision (e-14 or more).
> In fact, we use this to be able to verify that two independent modelers
> have created the same model.
>
> The fact that my model can't get the same results with two numerical
> evaluations of sympy code is concerning.
>
> My question in this email though, is how to improve the comparison code so
> that it doesn't crash for large expressions.
>
> Jason
> moorepants.info
> +01 530-601-9791
>
>
> On Thu, Oct 22, 2020 at 6:49 PM Oscar Gustafsson <
> oscar.gustafs...@gmail.com> wrote:
>
>> How much do they differ? I checked the issue and the latest result wasn't
>> that bad, although probably too much to blame floating-point errors.
>> However, if you push the numerical range into e.g. denormalized numbers I
>> guess it can simply be that. Not very likely though. Is it still that type
>> of numerical errors?
>>
>> BR Oscar
>>
>> Den tors 22 okt. 2020 09:37Jason Moore  skrev:
>>
>>> Howdy,
>>>
>>> I'm trying to track down a bug in this PR:
>>>
>>> https://github.com/pydy/pydy/pull/122
>>>
>>> I have quite large SymPy expressions which I evaluate with floating
>>> point numbers. I also cse those expressions and evaluate those with the
>>> same floating point numbers. I'm getting discrepancies in the results. I'd
>>> like to eliminate the possibility that the cse'd expressions are different
>>> from the original expression. I wrote this naive function to make the
>>> comparison:
>>>
>>> def compare_cse(expr):
>>>
>>> args = list(expr.free_symbols)
>>> args.remove(TIME)
>>> args += me.find_dynamicsymbols(expr)
>>> args = list(sm.ordered(args))
>>> vals = list(np.random.random(len(args)))
>>>
>>> eval_expr = sm.lambdify(args, expr)
>>> res = eval_expr(*vals)
>>>
>>> replacements, reduced_expr = sm.cse(expr)
>>> for repl in replacements:
>>> var = repl[0]
>>> sub_expr = repl[1]
>>> eval_sub_expr = sm.lambdify(args, sub_expr)
>>> vals.append(eval_sub_expr(*vals))
>>> args.append(var)
>>> eval_reduced_expr = sm.lambdify(args, reduced_expr[0])
>>> cse_res = eval_reduced_expr(*vals)
>>>
>>> np.testing.assert_allclose(res, cse_res)
>>>
>>> This seems to work for small expressions, but for the large expressions
>>> run this actually kills my Python interpreter after some minutes.
>>>
>>> Is there a better way to write this? I need some way to verify that
>>> numerical evaluation of a cse'd expression gives the same result as
>>> numerical evaluation of the expression that can actually complete in a
>>> reasonable amount of time.
>>>
>>> Jason
>>> moorepants.info
>>> +01 530-601-9791
>>>
>>> --
>>> 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 view this discussion on the web visit
>>> https://groups.google.com/d/msgid/sympy/CAP7f1Aivs53AcwyaoyO5%2Bzkn%2BPrX3TfoJJ5X85qQNZ6JRbo5JQ%40mail.gmail.com
>>> 
>>> .
>>>
>> --
>> 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 view this discussion on the web visit
>> https://groups.google.com/d/msgid/sympy/CAFjzj-LQNBFoF7Q4Rb%2BZG3J35-Et9b__SED-DMJvwUYqyKYADQ%40mail.gmail.com
>> 
>> .
>>
> --
> 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 view this discussion on the web visit
> https://groups.google.com/d/msgid/sympy/CAP7f1AgdkaC0D69s_ohka_AjDFvKCD%2BMoVKyFhs0Vn2b_DtFkw%40mail.gmail.com
>