Adding model with function only interface

pythonScripts/cModel.py

@@ -2,7 +2,8 @@ import numpy
 import json
 import os
 import scipy.interpolate
-
+#for partial function specializations
+import functools
 
 class model:
    def __init__(self):
@@ -24,113 +25,138 @@ class model:
       except KeyError:
          return 's'
 
-   def bind(self,sourceCompartment,targetCompartment,k,pc,scaleToVolume=0):
+   def bind(self,src,target,qName,pcName,useVolume=0):
       
 
-#establish a flow from one compartment to the other
-      cSrc=self.compartments[sourceCompartment]
-      cTarg=self.compartments[targetCompartment]
+      #establish a flow from source compartment to the target
 
-      vSPar=self.getVolumePar(sourceCompartment,scaleToVolume)
-      vTPar=self.getVolumePar(targetCompartment,scaleToVolume)
+      #get volume names
+      srcVName=self.getVolumePar(src,useVolume)
+      targetVName=self.getVolumePar(target,useVolume)
       
-      tu=self.getTimeUnit()
-      fk=get(tu,k)
-      fpc=get(tu,pc)
-      fvs=get(tu,vSPar)
-      fvt=get(tu,vTPar)
-
       #the source equation (where we subtract the current)
-      addValue(cSrc['targets'],sourceCompartment,-fk/fpc/fvs)
+      #in fact, this is the diagonal element
+      pSrc=self.couplingObject(-1,qName,pcName,srcVName)
+      #this includes derivatives and value!
+      self.addValueObject(src,src,pSrc)
+      
       #the target equation (where we add the current)
-      addValue(cTarg['targets'],sourceCompartment,fk/fpc/fvt)
-        
+      pTarget=self.couplingObject(1,qName,pcName,targetVName)
+      #equation is for target compartment, but binding for source
+      self.addValueObject(target,src,pTarget)
+   
+   def addValueObject(self,targetName,srcName,cObject):
+      #always binds equation id and a variable
+      targetList=self.compartments[targetName]['targets']
+      addValue(targetList,srcName,cObject["value"])
+      der=cObject["derivatives"]
+      for d in der:
+         targetSE=self.getSEJ_comp(d,targetName)
+         addValue(targetSE,srcName,der[d])
+
                     
-   def getDerivative(self,variable, sign, qPar, pcPar, vPar):
+   def couplingObject(self,sign, qParName, pcParName, vParName):
         
-      #for flow based transfer, k=Q/vS/pc, 
-      #and jacobian derivative is -Q/V/pc/pc
-      #for diffusion dominated transfer, k=k1/pc
-      
-      tu=self.getTimeUnit()
-      q=get(tu,qPar)
-      pc=get(tu,pcPar)
-      v=get(tu,vPar)
-
-      if variable=="partitionCoefficient":
-         return sign*(-1)*q/v/pc/pc
-      if variable=="diffusionCoefficient":
-         return sign/pc/v
-      if variable=="flow":
-         return sign/pc/v
-      if variable=="volume":
-         return sign*(-1)*q/pc/v/v
+      qPar=self.get(qParName)
+      pcPar=self.get(pcParName)
+      vPar=self.get(vParName)
+      q=qPar['value']
+      pc=pcPar['value']
+      v=vPar['value']
+      f=lambda t,q=q,pc=pc,v=v,s=sign:s*q(t)/v(t)/pc(t)
+
+      DPC=pcPar['derivatives']
+      dfdPC=lambda t,f=f,pc=pc:-f(t)/pc(t)
+      dPC={x:lambda t,df=dfdPC,D=DPC,x=x: df(t)*D[x](t) for x in DPC}
+      DQ=qPar['derivatives']
+      dfdQ=lambda t,f=f,q=q: f(t)/q(t)
+      dQ={x:lambda t,df=dfdQ,D=DQ,x=x: df(t)*D[x](t) for x in DQ}
+
+      DV=vPar['derivatives']
+      dfdV=lambda t,f=f,v=v: -f(t)/v(t)
+      dV={x:lambda t,df=dfdV,D=DV,x=x: df(t)*D[x](t) for x in DV}
+
+      #derivatives is the combination of the above
             
+      return derivedObject(f,[DPC,DQ,DV])
+
    def getVolumePar(self,compartment,useVolume=1):
+      #returnis volume name, if found and useVolume is directed, 
+      #or a standard parameter one
       if not useVolume:
