################################################################################################################################
##
## This code can be used in the open source CAS: SAGE (www.sagemath.org)
## to represent the McKean optimal stopping problem and the McKean 
## stochastic game, driven by a jump-diffusion with downwards directed
## exponential jumps. The code is based on the results derived in the 
## paper "Further calculations for the McKean stochastic game for a spectrally negative Levy process: from a point to an interval"
## written by Erik J. Baurdoux and Kees van Schaik.
##
## This is version 0.1 of the code, oct 2010.
##
## EXAMPLE:
## K=5; q=5/100
## JD=JumpDiff(sigma=3/10,llambda=2,theta=2,mu=1/10)  #create a jump-diffusion object
## JD.makeRiskNeutral(q)  #optional: redefine mu=q-sigma^2/2+llambda/(theta+1) so q=psi(1), i.e. risk neutral setting
## McKG=McKeanGame(JD,K,q)  #create a McKean stochastic game object
## delta=1/2
## McKG.getValue(delta,arg=None,xmin=-10,xmax=10)  #get the value function of the McKean game for delta=1/2, as a Piecewise object
## McKG.plot(delta,xmin=-10,xmax=10)  #get a plot of the McKean game for delta=1/2 (you don't need to call .getValue() first)
##
## TO DO:
## - fix Exponomial.simplify()
## - make case sigma=0 work properly
##
## Copyright 2010 Kees van Schaik (keesvanschaik@mail.com).
##
## This program is free software: you can redistribute it and/or modify
##    it under the terms of the GNU General Public License as published by
##    the Free Software Foundation, either version 2 of the License, or
##    (at your option) any later version.
##
##    This program is distributed in the hope that it will be useful,
##    but WITHOUT ANY WARRANTY; without even the implied warranty of
##    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
##    GNU General Public License for more details.
##
##    A full copy of the GNU General Public License can be found at <http://www.gnu.org/licenses/>.
################################################################################################################################

def findIdenticalInList(lst):
    """Returns a list of lists, each such list containing indices of
    elements in lst that are equal."""
    res=[]; indicesToSkip=[]
    for i in range(len(lst)):
        res_loc=[i]
        if not i in indicesToSkip:
            for j in range(i+1,len(lst)):
                if not j in indicesToSkip:
                    if lst[i]==lst[j]:
                        res_loc.append(j); indicesToSkip.append(j)
        if len(res_loc)>1:
            res.append(res_loc)
    return res

def flatten(lst):
    ###Only one level deep
    res=[]
    for el in lst:
        res+=el
    return res

def parseListsForPiecewise(intValList,functionsList,xmin,xmax):
    """Given a list of tuples (a,b), with assumed tpl[i][0]==tpl[i-1][1],
    and a list of functions, assumed to be of same length, returns a list
    of (tuple,function) pairs suitable to construct a Piecewise object
    with domain (xmin,xmax)."""
    if xmin is None or xmin<intValList[0][0]:
        xmin=intValList[0][0]
    if xmax is None or xmax>intValList[-1][1]:
        xmax=intValList[-1][1]
    piecewiseLst=[]; lowIndex=None; upIndex=None
    for i in range(len(intValList)):
        if not upIndex is None:
            break
        if xmin<intValList[i][1] and lowIndex is None:
            lowIndex=i
        if xmax<=intValList[i][1] and upIndex is None:
            upIndex=i
    #print 'lowIndex:',lowIndex
    #print 'upIndex:',upIndex
    if lowIndex==upIndex:
        piecewiseLst.append([(xmin,xmax),functionsList[lowIndex]])
    else:
        piecewiseLst.append([(xmin,intValList[lowIndex][1]),functionsList[lowIndex]])
        for i in range(lowIndex+1,upIndex):
            piecewiseLst.append([intValList[i],functionsList[i]])
        piecewiseLst.append([(intValList[upIndex][0],xmax),functionsList[upIndex]])
    return piecewiseLst
    

