Scaling derivative matrix with weights; implementation of static jacobi for combined system

pythonScripts/cModel.py

@@ -135,10 +135,16 @@ class model:
                   self.dM[self.lut[c]][self.lut[t]]=\
                      function.sumArray(arr)
             else:
-#just set once
+               #just set once
                self.fM[self.lut[c],self.lut[t]]=sum(arr)
-
-
+      #use fM to build static part of fJ
+      N=system.n*(system.m+1)
+      self.fJ=numpy.zeros((N,N))
+      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 system.jacobi(t) 
+    f
+    #
 
 #build SE part
       self.buildSE()
@@ -160,7 +166,9 @@ class model:
       #print(parList)
       self.m=len(parList)
       self.lutSE={c:i for (i,c) in zip(range(self.m),parList)}
-      
+      w=self.getWeights(self.lutSE)
+      w=numpy.sqrt(w)
+
       self.qSS={}
       self.SS=numpy.zeros((self.m,self.n,self.n))
       for parName in parList:
@@ -174,9 +182,9 @@ class model:
                #print('[{} {} {}] {}'.format(parName,compartment,t,targets[t]))
                arr=targets[t]
                if not function.contains(arr):
-                  self.SS[k,i,j]=sum(arr)
+                  self.SS[k,i,j]=w[k]*sum(arr)
                   continue
-               ft=function.sumArray(arr)
+               ft=function.sumArray(arr,w[k])
                try:
                   self.qSS[k][i][j]=ft
                except (KeyError,TypeError):
@@ -189,7 +197,13 @@ class model:
                      self.qSS[k][i][j]=ft
 
 
-                   
+                 #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):
+         self.fJ[i*self.n:(i+1)*self.n,i*self.n:(i+1)*self.n]=self.fM
+    
+   
 
 
    def inspect(self):
@@ -327,6 +341,17 @@ class model:
       ub=[f(t) for f in self.fu]
       return numpy.array(ub)
 
+   def jacobiFull(self,t):
+
+      #update jacobi created during build phase with time dependent values
+      for i in self.dM:
+         for j in self.dM[i]:
+            for k in range(system.m+1):
+               self.fJ[k*system.n+i,k*system.n+j]=self.dM[i][j](t)
+      return self.fJ
+
+
+
    def fSS(self,t):
       for k in self.qSS:
          for i in self.qSS[k]:
@@ -438,7 +463,9 @@ class model:
       
       s2out=numpy.zeros(se.shape[1:])
       se2=numpy.multiply(se,se)
-      return numpy.sqrt(numpy.dot(se2,self.getWeights(self.lutSE)))
+      #w=self.getWeights(self.lutSE)
+      w=numpy.ones((self.m))
+      return numpy.sqrt(numpy.dot(se2,w))
 
 
    def get(self,parName):

pythonScripts/function.py

@@ -208,8 +208,8 @@ def sumDict(dArray,x):
          continue
    return sumArray(fs)
 
-def sumArray(fs):
-   return lambda t,fs=fs: sum([to(f)(t) for f in fs])
+def sumArray(fs,w=1):
+   return lambda t,fs=fs,w=w: w*sum([to(f)(t) for f in fs])
 
 def to(a):
    if callable(a):

pythonScripts/ivp.py

@@ -59,18 +59,14 @@ def dfdySFull(t,S,system):
     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
+    return system.jacobi(t) 
 
 def solveSimultaneous(model,tmax,atol,rtol,method='LSODA',nt=201,\
       t0=0,y0=numpy.array([]), sIn=numpy.array([])):
 
    t_eval=numpy.linspace(t0,tmax,nt)
 
+
    s0=numpy.zeros((model.n,model.m+1))
 
    if sIn.size!=0:
@@ -81,6 +77,7 @@ def solveSimultaneous(model,tmax,atol,rtol,method='LSODA',nt=201,\
    #set initial condition
    s0[:,0]=y0
    s0=s0.ravel()
+
    sol=scipy.integrate.solve_ivp(dfdySFull,[t0, tmax],s0, args=(model,),\
                                  jac=jacobiSEFull, method=method, \
                                  atol=atol, rtol=rtol, t_eval=t_eval)
@@ -99,6 +96,7 @@ def solveSequential(model,tmax,atol,rtol,method='LSODA',nt=201,\
 
    if y0.size==0:
       y0=numpy.zeros(model.n)
+
    
    solIVP=scipy.integrate.solve_ivp(dfdy,[t0, tmax],y0, args=(model,), jac=jacobi,
                                   method=method, atol=atol, rtol=rtol, t_eval=t_eval)

pythonScripts/runSolver.py

@@ -55,7 +55,7 @@ def getScale(setup):
 
 
 
-def main(parFiles,jobDir):
+def main(parFiles,jobDir,startDir='NONE'):
    model=cModel.model()
    setupFile=parFiles[0]
    modelFile=parFiles[1]
@@ -65,7 +65,13 @@ def main(parFiles,jobDir):
    scale=getScale(setup)
    tmax=setup['tmax']*scale
    print('Using {}'.format(jobDir))
-   t0,y0,s0,lut,lutSE=getStartPoint(setup)
+   if startDir!='NONE':
+      print('Using solution from {}'.format(startDir))
+      t0,y0,s0,lut,lutSE=getStartPointFromDir(startDir)
+      print('t0={}'.format(t0))
+   else:
+      t0,y0,s0,lut,lutSE=getStartPoint(setup)
+
    #we should reorder S1 to match model.lutSE
    #S1 is nxm, so we should reorder columns
    sIn=numpy.zeros(s0.shape)
@@ -226,6 +232,9 @@ def loadSolutionFromDir(jobDir):
    qt=numpy.loadtxt(os.path.join(jobDir,'qt.txt'))
    return t,sol,se,sOut,qt,lut,lutSE
 
+def getStartPointFromDir(jobDir):
+   t,sol,se,sOut,qt,lut,lutSE=loadSolutionFromDir(jobDir)
+   return t[-1],sol[-1],sOut[-1],lut,lutSE
 
 def getStartPoint(setup):
    if setup['startFromRef']=='None':

setup/setupFast.json

@@ -1,6 +1,6 @@
 {"method":"LSODA",
-"atol":1e-8,
-"rtol":1e-8,
+"atol":1e-4,
+"rtol":1e-6,
 "mode":"IVPSimultaneous",
 "nt":201,
 "tmax":1,