-         return {"value":1}
+         return "one"
       try:
-         parVName=self.mod["volumes"][compartment]
-         parV=self.mod["parameters"][parVName]
+         return self.mod["volumes"][compartment]
+         #parV=self.mod["parameters"][parVName]
       except KeyError:
-         return {"value":1}
+         pass
 
-      parV["name"]=parVName
-      return parV
+      return "one"
 
-   def bindDerivative(self,parameterName,variable,src,target,q,pc,useVolume=1):
-      parName=parameterName
-      srcV=self.getVolumePar(src,useVolume)
-      doAdd=(variable!="volume" or calculateDerivative(srcV))
-         
-      if doAdd:
-         #the source equation (where we subtract the current)
-         if variable=="volume":
-             parName=srcV["name"]
-
-         csrc=self.getSEJ_comp(parName,src)
-         addValue(csrc,src,self.getDerivative(variable,-1,q,pc,srcV))
-
-      targetV=self.getVolumePar(target,useVolume)
-      doAdd=(variable!="volume" or calculateDerivative(targetV))
-
-      #the target equation (where we add the current)
-      if doAdd:
-         if variable=="volume":
-             parName=targetV["name"]
 
-         ctarget=self.getSEJ_comp(parName,target)
-         addValue(ctarget,src,self.getDerivative(variable,+1,q,pc,targetV))
-
-
-        
    def build(self):
       comps=self.compartments
       self.n=len(comps)
-      self.M=numpy.zeros((self.n,self.n))
       self.fu=[lambda t:0]*self.n
       self.lut={c:i for (i,c) in zip(range(self.n),comps.keys())}
+      self.fM={}
       for c in comps:
          comp=comps[c]
          if 'source' in comp:
              self.fu[self.lut[c]]=parseFunction(comp['source'])
                             
          for t in comp['targets']:
-             self.M[self.lut[c],self.lut[t]]=comp['targets'][t]
+            try:
+               self.fM[self.lut[c]][self.lut[t]]=\
+                  sumFunctionArray(comp['targets'][t])
+            except (KeyError,TypeError):
+               self.fM[self.lut[c]]={}
+               self.fM[self.lut[c]][self.lut[t]]=\
+                  sumFunctionArray(comp['targets'][t])
+
+
 #build SE part
-      self.m=len(self.seJ)
-      self.lutSE={c:i for (i,c) in zip(range(self.m),self.seJ.keys())}
-#MxNxN matrix
-      self.fSS=numpy.zeros((self.m,self.n,self.n))
-      for par in self.seJ:
-         sources=self.seJ[par]
+      self.buildSE()
+
+   def buildSE(self):
+      #check which parameterst to include
+      parList=[]
+
+      for parName in self.seJ:
+         #print(par)
+         par=self.mod['parameters'][parName]
+         usePar=calculateDerivative(par)
+         #print('[{}]: {}'.format(usePar,par))
+         if not usePar:
+            continue
+         parList.append(parName)
+      
+      #print(parList)
+      self.m=len(parList)
+      self.lutSE={c:i for (i,c) in zip(range(self.m),parList)}
+      
+      self.qSS={}
+      for parName in parList:
+         sources=self.seJ[parName]
          for compartment in sources:
-             targets=sources[compartment]
-             for t in targets:
-                k=self.lutSE[par]
-                i=self.lut[compartment]
-                j=self.lut[t]
-                self.fSS[k,i,j]=targets[t]
+            targets=sources[compartment]
+            for t in targets:
+               k=self.lutSE[parName]
+               i=self.lut[compartment]
+               j=self.lut[t]
+               #print('[{} {} {}] {}'.format(parName,compartment,t,targets[t]))
+               ft=sumFunctionArray(targets[t])
+               try:
+                  self.qSS[k][i][j]=ft
+               except (KeyError,TypeError):
+                  try:
+                     self.qSS[k][i]={}
+                     self.qSS[k][i][j]=ft
+                  except (KeyError,TypeError):
+                     self.qSS[k]={}
+                     self.qSS[k][i]={}
+                     self.qSS[k][i][j]=ft
+
 
                    
 
@@ -156,19 +182,19 @@ class model:
          fName=self.flows[f]
          fParName=self.mod['flows'][fName]
          fPar=pars[fParName]
