Opened 5 years ago

# substitute_function() can leave limits unevaluated

Reported by: Owned by: wonder major sage-6.4 symbolics N/A

### Description

Here is an example:

```┌────────────────────────────────────────────────────────────────────┐
│ Sage Version 6.4.1, Release Date: 2014-11-23                       │
│ Type "notebook()" for the browser-based notebook interface.        │
│ Type "help()" for help.                                            │
└────────────────────────────────────────────────────────────────────┘
sage: l = limit( function('f')(x), x=1 )
sage: l
limit(f(x), x, 1)
sage: ls = l.substitute_function( function('f'), (1+x).function(x) )
sage: ls
limit(x + 1, x, 1)
sage: simplify(ls)
limit(x + 1, x, 1)
sage: maxima(repr(ls))
2
```

This is a simplified case of a situation that's been biting me in an expression that's passed to odeint() after being computed: it breaks the integration because it fails to evaluate to a float expression.

### comment:1 Changed 5 years ago by wonder

Unfortunately that workaround doesn't work when exponentials are involved. The `repr` stage causes the constant `e` to be replaced by an unbound variable `e`:

```sage: l = limit( function('f')(x), x=0.1 )
sage: ls = l.substitute_function( function('f'), (1 - exp(x)).function(x) )
sage: ls
limit(-e^x + 1, x, 0.1)
sage: mls = maxima(repr(ls))
sage: mls
1-e^0.1
sage: N(mls)
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-22-518484fcab94> in <module>()
----> 1 N(mls)

/usr/local/sage/local/lib/python2.7/site-packages/sage/misc/functional.pyc in numerical_approx(x, prec, digits, algorithm)
1466             except (TypeError, ValueError):
1467                 pass
-> 1468         return sage.rings.complex_field.ComplexField(prec)(x)
1469
1470 n = numerical_approx

/usr/local/sage/local/lib/python2.7/site-packages/sage/rings/complex_field.pyc in __call__(self, x, im)
349         if im is not None:
350             x = x, im
--> 351         return Parent.__call__(self, x)
352
353     def _element_constructor_(self, x):

/usr/local/sage/local/lib/python2.7/site-packages/sage/structure/parent.so in sage.structure.parent.Parent.__call__ (build/cythonized/sage/structure/parent.c:9666)()

/usr/local/sage/local/lib/python2.7/site-packages/sage/structure/coerce_maps.so in sage.structure.coerce_maps.DefaultConvertMap_unique._call_ (build/cythonized/sage/structure/coerce_maps.c:4263)()

/usr/local/sage/local/lib/python2.7/site-packages/sage/structure/coerce_maps.so in sage.structure.coerce_maps.DefaultConvertMap_unique._call_ (build/cythonized/sage/structure/coerce_maps.c:4170)()

/usr/local/sage/local/lib/python2.7/site-packages/sage/rings/complex_field.pyc in _element_constructor_(self, x)
388                 pass
389             try:
--> 390                 return x._complex_mpfr_field_( self )
391             except AttributeError:
392                 pass

/usr/local/sage/local/lib/python2.7/site-packages/sage/interfaces/maxima_abstract.pyc in _complex_mpfr_field_(self, C)
1287             8.0751148893563733350506651837615871941533119425962889089783e-62 + 1.4142135623730950488016887242096980785696718753769480731767*I
1288         """
-> 1289         return C(self._sage_())
1290
1291     def _mpfr_(self, R):

/usr/local/sage/local/lib/python2.7/site-packages/sage/rings/complex_field.pyc in __call__(self, x, im)
349         if im is not None:
350             x = x, im
--> 351         return Parent.__call__(self, x)
352
353     def _element_constructor_(self, x):

/usr/local/sage/local/lib/python2.7/site-packages/sage/structure/parent.so in sage.structure.parent.Parent.__call__ (build/cythonized/sage/structure/parent.c:9666)()

/usr/local/sage/local/lib/python2.7/site-packages/sage/structure/coerce_maps.so in sage.structure.coerce_maps.DefaultConvertMap_unique._call_ (build/cythonized/sage/structure/coerce_maps.c:4263)()

/usr/local/sage/local/lib/python2.7/site-packages/sage/structure/coerce_maps.so in sage.structure.coerce_maps.DefaultConvertMap_unique._call_ (build/cythonized/sage/structure/coerce_maps.c:4170)()

/usr/local/sage/local/lib/python2.7/site-packages/sage/rings/complex_field.pyc in _element_constructor_(self, x)
388                 pass
389             try:
--> 390                 return x._complex_mpfr_field_( self )
391             except AttributeError:
392                 pass

/usr/local/sage/local/lib/python2.7/site-packages/sage/symbolic/expression.so in sage.symbolic.expression.Expression._complex_mpfr_field_ (build/cythonized/sage/symbolic/expression.cpp:8185)()

/usr/local/sage/local/lib/python2.7/site-packages/sage/symbolic/expression.so in sage.symbolic.expression.Expression._eval_self (build/cythonized/sage/symbolic/expression.cpp:7521)()

TypeError: Cannot evaluate symbolic expression to a numeric value.
sage: maxima(ls)
'limit(1-%e^_SAGE_VAR_x,_SAGE_VAR_x,0.1)
sage: maxima(SR(mls))
1-_SAGE_VAR__e^0.1
```

Here's what I'm trying out to work around this:

```sage: srls = SR(maxima(repr(ls)))
sage: srls
-_e^0.1 + 1
sage: srls.subs( { SR.symbol('_e'): e } )
-e^0.1 + 1
sage: N(srls.subs( { SR.symbol('_e'): e } ))
-0.105170918075648
```

