sympy / sympy

A computer algebra system written in pure Python
https://sympy.org/
Other
12.98k stars 4.43k forks source link

Derivatives of symbols with implicit dependencies #20096

Open ThePauliPrinciple opened 4 years ago

ThePauliPrinciple commented 4 years ago

Some symbols have implicit dependencies on other symbols that is not kept track of in sympy. For example, if one defines an indexed y[i]=x[i]**2, then a derivative of y[i] w.r.t. x[i] should exist. Currently sympy provides to option for this to the best of my ability. On can use Derivative, however, one can then not evaluate the derivative of something like y[i]=exp(x[i]) as .doit() will simply result in 0.

I suggest a new ImplicitDerivative, which is a subclass of Derivative and takes an "implicit" keywork argument. This specifies implicit dependencies such as above in a dictionairy, where the key is any sympy object and the value is a list of sympy objects on which it depends.

Please give any suggestions

Todo: using .doit() twice somehow fails and I am not sure how to fix it Some way to specify that a sympy object depends on any object it is derived to (e.g. by giving None instead of a list in the implicit dict and correct parsing) cleanup of code

from sympy import Derivative,Function,symbols,preorder_traversal
class ImplicitDerivative(Derivative):

    def __init__(self,*args,**kwargs):

        self.dummy_counter={}

        impdep=kwargs.get('implicits',{})
        self.impdep=impdep
        d_name='x{0}'
        i_name='f{0}'
        implicits=list(impdep.keys())
        implicit_dummies={ implicit:self.get_dummy(i_name,Function) for implicit in implicits  }

        dependencies=list(impdep.values())
        s=[] # I tried with a list comprehension with two for loops but it got into some kind of infinite loop, don't know why
        for dep in dependencies:               
            s.extend(dep)

        dependencies=s
        dependency_dummies={ dependency:self.get_dummy(d_name,symbols) for dependency in dependencies  }

        #for future: check if some implicits are also dependencies
        implicit_dummies={ imp:implicit_dummies[imp](*list(map(lambda d:dependency_dummies[d] ,dep))) for imp,dep in impdep.items()   }

        self.originals=implicits+dependencies
        self.dummies=[implicit_dummies[implicit] for implicit in implicits] + [dependency_dummies[dependency] for dependency in dependencies]

    def doit(self,**hints):
        variable_count=[ (expr.subs(zip(self.originals,self.dummies)),i) for expr,i in self.variable_count]

        subbed=Derivative(self.expr.subs(zip(self.originals,self.dummies)), *variable_count, evaluate=False)
        subbed_doit=subbed.doit()
        doit=subbed_doit.subs(zip(self.dummies,self.originals))
        #finally sub Derivatives for ImplicitDerivatives:
        deriv_subs = [ (arg,ImplicitDerivative(arg.expr,*arg.variable_count,implicits=self.impdep)) for arg in preorder_traversal(doit) if isinstance(arg,sym.Derivative)]
        final=doit.subs(deriv_subs)
        return final

    def get_dummy(self,name,call):
        di=self.dummy_counter.get(name,0)
        while self.expr.has(call(name.format(di))):
            di+=1
        self.dummy_counter[name]=di+1

        return call(name.format(di))
    def _eval_derivative(self,*args,**kwargs):
        print('hi')
        return super()._eval_derivatives(self,*args,**kwargs)
import sympy as sym
i=sym.Idx('i')
j=sym.Idx('j')
y=sym.IndexedBase('y')
x=sym.IndexedBase('x')
ddg=ImplicitDerivative(sym.tanh(y[i]),x[i],x[j],implicits={y[i]:[x[i],x[j]]})
ddg
ddg.doit()
ddg.doit().doit() #fails and becomes 0
oscarbenjamin commented 4 years ago

Rather than a different kind of derivative I think what we need is a way to generalise Function so that we can compose it with other kinds of symbols such as Indexed, MatrixSymbol etc. That is I want to declare a symbol that is both Indexed or IndexedBase and a Function at the same time. Or I want to declare a symbol that is both a MatrixSymbol and a Function representing a matrix that is a function of something else.

ThePauliPrinciple commented 4 years ago

Your suggestion would be nice way to solve it and this is what I tried first but it seems to me that it requires a complete redesign of both indexed and function, but would love to see a suggestion how to practically implement this.

W.r.t. other kind of derivative: the functionality in this class could also be added to the base derivative class, it only adds functionality through a keyword.

sylee957 commented 4 years ago

Although the classes like https://reference.wolfram.com/language/ref/Dt.html would be useful, the biggest concern I have is that we can end up incompatible classes of derivatives.

g-voelker commented 3 years ago

I am surprised that this is not implemented. Basically one can already generate an indexed function like this:

xx, nn = sympy.symbols("x n")
YY = sympy.IndexedBase(sympy.Function("y")(xx))

But then the derivative fails as a result calculated as follows is equal to zero

result = sympy.diff(YY[nn], xx)

The following calculates the derivative correctly but is no longer indexed.

result = sympy.diff(YY, xx)

That is reproducible and imho a bug in the current release.

Tested on Arch with python 3.8.5 and sympy 1.6.2

oscarbenjamin commented 3 years ago

one can already generate an indexed function

That's not really an indexed function. What you have is this:

In [31]: f = Function('f')                                                                                                        

In [32]: x = Symbol('x')                                                                                                          

In [33]: n = Symbol('n')                                                                                                          

In [34]: F = IndexedBase(f(x))                                                                                                    

In [35]: F                                                                                                                        
Out[35]: f(x)

In [36]: F[n]                                                                                                                     
Out[36]: f(x)[n]

Here F is an indexed with f(x) as its "symbol". The symbol is treated as an opaque parameter so the fact that it is a function of x does not mean that F itself is considered to be a function of x. The indexed will only recognise something that matches its symbol exactly:

In [46]: F[n].diff(x)                                                                                                             
Out[46]: 0

In [47]: F[n].diff(F[n])                                                                                                          
Out[47]: 1
g-voelker commented 3 years ago

Oh I see. Thanks for the explanation. Quite misleading then. I would suggest that before the new feature is introduced (if it will) a corresponding warning should be printed when attempting to index a function.