-         print('\t{}[{}]:{} [{}]'.format(f,fName,fParName,get(tu,fPar)))
+         print('\t{}[{}]:{} [{}]'.format(f,fName,fParName,self.get(fParName)))
 
       print('Volumes')
       for v in self.mod['volumes']:
          vParName=self.mod['volumes'][v]
          vPar=pars[vParName]
-         print('\t{}:{} [{}]'.format(v,vParName,get(tu,vPar)))
+         print('\t{}:{} [{}]'.format(v,vParName,self.get(vParName)))
 
       print('Partition coefficients')
       for pc in self.mod['partitionCoefficients']:
          pcParName=self.mod['partitionCoefficients'][pc]
          pcPar=pars[pcParName]
-         print('\t{}:{} [{}]'.format(pc,pcParName,get(tu,pcPar)))
+         print('\t{}:{} [{}]'.format(pc,pcParName,self.get(pcParName)))
 
 
       print('SE parameters')
@@ -184,13 +210,16 @@ class model:
                     
       with open(file,'r') as f:
          self.mod=json.load(f)
+
       for m in self.mod['compartments']:
          self.add_compartment(m)
+
+      self.add_default_parameters()
+      #standard parameters such as one,zero etc.
       for s in self.mod['sources']:
          #src=mod['sources'][s]
          self.add_source(s,self.mod['sources'][s])
       self.flows={}
-      print('timeUnit: {}'.format(self.getTimeUnit()))
       pars=self.mod['parameters']
       for f in self.mod['flows']:
          #skip comments
@@ -212,21 +241,12 @@ class model:
          comps=b.split('->')
          try:
              pcParName=self.mod['partitionCoefficients'][b]
-             pcPar=pars[pcParName]
          except KeyError:
-             pcPar={"value":1}
-         
+             pcParName="one"
+        
          kParName=self.mod['bindings']['diffusion'][b]
-         kPar=pars[kParName]
-         self.bind(comps[0],comps[1],kPar,pcPar,useVolume)
-         
-         #add derivatives calculateDerivative returns true
-         if calculateDerivative(kPar):
-             self.bindDerivative(kParName,"diffusionParameter",\
-               comps[0],comps[1],kPar,pcPar,useVolume)
-         if calculateDerivative(pcPar):
-             self.bindDerivative(pcParName,"partitionCoefficient",\
-             comps[0],comps[1],kPar,pcPar,useVolume)
+         #operate with names to allow for value/function/derived infrastructure
+         self.bind(comps[0],comps[1],kParName,pcParName,useVolume)
          
       for q in self.mod['bindings']['flow']:
          comps=q.split('->')
@@ -237,9 +257,8 @@ class model:
                  #get partition coefficient
                  try:
                      pcParName=self.mod['partitionCoefficients'][cs]
-                     pcPar=pars[pcParName]
                  except KeyError:
-                     pcPar={"value":1}
+                     pcParName="one"
                  
                  #get flow (direction could be reversed)
                  try:
@@ -248,24 +267,38 @@ class model:
                      qName=self.flows['{}:{}'.format(ct,cs)]
                  
                  flowParName=self.mod['flows'][qName]
-                 flowPar=pars[flowParName]
+                 #flowPar=pars[flowParName]
                  
-                 self.bind(cs,ct,flowPar,pcPar,1)
+                 self.bind(cs,ct,flowParName,pcParName,1)
                  
-                 if calculateDerivative(pcPar):
-                     self.bindDerivative(pcParName,"partitionCoefficient",\
-                        cs,ct,flowPar,pcPar)
-                 if calculateDerivative(flowPar):
-                     self.bindDerivative(flowParName,"flow",\
-                        cs,ct,flowPar,pcPar)
-                 self.bindDerivative("x","volume",cs,ct,flowPar,pcPar)
-                      
-             
       self.build()
+   
+   def add_default_parameters(self):
+      self.mod['parameters']['one']={'value':1}
+      self.mod['parameters']['zero']={'value':0}
+   
+
+   def M(self,t):
+      Mb=numpy.zeros((self.n,self.n))
+      for i in self.fM:
+         for j in self.fM[i]:
+            Mb[i,j]=self.fM[i][j](t)
+      #create an array and fill it with outputs of function at t
+      return Mb
 
    def u(self,t):
       ub=[f(t) for f in self.fu]
       return numpy.array(ub)