class Exponomial:
    def __init__(self,factors=[],exponents=[]):
        """Represents an expression of the form factors[0]*exp(exponents[0]*x)+...+factors[-1]*exp(exponents[-1]*x),
        where x is a free variable. Used here to handle certain operations on such expressions that Maxima 
        sometimes seems to fail to do succesfully."""
        if not len(factors)==len(exponents):
            raise Exception("In Exponomial.__init__(): provided factors and exponents are not of equal length")
        self.factors=factors; self.exponents=exponents
        
    def help_simplify(self):
        ###simplify() seems to throw occasional (maxima) error
        try:
            self.factors=[simplify(el) for el in self.factors]
        except:
            pass
        try:
            self.exponents=[simplify(el) for el in self.exponents]
        except:
            pass

    def n(self,digits=None):
        if digits is None:
            factsNew=[n(el) for el in self.factors]; exposNew=[n(el) for el in self.exponents]
        else:
            factsNew=[n(el,digits=digits) for el in self.factors]; exposNew=[n(el,digits=digits) for el in self.exponents]
        res=Exponomial(factors=factsNew,exponents=exposNew)
        #res.simplify()
        return res
        #for i in range(len(self.exponents)):
        #    if digits is None:
        #        res+=n(self.factors[i])*exp(n(self.exponents[i])*x)
        #    else:
        #        res+=n(self.factors[i],digits=digits)*exp(n(self.exponents[i],digits=digits)*x)
        #return res

    def stripExponentsWithCriterium(self,criteriumFunc):
        toRemove=[]
        for i in range(len(self.exponents)):
            if criteriumFunc(self.exponents[i]):
                toRemove.append(i)
        if len(toRemove)>0:
            self.exponents=[self.exponents[i] for i in range(len(self.exponents)) if not i in toRemove]
            self.factors=[self.factors[i] for i in range(len(self.factors)) if not i in toRemove]
        
    def simplify(self):
        return ###For now, there seems to be some problem...
        print 'Simplfy called on ',self
        print '  .n():',self.n()
        print '  simplify 1: len=',self.__len__()
        self.help_simplify()
        print '  simplify 2: len=',self.__len__()
        ###Add terms with identical exponents:
        identLst=findIdenticalInList(self.exponents)
        factorsNew=[]; exponentsNew=[]
        for el in identLst:
            exponentsNew.append(self.exponents[el[0]])
            factor=0
            for el2 in el:
                factor+=self.factors[el2]
            factorsNew.append(factor)
        flatIdentLst=flatten(identLst)
        for i in range(len(self.exponents)):
            if not i in flatIdentLst:
                factorsNew.append(self.factors[i]); exponentsNew.append(self.exponents[i])
        self.factors=factorsNew; self.exponents=exponentsNew
        print '  simplify 3: len=',self.__len__()
        self.help_simplify()
        print '  simplify 4: len=',self.__len__()
        ###Remove terms with factors that are 0:
        toRemove=[]
        for i in range(len(self.factors)):
            if self.factors[i]==0:
                toRemove.append(i)
        cnt=0
        for i in range(len(toRemove)):
            self.factors.pop(i-cnt); self.exponents.pop(i-cnt); cnt=cnt-1
        print '  simplify 5: len=',self.__len__()
            
    def linTransform(self,a,b):
        """In function form: if self is c1*e^(c2*x) then after this transformation we
        have that self equals c1*e^(c2*(a*x+b))=(c1*e^(c2*b))*e^(c2*a*x)."""
        for i in range(len(self.factors)):
            self.factors[i]=self.factors[i]*exp(self.exponents[i]*b)
            self.exponents[i]=self.exponents[i]*a
        self.simplify()
        
    def asExpression(self):
        tmp=0
        for i in range(len(self.factors)):
            tmp+=self.factors[i]*exp(self.exponents[i]*x)
        return tmp
    
    def asFunction(self):
        return self.asExpression().function(x)
    
    def __call__(self,arg):
        return self.asExpression().subs(x=arg)
    
    def __len__(self):
        return len(self.exponents)
    
    def __repr__(self):
        return "Exponomial: "+str(self.asFunction())
    
    def __add__(self,other):
        new=Exponomial(factors=self.factors+other.factors,exponents=self.exponents+other.exponents)
        new.simplify()
        return new

    def __mul__(self,other):
        if isinstance(other,Exponomial):
            factorsNew=[]; exposNew=[]
            for i in range(len(self.factors)):
                for j in range(len(other.factors)):
                    factorsNew.append(self.factors[i]*other.factors[j])
                    exposNew.append(self.exponents[i]+other.exponents[j])
            new=Exponomial(factors=factorsNew,exponents=exposNew)
        else:
            new=Exponomial(factors=[other*el for el in self.factors],exponents=self.exponents)
        new.simplify()
        return new     
    
    def __rmul__(self,other):
        return self.__mul__(other)
    
    def integrate(self,a=None,b=None):
        antiDer=0
        for i in range(len(self.exponents)):
            if self.exponents[i]==0:
                antiDer+=self.factors[i]*x
            else:
                antiDer+=self.factors[i]*exp(self.exponents[i]*x)/self.exponents[i]
        antiDer=antiDer.function(x)
        if a is None and b is None:
            return antiDer
        else:
            return antiDer(b)-antiDer(a)
        
    def runningIntegral(self,a):
        """Returns \int_a^x self(z) dz as an Exponomial."""
        Cst=0; factsNew=[]; exposNew=[]
        for i in range(len(self.exponents)):
            if self.exponents[i]==0:
                raise Exception("In Exponomial.runningIntegral(): one of the exponents is 0, hence the running integral can not be expressed as Exponomial")
            else:
                exposNew.append(self.exponents[i])
                factsNew.append(self.factors[i]/self.exponents[i])
                Cst+=-self.factors[i]*exp(self.exponents[i]*a)/self.exponents[i]
        return Exponomial(factors=factsNew+[Cst,],exponents=exposNew+[0,])
    
    def derivative(self):
        return Exponomial(factors=map(lambda x,y: x*y,self.factors,self.exponents),exponents=self.exponents)

