Browse Source

First remote version

Andrej 2 years ago
parent
commit
523813258e

+ 3 - 49
models/humanHG.json

@@ -51,7 +51,7 @@
          "redBloodCells->plasma":"kRBC",
          "kidney->urine":"kU",
          "slowlyPerfused->hair":"kH",
-         "hair->slowlyPerfused":"kH",
+         "hair->venous":"kH",
          "brainBlood->brain":"kBR",
          "liver->intestine":"kB",
          "liver->inorganicMercury":"kI",
@@ -70,53 +70,7 @@
        "liver":"liverPC",
        "gut":"gutPC",
        "brainBlood->brain":"brainPC",
-       "hair->slowlyPerfused":"hairPC",
+       "hair->venous":"hairPC",
        "redBloodCells->plasma":"rbcPC"},
-     "parameters":{
-        "kidneyPC":{"value":4.0},
-        "richlyPerfusedPC":{"value":1},
-        "fatPC":{"value":0.15},
-        "slowlyPerfusedPC":{"value":2},
-        "brainBloodPC":{"value":1},
-        "brainPC":{"value":3},
-        "hairPC":{"value":248.7},
-        "placentaPC":{"value":2},
-	     "liverPC":{"value":5},
-        "gutPC":{"value":1},
-        "rbcPC":{"value":12},
-        "flowNotes":"assuming 70kg individual, fetus at 0.1kg",
-        "plasmaFlow":{"value":8.067,"unit":"l/min"},
-        "kidneyFlow":{"value":1.412,"unit":"l/min"},
-        "richlyPerfusedFlow":{"value":1.484,"unit":"l/min"},
-        "fatFlow":{"value":0.419,"unit":"l/min"},
-        "slowlyPerfusedFlow":{"value":2.01,"unit":"l/min"},
-        "brainBloodFlow":{"value":0.920,"unit":"l/min"},
-        "placentaFlow":{"value":0.173,"unit":"l/min"},
-        "liverInFlow":{"value":0.391,"unit":"l/min"},
-        "liverOutFlow":{"value":1.928,"unit":"l/min"},
-        "gutFlow":{"value":1.538,"unit":"l/min"},
-        "plasmaVolume":{"value":1.64,"unit":"kg"},
-        "redBloodCellsVolume":{"value":1.64,"unit":"kg"},
-        "kidneyVolume":{"value":0.28,"unit":"kg"},
-        "richlyPerfusedVolume":{"value":7,"unit":"kg"},
-        "fatVolume":{"value":19.11,"unit":"kg"},
-        "hairVolume":{"value":0.14,"unit":"kg"},
-        "slowlyPerfusedVolume":{"value":24.5,"unit":"kg"},
-        "brainBloodVolume":{"value":0.49,"unit":"kg"},
-        "brainVolume":{"value":1.4,"unit":"kg"},
-        "placentaVolume":{"value":0.0431,"unit":"kg"},
-        "liverVolume":{"value":1.82,"unit":"kg"},
-        "gutVolume":{"value":1.19,"unit":"kg"},
-        "intestineVolume":{"value":0.98,"unit":"kg"},
-        "kRBC":{"value":0.60,"unit":"l/min"},
-        "kU":{"value":0,"unit":"l/min"},
-        "commentkH":"assume pc-times larger input flow, 2.823e-6*248.7",
-        "kH":{"value":7.021e-4,"unit":"l/min"},
-        "kBR":{"value":4.033e-3,"unit":"l/min"},
-        "kI":{"value":4.033e-6,"unit":"l/min"},
-        "kB":{"value":4.033e-5,"unit":"l/min"},
-        "kD":{"value":4.033e-5,"unit":"l/min"},
-        "kF":{"value":8.066e-5,"unit":"l/min"},
-        "kF1":{"value":0,"unit":"l/min"},
-        "kR":{"value":2.016e-3,"unit":"l/min"}}
+     "parameters":"humanHG_parameters1.json"
 }

+ 48 - 0
models/humanHG_parameters.json

@@ -0,0 +1,48 @@
+{"parameters":{
+        "kidneyPC":{"value":4.0},
+        "richlyPerfusedPC":{"value":1},
+        "fatPC":{"value":0.15},
+        "slowlyPerfusedPC":{"value":2},
+        "brainBloodPC":{"value":1},
+        "brainPC":{"value":3},
+        "hairPC":{"value":248.7},
+        "placentaPC":{"value":2},
+	     "liverPC":{"value":5},
+        "gutPC":{"value":1},
+        "rbcPC":{"value":12},
+        "flowNotes":"assuming 70kg individual, fetus at 0.1kg",
+        "plasmaFlow":{"value":8.067,"unit":"l/min"},
+        "kidneyFlow":{"value":1.412,"unit":"l/min"},
+        "richlyPerfusedFlow":{"value":1.484,"unit":"l/min"},
+        "fatFlow":{"value":0.419,"unit":"l/min"},
+        "slowlyPerfusedFlow":{"value":2.01,"unit":"l/min"},
+        "brainBloodFlow":{"value":0.920,"unit":"l/min"},
+        "placentaFlow":{"value":0.173,"unit":"l/min"},
+        "liverInFlow":{"value":0.391,"unit":"l/min"},
+        "liverOutFlow":{"value":1.928,"unit":"l/min"},
+        "gutFlow":{"value":1.538,"unit":"l/min"},
+        "plasmaVolume":{"value":1.64,"unit":"kg"},
+        "redBloodCellsVolume":{"value":1.64,"unit":"kg"},
+        "kidneyVolume":{"value":0.28,"unit":"kg"},
+        "richlyPerfusedVolume":{"value":7,"unit":"kg"},
+        "fatVolume":{"value":19.11,"unit":"kg"},
+        "hairVolume":{"value":0.140,"unit":"kg","dist":"lognormal","cv":0.5},
+        "slowlyPerfusedVolume":{"value":24.5,"unit":"kg"},
+        "brainBloodVolume":{"value":0.49,"unit":"kg"},
+        "brainVolume":{"value":1.4,"unit":"kg"},
+        "placentaVolume":{"value":0.0431,"unit":"kg"},
+        "liverVolume":{"value":1.82,"unit":"kg"},
+        "gutVolume":{"value":1.19,"unit":"kg"},
+        "intestineVolume":{"value":0.98,"unit":"kg"},
+        "kRBC":{"value":0.60,"unit":"l/min"},
+        "kU":{"value":0,"unit":"l/min"},
+        "commentkH":"assume pc-times larger input flow, 2.823e-6*248.7",
+        "kH":{"value":2.823e-6,"unit":"l/min"},
+        "kBR":{"value":4.033e-3,"unit":"l/min"},
+        "kI":{"value":4.033e-6,"unit":"l/min"},
+        "kB":{"value":4.033e-5,"unit":"l/min"},
+        "kD":{"value":4.033e-5,"unit":"l/min"},
+        "kF":{"value":8.066e-5,"unit":"l/min"},
+        "kF1":{"value":0,"unit":"l/min"},
+        "kR":{"value":2.016e-3,"unit":"l/min"}}
+}