+
+   def fSS(self,t):
+      fSS=numpy.zeros((self.m,self.n,self.n))
+      for k in self.qSS:
+         for i in self.qSS[k]:
+            for j in self.qSS[k][i]:
+               #print('[{},{},{}] {}'.format(k,i,j,self.qSS[k][i][j]))
+               fSS[k,i,j]=(self.qSS[k][i][j])(t)
+      return fSS
+ 
                      
    def fSY(self,y,t):
       #M number of sensitivity parameters
@@ -274,7 +307,7 @@ class model:
 
       #assume a tabulated solution y(t) at t spaced intervals
 
-      qS=self.fSS.dot(y)
+      qS=self.fSS(t).dot(y)
       #qS is MxN
       #but NxM is expected, so do a transpose
 
@@ -293,7 +326,7 @@ class model:
         
    #assume a tabulated solution y(t) at t spaced intervals
         
-      qS=self.fSS.dot(self.getY(t))
+      qS=self.fSS(t).dot(self.getY(t))
       return numpy.transpose(qS)
                      
    def getSEJ(self,parName):
@@ -351,7 +384,186 @@ class model:
          se2=numpy.multiply(fse,fse)
          s2out+=se2*s2
       return numpy.sqrt(s2out)
-        
+
+   def get(self,parName):
+      par=self.mod['parameters'][parName]
+      par['name']=parName
+      if "value" in par:
+         return self.getValue(par)
+      if "function" in par:
+         return self.getFunction(par)
+      if "derived" in par:
+         return self.getDerived(par)
+
+   def getValue(self,par):
+
+      v=par["value"]
+      parName=par['name']
+      #convert to seconds
+      try:
+         parUnits=par['unit'].split('/')
+      except (KeyError,IndexError):
+         #no unit given
+         return valueObject(v,parName)
+     
+      timeUnit=self.getTimeUnit()
+
+      try:
+         if parUnits[1]==timeUnit:
+            return valueObject(v,parName)
+      except IndexError:
+         #no / in unit name
+         return valueObject(v,parName)
+
+      if parUnits[1]=='min' and timeUnit=='s':
+         return valueObject(lambda t,v=v: v(t)/60,parName)
+      
+      if parUnits[1]=='s' and timeUnit=='min':
+         return valueObject(lambda t,v=v: 60*v(t),parName)
+
+      #no idea what to do
+      return valueObject(v,parName)
+
+   def getFunction(self,par):
+      fName=par['function']
+      #print('[{}]: getFunction({})'.format(par['name'],par['function']))
+      df=self.mod['functions'][fName]
+      if df['type']=='linearGrowth':
+         skip=['type']
+         par1={x:self.get(df[x]) for x in df if x not in skip}
+         #print(par1)
+         return linearGrowth(par1)
+
+   def getDerived(self,par):
+      dName=par['derived']
+      d=self.mod['derivedParameters'][dName]
+      #print('Derived [{}]: type {}'.format(dName,d['type']))
+      if d['type']=='product':
+         #print('{}*{}'.format(d['a'],d['b']))
+         pA=self.get(d['a'])
+         a=pA['value']
+         DA=pA['derivatives']
+         pB=self.get(d['b'])
+         b=pB['value']
+         DB=pB['derivatives']
+         #even more generic -> df/dp=[df/dA*dA/dp+df/dB*dB/dp]
+         f=lambda t,a=a,b=b,:a(t)*b(t)
+         dfdA=lambda t,b=b: b(t)
+         dfdB=lambda t,a=a: a(t)
+         dA={x:lambda t,df=dfdA,D=DA,x=x: df(t)*D[x](t) for x in DA}
+         dB={x:lambda t,df=dfdB,D=DB,x=x: df(t)*D[x](t) for x in DB}
+         return derivedObject(f,[dA,dB])
+
+      if d['type']=='power':
+         #print('{}^{}'.format(d['a'],d['n']))
+         pA=self.get(d['a'])
+         a=pA['value']
+         DA=pA['derivatives']
+         pN=self.get(d['n'])
+         n=pN['value']
+         DN=pN['derivatives']
+         f=lambda t,a=a,n=n:numpy.power(a(t),n(t))
+         dfdA=lambda t,n=n,f=f,a=a:n(t)*f(t)/a(t)
+         dfdN=lambda t,a=a,f=f:numpy.log(a(t))*f(t)
+         dA={x:lambda t,df=dfdA,D=DA,x=x: df(t)*D[x](t) for x in DA}
+         dN={x:lambda t,df=dfdN,D=DN,x=x: df(t)*D[x](t) for x in DN}
+         return derivedObject(f,[dA,dN])
+      if d['type']=='ratio':
+         #print('{}/{}'.format(d['a'],d['b']))
+         pA=self.get(d['a'])
+         a=pA['value']
+         DA=pA['derivatives']
+         pB=self.get(d['b'])
+         b=pB['value']
+         DB=pB['derivatives']
+         #even more generic -> df/dp=[df/dA*dA/dp+df/dB*dB/dp]
+         f=lambda t,a=a,b=b,:a(t)/b(t)
+         dfdA=lambda t,f=f,a=a: f(t)/a(t)
+         dfdB=lambda t,f=f,b=b: -f(t)/b(t)
+         dA={x:lambda t,df=dfdA,D=DA,x=x: df(t)*D[x](t) for x in DA}
+         dB={x:lambda t,df=dfdB,D=DB,x=x: df(t)*D[x](t) for x in DB}
+         return derivedObject(f,[dA,dB])
+      if d['type']=='sum':
+         #print('{}+{}'.format(d['a'],d['b']))
+         pA=self.get(d['a'])
+         a=pA['value']
+         DA=pA['derivatives']
+         pB=self.get(d['b'])
+         b=pB['value']
+         DB=pB['derivatives']
+         #even more generic -> df/dp=[df/dA*dA/dp+df/dB*dB/dp]
+         f=lambda t,a=a,b=b,:a(t)+b(t)
+         dfdA=lambda t: 1
+         dfdB=lambda t: 1
+         dA={x:lambda t,df=dfdA,D=DA,x=x: df(t)*D[x](t) for x in DA}
+         dB={x:lambda t,df=dfdB,D=DB,x=x: df(t)*D[x](t) for x in DB}
+         return derivedObject(f,[dA,dB])
+
+
+
+
+
+
+def calculateDerivative(par):
+   #add derivatives if dist(short for distribution) is specified
+   return "dist" in par    
+
+   
+
+def sumValues(dArray,x):
+   s=0
+   for a in dArray:
+      try:
+         s+=a[x]
+      except KeyError:
+         continue
+   return s
+         
+def sumFunctions(dArray,x):
+   #print('sumFunctions for {}'.format(x))
+   fs=[]
+   for a in dArray:
+      #print('\tChecking {}'.format(a))
+      try:
+         #check if a[x] is a function with a dummy parameter
+         v=a[x](1e7)
+         fs.append(a[x])
+         #print('\tAdding [{}]'.format(v))
+      except KeyError:
+         #print('\t{} not in table'.format(x))
+         continue
+      except TypeError:
+         #print('\tConstructing lambda for {}'.format(x))
+         fs.append(lambda t,a=a,x=x:a[x])
+   return sumFunctionArray(fs)
+
+def sumFunctionArray(fs):
+   return lambda t,fs=fs: sum([f(t) for f in fs])
+
+def valueObject(v,parName):
+   #convert everything to functions
+   d0={parName:lambda t:1}
+   try:
+      v(3)
+      return {'value':v,'derivatives':d0}
+   except TypeError:
+      pass
+   return {'value':lambda t,v=v:v,'derivatives':d0}
+
+def derivedObject(v,dArray):
+   try:
+      v(3)
+      o={'value':v}
+   except TypeError:
+      o={'value':lambda t,v=v:v}
+   allKeys=[]
+   for x in dArray:
+      allKeys.extend(x.keys()) 
+   allKeys=list(set(allKeys))
+   o['derivatives']={x:sumFunctions(dArray,x) for x in allKeys}
+   return o
+ 
+
 def splitVector(v):
    if v.find('(')<0:
       return [v]