class JumpDiff:
    """Represents a jump-diffusion with downwards directed exponential jumps, i.e. a process X with
    X_t = sigma*W_t + mu*t - \sum_{i=1}^N_t Y_i,
    with W a Brownian motion, mu in R, N a Poisson process with rate llambda, (Y_i) a sequence of
    iid rv's with parameter theta."""
    def __init__(self,sigma,llambda,theta,mu=None):
        self.sigma=sigma; self.mu=mu; self.llambda=llambda; self.theta=theta; self.riskNeutral=False
        if not (sigma>=0 and llambda>0 and theta>=0):
            raise Exception(r"In JumpDiff.__init__(): please provide parameters satisfying \sigma>=0, \lambda>0, \theta>0")
        
    def makeRiskNeutral(self,q):
        self.riskNeutral=True
        self.mu=q-self.sigma^2/2+self.llambda/(self.theta+1)
        
    def isSubordinator(self):
        if self.sigma==0 and self.mu<=0:
            return True
        return False
        
    def isRiskNeutralFor(self,q):
        func=self.getLaplacian()
        if q==func(1):
            return True
        return False
    
    def isConsistentWith(self,q):
        func=self.getLaplacian()
        if q>0 and func(1)>=0 and func(1)<=q:
            return (True,)
        return (False,func(1))
        
    def getLaplacian(self):
        return (self.sigma^2*x^2/2+self.mu*x-self.llambda+self.llambda*self.theta/(self.theta+x)).function(x)
    
    def getInvLaplacianAt(self,arg):
        func=self.getLaplacian()-arg
        lb=0; ub=1
        while True:
            if func(ub)>0:
                break
            elif func(ub)==0:
                return ub
            else:  ###func(ub)<0
                lb=ub; ub+=1
        return find_root(func,lb,ub)
    