+ 99 - 0
models/humanHG_parameters1.json

@@ -0,0 +1,99 @@
+{"parameters":{
+        "bodyWeight":{"value":70,"units":"kg"},
+        "bodyWeight0p75":{"derived":"bwPower0p75"},
+        "kidneyPC":{"value":4.0},
+        "richlyPerfusedPC":{"value":1},
+        "fatPC":{"value":0.15},
+        "slowlyPerfusedPC":{"value":2},
+        "brainBloodPC":{"value":1},
+        "brainPC":{"value":3},
+        "hairPCBlood":{"value":248.7,"dist":"lognormal","cv":0.7},
+        "hairPC":{"derived":"hairPC"},
+        "placentaPC":{"value":2},
+	     "liverPC":{"value":5},
+        "gutPC":{"value":1},
+        "rbcPC":{"value":12},
+        "flowNotes":"assuming 70kg individual, fetus at 0.1kg",
+        "plasmaFlow":{"value":8.067,"unit":"l/min"},
+        "plasmaFlowScaled":{"value":0.333,"unit":"l/min"},
+        "kidneyFlow":{"value":1.412,"unit":"l/min"},
+        "richlyPerfusedFlow":{"value":1.484,"unit":"l/min"},
+        "fatFlow":{"value":0.419,"unit":"l/min"},
+        "slowlyPerfusedFlow":{"value":2.01,"unit":"l/min"},
+        "brainBloodFlow":{"value":0.920,"unit":"l/min"},
+        "placentaFlow":{"value":0.173,"unit":"l/min"},
+        "liverInFlow":{"value":0.391,"unit":"l/min"},
+        "liverOutFlow":{"value":1.928,"unit":"l/min"},
+        "gutFlow":{"value":1.538,"unit":"l/min"},
+        "plasmaVolume":{"value":1.64,"unit":"kg"},
+        "redBloodCellsVolume":{"value":1.64,"unit":"kg"},
+        "kidneyVolume":{"value":0.28,"unit":"kg"},
+        "richlyPerfusedVolume":{"value":7,"unit":"kg"},
+        "fatVolume":{"value":19.11,"unit":"kg"},
+        "hairVolume":{"function":"hairGrowth"},
+        "hairGrowthDuration":{"derived":"hairGrowthDuration"},
+        "hairVolumeMax":{"derived":"hairVolumeMax"},
+        "hairVolumeFraction":{"value":0.002,"dist":"lognormal","cv":0.5},
+        "slowlyPerfusedVolume":{"value":24.5,"unit":"kg"},
+        "brainBloodVolume":{"value":0.49,"unit":"kg"},
+        "brainVolume":{"value":1.4,"unit":"kg"},
+        "placentaVolume":{"value":0.0431,"unit":"kg"},
+        "liverVolume":{"value":1.82,"unit":"kg"},
+        "gutVolume":{"value":1.19,"unit":"kg"},
+        "intestineVolume":{"value":0.98,"unit":"kg"},
+        "kRBC":{"value":0.60,"unit":"l/min"},
+        "kU":{"value":0,"unit":"l/min"},
+        "commentkH":"assume pc-times larger input flow, 2.823e-6*248.7",
+        "kH":{"value":2.823e-6,"unit":"l/min"},
+        "kBR":{"value":4.033e-3,"unit":"l/min"},
+        "kI":{"value":4.033e-6,"unit":"l/min"},
+        "kB":{"value":4.033e-5,"unit":"l/min"},
+        "kD":{"value":4.033e-5,"unit":"l/min"},
+        "kF":{"value":8.066e-5,"unit":"l/min"},
+        "kF1":{"value":0,"unit":"l/min"},
+        "kR":{"value":2.016e-3,"unit":"l/min"},
+        "zero":{"value":0},
+        "threeQuarters":{"value":0.75},
+        "heavisideTransition":{"value":100},
+        "hairGrowthStart":{"value":0e6,"unit":"min"},
+        "hairGrowthStop":{"derived":"hairGrowthStop"},
+        "hairGrowth":{"derived":"hairGrowth"}},
+    "functions":{
+       "commentHairGrowth":"2yr hair growth to a total mass of 140g",
+       "hairGrowth":{
+           "type":"linearGrowth",
+           "t1":"hairGrowthStart",
+           "W1":"heavisideTransition",
+           "t2":"hairGrowthStop",
+           "W2":"heavisideTransition",
+           "C":"hairVolumeMax"}},
+    "derivedParameters":{
+       "hairVolumeMax":{
+           "type":"product",
+           "a":"hairVolumeFraction",
+           "b":"bodyWeight"},
+        "hairGrowthStop":{
+           "type":"sum",
+           "a":"hairGrowthStart",
+           "b":"hairGrowthDuration"},
+        "hairGrowthDuration":{
+            "type":"ratio",
+            "a":"hairVolumeMax",
+            "b":"hairGrowth"},
+        "hairGrowth":{
+            "type":"ratio",
+            "a":"kH",
+            "b":"hairPC"},
+       "bwPower0p75":{
+          "type":"power",
+          "a":"bodyWeight",
+          "n":"threeQuarters"},
+       "plasmaFlow":{
+          "type":"product",
+          "a":"plasmaFlowScaled",
+          "b":"bodyWeight0p75"},
+       "hairPC":{
+          "type":"ratio",
+          "a":"hairPCBlood",
+          "b":"slowlyPerfusedPC"}}
+}