I've tried poking around the maxima_lib module for a more straightforward way to ask Maxima to evaluate the limit without the data loss involved in the `repr` step, but so far I haven't found my way. Advice would be welcome. Thanks for this wonderfully useful package and for all the bug fixes!

### comment:2 Changed 5 years ago by nbruin

It seems that the inert form of "limit" isn't easily transformed back into one that is considered for evaluation, so once the limit is sitting there as an unevaluated one, substituting into it doesn't give you a limit that is attempted to be evaluated. Perhaps we need a "simplify_limit" that substitutes all the "%LIMIT" symbols for "\$LIMIT" in a maxima expression:

```sage: l = limit( function('f')(x), x=1 )
sage: ls = l.substitute_function( function('f'), (1+x).function(x) )
sage: from sage.libs.ecl import EclObject
sage: from sage.interfaces.maxima_lib import max_limit
sage:
sage: A=maxima_calculus(ls).ecl()
sage: B=EclObject([max_limit]).cons(A.cdr()) #the SIMP doesn't matter
sage: A
<ECL: ((%LIMIT SIMP) ((MPLUS SIMP) 1 |\$_SAGE_VAR_x|) |\$_SAGE_VAR_x| 1)>
sage: B
<ECL: ((\$LIMIT) ((MPLUS SIMP) 1 |\$_SAGE_VAR_x|) |\$_SAGE_VAR_x| 1)>
sage: maxima_calculus(B)
2
```

Alternatively, perhaps we should change the translation of the "limit" operator to be "\$LIMIT" instead of "%LIMIT" (at the cost of re-evaluating any limits when they pass from sage to maxima). I don't know where exactly that conversion happens, because it does seem to be partially in place already:

```sage: maxima_calculus(ls).ecl()
<ECL: ((%LIMIT SIMP) ((MPLUS SIMP) 1 |\$_SAGE_VAR_x|) |\$_SAGE_VAR_x| 1)>
sage: maxima_calculus(ls.operator()).ecl()
<ECL: \$LIMIT>
```

I suspect that this is the culprit:

```sage: L=ls.operator()
sage: L._maxima_lib_init_()
"'limit"
```

(note the explicit "'" there: it's the inert limit). We also have See `sage.symbolic.function_factory.function_factory` which defines the class `NewSymbolicFunction` (is that really a good idea to have that class closed over?) with the relevant `_maxima_init_`.

```sage: sage.interfaces.maxima_lib.sage_op_dict[L]
<ECL: %LIMIT>
```

which gets picked up because the first encounter of the limit symbol is the inert "%limit" on the maxima side, which primes the dict with that translation. This is easily overridden if desired. However, "simplify" and friends don't use `sr_to_max` yet, so it wouldn't make much of a difference here.

In short, we may need to instantiate our own "limit" operator in ginac to accommodate the different translation options, unless we're happy with a "simplify_limit" routine instead.

### comment:3 Changed 5 years ago by wonder

Yes, thank you. I was just converging on the same conclusions about `'limit` and using a python function to force re-evaluation.

I've been assuming that when I see "limit( ..., x, 0.1 )" in a result, it's going to be passed back to maxima in a form that'll re-evaluate the limit. I'd vote for implementing it that way rather than requiring the user to do "simplify_limit" manually at the appropriate times. But I'm sure other users are doing infinite sums and such things that are better not to re-evaluate...

### comment:4 follow-up: ↓ 5 Changed 5 years ago by wonder

Hmm, I didn't realize you could use python functions in `substitute_function`. Just noticed that when looking through `expression_conversions.py`. Given that, this seems to work:

```sage: l = limit( function('f')(x), x=0.1 )
sage: ls = l.substitute_function( function('f'), (1 - exp(-x)).function(x) )
sage: ls
limit(-e^(-x) + 1, x, 0.1)
sage: def eval_limit( ex, var, val ):
kv = { str(var):val }
return limit( ex, **kv )
....:
sage: ls.substitute_function( l.operator(), eval_limit )
0.09516258196404048
```

:-)

Note using `sage.symbolic.function_factory.function('limit')` in place of `l.operator()` does not work.

Update: here it is in `simplify_limits` form:

```limop = limit( SR('f(x)'), x=0 ).operator()
def simplify_limits( expr ):
def eval_limit( expr, var, val ):
kv = { str(var):val }
return limit( expr, **kv )
return expr.substitute_function( limop, eval_limit )
```
Last edited 5 years ago by wonder (previous) (diff)

### comment:5 in reply to: ↑ 4 ; follow-up: ↓ 6 Changed 5 years ago by nbruin

```limop = limit( SR('f(x)'), x=0 ).operator()
def simplify_limits( expr ):
def eval_limit( expr, var, val ):
kv = { str(var):val }
return limit( expr, **kv )
return expr.substitute_function( limop, eval_limit )
```

Smart solution. There is an `eval_limit` with the right interface already available, though:

```limop = limit( SR('f(x)'), x=0 ).operator()
def simplify_limits( expr ):
return expr.substitute_function( limop, maxima_calculus.sr_limit )
```

should do the trick. The fact that

```sage: limsym=sage.symbolic.function_factory.function('limit')
sage: limsym == limop
False
```

indicates that we really should be making a symbol for "limit" beforehand (and equip it with the right translations)

Version 0, edited 5 years ago by nbruin (next)

### comment:6 in reply to: ↑ 5 Changed 5 years ago by wonder

```    return expr.substitute_function( limop, maxima_calculus.sr_limit )