class McKeanGame:
    """An object that represents a stochastic game with lower payoff process L_t = (K-exp(X_t))^+ and upper payoff
    process U_t = (K-exp(X_t))^+ + delta, where X is a jump-diffusion represented by a JumpDiff object. The discount rate
    is q. A saddle point is (\sigma^*,\tau^*), where \sigma^*=\inf \{ t >= 0 | X_t <= x^* \} and \tau^* is either
    infinite a.s. (when self.getMode() returns 'optStopp'); is \inf \{ t >= 0 | X_t = \log(K) \} (when self.getMode() 
    returns 'game_point'); is \inf \{ t >= 0 | X_t \in [\log(K),y^*] \} for some y^*>\log(K) (when self.getMode() returns
    'game_interval'). self.getValue() returns the value function V (as a function of the starting point of X) of
    this game."""
    def __init__(self,JumpDiffObj,K,q):
        self.JumpDiff=JumpDiffObj; self.K=K; self.q=q
        assert (not self.JumpDiff.isSubordinator()), r"In McKeanGame.__init__(): the provided JumpDiff object is a subordinator"
        testCond=self.JumpDiff.isConsistentWith(q)
        assert testCond[0], r"In McKeanGame.__init__(): the condition $q>0$ and $0 \leq \psi(1) \leq q$ fails to hold, with n(\psi(1))="+str(n(testCond[1]))
        self.McKeanOptStopp=McKeanOptStopp(JumpDiffObj,K,q)
        self.deltaBar=self.McKeanOptStopp.getValue(arg=log(K))
        self.xStar={}; self.yStar={}; self.deltaZero=None; self.valFunc={}
        self.maxInternalError=10^(-10)
        
    def getxStar(self,delta,doCache=True):
        if delta>=self.deltaBar:
            return self.McKeanOptStopp.kStar
        if doCache:
            if str(delta) in self.xStar.keys():
                return self.xStar[str(delta)]            
        ###else:  Using equation (7) from paper:
        func=self.JumpDiff.getLaplacian()
        Z0=ScaleFuncZ(self.JumpDiff,0,self.q).toFunction()
        Z1=ScaleFuncZ(self.JumpDiff,1,self.q-func(1)).toFunction()       
        rootFunc=(Z0(log(self.K)-x)-Z1(log(self.K)-x)-delta/self.K).function(x)
        rb=log(self.K); lb=log(self.K)-1
        sgnR=sgn(rootFunc(rb)); lb=log(self.K)-1
        while sgn(rootFunc(lb))==sgnR:
            rb=lb; lb=lb-1
        rt=find_root(rootFunc,lb,rb)
        if doCache:
            self.xStar[str(delta)]=rt
        return rt
    
    def getyStar(self,delta,doCache=True):
        if delta>self.deltaBar:
            return None
        if doCache:
            if str(delta) in self.yStar.keys():
                return self.yStar[str(delta)]
        ###We don't need to compute deltaZero first (takes much time), just compute
        ###y^* according to the equation below, and return max(log(K),y^*).
        #if self.deltaZero is None:
        #    self.deltaZero=self.getDeltaZero()
        #if delta>=self.deltaZero:
        #    return log(self.K)
        K=self.K; q=self.q; theta=self.JumpDiff.theta; llambda=self.JumpDiff.llambda
        xStar=self.getxStar(delta)
        func=self.JumpDiff.getLaplacian()
        Z0=ScaleFuncZ(self.JumpDiff,0,self.q).toExponomial()
        Z1=ScaleFuncZ(self.JumpDiff,1,self.q-func(1)).toExponomial()
        int1=(Z0*Exponomial(factors=[1,],exponents=[theta,])).integrate(a=0,b=log(K)-xStar)
        int2=(Z1*Exponomial(factors=[1,],exponents=[theta+1,])).integrate(a=0,b=log(K)-xStar)
        int=exp(theta*xStar)*K*int1-exp(xStar*(theta+1))*int2
        lhs=int+K*exp(theta*xStar)/theta-exp((theta+1)*xStar)/(theta+1)-delta*K^theta/theta 
        invLapl=self.JumpDiff.getInvLaplacianAt(q)
        rhs=delta*q*(theta+invLapl)/(llambda*theta*invLapl)
        yStar=-log(rhs/lhs)/theta
        if yStar<=log(K):
            yStar=log(K)
        #assert yStar>log(K), "In McKeanGame.getyStar(): y^*="+str(n(yStar))+", log(K)="+str(n(log(K)))+", hence assertion y*>log(K) failed"
        if doCache:
            self.yStar[str(delta)]=yStar
        return yStar
        
    def getMode(self,delta):
        if delta>self.deltaBar:
            return 'optStopp'
        if self.getyStar(delta)==log(self.K):  ###avoid having to compute deltaZero
            return 'game_point'
        return 'game_interval'
        #if self.deltaZero is None:
        #    self.deltaZero=self.getDeltaZero()
        #if delta>=self.deltaZero:
        #    return 'game_point'
        #return 'game_interval'       
    
    def _help_deltaZero_equation(self,delta):
        """This returns the lhs of the formula on page P7 of the notes, which is 
        a simplified version of the formula in Thm 8 (i)."""
        theta=self.JumpDiff.theta; llambda=self.JumpDiff.llambda; K=self.K; q=self.q
        invLapl=self.JumpDiff.getInvLaplacianAt(q)
        xStar=self.getxStar(delta,doCache=True)
        func=self.JumpDiff.getLaplacian()
        Z0=ScaleFuncZ(self.JumpDiff,0,self.q).toExponomial()
        Z1=ScaleFuncZ(self.JumpDiff,1,self.q-func(1)).toExponomial()
        int1=(Z0*Exponomial(factors=[1,],exponents=[theta,])).integrate(a=0,b=log(K)-xStar)
        int2=(Z1*Exponomial(factors=[1,],exponents=[theta+1,])).integrate(a=0,b=log(K)-xStar)
        int=exp(theta*xStar)*K*int1-exp(xStar*(theta+1))*int2 
        tmp=q*(invLapl+theta)/(llambda*theta*invLapl)+1/theta
        return int+K*exp(theta*xStar)/theta-exp((theta+1)*xStar)/(theta+1)-delta*K^theta*tmp 
    
    def getDeltaZero(self,showProgress=False):
        """Finds the solution to lhs=0, where lhs is as returned by self._help_deltaZero_equation()."""
        if not self.deltaZero is None:
            return self.deltaZero
        if showProgress:
            print "Entering self.getDeltaZero()..."
        func=self._help_deltaZero_equation
        if showProgress:
            print " ...finding lb and ub"; t=cputime()
        lb=self.deltaBar/2; ub=self.deltaBar/2
        while sgn(func(lb))==sgn(func(ub)):
            lb=lb/2; ub=(self.deltaBar+ub)/2
        if showProgress:
            print " ...took "+str(cputime()-t)+" secs"
            print " ...actual root finding"; t=cputime()
        rt=find_root(func,lb,ub)
        if showProgress:
            print " ...took "+str(cputime()-t)+" secs"
        self.deltaZero=rt
        return rt
    
    def _getValueFuncLeft_asExponomial(self,xStar,delta):
        func=self.JumpDiff.getLaplacian()
        Z0=ScaleFuncZ(self.JumpDiff,0,self.q)
        Z1=ScaleFuncZ(self.JumpDiff,1,self.q-func(1))
        Z0.shiftBy(-xStar); Z1.shiftBy(-xStar)
        res=self.K*Z0.toExponomial()+Exponomial(factors=[-1,],exponents=[1,])*Z1.toExponomial()
        assert abs(res(xStar)-(self.K-exp(xStar)))<self.maxInternalError, r"In McKeanGame._getValueFuncLeft_asExponomial: pasting problem: |V(x^*+)-V(x^*-)| is larger than internal tolerance bound "+str(self.maxInternalError)
        assert abs(res(log(self.K))-delta)<self.maxInternalError, r"In McKeanGame._getValueFuncLeft_asExponomial: pasting problem: |V(log(K))-\delta| is larger than internal tolerance bound "+str(self.maxInternalError)
        if self.JumpDiff.sigma>0:
            resDer=res.derivative()
            assert abs(resDer(xStar)-(-exp(xStar)))<self.maxInternalError, r"In McKeanGame._getValueFuncLeft_asExponomial: pasting problem: |V'(x^*+)-V'(x^*-)| is larger than internal tolerance bound "+str(self.maxInternalError)
        return res
        #return (self.K*Z0.toFunction()-exp(x)*Z1.toFunction()).function(x)        
    
    def _getValueFuncLeft(self,xStar,delta):
        return self._getValueFuncLeft_asExponomial(xStar,delta).asFunction()
        #func=self.JumpDiff.getLaplacian()
        #Z0=ScaleFuncZ(self.JumpDiff,0,self.q)
        #Z1=ScaleFuncZ(self.JumpDiff,1,self.q-func(1))
        #Z0.shiftBy(-xStar); Z1.shiftBy(-xStar)
        #return (self.K*Z0.toFunction()-exp(x)*Z1.toFunction()).function(x)
        
    def _getValueFuncRight_interval_asExponomial(self,xStar,yStar,delta):
        func=self.JumpDiff.getLaplacian(); funcDer=func.derivative; invLapl=self.JumpDiff.getInvLaplacianAt(self.q)
        theta=self.JumpDiff.theta; llambda=self.JumpDiff.llambda
        ###Compute constant ( Cst = \lambda*\theta*e^{\theta*y^*}*\int_{t<0} (w_{\delta}(t+y^*)-\delta)*e^{\theta t} dt ):
        expo=self._getValueFuncLeft_asExponomial(xStar,delta)*Exponomial(factors=[1,],exponents=[theta])
        Cst=expo.integrate(a=xStar,b=log(self.K))
        Cst+=self.K*exp(theta*xStar)/theta-exp((theta+1)*xStar)/(theta+1)-delta*self.K^theta/theta
        Cst=Cst*exp(-theta*yStar)
        Cst=Cst*llambda*theta*exp(theta*yStar)
        ###Compute F  ( F(x) = \int_{z=0}^{x-y^*} W^{(q)}(z)*e^{\theta*z} dz ):
        expo=ScaleFuncW(self.JumpDiff,0,self.q).toExponomial()*Exponomial(factors=[1,],exponents=[theta])
        expo=expo.runningIntegral(0)
        #print 'expo:',expo.n()
        #print 'len(expo):',len(expo)
        #print 'n(expo(0)):',n(expo(0))
        expo.linTransform(1,-yStar)
        F=expo
        #print 'F:',F.n()
        #print 'len(F):',len(F)
        #print 'F(yStar):',F(yStar)
        #print 'n(F(yStar)):',n(F(yStar))
        ###End result  ( = x \mapsto \delta*Z^{(q)}(x-y^*) - Cst*e^{-\theta*x}*F(x) ):
        res=ScaleFuncZ(self.JumpDiff,0,self.q).toExponomial()
        res.linTransform(1,-yStar)
        res=delta*res+Exponomial(factors=[-Cst,],exponents=[-theta,])*F
        res.stripExponentsWithCriterium(lambda x: x>0)  ###In theory no positive exponents should be present, this is to ensure any such terms present due to numerical noise get removed
        assert abs(res(yStar)-delta)<self.maxInternalError, r"In McKeanGame._getValueFuncRight_interval_asExponomial: pasting problem: |V(y^*)-\delta| is larger than internal tolerance bound "+str(self.maxInternalError)
        if self.JumpDiff.sigma>0:
            resDer=res.derivative()
            assert abs(resDer(yStar))<self.maxInternalError, r"In McKeanGame._getValueFuncRight_interval_asExponomial: pasting problem: |V'(y^*)| is larger than internal tolerance bound "+str(self.maxInternalError)
        return res
    
    def _getValueFuncRight_point(self,xStar,delta):
        func=self.JumpDiff.getLaplacian(); funcDer=func.derivative; invLapl=self.JumpDiff.getInvLaplacianAt(self.q)
        W0=ScaleFuncW(self.JumpDiff,0,self.q)
        W0.shiftBy(-log(self.K))
        if self.JumpDiff.isRiskNeutralFor(self.q):
            alpha=exp(xStar)*funcDer(1)-self.K*func(1)
        else:
            alpha=exp(xStar)*(self.q-func(1))/(invLapl-1)-self.q*self.K/invLapl
        res=(alpha*exp(invLapl*(log(self.K)-xStar))*W0.toFunction()).function(x)
        assert abs(res(log(self.K))-delta)<self.maxInternalError, r"In McKeanGame._getValueFuncRight_point: pasting problem: |V(log(K))-\delta| is larger than internal tolerance bound "+str(self.maxInternalError)
        return res
    
    def getValue(self,delta,arg=None,xmin=None,xmax=None):
        """Returns the value of the game. If arg is not None, the value function evaluated at argument arg
        is returned. If arg is None, the value function is returned as a Piecewise object with domain (xmin,xmax)."""
        mode=self.getMode(delta) ###is either 'optStopp', 'game_point' or 'game_interval'
        if mode=='optStopp':
            if not arg is None:
                return self.McKeanOptStopp.getValue(arg=arg)
            return self.McKeanOptStopp.getValue(xmin=xmin,xmax=xmax)
        ###First part: produce lists intVals and functions covering whole domain:
        xStar=self.getxStar(delta)
        f1=(self.K-exp(x)).function(x)
        f2=self._getValueFuncLeft(xStar,delta)
        intVals=[(-infinity,xStar),(xStar,log(self.K))]
        funcLst=[f1,f2]
        if mode=='game_point':
            intVals.append((log(self.K),infinity))
            funcLst.append(f2+self._getValueFuncRight_point(xStar,delta))
        if mode=='game_interval':
            yStar=self.getyStar(delta)
            intVals.append((log(self.K),yStar))
            funcLst.append((delta+0*x).function(x))
            intVals.append((yStar,infinity))
            f4=self._getValueFuncRight_interval_asExponomial(xStar,yStar,delta).asFunction()
            funcLst.append(f4)
            #_getValueFuncRight_interval_asExponomial(self,xStar,yStar,delta)
        ###Second part: return either value or piecewise function:
        if not arg is None:
            for i in range(len(intVals)):
                if arg<intVals[i][1]: ###arbitrary choice between < and <=
                    return funcLst[i](arg)
        else:
            return Piecewise(parseListsForPiecewise(intVals,funcLst,xmin,xmax))
        
    def plot(self,delta,xmin,xmax,withLowerPayoff=True,withUpperPayoff=True):
        plots=[]
        if withLowerPayoff:
            intVals=[(-infinity,log(self.K)),(log(self.K),infinity)]
            funcLst=[(self.K-exp(x)).function(x),(0*x).function(x)]
            g=Piecewise(parseListsForPiecewise(intVals,funcLst,xmin,xmax))
            plots.append(g.plot(rgbcolor='black'))
        if withUpperPayoff:
            intVals=[(-infinity,log(self.K)),(log(self.K),infinity)]
            funcLst=[(self.K-exp(x)+delta).function(x),(0*x+delta).function(x)]
            g=Piecewise(parseListsForPiecewise(intVals,funcLst,xmin,xmax))
            plots.append(g.plot(rgbcolor='black'))
        f=self.getValue(delta,None,xmin,xmax)
        plots.append(f.plot(rgbcolor='red'))
        return sum(plots)

    