@@ -361,22 +573,33 @@ def parseFunction(formula):
    if formula['name']=='exponential':
       c0=formula['constant']
       k=formula['k']
-      return lambda t:c0*numpy.exp(k*t)
+      return lambda t,c=c0,k=k:c*numpy.exp(k*t)
    if formula['name']=='constant':
       c0=formula['value']
-      return lambda t:c0
+      return lambda t,c0=c0:c0
    if formula['name']=='Heavyside':
       t0=formula['limit']
       v=formula['value']
-      return lambda t:v if t<t0 else 0
+      return lambda t,v=v,t0=t0:v if t<t0 else 0
    return lambda t:1
 
 def addValue(qdict,compName,v):
-    #add real number v to attribute compName of dictionary qdict, check if compName exists and handle the potential error
+   #add function to compName of dictionary qdict, 
+   #check if compName exists and handle the potential error
+   #lambda functions can't be summed directly, so qdict is a list
+   #that will be merged at matrix generation time
    try:
-      qdict[compName]+=v
+      qdict[compName].append(v)
    except KeyError:
-      qdict[compName]=v
+      qdict[compName]=[v]
+
+   #also add derivatives
+   #
+   #   for d in dTarget:
+   #      ctarget=self.getSEJ_comp(d,target)
+   #      addValue(ctarget,target,dTarget[d])
+
+
 
 
 def get(timeUnit,par):
