substitute_function() can leave limits unevaluated

[wonder]

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.

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

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

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

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

[nbruin]

Replying to wonder:

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)

CORRECTION: The symbol to use is sage.calculus.calculus._limit and the reason for the inequality is an extra attribute for proper latex representation of limits:

sage: limsym=sage.symbolic.function_factory.function('limit',print_latex_func=sage.calculus.calculus._limit_latex_)
sage: limsym == limop
True

[wonder]

Replying to nbruin:

    return expr.substitute_function( limop, maxima_calculus.sr_limit )

Thanks! Yes, that also works.