class McKeanOptStopp:
    """An object that represents an optimal stopping problem with payoff process L_t = (K-exp(X_t))^+, where X is a 
    jump-diffusion represented by a JumpDiff object. The discount rate is q. The optimal stopping time is given by
    \inf \{ t >= 0 | X_t <= k^* \}. self.getValue() returns the value function V (as a function of the starting point of X) of
    this optimal stopping problem."""
    def __init__(self,JumpDiffObj,K,q):
        self.JumpDiff=JumpDiffObj; self.K=K; self.q=q
        assert (not self.JumpDiff.isSubordinator()), r"In McKeanOptStopp.__init__(): the provided JumpDiff object is a subordinator"
        testCond=self.JumpDiff.isConsistentWith(q)
        self.maxInternalError=10^(-10)
        assert testCond[0], r"In McKeanOptStopp.__init__(): the condition $q>0$ and $0 \leq \psi(1) \leq q$ fails to hold, with n(\psi(1))="+str(n(testCond[1]))
        if JumpDiffObj.isRiskNeutralFor(q):
            func=JumpDiffObj.getLaplacian()
            funcPr=func.derivative()
            self.kStar=log(K*func(1)/funcPr(1))
        else:
            func=JumpDiffObj.getLaplacian()
            tmp=(JumpDiffObj.getInvLaplacianAt(q)-1)/(q-func(1))
            self.kStar=log(K*q*tmp/JumpDiffObj.getInvLaplacianAt(q))
        assert self.kStar<log(K), r"In McKeanOptStopp.__init__(): assertion k*<log(K) failed"

    def _getValueRight_asExponomial(self):
        func=self.JumpDiff.getLaplacian()
        Z0=ScaleFuncZ(self.JumpDiff,0,self.q)
        Z0.shiftBy(-self.kStar)
        #print 'Z0(k^*)=',n(Z0(self.kStar))
        Z1=ScaleFuncZ(self.JumpDiff,1,self.q-func(1))
        Z1.shiftBy(-self.kStar)
        #print 'Z1(k^*)=',n(Z1(self.kStar))
        res=self.K*Z0.toExponomial()+Exponomial(factors=[-1,],exponents=[1,])*Z1.toExponomial()
        res.stripExponentsWithCriterium(lambda x: x>0) ###remove any terms with positive exponents as they cancel out in theory, possibly present here only due to numerical noise
        res.simplify()
        print res.n()
        #print 'res(k^*)=',n(res(self.kStar))
        assert abs(res(self.kStar)-(self.K-exp(self.kStar)))<self.maxInternalError, r"In McKeanOptStopp._getValueRight_asExponomial: pasting problem: |V(k^*+)-V(k^*-)| is larger than internal tolerance bound "+str(self.maxInternalError)
        if self.JumpDiff.sigma>0:
            resDer=res.derivative()
            assert abs(resDer(self.kStar)-(-exp(self.kStar)))<self.maxInternalError, r"In McKeanOptStopp._getValueRight_asExponomial: pasting problem: |V'(k^*+)-V'(k^*-)| is larger than internal tolerance bound "+str(self.maxInternalError)
        return res
        
    def getValue(self,arg=None,xmin=None,xmax=None):
        """Returns the value of the optimal stopping problem. If arg is not None, the value function 
        evaluated at argument arg is returned. If arg is None, the value function is returned as a 
        Piecewise object with domain (xmin,xmax)."""
        if not arg is None:
            if arg<=self.kStar:
                return self.K-exp(arg)
            func=self.JumpDiff.getLaplacian()
            Z0=ScaleFuncZ(self.JumpDiff,0,self.q)
            Z1=ScaleFuncZ(self.JumpDiff,1,self.q-func(1))
            return self.K*Z0(arg-self.kStar)-exp(arg)*Z1(arg-self.kStar)
        f1=(self.K-exp(x)).function(x)
        f2=self._getValueRight_asExponomial().asFunction()
        intVals=[(-infinity,self.kStar),(self.kStar,infinity)]
        funcLst=[f1,f2]
        return Piecewise(parseListsForPiecewise(intVals,funcLst,xmin,xmax))

    def plot(self,xmin,xmax,withPayoff=True):
        plots=[]
        if withPayoff:
            intVals=[(-infinity,log(self.K)),(log(self.K),infinity)]
            funcLst=[(self.K-exp(x)).function(x),(0*x).function(x)]
            g=Piecewise(parseListsForPiecewise(intVals,funcLst,xmin,xmax))
            plots.append(g.plot(rgbcolor='black'))
        f=self.getValue(None,xmin,xmax)
        plots.append(f.plot(rgbcolor='red'))
        return sum(plots)
    
