Converting sources to inputs, adjusting to account for compartment mass at injection, re-qualifying inputs as parameters, administrative move of parameter combinations (multiply, divide, add, power) to function.py

pythonScripts/cModel.py

@@ -14,8 +14,8 @@ class model:
       self.seJ={}
       self.scaled=[]
         
-   def add_source(self,compartmentName,formula):
-      self.compartments[compartmentName]['source']=formula
+   def add_input(self,compartmentName,parameterName):
+      self.compartments[compartmentName]['input']=parameterName
 
    def add_compartment(self,compartmentName):
       self.compartments[compartmentName]={}
@@ -91,7 +91,7 @@ class model:
          return function.Object(f,[dPC,dQ,dV])
       else:
          f=sign*q/v/pc
-         return derivedObject(sign*q/v/pc,\
+         return function.derivedObject(sign*q/v/pc,\
             [{'df':-f/pc,'D':DPC},\
             {'df':sign/v/pc,'D':DQ},\
             {'df':-f/v,'D':DV}])
@@ -113,14 +113,30 @@ class model:
    def build(self):
       comps=self.compartments
       self.n=len(comps)
-      self.fu=[lambda t:0]*self.n
+      #numeric representation of the input
+      self.fu=numpy.zeros((self.n))
+      #dictionary that holds potential input function objects
+      self.du={}
       self.lut={c:i for (i,c) in zip(range(self.n),comps.keys())}
       self.dM={}
       self.fM=numpy.zeros((self.n,self.n))
+      self.inputDerivatives={}
+      self.uTotal=[]
       for c in comps:
          comp=comps[c]
-         if 'source' in comp:
-             self.fu[self.lut[c]]=parseFunction(comp['source'])
+         if 'input' in comp:
+            qs=self.get(comp['input'])
+            self.uTotal.append(qs["value"])
+            qV=self.getVolumePar(c)
+            #input is a quotient (amount of exogen per unit time per volume(mass) of input compartment)
+            qs1=function.ratio(qs,self.get(qV))
+            if function.isFunction(qs1):
+               self.du[self.lut[c]]=qs1
+            else:
+               self.fu[self.lut[c]]=qs1['value']
+               #let buildSE know we have to include this derivatives
+               self.inputDerivatives[c]=qs1['derivatives']
+              
                             
          for t in comp['targets']:
             arr=comp['targets'][t]
@@ -136,7 +152,12 @@ class model:
             else:
                #just set once
                self.fM[self.lut[c],self.lut[t]]=sum(arr)
-
+      #generate total from self.uTotal
+      #ignore derivatives; uTotal is just a scaling shorthand
+      if function.contains(self.uTotal):
+         self.du[self.lut['total']]=function.Object(sumArray(self.uTotal),[])
+      else:
+         self.fu[self.lut['total']]=sum(self.uTotal)
       #build SE part
       self.buildSE()
 
@@ -144,7 +165,19 @@ class model:
       #check which parameterst to include
       parList=[]
       pars=self.parSetup['parameters']
-      for parName in self.seJ:
+      #add derivatives to jacobi terms
+      parCandidates=list(self.seJ.keys())
+      #add derivatives of input terms
+      for x in self.inputDerivatives:
+         D=self.inputDerivatives[x]
+         parCandidates.extend(list(D.keys()))
+      for x in self.du:
+         D=self.du[x]['derivatives']
+         parCandidates.extend(list(D.keys()))
+      #get rid of duplicates
+      parCandidates=list(set(parCandidates))
+
+      for parName in parCandidates:
          #print(par)
          
          par=pars[parName]
@@ -153,7 +186,7 @@ class model:
          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)}
@@ -162,8 +195,12 @@ class model:
 
       self.qSS={}
       self.SS=numpy.zeros((self.m,self.n,self.n))
+      #elements of SS will be w_p*dM_i,j/dp
       for parName in parList:
-         sources=self.seJ[parName]
+         try:
+            sources=self.seJ[parName]
+         except KeyError:
+            continue
          for compartment in sources:
             targets=sources[compartment]
             for t in targets:
@@ -188,7 +225,17 @@ class model:
                      self.qSS[k][i][j]=ft
 
 
-                 #use fM to build static part of fJ
+      #derivatives of inputs
+      #time dependent derivatives are handled in self.Su(t)
+      self.fSu=numpy.zeros((self.m,self.n))
+      for x in self.inputDerivatives:
+         D=self.inputDerivatives[x]
+         for p in D:
+            if p in parList:
+               k=self.lutSE[p]
+               self.fSu[self.lutSE[p],self.lut[x]]=D[p]*w[k]
+
+      #use fM to build static part of fJ
       N=self.n*(self.m+1)
       self.fJ=numpy.zeros((N,N))
       for i in range(self.m+1):
@@ -208,8 +255,8 @@ class model:
       for c in comps:
          print('{}/{}:'.format(c,self.lut[c]))
          comp=comps[c]
-         if 'source' in comp:
-             print('\tsource\n\t\t{}'.format(comp['source']))
+         if 'input' in comp:
+             print('\tinput\n\t\t{}'.format(comp['input']))
          print('\ttargets')
          for t in comp['targets']:
              print('\t\t{}[{},{}]: {}'.format(t,self.lut[c],self.lut[t],\
@@ -232,7 +279,7 @@ class model:
          pcParName=self.mod['partitionCoefficients'][pc]
          pcPar=pars[pcParName]
          print('\t{}:{} [{}]'.format(pc,pcParName,self.get(pcParName)))
-   
+
    def inspectSE(self):
 
       print('SE parameters')
@@ -262,9 +309,9 @@ class model:
 
       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])
+      for s in self.mod['inputs']:
+         #src=mod['inputs'][s]
+         self.add_input(s,self.mod['inputs'][s])
       self.flows={}
       #pars=self.mod['parameters']
       pars=self.parSetup['parameters']
@@ -361,9 +408,19 @@ class model:
             self.fM[i,j]=k[j]/(y[it]+eps)
 
    def u(self,t):
-      ub=[f(t) for f in self.fu]
-      ub[self.lut['total']]=sum(ub)
-      return numpy.array(ub)
+      for x in self.du:
+         self.fu[x]=self.du[x]['value'](t)
+      #this should be done previously
+      return self.fu
+
+   def Su(self,t):
+      #add time dependent values
+      for x in self.du:
+         D=self.du[x]['derivatives']
+         for p in D:
+            k=self.lutSE[p]
+            self.fSu[self.lutSE[p],self.lut[x]]=w[k]*D[p](t)
+      return self.fSu
 
    def jacobiFull(self,t):
 
@@ -481,6 +538,7 @@ class model:
       else:
          #for Gaussian, cv is sigma/value; get sigma by value multiplication
          return par["cv"]*par["cv"]*v*v
+
       
    def getMax(lutSE):
       fm=-1
@@ -497,6 +555,17 @@ class model:
          wts[j]=self.getWeight(parName)
       return wts
 
+      
+   def getVolumes(self):
+      m=numpy.zeros((len(self.lut)))
+      for p in self.lut:
+         m[self.lut[p]]=self.getVolume(p)
+      return m
+   
+   def getVolume(self,p):
+      pV=self.getVolumePar(p)
+      return self.get(pV)['value']
+
    def getDerivatives(self,se,i):
       #return latest point derivatives
       fse=se[-1][i]
@@ -580,142 +649,31 @@ class model:
       d=self.parSetup['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]
-         if any(['function' in pA,'function' in pB]):
-            fa=function.to(a)
-            fb=function.to(b)
-            f=lambda t,a=fa,b=fb:a(t)*b(t)
-            dfdA=lambda t,b=fb: b(t)
-            dfdB=lambda t,a=fa: a(t)
-            dA=function.generate(dfdA,DA)
-            dB=function.generate(dfdB,DB)
-            return function.Object(f,[dA,dB])
-         else:
-            return derivedObject(a*b,[{'df':b,'D':DA},{'df':a,'D':DB}])
+         return function.product(self.get(d['a']),self.get(d['b']))
 
       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']
-         if any(['function' in pA,'function' in pN]):
-            fa=function.to(a)
-            fn=function.to(n)
-            f=lambda t,a=fa,n=fn:numpy.power(a(t),n(t))
-            dfdA=lambda t,n=fn,f=f,a=fa:n(t)*f(t)/a(t)
-            dfdN=lambda t,a=fa,f=f:numpy.log(a(t))*f(t)
-            dA=function.generate(dfdA,DA)
-            dN=function.generate(dfdN,DN)
-            return function.Object(f,[dA,dN])
-         else:
-            f=numpy.power(a,n)
-            return derivedObject(f,[{'df':n*f/a,'D':DA},{'df':f*numpy.log(a),'D':DN}])
+         return function.power(self.get(d['a']),self.get(d['n']))
       
       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]