@@ -403,6 +626,59 @@ def get(timeUnit,par):
    #no idea what to do
    return v
 
-def calculateDerivative(par):
-   #add derivatives if dist(short for distribution) is specified
-   return "dist" in par    
+def logistic(par):
+   t0=par['t0']['value']
+   W=par['W']['value']
+   L=lambda t,t0=t0,W=W:1/(1+numpy.exp(-(t-t0(t))/W(t)))
+   #D1[x] could be a function !!!
+   Dt0=par['t0']['derivatives']
+   dt0={x:lambda t,L=L,W=W,D=Dt0,x=x: L(t)*(1-L(t))/W(t)*D[x](t) for x in Dt0}
+   #D2[x] could be a function !!!
+   DW=par['W']['derivatives']
+   dW={x:lambda t,L=L,t0=t0,W=W,D=DW,x=x: L(t)*(1-L(t))*(t-t0(t))/W(t)/W(t)*D[x](t) for x in DW}
+   #merge
+   return derivedObject(L,[dt0,dW])
+
+def linearGrowth(par):
+   #make sure we never return zero or negative
+   C=par['C']['value']
+   t1=par['t1']['value']
+   t2=par['t2']['value']
+   eps=1e-8*C(0)
+   pL1=logistic({'t0':par['t1'],'W':par['W1']})
+   pL2=logistic({'t0':par['t2'],'W':par['W2']})
+   L1=pL1['value']
+   L2=pL2['value']
+   
+   deltaW=lambda t,t1=t1,t2=t2:t2(t)-t1(t)
+   deltaT=lambda t,t1=t1:t-t1(t)
+   fq=lambda t,dt=deltaT,w=deltaW:dt(t)/w(t)
+   f1=lambda t,eps=eps,q=fq,L1=L1,L2=L2,C=C:C(t)*q(t)*L1(t)*(1-L2(t))
+   f2=lambda t,L1=L1,L2=L2,C=C: C(t)*L1(t)*L2(t)
+   
+   f=lambda t,eps=eps,f1=f1,f2=f2:eps+f1(t)+f2(t)
+
+   Dt1=par['t1']['derivatives']
+   dfdt1=lambda t,q=fq,C=C,w=deltaW,L1=L1,L2=L2:\
+      -C(t)*L1(t)*(1-L2(t))*(1+q(t))/w(t)
+   dt1={x:lambda t,df=dfdt1,D=Dt1,x=x:df(t)*D[x](t) for x in Dt1}
+   
+   Dt2=par['t2']['derivatives']
+   dfdt2=lambda t,f1=f1,w=deltaW:-f1(t)/w(t)
+   dt2={x:lambda t,df=dfdt2,D=Dt2,x=x:df(t)*D[x](t) for x in Dt2}
+
+   DL1=pL1['derivatives']
+   dfdL1=lambda t,C=C,L2=L2,q=fq:C(t)*(q(t)*(1-L2(t))+L2(t))
+   dL1={x:lambda t,df=dfdL1,D=DL1,x=x:df(t)*D[x](t) for x in DL1}
+
+   DL2=pL2['derivatives']
+   dfdL2=lambda t,C=C,L1=L1,q=fq:C(t)*L1(t)*(1-q(t))
+   dL2={x:lambda t,df=dfdL2,D=DL2,x=x: df(t)*D[x](t) for x in DL2}
+
+   DC=par['C']['derivatives']
+   dfdC=lambda t,L1=L1,L2=L2,q=fq:L1(t)*(q(t)*(1-L2(t))+L2(t))
+   dC={x:lambda t,df=dfdC,D=DC,x=x:df(t)*D[x](t) for x in DC}
+
+   return derivedObject(f,[dt1,dt2,dL1,dL2,dC])
+
+