class ScaleFunc:
    """A super class representing a scale function (W or Z) for a jump-diffusion X as 
    represented by a JumpDiff object.""" 
    def __init__(self,JumpDiffObj,c,r):
        if JumpDiffObj.sigma==0:
            raise Exception("In ScaleFunc.__init__(): the case sigma=0 is not yet (correctly) implemented")
        if c<0:
            raise Exception("In ScaleFunc.__init__(): parameter c can not be negative")
        if r<0:
            raise Exception("In ScaleFunc.__init__(): parameter r can not be negative")
        if JumpDiffObj.isSubordinator():
            raise Exception("In ScaleFunc.__init__(): the process object JumpDiffObj seems to be a subordinator, hence scale function not defined")
        self.JumpDiffObj=JumpDiffObj; self.c=c; self.r=r
        #self.betas=[0,]+self.getBetas()
        #self.constants=[0,]+self.getConstants() 
        self.betas=None ###will be overriden by children
        self.constants=None ###will be overriden by children
        self.asExpression=None ###will be overriden by children
        self.valForNegativeArguments=None ###will be overriden by children
        self.shiftAmount=0 ###Shifts the function by this amount, i.e. x -> ScaleFunc(x+self.shiftBy)
        
    def shiftBy(self,shift):
        self.shiftAmount=shift
        self.asExpression=self.asExpression.subs(x=x+self.shiftAmount)
        
    def __call__(self,arg):
        if arg<-self.shiftAmount:
            return self.valForNegativeArguments
        return self.asExpression.subs(x=arg)
    
    def toPiecewiseFunction(self,xmin,xmax):
        if xmin>-self.shiftAmount or xmax<-self.shiftAmount:
            raise Exception("In ScaleFunc.toPiecewiseFunction(): constraint xmin<=-self.shiftAmount and xmax>=-self.shiftAmount not satisfied")
        return Piecewise([[(xmin,-self.shiftAmount),(self.valForNegativeArguments+0*x).function(x)],[(-self.shiftAmount,xmax),self.toFunction()]])
    
    def toFunction(self):
        """CAREFUL: the returned formula only holds for x>-self.shiftAmount!"""
        return self.asExpression.function(x)
        
    def getConstants(self):
        if self.betas is None:
            self.betas=[0,]+self.getBetas()
        res=[]; theta=self.JumpDiffObj.theta; sigma=self.JumpDiffObj.sigma; mu=self.JumpDiffObj.mu
        if sigma>0:
            for i in range(1,4):
                tmp=1
                for j in range(1,4):
                    if not j==i:
                        tmp=tmp*(self.betas[j]-self.betas[i])
                res.append(2*(theta+self.c+self.betas[i])/(sigma^2*tmp))
            return res
        else:
            return [1/((self.betas[1]-self.betas[2])*mu),1/((self.betas[2]-self.betas[1])*mu)]
        
    def getBetas(self):
        if self.JumpDiffObj.sigma>0:
            theta=self.JumpDiffObj.theta; sigma=self.JumpDiffObj.sigma; llambda=self.JumpDiffObj.llambda; mu=self.JumpDiffObj.mu
            res=[]
            if self.r>0:
                if self.JumpDiffObj.isRiskNeutralFor(self.r):
                    tmp=theta/2+self.r/sigma^2+llambda/(sigma^2*(theta+1))
                    tmp2=tmp^2-2*self.r*theta/sigma^2
                    res=[-tmp-sqrt(tmp2),-tmp+sqrt(tmp2),1]
                else:
                    laplac=self.JumpDiffObj.getLaplacian()
                    rootFunc=(laplac(x+self.c)-laplac(self.c)-self.r).function(x)
                    ###\beta_1:
                    lb=-theta-self.c-1; ub=-theta-self.c-1
                    while sgn(rootFunc(ub))>=0:
                        ub=(-theta-self.c+ub)/2
                    while sgn(rootFunc(lb))<=0:
                        lb=lb-1
                    res.append(find_root(rootFunc,lb,ub))
                    ###\beta_2:
                    lb=(-theta-self.c)/2; ub=0
                    while sgn(rootFunc(lb))<=0:
                        lb=(-theta-self.c+lb)/2
                    res.append(find_root(rootFunc,lb,ub))
                    ###\beta_3:
                    lb=0; ub=1
                    while sgn(rootFunc(ub))<=0:
                        ub=ub+1
                    res.append(find_root(rootFunc,lb,ub))   
                assert (res[0]<-theta-self.c and res[1]>-theta-self.c and res[1]<0 and res[2]>0), r"In ScaleFunc.getBetas(), case r>0: assertion \beta_1<-theta-c<\beta_2<0<\beta_3 failed"
                return res      
            ###else: cases where r=0:
            laplac=self.JumpDiffObj.getLaplacian()
            res.append(0); rootFunc=(laplac(x+self.c)-laplac(self.c)).function(x)
            aa=sigma^2/2; bb=sigma^2*theta/2+3*sigma^2*self.c/2+mu
            cc=(sigma^2*self.c+mu)*(theta+self.c)-llambda*theta/(theta+self.c)
            tmp=bb^2-4*aa*cc
            if tmp<=0:
                raise Exception("In ScaleFunc.getBetas(): quadratic equation for finding last two roots appears to have less than two (real) solutions")
            res.append((-bb-sqrt(tmp))/(2*aa)); res.append((-bb+sqrt(tmp))/(2*aa))
            res.sort()
            assert (res[0]<-theta-self.c and res[1]>-theta-self.c and res[2]>=res[1]), r"In ScaleFunc.getBetas(), case r=0: assertion \beta_1<-theta-c<\beta_2<=\beta_3 failed"
            assert (rootFunc(res[0])==0 and rootFunc(res[1])==0 and rootFunc(res[2])==0), r"In ScaleFunc.getBetas(), case r=0: assertion \psi_c(\beta_i)=0 for i=1,2,3 failed"
            return res
        else:  ###case self.JumpDiffObj.sigma==0:
            laplac=self.JumpDiffObj.getLaplacian(); rootFunc=(laplac(x+self.c)-laplac(self.c)-self.r).function(x)
            theta=self.JumpDiffObj.theta; llambda=self.JumpDiffObj.llambda; mu=self.JumpDiffObj.mu
            aa=mu; bb=mu*(self.c+theta)-llambda*theta/(theta+self.c)-self.r; cc=-self.r*(self.c+theta)
            tmp=bb^2-4*aa*cc
            #print 'rootFunc:',rootFunc
            #print 'aa:',n(aa)
            #print 'bb:',n(bb)
            #print 'cc:',n(cc)
            assert tmp>0, r"In ScaleFunc.getBetas(): quadratic equation for finding the two roots appears to have less than two (real) solutions"
            res=[((-bb-sqrt(tmp))/(2*aa)),((-bb+sqrt(tmp))/(2*aa))]
            assert (rootFunc(res[0])==0 and rootFunc(res[1])==0), r"In ScaleFunc.getBetas(): assertion \psi_c(\beta_i)=0 for i=1,2 failed"
            return res
            