+ 80 - 0
models/humanHG_parameters2.json

@@ -0,0 +1,80 @@
+{"parameters":{
+        "bodyWeight":{"value":70,"units":"kg"},
+        "bodyWeight0p75":{"derived":"bwPower0p75"},
+        "kidneyPC":{"value":4.0},
+        "richlyPerfusedPC":{"value":1},
+        "fatPC":{"value":0.15},
+        "slowlyPerfusedPC":{"value":2},
+        "brainBloodPC":{"value":1},
+        "brainPC":{"value":3},
+        "hairPC":{"value":248.7},
+        "placentaPC":{"value":2},
+	     "liverPC":{"value":5},
+        "gutPC":{"value":1},
+        "rbcPC":{"value":12},
+        "flowNotes":"assuming 70kg individual, fetus at 0.1kg",
+        "plasmaFlow":{"value":8.067,"unit":"l/min"},
+        "plasmaFlowScaled":{"value":0.333,"unit":"l/min"},
+        "kidneyFlow":{"value":1.412,"unit":"l/min"},
+        "richlyPerfusedFlow":{"value":1.484,"unit":"l/min"},
+        "fatFlow":{"value":0.419,"unit":"l/min"},
+        "slowlyPerfusedFlow":{"value":2.01,"unit":"l/min"},
+        "brainBloodFlow":{"value":0.920,"unit":"l/min"},
+        "placentaFlow":{"value":0.173,"unit":"l/min"},
+        "liverInFlow":{"value":0.391,"unit":"l/min"},
+        "liverOutFlow":{"value":1.928,"unit":"l/min"},
+        "gutFlow":{"value":1.538,"unit":"l/min"},
+        "plasmaVolume":{"value":1.64,"unit":"kg"},
+        "redBloodCellsVolume":{"value":1.64,"unit":"kg"},
+        "kidneyVolume":{"value":0.28,"unit":"kg"},
+        "richlyPerfusedVolume":{"value":7,"unit":"kg"},
+        "fatVolume":{"value":19.11,"unit":"kg"},
+        "hairVolume":{"function":"hairGrowth"},
+        "hairGrowth":{"value":5e-8,"unit":"kg/min"},
+        "hairGrowthDuration":{"value":3e6,"dist":"lognormal","cv":0.5},
+        "slowlyPerfusedVolume":{"value":24.5,"unit":"kg"},
+        "brainBloodVolume":{"value":0.49,"unit":"kg"},
+        "brainVolume":{"value":1.4,"unit":"kg"},
+        "placentaVolume":{"value":0.0431,"unit":"kg"},
+        "liverVolume":{"value":1.82,"unit":"kg"},
+        "gutVolume":{"value":1.19,"unit":"kg"},
+        "intestineVolume":{"value":0.98,"unit":"kg"},
+        "kRBC":{"value":0.60,"unit":"l/min"},
+        "kU":{"value":0,"unit":"l/min"},
+        "commentkH":"assume pc-times larger input flow, 2.823e-6*248.7",
+        "kH":{"value":2.823e-6,"unit":"l/min"},
+        "kBR":{"value":4.033e-3,"unit":"l/min"},
+        "kI":{"value":4.033e-6,"unit":"l/min"},
+        "kB":{"value":4.033e-5,"unit":"l/min"},
+        "kD":{"value":4.033e-5,"unit":"l/min"},
+        "kF":{"value":8.066e-5,"unit":"l/min"},
+        "kF1":{"value":0,"unit":"l/min"},
+        "kR":{"value":2.016e-3,"unit":"l/min"},
+        "zero":{"value":0},
+        "threeQuarters":{"value":0.75},
+        "heavisideTransition":{"value":100},
+        "hairGrowthStart":{"value":0e6,"unit":"min"},
+        "hairGrowthStop":{"derived":"hairGrowthStop"}},
+    "functions":{
+       "commentHairGrowth":"2yr hair growth to a total mass of 140g",
+       "hairGrowth":{
+           "type":"linearGrowthFixedSlope",
+           "t1":"hairGrowthStart",
+           "W1":"heavisideTransition",
+           "t2":"hairGrowthStop",
+           "W2":"heavisideTransition",
+           "alpha":"hairGrowth"}},
+    "derivedParameters":{
+        "hairGrowthStop":{
+           "type":"sum",
+           "a":"hairGrowthStart",
+           "b":"hairGrowthDuration"},
+       "bwPower0p75":{
+          "type":"power",
+          "a":"bodyWeight",
+          "n":"threeQuarters"},
+       "plasmaFlow":{
+          "type":"product",
+          "a":"plasmaFlowScaled",
+          "b":"bodyWeight0p75"}}
+}