-         if any(['function' in pA,'function' in pB]):
-            fa=function.to(a)
-            fb=function.to(b)
-            f=lambda t,a=fa,b=fb,:a(t)/b(t)
-            dfdA=lambda t,f=f,a=fa: f(t)/a(t)
-            dfdB=lambda t,f=f,b=fb: -f(t)/b(t)
-            dA=function.generate(dfdA,DA)
-            dB=function.generate(dfdB,DB)
-            return function.Object(f,[dA,dB])
-         else:
-            return derivedObject(a/b,[{'df':1/b,'D':DA},{'df':-a/b/b,'D':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]
-         if any(['function' in pA,'function' in pB]):
-            fa=function.to(a)
-            fb=function.to(b)
-            f=lambda t,a=fa,b=fb,:a(t)+b(t)
-            dfdA=lambda t: 1
-            dfdB=lambda t: 1
-            dA=function.generate(dfdA,DA)
-            dB=function.generate(dfdB,DB)
-            return function.Object(f,[dA,dB])
-         else:
-            return derivedObject(a+b,[{'df':1,'D':DA},{'df':1,'D':DB}])
+         return function.ratio(pA=self.get(d['a']),pB=self.get(d['b']))
 
+      if d['type']=='sum':
+         return function.add(pA=self.get(d['a']),pB=self.get(d['b']))
 
 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 valueObject(v,parName):
    #convert everything to functions
    d0={parName:1}
    return {'value':v,'derivatives':{parName:1}}
 
-def derivedObject(f,ders):
-   o={'value':f}
-   DD=[]
-   for d in ders:
-      df=d['df']
-      D=d['D']
-      DD.append({x:df*D[x] for x in D})
-   allKeys=[]
-   for x in DD:
-      allKeys.extend(x.keys()) 
-   allKeys=list(set(allKeys))
-   o['derivatives']={x:sumValues(DD,x) for x in allKeys}
-   return o
- 
-
 def splitVector(v):
    if v.find('(')<0:
       return [v]
    return v[1:-1].split(',')
 
-def parseFunction(formula):
-   if formula['name']=='exponential':
-      c0=formula['constant']
-      k=formula['k']
-      return lambda t,c=c0,k=k:c*numpy.exp(k*t)
-   if formula['name']=='constant':
-      c0=formula['value']
-      return lambda t,c0=c0:c0
-   if formula['name']=='Heavyside':
-      t0=formula['limit']
-      v=formula['value']
-      return lambda t,v=v,t0=t0:v if t<t0 else 0
-   return lambda t:1
-
 def addValue(qdict,compName,v):
    #add function to compName of dictionary qdict, 
    #check if compName exists and handle the potential error

pythonScripts/function.py

@@ -194,6 +194,21 @@ def Object(v,dArray):
    o['derivatives']={x:sumDict(dArray,x) for x in allKeys}
    return o
 
+def derivedObject(f,ders):
+   o={'value':f}
+   DD=[]
+   for d in ders:
+      df=d['df']
+      D=d['D']
+      DD.append({x:df*D[x] for x in D})
+   allKeys=[]
+   for x in DD:
+      allKeys.extend(x.keys()) 
+   allKeys=list(set(allKeys))
+   o['derivatives']={x:sumValues(DD,x) for x in allKeys}
+   return o
+ 
+
 def sumDict(dArray,x):
    #print('sumFunctions for {}'.format(x))
    fs=[]
@@ -211,6 +226,20 @@ def sumDict(dArray,x):
 def sumArray(fs,w=1):
    return lambda t,fs=fs,w=w: w*sum([to(f)(t) for f in fs])
 
+def sumValues(dArray,x):
+   s=0
+   for a in dArray:
+      try:
+         s+=a[x]
+      except KeyError:
+         continue
+   return s
+
+def isFunction(vObject):
+   if callable(vObject['value']):
+      return True
+   return False
+
 def to(a):
    if callable(a):
       return a
@@ -230,3 +259,83 @@ def multiply(t,a,b):
 
 def generate(df,D):
    return {x:functools.partial(multiply,a=D[x],b=df) for x in D}
+
+def product(pA,pB):
+   #return a product of two values with derivatives
+   #print('{}*{}'.format(d['a'],d['b']))
+   a=pA['value']
+   DA=pA['derivatives']
+   b=pB['value']
+   DB=pB['derivatives']
+   if any(['function' in pA,'function' in pB]):
+      fa=to(a)
+      fb=to(b)
+      f=lambda t,a=fa,b=fb,:a(t)*b(t)
+      dfdA=lambda t,b=fb: b(t)
+      dfdB=lambda t,a=fa: a(t)
+      dA=generate(dfdA,DA)
+      dB=generate(dfdB,DB)
+      return Object(f,[dA,dB])
+   else:
+      return derivedObject(a*b,[{'df':b,'D':DA},{'df':a,'D':DB}])
+
+def power(pA,pN):
+   #return pA to the power of pN, numpy.power(pA,pN), with derivatives
+   a=pA['value']
+   DA=pA['derivatives']
+   n=pN['value']
+   DN=pN['derivatives']
+   if any(['function' in pA,'function' in pN]):
+      fa=to(a)
+      fn=to(n)
+      f=lambda t,a=fa,n=fn:numpy.power(a(t),n(t))
+      dfdA=lambda t,n=fn,f=f,a=fa:n(t)*f(t)/a(t)
+      dfdN=lambda t,a=fa,f=f:numpy.log(a(t))*f(t)
+      dA=generate(dfdA,DA)
+      dN=generate(dfdN,DN)
+      return Object(f,[dA,dN])
+   else:
+      f=numpy.power(a,n)
+      return derivedObject(f,[{'df':n*f/a,'D':DA},{'df':f*numpy.log(a),'D':DN}])
+
+def ratio(pA,pB):
+   #return ratio pA/pB with all derivatives
+   #print('{}/{}'.format(d['a'],d['b']))
+   a=pA['value']
+   DA=pA['derivatives']
+   b=pB['value']
+   DB=pB['derivatives']
+   if any(['function' in pA,'function' in pB]):
+      fa=to(a)
+      fb=to(b)
+      f=lambda t,a=fa,b=fb,:a(t)/b(t)
+      dfdA=lambda t,f=f,a=fa: f(t)/a(t)
+      dfdB=lambda t,f=f,b=fb: -f(t)/b(t)
+      dA=generate(dfdA,DA)
+      dB=generate(dfdB,DB)
+      return Object(f,[dA,dB])
+   else:
+      return derivedObject(a/b,[{'df':1/b,'D':DA},{'df':-a/b/b,'D':DB}])
+
+def add(pA,pB):
+   #return sum of pA and pB with derivatives
+#print('{}+{}'.format(d['a'],d['b']))
+   a=pA['value']
+   DA=pA['derivatives']
+   b=pB['value']
+   DB=pB['derivatives']
+#even more generic -> df/dp=[df/dA*dA/dp+df/dB*dB/dp]
+   if any(['function' in pA,'function' in pB]):
+      fa=to(a)
+      fb=to(b)
+      f=lambda t,a=fa,b=fb,:a(t)+b(t)
+      dfdA=lambda t: 1
+      dfdB=lambda t: 1
+      dA=generate(dfdA,DA)
+      dB=generate(dfdB,DB)
+      return Object(f,[dA,dB])
+   else:
+      return derivedObject(a+b,[{'df':1,'D':DA},{'df':1,'D':DB}])
+
+
+