class ScaleFuncW(ScaleFunc):
    """A class representing the scale function W for a jump-diffusion X as 
    represented by a JumpDiff object.""" 
    def __init__(self,JumpDiffObj,c,r):
        ScaleFunc.__init__(self,JumpDiffObj,c,r)
        #assert self.JumpDiffObj.sigma>0, "In ScaleFuncW.__init__(): not (yet) implemented for \sigma=0"
        self.betas=[0,]+self.getBetas()
        if self.JumpDiffObj.sigma==0 and self.betas[2]==0 and self.betas[3]==0:
            theta=self.JumpDiffObj.theta; sigma=self.JumpDiffObj.sigma
            expr=2*((1-self.c-theta)*exp(self.betas[1]*x)-(theta+self.c)*x+theta+self.c-1)/(sigma^2*self.betas[1])
        else:
            self.constants=[0,]+self.getConstants()
            expr=0
            for i in range(1,len(self.constants)):
                expr+=self.constants[i]*exp(self.betas[i]*x)
        self.asExpression=expr
        self.valForNegativeArguments=0
    
    def toExponomial(self):
        """CAREFUL: the returned formula only holds for x>-self.shiftAmount!"""
        if self.r==0:
            raise Exception("In ScaleFuncW.asExponomial(): for r=0 the scale function cannot be expressed as an exponomial")
        expo=Exponomial(factors=self.constants[1:],exponents=self.betas[1:])
        expo.linTransform(1,self.shiftAmount)
        return expo
    