+ 138 - 165
pythonScripts/cModel.py

@@ -4,6 +4,9 @@ import os
 import scipy.interpolate
 #for partial function specializations
 import functools
+import function
+import importlib
+importlib.reload(function)
 
 class model:
    def __init__(self):
@@ -63,22 +66,30 @@ class model:
       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}
 
+      if any(['function' in qPar,'function' in pcPar, 'function' in vPar]):
+         fq=function.to(q)
+         fpc=function.to(pc)
+         fv=function.to(v)
+         f=lambda t,q=fq,pc=fpc,v=fv,s=sign:s*q(t)/v(t)/pc(t)
+         dfdPC=lambda t,f=f,pc=fpc:-f(t)/pc(t)
+         dPC=function.generate(dfdPC,DPC)
+         dfdQ=lambda t,f=f,q=fq: f(t)/q(t)
+         dQ={x:lambda t,df=dfdQ,D=DQ,x=x: df(t)*D[x](t) for x in DQ}
+         dfdV=lambda t,f=f,v=fv: -f(t)/v(t)
+         dV={x:lambda t,df=dfdV,D=DV,x=x: df(t)*D[x](t) for x in DV}
+         return function.Object(f,[dPC,dQ,dV])
+      else:
+         f=sign*q/v/pc
+         return derivedObject(sign*q/v/pc,\
+            [{'df':-f/pc,'D':DPC},\
+            {'df':sign/v/pc,'D':DQ},\
+            {'df':-f/v,'D':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, 
@@ -99,20 +110,28 @@ class model:
       self.n=len(comps)
       self.fu=[lambda t:0]*self.n
       self.lut={c:i for (i,c) in zip(range(self.n),comps.keys())}
-      self.fM={}
+      self.dM={}
+      self.fM=numpy.zeros((self.n,self.n))
       for c in comps:
          comp=comps[c]
          if 'source' in comp:
              self.fu[self.lut[c]]=parseFunction(comp['source'])
                             
          for t in comp['targets']:
-            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])
+            arr=comp['targets'][t]
+            
+            if function.contains(arr):
+               try:
+                  self.dM[self.lut[c]][self.lut[t]]=\
+                     function.sumArray(arr)
+               except (KeyError,TypeError):
+                  self.dM[self.lut[c]]={}
+                  self.dM[self.lut[c]][self.lut[t]]=\
+                     function.sumArray(arr)
+            else:
+#just set once
+               self.fM[self.lut[c],self.lut[t]]=sum(arr)
+
 
 
 #build SE part
@@ -121,10 +140,11 @@ class model:
    def buildSE(self):
       #check which parameterst to include
       parList=[]
-
+      pars=self.parSetup['parameters']
       for parName in self.seJ:
          #print(par)
-         par=self.mod['parameters'][parName]
+         
+         par=pars[parName]
          usePar=calculateDerivative(par)
          #print('[{}]: {}'.format(usePar,par))
          if not usePar:
@@ -136,6 +156,7 @@ class model:
       self.lutSE={c:i for (i,c) in zip(range(self.m),parList)}
       
       self.qSS={}
+      self.SS=numpy.zeros((self.m,self.n,self.n))
       for parName in parList:
          sources=self.seJ[parName]
          for compartment in sources:
@@ -145,7 +166,11 @@ class model:
                i=self.lut[compartment]
                j=self.lut[t]
                #print('[{} {} {}] {}'.format(parName,compartment,t,targets[t]))
-               ft=sumFunctionArray(targets[t])
+               arr=targets[t]
+               if not function.contains(arr):
+                  self.SS[k,i,j]=sum(arr)
+                  continue
+               ft=function.sumArray(arr)
                try:
                   self.qSS[k][i][j]=ft
                except (KeyError,TypeError):
@@ -163,7 +188,8 @@ class model:
 
    def inspect(self):
       comps=self.compartments
-      pars=self.mod['parameters']
+      pars=self.parSetup['parameters']
+      #pars=self.mod['parameters']
       
       tu=self.getTimeUnit()
       print('Time unit: {}'.format(tu))
@@ -206,11 +232,15 @@ class model:
              for t in targets:
                  print('\t SE bind {}/{}:{}'.format(compartment,t,targets[t]))
     
-   def parse(self,file):
+   def parse(self,setupFile,parameterFile):
                     
-      with open(file,'r') as f:
+      with open(setupFile,'r') as f:
          self.mod=json.load(f)
 
+      
+      with open(parameterFile,'r') as f:
+         self.parSetup=json.load(f)
+
       for m in self.mod['compartments']:
          self.add_compartment(m)
 
@@ -220,7 +250,8 @@ class model:
          #src=mod['sources'][s]
          self.add_source(s,self.mod['sources'][s])
       self.flows={}
-      pars=self.mod['parameters']
+      #pars=self.mod['parameters']
+      pars=self.parSetup['parameters']
       for f in self.mod['flows']:
          #skip comments
          if f.find(':')<0:
@@ -274,30 +305,29 @@ class model:
       self.build()
    
    def add_default_parameters(self):
-      self.mod['parameters']['one']={'value':1}
-      self.mod['parameters']['zero']={'value':0}
+      pars=self.parSetup['parameters']
+      pars['one']={'value':1}
+      pars['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)
+      for i in self.dM:
+         for j in self.dM[i]:
+            self.fM[i,j]=self.dM[i][j](t)
       #create an array and fill it with outputs of function at t
-      return Mb
+      return self.fM
 
    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
+               self.SS[k,i,j]=(self.qSS[k][i][j])(t)
+      return self.SS
  
                      
    def fSY(self,y,t):
@@ -362,7 +392,8 @@ class model:
     
    def calculateUncertainty(self,sol,se):
       s2out=numpy.zeros(sol.shape)
-      pars=self.mod['parameters']
+      #pars=self.mod['parameters']
+      pars=self.parSetup['parameters']
       for parName in pars:
          par=pars[parName]
          if not calculateDerivative(par):
@@ -379,14 +410,15 @@ class model:
             s2=par["cv"]*par["cv"]*v*v
          j=self.lutSE[parName]
          fse=se[:,:,j]
-         print('Calculating for {}/{}:{}'.format(parName,j,fse.shape))
+         print('Calculating for {}/{}:{} [{}]'.format(parName,j,fse.shape,s2))
          #se2 is nt x n
          se2=numpy.multiply(fse,fse)
          s2out+=se2*s2
       return numpy.sqrt(s2out)
 
    def get(self,parName):
-      par=self.mod['parameters'][parName]
+      pars=self.parSetup['parameters']
+      par=pars[parName]
       par['name']=parName
       if "value" in par:
          return self.getValue(par)
@@ -394,6 +426,7 @@ class model:
          return self.getFunction(par)
       if "derived" in par:
          return self.getDerived(par)
+      print('Paramter {} not found!'.format(parName))
 
    def getValue(self,par):
 
@@ -427,16 +460,20 @@ class model:
    def getFunction(self,par):
       fName=par['function']
       #print('[{}]: getFunction({})'.format(par['name'],par['function']))
-      df=self.mod['functions'][fName]
+      df=self.parSetup['functions'][fName]
+      skip=['type']
+      par1={x:self.get(df[x]) for x in df if x not in skip}
       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)
+         return function.linearGrowth(par1)
+      if df['type']=='linearGrowthFixedSlope':
+         return function.linearGrowthFixedSlope(par1)
+         
+      print('Function {} not found!'.format(fName))
 
    def getDerived(self,par):
       dName=par['derived']
-      d=self.mod['derivedParameters'][dName]
+      d=self.parSetup['derivedParameters'][dName]
       #print('Derived [{}]: type {}'.format(dName,d['type']))
       if d['type']=='product':
          #print('{}*{}'.format(d['a'],d['b']))
@@ -447,12 +484,17 @@ class model:
          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 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={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 function.Object(f,[dA,dB])
+         else:
+            return derivedObject(a*b,[{'df':b,'D':DA},{'df':a,'D':DB}])
 
       if d['type']=='power':
          #print('{}^{}'.format(d['a'],d['n']))
@@ -462,12 +504,19 @@ class model:
          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 any(['function' in pA,'function' in pB]):
+            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={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 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})
+      
       if d['type']=='ratio':
          #print('{}/{}'.format(d['a'],d['b']))
          pA=self.get(d['a'])
@@ -477,12 +526,17 @@ class model:
          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 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={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 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'])
@@ -492,24 +546,23 @@ class model:
          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])