class ScaleFuncZ(ScaleFunc):
    """A class representing the scale function Z for a jump-diffusion X as 
    represented by a JumpDiff object.""" 
    def __init__(self,JumpDiffObj,c,r):
        ScaleFunc.__init__(self,JumpDiffObj,c,r)
        #assert self.JumpDiffObj.sigma>0, "In ScaleFuncZ.__init__(): not (yet) implemented for \sigma=0"
        self.betas=[0,]+self.getBetas(); 
        if self.r>0:
            self.constants=[0,]+self.getConstants(); expr=0
            for i in range(1,len(self.constants)):
                expr+=self.constants[i]*exp(self.betas[i]*x)/self.betas[i]
            expr=self.r*expr
            if self.JumpDiffObj.sigma==0:
                expr=expr+1-self.r*sum([self.constants[i]/self.betas[i] for i in range(1,len(self.constants))])
        else:
            expr=1+0*x
        self.asExpression=expr    
        self.valForNegativeArguments=1
        
    def toExponomial(self):
        """CAREFUL: the returned formula only holds for x>-self.shiftAmount!"""
        if self.r>0:
            facts=[self.r*self.constants[i]/self.betas[i] for i in range(1,len(self.constants))]
            expos=self.betas[1:]
            if self.JumpDiffObj.sigma==0:
                facts.append(1-self.r*sum([self.constants[i]/self.betas[i] for i in range(1,len(self.constants))]))
                expos.append(0)
            expo=Exponomial(factors=facts,exponents=expos)
        else:
            expo=Exponomial(factors=[1,],exponents=[0,])
        expo.linTransform(1,self.shiftAmount)
        return expo