-
-
-
-
+         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={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 function.Object(f,[dA,dB])
+         else:
+            return derivedObject(a+b,[{'df':1,'D':DA},{'df':1,'D':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:
@@ -518,49 +571,25 @@ def sumValues(dArray,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}
+   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 dArray:
+   for x in DD:
       allKeys.extend(x.keys()) 
    allKeys=list(set(allKeys))
-   o['derivatives']={x:sumFunctions(dArray,x) for x in allKeys}
+   o['derivatives']={x:sumValues(DD,x) for x in allKeys}
    return o
  
 
@@ -625,60 +654,4 @@ def get(timeUnit,par):
 
    #no idea what to do
    return v
-
-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])
-
      

File diff suppressed because it is too large
+ 32 - 171
pythonScripts/compartmentModel.ipynb


+ 232 - 0
pythonScripts/function.py

@@ -0,0 +1,232 @@
+import numpy
+import functools
+
+def fval(t,W):
+   if callable(W):
+      return W(t)
+   return W
+
+def L(t,t0,W):
+   ft0=fval(t,t0)
+   fW=fval(t,W)
+   x=(t-ft0)/fW
+   if x<-18:
+      return 0
+   if x>18:
+      return 1
+   return 1/(1+numpy.exp(-x))
+def dL(t,t0,W):
+   return L(t,t0,W)*(1-L(t,t0,W))
+
+def dLdW(t,t0,W):
+   fW=fval(t,W)
+   ft0=fval(t,t0)
+   return dL(t,t0,W)*(t-ft0)/fW/fW
+
+def dLdt0(t,t0,W):
+   fW=fval(t,W)
+   return dL(t,t0,W)/fW
+
+def logisticPar(par):
+   t0=par['t0']['value']
+   W=par['W']['value']
+   f=functools.partial(L,t0=t0,W=W)
+   
+   #D1[x] could be a function !!!
+   Dt0=par['t0']['derivatives']
+   dfdt0=functools.partial(dLdt0,t0=t0,W=W)
+   dt0={x:lambda t,df=dfdt0,D=Dt0,x=x: df(t)*to(D[x])(t) for x in Dt0}
+   #D2[x] could be a function !!!
+   DW=par['W']['derivatives']
+   dfdW=functools.partial(dLdW,t0=t0,W=W)
+   dW={x:lambda t,df=dfdW,D=DW,x=x: df(t)*to(D[x])(t) for x in DW}
+   #merge
+   return Object(f,[dt0,dW])
+
+def deltaW(t,t1,t2):
+   ft1=fval(t,t1)
+   ft2=fval(t,t2)
+   return ft2-ft1
+
+def deltaT(t,t1):
+   ft1=fval(t,t1)
+   return t-t1
+
+def fq(t,t1,t2):
+   return deltaT(t,t1)/deltaW(t,t1,t2)
+
+def f1(t,C,t1,t2,L1,L2):
+   fC=fval(t,C)
+   return fC*fq(t,t1,t2)*L1(t)*(1-L2(t))
+
+def f2(t,C,L1,L2):
+   fC=fval(t,C)
+   return fC*L1(t)*L2(t)
+
+def lg(t,C,t1,t2,eps,L1,L2):
+   fC=fval(t,C)
+   return eps*fC+f1(t,C,t1,t2,L1,L2)+f2(t,C,L1,L2)
+
+def dlgdt1(t,C,t1,t2,L1,L2):
+   fC=fval(t,C)
+   return -fC*L1(t)*(1-L2(t))*(1+fq(t,t1,t2))/deltaW(t,t1,t2)
+
+def dlgdt2(t,C,t1,t2,L1,L2):
+   return -f1(t,C,t1,t2,L1,L2)/deltaW(t,t1,t2)
+
+def dlgdL1(t,C,t1,t2,L1,L2):
+   fC=fval(t,C)
+   fL2=L2(t)
+   return fC*(fq(t,t1,t2)*(1-fL2)+fL2)
+
+def dlgdL2(t,C,t1,t2,L1,L2):
+   fC=fval(t,C)
+   return fC*L1(t)*(1-fq(t,t1,t2))
+
+def dlgdC(t,C,t1,t2,L1,L2):
+   fL2=L2(t)
+   return L1(t)*(fq(t,t1,t2)*(1-fL2)+fL2)
+
+def lgFixedSlope(t,alpha,t1,t2,eps,L1,L2):
+   fA=fval(t,alpha)
+   ft1=fval(t,t1)
+   dW=deltaW(t,t1,t2)
+   return fA*(dW*eps+L1(t)*((t-ft1)*(1-L2(t))+dW*L2(t)))
+
+def dlgFixedSlopedt1(t,alpha,L1,L2):
+   fA=fval(t,alpha)
+   return -fA*L1(t)
+
+def dlgFixedSlopedt2(t,alpha,L1,L2):
+   fA=fval(t,alpha)
+   return -fA*L1(t)*L2(t)
+
+def dlgFixedSlopedL1(t,alpha,t1,t2,L1,L2):
+   return lgFixedSlope(t,alpha,t1,t2,L1,L2)/L1(t)
+
+def dlgFixedSlopedL2(t,alpha,t1,t2,L1,L2):
+   fA=fval(t,alpha)
+   ft2=fval(t,t2)
+   return fA*L1(t)*(ft2-t)
+
+def dlgFixedSlopedalpha(t,alpha,t1,t2,L1,L2):
+   fA=fval(t,alpha)
+   return lgFixedSlope(t,alpha,t1,t2,L1,L2)/fA
+
+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
+
+   pL1=logisticPar({'t0':par['t1'],'W':par['W1']})
+   pL2=logisticPar({'t0':par['t2'],'W':par['W2']})
+   L1=pL1['value']
+   L2=pL2['value']
+   f=functools.partial(lg,C=C,t1=t1,t2=t2,eps=eps,L1=L1,L2=L2)
+
+   Dt1=par['t1']['derivatives']
+   dfdt1=functools.partial(dlgdt1,C=C,t1=t1,t2=t2,L1=L1,L2=L2)
+   dt1={x:lambda t,df=dfdt1,D=Dt1,x=x:df(t)*to(D[x])(t) for x in Dt1}
+   
+   Dt2=par['t2']['derivatives']
+   dfdt2=functools.partial(dlgdt2,C=C,t1=t1,t2=t2,L1=L1,L2=L2)
+   dt2={x:lambda t,df=dfdt2,D=Dt2,x=x:df(t)*to(D[x])(t) for x in Dt2}
+
+   DL1=pL1['derivatives']
+   dfdL1=functools.partial(dlgdL1,C=C,t1=t1,t2=t2,L1=L1,L2=L2)
+   dL1={x:lambda t,df=dfdL1,D=DL1,x=x:df(t)*D[x](t) for x in DL1}
+
+   DL2=pL2['derivatives']
+   dfdL2=functools.partial(dlgdL2,C=C,t1=t1,t2=t2,L1=L1,L2=L2)
+   dL2={x:lambda t,df=dfdL2,D=DL2,x=x: df(t)*D[x](t) for x in DL2}
+
+   DC=par['C']['derivatives']
+   dfdC=functools.partial(dlgdC,C=C,t1=t1,t2=t2,L1=L1,L2=L2)
+   dC={x:lambda t,df=dfdC,D=DC,x=x:df(t)*to(D[x])(t) for x in DC}
+
+   return Object(f,[dt1,dt2,dL1,dL2,dC])
+
+def linearGrowthFixedSlope(par):
+   #make sure we never return zero or negative
+   alpha=par['alpha']['value']
+   t1=par['t1']['value']
+   t2=par['t2']['value']
+   eps=1e-8
+
+   pL1=logisticPar({'t0':par['t1'],'W':par['W1']})
+   pL2=logisticPar({'t0':par['t2'],'W':par['W2']})
+   L1=pL1['value']
+   L2=pL2['value']
+   f=functools.partial(lgFixedSlope,alpha=alpha,t1=t1,t2=t2,eps=eps,L1=L1,L2=L2)
+
+   Dt1=par['t1']['derivatives']
+   dfdt1=functools.partial(dlgFixedSlopedt1,alpha=alpha,L1=L1,L2=L2)
+   dt1={x:lambda t,df=dfdt1,D=Dt1,x=x:df(t)*to(D[x])(t) for x in Dt1}
+   
+   Dt2=par['t2']['derivatives']
+   dfdt2=functools.partial(dlgFixedSlopedt2,alpha=alpha,L1=L1,L2=L2)
+   dt2={x:lambda t,df=dfdt2,D=Dt2,x=x:df(t)*to(D[x])(t) for x in Dt2}
+
+   DL1=pL1['derivatives']
+   dfdL1=functools.partial(dlgFixedSlopedL1,alpha=alpha,t1=t1,t2=t2,L1=L1,L2=L2)
+   dL1={x:lambda t,df=dfdL1,D=DL1,x=x:df(t)*D[x](t) for x in DL1}
+
+   DL2=pL2['derivatives']
+   dfdL2=functools.partial(dlgFixedSlopedL2,alpha=alpha,t1=t1,t2=t2,L1=L1,L2=L2)
+   dL2={x:lambda t,df=dfdL2,D=DL2,x=x: df(t)*D[x](t) for x in DL2}
+
+   DA=par['alpha']['derivatives']
+   dfdA=functools.partial(dlgFixedSlopedalpha,alpha=alpha,t1=t1,t2=t2,L1=L1,L2=L2)
+   dA={x:lambda t,df=dfdA,D=DA,x=x:df(t)*to(D[x])(t) for x in DA}
+
+   return Object(f,[dt1,dt2,dL1,dL2,dA])
+
+
+
+def Object(v,dArray):
+   o={'value':v,'function':True}
+   allKeys=[]
+   for x in dArray:
+      allKeys.extend(x.keys()) 
+   allKeys=list(set(allKeys))
+   o['derivatives']={x:sumDict(dArray,x) for x in allKeys}
+   return o
+
+def sumDict(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
+         fs.append(a[x])
+         #print('\tAdding [{}]'.format(v))
+      except KeyError:
+         #print('\t{} not in table'.format(x))
+         continue
+   return sumArray(fs)
+
+def sumArray(fs):
+   return lambda t,fs=fs: sum([to(f)(t) for f in fs])
+
+def to(a):
+   if callable(a):
+      return a
+   return lambda t,a=a:a
+ 
+def contains(arr):
+   for a in arr:
+      if callable(a):
+         return True
+   return False
+
+def multiply(t,a,b):
+   #abstract actual type of a and b
+   fa=fval(t,a)
+   fb=fval(t,b)
+   return fa*fb
+
+def generate(df,D):
+   return {x:functools.partial(multiply,a=D[x],b=df) for x in D}

+ 133 - 0
pythonScripts/ivp.py

@@ -0,0 +1,133 @@
+import numpy
+import scipy.integrate
+#scipy solve_ivp interface using system of n equations
+#infrastructure for sensitivity analysis of y with respect
+#to model parameters p also integrated
+#
+#General form of equation:
+#  dy/dt=M*y+u
+#where y is an n-vector and M is nxn Jacobi matrix 
+#
+#implementation of system has to provide attributes:
+# - n, number of variables
+# - m, number of parameters
+# - M(t), nxn, a time dependent Jacobi matrix,
+# - u(t), 1xn, a time dependent source term
+# - fS(t) mxn, a time dependent array of derivatives 
+#           fS(t)[p,:] = d(M*y)/dp (t)
+#   where y is internal interpolated y-only ivp solution
+# - fSY(t,y), mxn, a time dependent array of derivatives:
+#           fS(t)[p,:] = d(M*y)/dp (t)
+#   where current y is supplied at each step
+
+def dfdy(t,y,system):
+    dfdy=system.M(t).dot(y)+system.u(t)
+    return dfdy
+
+def jacobi(t,y,system):
+    return system.M(t)
+
+#SE post calculation
+def dfdyS(t,S,system):
+    #unwrap S to NxM where M is number of parameters
+    mS=numpy.reshape(S,(system.n,system.m))
+    mOut=system.M(t).dot(mS)+system.fS(t)
+    return numpy.ravel(mOut)
+
+def jacobiSE(t,S,system):
+    N=system.n*(system.m)
+    fJ=numpy.zeros((N,N))
+    #print('fJ shape {}'.format(fJ.shape))
+    for i in range(system.m):
+        fJ[i*system.n:(i+1)*system.n,i*system.n:(i+1)*system.n]=system.M(t)
+    return fJ
+
+#SE simultaeneous calculation
+def dfdySFull(t,S,system):
+    #unwrap S to NxM where M is number of parameters
+    mS=numpy.reshape(S,(system.n,system.m+1))
+    #system.fS(y,t) is NxM matrix where M are parameters
+    y=mS[:,0]
+    mOut=system.M(t).dot(mS)+system.fSY(y,t)
+    return numpy.ravel(mOut)
+
+def jacobiSEFull(t,S,system):
+    N=system.n*(system.m+1)
+    fJ=numpy.zeros((N,N))
+    #print('fJ shape {}'.format(fJ.shape))
+    for i in range(system.m+1):
+        fJ[i*system.n:(i+1)*system.n,i*system.n:(i+1)*system.n]=system.M(t)
+    return fJ
+
+def solveSimultaneous(sys,tmax,atol,rtol,method='LSODA'):
+    S1=numpy.zeros((sys.n,sys.m+1))
+    y0=numpy.zeros(sys.n)
+    #set initial condition
+    S1[:,0]=y0
+    S1=S1.ravel()
+    sol=scipy.integrate.solve_ivp(dfdySFull,[0, tmax],S1, args=(sys,), jac=jacobiSEFull,
+                                  method=method, atol=atol, rtol=rtol)
+    t=sol.t
+    sFull=numpy.reshape(numpy.transpose(sol.y),(len(t),sys.n,sys.m+1))
+    s1=sFull[:,:,1:]
+    ysol=sFull[:,:,0]
+    se=sys.calculateUncertainty(ysol,s1)
+    print('Done simultaneous LSODA SE')
+    return t,ysol,se
+
+def solveSequential(sys,tmax,atol,rtol,method='LSODA'):
+   y0=numpy.zeros(sys.n)
+   solIVP=scipy.integrate.solve_ivp(dfdy,[0, tmax],y0, args=(sys,), jac=jacobi,
+                                  method=method, atol=atol, rtol=rtol)
+#y is n x nt (odeint nt x n)
+   sol=numpy.transpose(solIVP.y)
+   t=solIVP.t
+   print('shape (y) {}'.format(sol.shape))
+   sys.setY(t,sol)
+   S0=numpy.zeros((sys.n,sys.m))
+   S0=S0.ravel()
+    
+   solIVPSE=scipy.integrate.solve_ivp(dfdyS,[0, tmax],S0, args=(sys,), jac=jacobiSE,
+                               method=method, atol=atol, rtol=rtol)
+   sraw=numpy.reshape(numpy.transpose(solIVPSE.y),(len(solIVPSE.t),sys.n,sys.m))
+#interpolate on t
+   s1=numpy.zeros((len(t),sys.n,sys.m))
+   for i in range(sys.n):
+      for j in range(sys.m):
+         tck = scipy.interpolate.splrep(solIVPSE.t, sraw[:,i,j], s=0)
+         s1[:,i,j]=scipy.interpolate.splev(t, tck, der=0)
+
+   se=sys.calculateUncertainty(sol,s1)
+   return t,sol,se
+
+def solveSimultaneousOdeint(sys,tmax,nt=201):
+   t = numpy.linspace(0,tmax, nt)
+   y0=numpy.zeros(sys.n)
+   S1=numpy.zeros((sys.n,sys.m+1))
+#set initial condition
+   S1[:,0]=y0
+   S1=S1.ravel()
+   solSE1=scipy.integrate.odeint(dfdySFull, S1, t, args=(sys,),Dfun=jacobiSEFull,tfirst=True)
+   sFull=numpy.reshape(solSE1,(len(t),sys.n,sys.m+1))
+   s1=sFull[:,:,1:]
+   sol=sFull[:,:,0]
+   se=sys.calculateUncertainty(sol,s1)
+   print('Done simultaneous SE')
+   return t,sol,se
+
+def solveSequentialOdeint(sys,tmax,nt=201):
+   t = numpy.linspace(0,tmax, nt)
+   y0=numpy.zeros(sys.n)
+   sol = scipy.integrate.odeint(dfdy, y0=y0, t=t, args=(sys,),Dfun=jacobi,tfirst=True)
+   print('shape (y) {}'.format(sol.shape))
+
+   sys.setY(t,sol)
+   S0=numpy.zeros((sys.n,sys.m))
+   S0=S0.ravel()
+   solSE=scipy.integrate.odeint(dfdyS, S0, t, args=(sys,),Dfun=jacobiSE,tfirst=True)
+   s1=numpy.reshape(solSE,(len(t),sys.n,sys.m))
+   se=sys.calculateUncertainty(sol,s1)
+   print('Done sequential SE')
+   return t,sol,se
+
+

+ 66 - 0
pythonScripts/runSolver.py

@@ -0,0 +1,66 @@
+import ivp
+import os
+import cModel
+import sys
+import json
+import numpy
+import time
+
+def get(setup,par):
+   defaultValues={\
+      'method':'LSODA',\
+      'atol':1e-4,\
+      'rtol':1e-4,\
+      'tmax':1,\
+      'mode':'IVPSimultaneous',\
+      'nt':201,
+      'tUnit':'day'}
+   try:
+      return setup[par]
+   except KeyError:
+      pass
+   return defaultValues[par]
+
+def main(parFiles,jobDir):
+   sys=cModel.model()
+   setupFile=parFiles[1]
+   parameterFile=parFiles[2]
+   sys.parse(setupFile,parameterFile)
+   with open(parFiles[0],'r') as f:
+      setup=json.load(f)
+   tmax=get(setup,'tmax')
+   atol=get(setup,'atol')
+   rtol=get(setup,'rtol')
+   mode=get(setup,'mode')
+   method=get(setup,'method')
+   nt=get(setup,'nt')
+   tUnit=get(setup,'tUnit')
+   if tUnit=='min':
+      pass
+   if tUnit=='hour':
+      tmax*=60
+   if tUnit=='day':
+      tmax*=24*60
+   if tUnit=='month':
+      tmax*=24*60*30
+   if tUnit=='year':
+      tmax*=24*60*30*365
+
+   start_time=time.time()
+   if mode=='SequentialOdeint':
+      t,sol,se=ivp.solveSequentialOdeint(sys,tmax,nt)
+   if mode=='SimultaneousOdeint':
+      t,sol,se=ivp.solveSimultaneousOdeint(sys,tmax,nt)
+
+   if mode=='IVP':
+      t,sol,se=ivp.solveSequential(sys,tmax,atol=atol,rtol=rtol,method=method)
+        
+   if mode=='IVPSimultaneous':
+      t,sol,se=ivp.solveSimultaneous(sys,tmax,atol=atol,rtol=rtol,method=method)
+
+   end_time=time.time()
+   print('Time: {:.3f} s'.format(end_time-start_time))
+   #store
+   numpy.savetxt(os.path.join(jobDir,'t.txt'),t)
+   numpy.savetxt(os.path.join(jobDir,'sol.txt'),sol)
+   numpy.savetxt(os.path.join(jobDir,'se.txt'),se)

Some files were not shown because too many files changed in this diff