Updates to run with solve matrix

models/cDiazepam_parameters.json

@@ -1,10 +1,10 @@
 {
      "parameters":{
         "commentsInput":"X",
-        "_venousInput":{"value":0.010,"commentUnit":"sameUnitAsVolumes"},
-        "venousInput":{"function":"venousInput"},
-        "inputA":{"value":0.25e6,"commentUnit":"AU;reported concentrations are in AU/volumeUnit","dist":"normal","cv":0},
-        "inputTau":{"value":1,"commentUnit":"sameAsTimeUnit","dist":"normal","cv":0},
+        "venousInput":{"value":0.3e6,"commentUnit":"sameUnitAsVolumes"},
+        "_venousInput":{"function":"venousInput"},
+        "inputA":{"value":0.01,"commentUnit":"sameAsVolumes"},
+        "inputTau":{"value":1,"commentUnit":"sameAsTimeUnit"},
         "commentPC":"X",
         "adiposePC":{"value":23.1,"dist":"lognormal","cv":0.3},
 	     "brainPC":{"value":2.2,"dist":"lognormal","cv":0.3},
@@ -46,7 +46,7 @@
         "stomachVolume":{"value":1.1,"unit":"g"},
         "testesVolume":{"value":2.5,"unit":"g"},
         "remainderVolume":{"value":15.8,"unit":"g"},
-        "venousVolume":{"value":13.6,"unit":"g","dist":"normal","cv":0},
+        "venousVolume":{"value":13.6,"unit":"g"},
         "arterialVolume":{"value":6.8,"unit":"g"},
         "liverToExcrementNotes":"CL=422 ml/min,PC=9.1,f=0.15,V=11g->k1=0.0105/s",
         "k1":{"value":0.0105,"unit":"1/s"}},

pythonScripts/cModel.py

@@ -28,6 +28,51 @@ class model:
          return self.mod['timeUnit']
       except KeyError:
          return 's'
+   def getParameters(self):
+        """Returns the model parameters."""
+        return self.parSetup['parameters']
+
+
+   def updateParameters(self, param_updates, jsonFile, newFile):
+      """
+      Updates the parameters directly within the pars structure and writes them to a new file,
+      keeping the original structure intact.
+
+      Parameters:
+         param_updates (dict): Dictionary where keys are parameter names and values are new parameter values.
+         jsonFile (str): Path to the original JSON file (for structure reference).
+         newFile (str): Path to the new file where updated parameters will be saved.
+      """
+      # Update the parameters
+      pars = self.parSetup['parameters']
+      for parName, newValue in param_updates.items():
+         if parName in pars:
+               if "value" in pars[parName]:
+                  pars[parName]["value"] = newValue
+                  print(f"Parameter '{parName}' updated to {newValue}.")
+               else:
+                  print(f"Parameter '{parName}' does not have a 'value' field.")
+         else:
+               print(f"Parameter '{parName}' not found in the model. Adding it.")
+               pars[parName] = {"value": newValue}
+
+      # Read the original file to get the structure
+      with open(jsonFile, 'r') as file:
+         data = json.load(file)
+
+      # Merge updated parameters into the original data structure
+      data['parameters'] = self.parSetup['parameters']
+
+      # Save the updated data to a new file (not overwriting the original file)
+      with open(newFile, 'w') as file:
+         json.dump(data, file, indent=4)
+      
+      # Get relative path from the 'rezultati' directory onwards
+      relative_path = os.path.relpath(newFile, start='/home/jakob/Documents/Sola/IJS/PBPK_ociscen')
+
+      print(f"Parameters have been written to the new file: {relative_path}")
+
+
 
    def bind(self,src,target,qName,pcName):
       
@@ -91,6 +136,7 @@ class model:
          return function.Object(f,[dPC,dQ,dV])
       else:
          f=sign*q/v/pc
+         #print('Adding: {}->{} [{}] {}*{}/{}/{}: {}'.format(qParName,vParName,pcParName,sign,q,v,pc,f))
          return function.derivedObject(sign*q/v/pc,\
             [{'df':-f/pc,'D':DPC},\
             {'df':sign/v/pc,'D':DQ},\
@@ -152,6 +198,7 @@ class model:
             else:
                #just set once
                self.fM[self.lut[c],self.lut[t]]=sum(arr)
+               #print('Setting: {}->{}: {}'.format(c,t,self.fM[self.lut[c],self.lut[t]]))
       #generate total from self.uTotal
       #ignore derivatives; uTotal is just a scaling shorthand
       if function.contains(self.uTotal):
@@ -310,6 +357,22 @@ class model:
          self.scaled.append(m)
 
       self.add_default_parameters()
+      self.applyValues()
+
+   def clearCompartments(self):
+      for c in self.compartments:
+         for t in self.compartments[c]:
+            self.compartments[c][t]={}
+
+       #clear SE compartments
+      for parName in self.seJ:
+         for c in self.seJ[parName]:
+            self.seJ[parName][c]={}
+
+
+   def applyValues(self):
+      self.clearCompartments()
+      #print(self.compartments)
       #standard parameters such as one,zero etc.
       for s in self.mod['inputs']:
          #src=mod['inputs'][s]
@@ -543,7 +606,10 @@ class model:
          #return sln2*v*v
       #else:
          #for Gaussian, cv is sigma/value; get sigma by value multiplication
-      return par["cv"]*par["cv"]*v*v
+      try:
+         return par["cv"]*par["cv"]*v*v
+      except KeyError:
+         return 0
 
       
    def getMax(lutSE):
@@ -580,15 +646,55 @@ class model:
 
   
 
-   def calculateUncertainty(self,se):
+   def calculateUncertainty(self,s):
       
-      s2out=numpy.zeros(se.shape[1:])
-      se2=numpy.multiply(se,se)
+      #s2out=numpy.zeros(s1.shape[1:])
+      s2=numpy.multiply(s,s)
       #w=self.getWeights(self.lutSE)
       w=numpy.ones((self.m))
-      return numpy.sqrt(numpy.dot(se2,w))
+      return numpy.sqrt(numpy.dot(s2,w))
+
+   def getParameter(self,parName):
+#return parameter (this actually returns a reference, so changing it changes the value of internal parameter)
+      try:
+         return self.parSetup['parameters'][parName]
+      except:
+         print(f'Failed to find parameter {parName}')
 
+      return None
 
+   def setParameter(self,parName,parObj):
+#replace parameter with a new object
+#typically change distribution or cv along with the value
+      self.parSetup['parameters'][parName]=parObj
+
+   def setValue(self, parName, parValue):
+#change a single parameter parName to value parValue
+#should run applyValues after all values are set
+      pars=self.parSetup['parameters']
+      try:
+         par=pars[parName]
+      except:
+         print(f'Failed to find parameter {parName}')
+         return False
+      try:
+         par['value']=parValue
+      except KeyError:
+         print(f'Failed to set value for {parName}')
+         return False
+      return True
+
+
+
+   def setValues(self,parNames,parValues):
+#change a set of parameters in list parNames to 
+#accordingly ordered set of values parValues
+#also recalculates the matrix
+
+      for p,v in zip(parNames,parValues):
+         self.setValue(p,v)
+      self.applyValues()
+   
    def get(self,parName):
       pars=self.parSetup['parameters']
       par=pars[parName]
@@ -724,4 +830,3 @@ def get(timeUnit,par):
 
    #no idea what to do
    return v
-     

pythonScripts/runSolver.py

@@ -97,35 +97,58 @@ def getFiles(setup, getAll=False):
       inFiles.append({x:'{}{}.txt'.format(x,l) for x in readable})
    return inFiles
 
-
-def main(parFiles,jobDir,startDir='NONE'):
-   model=cModel.model()
-   setupFile=parFiles[0]
-   modelFile=parFiles[1]
-   parameterFile=parFiles[2]
-   model.parse(modelFile,parameterFile)
-   setup=parseSetup(setupFile)
-   scale=getScale(setup)
-   tmax=float(setup['tmax'])*scale
-   print('Using {}'.format(jobDir))
+def findStartPoint(startDir,setup):
    if startDir!='NONE':
       print('Using solution from {}'.format(startDir))
       #t0,y0,s0,lut,lutSE=getStartPointFromDir(startDir)
       startPoint=getStartPointFromDir(startDir)
       print('t0={}'.format(startPoint['t0']))
-   else:
-      try:
-         print('Using solution from {}'.format(setup['startFromRef']))
-      except KeyError:
-         pass
-      startPoint=getStartPoint(setup)
-      print('t0={}'.format(startPoint['t0']))
+      return startPoint
+   try:
+      print('Using solution from {}'.format(setup['startFromRef']))
+   except KeyError:
+      pass
+   startPoint=getStartPoint(setup)
+   print('t0={}'.format(startPoint['t0']))
+   return startPoint
+
+def checkForExistingSolution(model,setup,step,jobDir):
+
+   useStored=setup.get('useStored',True)
+   if not useStored:
+      return {'status':False}
+
+   #check if solution exists
+   readable=['t','sol','se','qt','sOut']
+   l=step['label']
+   files=[os.path.join(jobDir,'{}{}.txt'.format(x,l)) for x in readable]
+   filesPresent=all( [os.path.isfile(x) for x in files])
+   if not filesPresent:
+      return {'status':False}
+
+   #set t0,y0,sIn
+   latestFiles={x:'{}{}.txt'.format(x,l) for x in readable}
+   data=loadSolutionFromFiles(setup,jobDir,[latestFiles])
+   copyFields=['t','se','sol','qt']
+   out={x:data[x] for x in copyFields}
+   out.update({'status':True})
+   s0=data['sOut']#Nxnxm
+   se=data['se']
+#could be that sOut is ordered differently than in model
+   out['s1']=numpy.zeros(s0.shape)
+   for x in model.lutSE:
+      out['s1'][:,:,model.lutSE[x]]=s0[:,:,data['lutSE'][x]]
+   return out
+      
 
+def solve(model,setup,startPoint,jobDir):
+
+   #only uses values of parameters from memory
    #we should reorder S1 to match model.lutSE
-   #S1 is nxm, so we should reorder columns
    t0=startPoint['t0']
    y0=startPoint['y0']
    s0=startPoint['s0']
+   #S1 is nxm, so we should reorder columns
    sIn=numpy.zeros(s0.shape)
    if s0.size>0:
       for x in model.lutSE:
@@ -143,27 +166,30 @@ def main(parFiles,jobDir,startDir='NONE'):
       pass
    setup['t0']=t0 
    strides=getStrides(setup)
+   full_t=[t0]
+   full_sol=[y0]
+   full_sOut=[sIn]
+   full_qt=[t0]
    for step in strides:
-
-      #check if solution exists
-      readable=['t','sol','se','qt','sOut']
-      l=step['label']
-      files=[os.path.join(jobDir,'{}{}.txt'.format(x,l)) for x in readable]
-      filesPresent=all( [os.path.isfile(x) for x in files])
-      if filesPresent:
-         #set t0,y0,sIn
-         latestFiles={x:'{}{}.txt'.format(x,l) for x in readable}
-         data=loadSolutionFromFiles(setup,jobDir,[latestFiles])
-         startPoint=startPointObject(data)
-         t0=startPoint['t0']
-         y0=startPoint['y0']
-         s0=startPoint['s0']
-         for x in model.lutSE:
-            sIn[:,model.lutSE[x]]=s0[:,startPoint['lutSE'][x]]
-         print('Completed step {}: {}'.format(l,filesPresent))
+      obj=checkForExistingSolution(model,setup,step,jobDir)
+      if obj['status']:
+         print('Completed step {}'.format(step['label']))
+         t=obj['t']
+         sOut=obj['s1']
+         sol=obj['sol']
+         se=obj['se']
+         qt=obj['qt']
+         full_qt.extend(qt[1:])
+         full_sOut.extend(sOut[1:])
+         full_t.extend(t[1:])
+         full_sol.extend(sol[1:])
+#set starting values for next step
+         t0=obj['t'][-1]
+         y0=obj['sol'][-1]
+         sIn=obj['s1'][-1]
          continue
       
-      print('Step required {}'.format(l))
+      print('Step required {}'.format(step['label']))
 
       t1=t0+step['length']
       start_time=time.time()
@@ -195,7 +221,8 @@ def main(parFiles,jobDir,startDir='NONE'):
       #interpolate on s1
       #s1 has shape ft x n x m
       #where ft is LSODE driven number of points
-      qt,sOut=interpolate(setup,model,t,s1,t0,t1)
+      qt=numpy.linspace(t0,t1,setup['nt'])
+      sOut=interpolate3d(qt,t,s1)
       print('Time: {:.3f} s'.format(end_time-start_time))
       y0=sol[-1]
       t0=t[-1]
@@ -203,15 +230,59 @@ def main(parFiles,jobDir,startDir='NONE'):
 
       #store
       writeable={'t':t,'sol':sol,'se':se,'qt':qt,'sOut':sOut}
-      writeable={os.path.join(jobDir,'{}{}.txt'.format(x,l)):writeable[x]\
+      writeable={os.path.join(jobDir,'{}{}.txt'.format(x,step['label'])):writeable[x]\
          for x in writeable}
       for x in writeable:
          if x.find('sOut')>-1:
             write3D(x,writeable[x],model.lut,model.lutSE)
          else:
             numpy.savetxt(x,writeable[x])
-      print('Completed step {}'.format(l))
+
+      full_t.extend(t[1:])
+      full_sol.extend(sol[1:])
+      full_qt.extend(qt[1:])
+      full_sOut.extend(sOut[1:])
+      print('Completed step {}'.format(step['label']))
+   
    print('Completed run')
+   #interpolate on this
+   setFirstElementToCorrectShape(full_sOut)
+   setFirstElementToCorrectShape(full_sol)
+
+   full_sOut=numpy.array(full_sOut)
+   full_qt=numpy.array(full_qt)
+   full_sol=numpy.array(full_sol)
+   print(full_sOut.shape)
+   qtFull=numpy.linspace(full_t[0],full_t[-1],setup['nt'])
+   solFull=interpolate2d(qtFull,full_t,full_sol)
+   sOutFull=interpolate3d(qtFull,full_qt,full_sOut)
+   return qtFull,solFull,sOutFull
+
+def setFirstElementToCorrectShape(array):
+   if not array[0].size:
+      array[0]=numpy.zeros(array[-1].shape)
+   
+
+
+def main(parFiles,jobDir,startDir='NONE'):
+   print(parFiles)
+   model=cModel.model()
+   setupFile=parFiles[0]
+   modelFile=parFiles[1]
+   parameterFile=parFiles[2]
+   model.parse(modelFile,parameterFile)
+   setup=parseSetup(setupFile)
+   scale=getScale(setup)
+   tmax=float(setup['tmax'])*scale
+   print('Using {}'.format(jobDir))
+
+   # set start point
+   startPoint=findStartPoint(startDir,setup)
+   t,sol,sOut=solve(model,setup,startPoint,jobDir)
+   print(t)
+   writeSetup(jobDir,setup,modelFile,parameterFile)
+
+def writeSetup(jobDir,setup,modelFile,parameterFile):
    #update setup with new content (particularly t0)
    setupOut=os.path.join(jobDir,'setup.json')
    with open(setupOut,'w+') as f:
@@ -227,18 +298,30 @@ def main(parFiles,jobDir,startDir='NONE'):
       print('Written {}'.format(xOut))
 
 
-def interpolate(setup,model,t,s1,t0,tmax):
+def interpolate3d(qt,t,s1):
    #interpolate on s1
    #s1 has shape ft x n x m
    #where ft is LSODE driven number of points
-   sOut=numpy.zeros((setup['nt'],model.n,model.m))
-   qt=numpy.linspace(t0,tmax,setup['nt'])
-   for i in range(model.n):
-      for j in range(model.m):
+   fN=s1.shape[1]
+   try:
+      fM=s1.shape[2]
+   except IndexError:
+      return interpolate2d(qt,t,s1)
+
+   sOut=numpy.zeros((len(qt),fN,fM))
+   for i in range(fN):
+      for j in range(fM):
          y=s1[:,i,j]
          f = scipy.interpolate.interp1d(t, y, kind='cubic')
          sOut[:,i,j]=f(qt)
-   return qt,sOut
+   return sOut
+
+def interpolate2d(qt,t,y):
+   yOut=numpy.zeros((len(qt),y.shape[1]))
+   for i in range(y.shape[1]):
+      f=scipy.interpolate.interp1d(t,y[:,i],kind='cubic')
+      yOut[:,i]=f(qt)
+   return yOut
 
 def invertLUT(lut):
    fcode=[None]*len(lut)
@@ -436,6 +519,12 @@ def getStartPoint(setup):
    data=loadSolutionFromRef(setup)
    return startPointObject(data)
 
+def storeStartPointObject(fname,obj):
+   numpy.savez(fname,startobj=obj)
+
+def loadStartPointObject(fname):
+   return numpy.load(fname,allow_pickle=True)['startobj'].item()
+
 if __name__=="__main__":
     parFiles=sys.argv[1].split(';')
     jobDir=sys.argv[2]

pythonScripts/solveMatrix.py

@@ -125,202 +125,13 @@ def removeImMatrix(A,eps):
    f=removeImVector(f,eps)
    return f.reshape(A.shape)
  
-def solveMatrix(model,tmax,nt=201,t0=0,y0=numpy.array([]), sIn=numpy.array([]),method='full'):
-   if method=='full':
-      obj=_solveMatrix(model,tmax,nt,t0,y0,sIn)
-   if method=='separately':
-      obj=_solveMatrixSeparately(model,tmax,nt,t0,y0,sIn)
-
+def solveMatrix(model,tmax,nt=201,t0=0,y0=numpy.array([]), sIn=numpy.array([]),method='fraction'):
+#method=fraction sets output of excrement compartments to fractions of total ingested chemical
+#method=derivate sets output to amount excreted per unit time (in native units, ng/min)
+   obj=_solveMatrix(model,tmax,nt,t0,y0,sIn,method)
    return obj['t'],obj['sol'],obj['se'],obj['s1']
 
-def _solveMatrix(model,tmax,nt=201,t0=0,y0=numpy.array([]), sIn=numpy.array([])):
-   t=numpy.linspace(0,tmax-t0,nt)
-
-#solve system Q at time points t for initial values y0
-#SS is gradient of system nxn matrix M with respect to m parameters p 
-   SS=model.fSS(0)
-#large (M+1)*N by (M+1)*N matrix
-   m=len(model.lutSE)
-   A=removeScaled(model,model.M(0))
-#adjusted n where only non-terminal(scaled) compartments are counted
-   n=A.shape[0]
-   SM=numpy.zeros(((m+1)*n,(m+1)*n))
-   for k in range(m):
-#A is on diagonal
-      SM[k*n:(k+1)*n,k*n:(k+1)*n]=A
-#on the far right, gradient matrix SS sits
-      SM[k*n:(k+1)*n,-n:]=removeScaled(model,SS[k])
-#for original solutions, matrix A is used
-      SM[-n:,-n:]=A
-#matrix SM is now of full rank and can be diagonalized
-#print('{}/{}'.format(numpy.linalg.matrix_rank(SM),(m+1)*n))
-
-#RHS
-#SU has shape mxn
-#we have to remove scaled columns
-#Add plain rhs
-   SU=removeScaled(model,model.Su(0),mode='columnsOnly').flatten()
-   u=removeScaled(model,model.u(0),mode='vector')
-   print('{} {}'.format(SU.shape,u.shape))
-   SU=numpy.concatenate((SU,u))
-   print(SU.shape)
-#diagonalize matrix
-   lam,P=numpy.linalg.eig(SM)
-#no effect down to eps=1e-2
-   #eps=1e-4
-   #print(f'Using eps={eps}')
-   #lam=removeImVector(lam,eps)
-   #P=removeImMatrix(P,eps)
-   P1=numpy.linalg.inv(P)
-   #P1=removeImMatrix(P1,eps)
-   D=numpy.diag(lam)
-
-#verify
-#print(numpy.max(P @ D @ P1-M))
-
-#shift input to diagonal system
-   v=P1 @ SU
-   #also no effect
-   #maxV=numpy.max(numpy.abs(v))
-   #v=v*(numpy.abs(v)>eps*maxV)
-   #v=removeImVector(v,eps)
-
-   #print(f'v={v}')
-   #print(f'lam={lam}')
-
-#set initial values if missing
-   if not y0.size:
-      y0p=numpy.zeros(n)
-   else:
-#should move this to diagonal space
-      y0p=removeScaled(model,y0,mode='vector')
-
-   if not sIn.size:
-      s0=numpy.zeros(n*m)
-   else:
-#sIn is n by m matrix, so transpose before reshaping
-      s0=numpy.transpose(sIn)
-      s0=removeScaled(model,s0,mode='columnsOnly')
-      s0=s0.flatten()
-   print('s0 {} y0p {}'.format(s0.shape,y0p.shape))
-   S0=numpy.concatenate((s0,y0p))
-   S0= P1 @ S0   
-
-#present results as a (n*(m+1) x N array)
-#n number of variables/compartments
-#N number of time points
-   W=numpy.outer(lam,t)
-   tE=numpy.ones(t.shape)
-   V0=numpy.outer(S0,tE)
-   V=numpy.outer(v/lam,tE)
-   y1=(V0+V)*numpy.exp(W)-V
-
-#convert back to true system
-   y=P @ y1
-
-#reinsert scaled into array
-#start with y of shape (m+1)*n by N where N is number of time-points and n are only non-scaled compartments
-#should end up with y of shape (m+1)*n by N, where n is now number of ALL compartments
-   originalMap=model.lut
-   updatedMap=skipRedundant(model)
-   n0=len(originalMap)
-   n=len(updatedMap)
-   N=y.shape[1]
-   sOut=numpy.zeros((N,n0,m+1))
-   for k in range(m+1):
-      for c in updatedMap:
-         i=updatedMap[c]
-         i0=originalMap[c]
-         sOut[:,i0,k]=y[k*n+i,:]
-
-#print('Shape: {}, n={} N={} m={}'.format(yout.shape,n0,N,m))
-#equivalent of Sout (N,n,m)
-#return numpy.transpose(yout).reshape(N,n0,m+1)
-   sOut=calculateScaled(model,t,sOut,sIn,y0)
-   sol=sOut[:,:,-1]
-   s1=sOut[:,:,:-1]
-   se=getSE(s1)
-   return {'t':t+t0,'sol':sol,'se':se,'s1':s1,'P':P,'D':D,'SM':SM}
-
-
-
-
-def calculateScaled(model,t,sOut,sIn=numpy.zeros(0),y0=numpy.zeros(0)):
-
-#update sOut w/ calculated scaled values
-#sIn,yIn are initial values
-#sIn is of shape (n,m)
-#yIn is (n,)
-#sOut is (N,n,m+1)
-   lutSE=model.lutSE
-   lut=model.lut
-
-#if reading a model, keep the read lut
-#lutSE=Q['lutSE']
-#lut=Q['lut']
-
-   m=len(lutSE)
-   n=len(lut)
-   N=len(t)
-   
-   if not sIn.size:
-      sIn=numpy.zeros((n,m))
-
-   if not y0.size:
-      y0=numpy.zeros(n)
-
-   #add column for initial values
-   sIn=numpy.c_[sIn,y0]
-
-
-#reshape to n*(m+1) by N
-   y=numpy.zeros((n*(m+1),N))
-
-   for k in range(m+1):
-      for i in range(n):
-         y[n*k+i,:]=sOut[:,i,k]
-
-#full version of SM
-
-   SS=model.fSS(0)
-   A=model.M(0)
-   SM=numpy.zeros(((m+1)*n,(m+1)*n))
-   for k in range(m):
-#A is on diagonal
-      SM[k*n:(k+1)*n,k*n:(k+1)*n]=A
-#on the far right, gradient matrix SS sits
-      SM[k*n:(k+1)*n,-n:]=model.SS[k]
-
-   SM[-n:,-n:]=A
-
-#full version of RHS
-   SU=model.Su(0).flatten()
-   u=model.u(0)
-#n*(m+1) by N
-   SU=numpy.outer(numpy.concatenate((SU,u)),numpy.ones(t.shape))
-
-#integral, n*(m+1) by N
-   fI=integrate(t,y)
-   fU=integrate(t,SU)
-
-#apply couplings
-   fI= SM @ fI
-
-#update values; scale compartments to total exposure
-   iT=lut['total']
-   fTotal=y0[iT]+fU[m*n+iT,:]
-   for k in range(m+1):
-      for c in model.scaled:
-         i=lut[c]
-#sOut[1:,i,k]=fI[k*n+i,:]+fU[k*n+i,:]
-         sOut[1:,i,k]=(sIn[i,k]*y0[iT]+fI[k*n+i,:]+fU[k*n+i,:])/fTotal
-         sOut[0,i,k]=sIn[i,k]
-#set total for solution only
-   sOut[1:,iT,-1]=fTotal
-   sOut[0,iT,-1]=y0[iT]
-   return sOut
-
-def _solveMatrixSeparately(model,tmax,nt=201,t0=0,y0=numpy.array([]), sIn=numpy.array([])):
+def _solveMatrix(model,tmax,nt=201,t0=0,y0=numpy.array([]), sIn=numpy.array([]),method='fraction'):
 #sIn (n,m)
 #y0 (n,)
    t=numpy.linspace(0,tmax-t0,nt)
@@ -342,8 +153,7 @@ def _solveMatrixSeparately(model,tmax,nt=201,t0=0,y0=numpy.array([]), sIn=numpy.
   
    sOut=numpy.zeros([N,n,m+1])
 #y is (N,n)
-   #y=solveUnit(model,t,y0,u0,model.u(0))   
-   y=solver.solve(y0,model.u(0))   
+   y=solver.solve(y0,model.u(0),method=method)   
    sOut[:,:,-1]=y
 #(N,n,1)
    y1=numpy.expand_dims(y,axis=2)
@@ -364,12 +174,9 @@ def _solveMatrixSeparately(model,tmax,nt=201,t0=0,y0=numpy.array([]), sIn=numpy.
 #yIn is (n,)
       yIn=sIn[:,k]
 #return in (N,n)
-#      sOut[:,:,k]=solveUnit(model,t,yIn,u0,u,rhs)
-      #sOut[:,:,k]=solver.solve(sIn[:,k],SU[k],rhs)
-      sOut[:,:,k]=solver.solveS(SS[k],y1,y0,u=u,z0=sIn[:,k],r=SU[k])
+      sOut[:,:,k]=solver.solveS(SS[k],y,y0,u=u,z0=sIn[:,k],r=SU[k],method=method)
 
 
-#      break
       
    sol=sOut[:,:,-1]
    s1=sOut[:,:,:-1]
@@ -391,8 +198,10 @@ class solveUnitObj:
       self.n=len(self.updatedMap)
  
       self.N=len(t)
+
+      self.M=self.model.M(0)
       
-      self.A=removeScaled(model,model.M(0))
+      self.A=removeScaled(model,self.M)
 
 #diagonalize matrix
       self.lam,self.P=numpy.linalg.eig(self.A)
@@ -448,12 +257,11 @@ class solveUnitObj:
 #
 
 
-   def solve(self,y0=numpy.zeros(0),u=numpy.zeros(0),rhs=numpy.zeros(0)):
+   def solve(self,y0=numpy.zeros(0),u=numpy.zeros(0),method='fraction'):
 
-#solve system Q at time points t for initial values y0
-#u is the time independent nohomogeneous part
-#rhs is tabulated non-homogeneous part (N,n)
-#u0 is initial accumulated xenobiote
+#solve system
+#y0 (n,) initial values
+#u (n,) the time independent input
 
 
 #use Green function theorem, since impulse response is simple
@@ -479,66 +287,35 @@ class solveUnitObj:
          y0=numpy.zeros(u.shape)
       y0p=removeScaled(self.model,y0,mode='vector')
 
-#overload w/ rhs
-#RHS has a time component
-#typically of shape n by N
-      if rhs.size:
-#to N by n shape
-#      rhs=rhs.transpose()
-#remove trivial components
-         rhsp=removeScaled(self.model,rhs,mode='columnsOnly')
-#to (N,n,1)
-         rhsp=numpy.expand_dims(rhsp,axis=2) 
-      else:
-         rhsp=numpy.zeros(0)
-   #print(numpy.linalg.matrix_rank(A))
-
-
-
-#(n,n) @ (N,n,n) @ (n,n) @ (n,)
-#y1 is (N,n)
-#sum works for (N,n)+(n,)
-   #print(up)
-   #print(P1 @ up)
-   #print( (W-E)*invL @ P1 @ up)
+#(N,n)=(n,n) @ (N,n,n) @ (n,n) @ (n,) + (n,n) @ (N,n,n) @ (n,n) @ (n,)
       y= self.P @ self.W @ self.P1 @ y0p + self.P @ self.Q @ self.P1 @ up
 
-      if rhsp.size:
-
-#time dependent component
-#(a,b,c,e)=(a,b,c,d) @ (b,d,e)
-#(N,N,n,1)=(n,n) @ (N,N,n,n) @ (n,n) @ (N,n,1)
-#Y1 is (N,N,n,1)
-         Y1=numpy.real(self.P @ self.W1 @ self.P1 @ rhsp)
-#to apply Green's theorem, integral along axis=1 should be performed
-#(N,n,1)
-      
-         Y2=integrateFull(self.t,Y1,axis=1)
-#to remove last dimension, squeeze at the end
-#(N,n), same as y1
-         y+=Y2.squeeze()
-     
 #back to original system, add trivial compartments
+#(N,n)
       yOut=numpy.zeros((self.N,self.n0))
-   
+
+#shift compartments from modified index i (from y) to original index i0 (to yOut)
       for c in self.updatedMap:
          i=self.updatedMap[c]
          i0=self.originalMap[c]
          yOut[:,i0]=y[:,i]
   
-      #return yOut
-#here we need full sized y0,u,rhs, so make sure they are not overwritten by reduced versions
-      return getScaled(self.model,self.t,yOut,self.u0,y0,u,rhs)
+#here we need original, full sized y0,u,rhs, so make sure they are not overwritten by reduced versions
+#content determined by method
+      if method=='fraction':
+         return getScaled(self.model,self.t,yOut,self.u0,y0,u)
+      if method=='derivative':
+         return setScaledDerivative(self.model,yOut,u)
 
    def solveS(self,S,y,y0=numpy.zeros(0),u=numpy.zeros(0),\
-      z0=numpy.zeros(0),r=numpy.zeros(0)):
+      z0=numpy.zeros(0),r=numpy.zeros(0),method='fraction'):
 
 #solve gradient system (S) at time points t 
-#y (N,n,1) is solution of the system 
+#y (N,n) is solution of the system 
 #y0 (n,) initial values of solution  
-#u (n,)is the non-homogenous (but constant) part of equation
-#z0 (n,)initial values of gradients
-#r (n,) is the time independent gradient of nonhomogenous part
+#u (n,) is the non-homogenous (but constant) part of equation
+#z0 (n,) initial values of gradients
+#r (n,) is the time independent gradient of nonhomogenous part / input (grad u)
 
 #use Green function theorem, since impulse response is simple
 #yGREEN=exp(lam*(t-t0) and integrate it for exponentail solutions 
@@ -581,25 +358,22 @@ class solveUnitObj:
 #(N,n,n)
       HS= numpy.einsum('a,bc',self.tE,Spp) * self.H
 
-#z=P @ W @ P1 @ z0 + P @ Q @ P1 r + 
-# + P @ (S * H) @ P1 @ y0 + P @ (S * H) @ D1 @ P1 @ u +P @ S @ Q @ D1 @ P1 @ u
-#(n,n) @ (N,n,n) @ (n,n) @ (n,)
-#y1 is (N,n)
-#sum works for (N,n)+(n,)
+#(N,n) = (n,n) @ (N,n,n) @ (n,n) @ (n,)
 
       z1 = self.P @ self.W @ self.P1 @ z0p 
       z2 = self.P @ self.Q @ self.P1 @ rp
 
 #intermediate variable for speed-up
-#(n,n) @ (n,n) @(n,)
+#(n,) = (n,n) @ (n,n) @(n,)
       _f=self.D1 @ self.P1 @ up
 
       #combine above for speed-up
-#(n,n) @ (N,n,n) @ ( (n,n) @ (n,) + (n,) )
+#(N,n) = (n,n) @ (N,n,n) @ ( (n,n) @ (n,) + (n,) )
       z3 = self.P @ HS @ (self.P1 @ y0p + _f)
-#(n,n) @ (n,n) @ (N,n,n) @ (n,)
+#(N,n) = (n,n) @ (n,n) @ (N,n,n) @ (n,)
       z4 = self.P @ self.Q @ Spp @ _f
      
+#(N,n)
       z=z1+z2+(z3-z4)
 
 #back to original system, add trivial compartments
@@ -610,134 +384,18 @@ class solveUnitObj:
          i0=self.originalMap[c]
          sOut[:,i0]=z[:,i]
   
-      #return yOut
 #here we need full sized y0,u,rhs, so make sure they are not overwritten by reduced versions
-      #y = self.P @ self. W @ self.P1 @ y + self.P @ self.Q @ self.P1 @ u
-      #(N,n0) = (n0,n0) @ (N,n0,1)
-      rhs = (S @ y).squeeze()
-      return getScaled(self.model,self.t,sOut,self.u0,z0,r,rhs)
+#(N,n) = (n,n) @ (N,n,1)
+      rhs = (S @ numpy.expand_dims(y,axis=2)).squeeze()
+      if method=='fraction':
+         return setScaled(self.model,self.t,sOut,self.u0,z0,r,rhs)
+      if method=='derivative':
+         return setScaledDerivative(self.model,sOut,r,rhs)
 
 
  
 
-def solveUnit(model,t,y0=numpy.zeros(0),u0=0,u=numpy.zeros(0),rhs=numpy.zeros(0)):
-#solve system Q at time points t for initial values y0
-#u is the time independent nohomogeneous part
-#rhs is tabulated non-homogeneous part (N,n)
-#u0 is initial accumulated xenobiote
-
-#use Green function theorem, since impulse response is simple
-#yGREEN=exp(lam*(t-t0)
-   
-#some dimensions
-   originalMap=model.lut
-   updatedMap=skipRedundant(model)
-   n0=len(originalMap)
-   n=len(updatedMap)
-   N=len(t)
-
-#matrix solution works for reduced space where all
-#non-trivial components are reduced
-
-#typical variables are denoted with trailing p (for prime)
-#and have a smaller size (fewer compartments) than
-#full variables
-
-#get the matrix and remove trivial components
-   A=removeScaled(model,model.M(0))
-#remove trivial components from non-homogeneous parts
-   if not u.size:
-      u=numpy.zeros(n0)
-
-   up=removeScaled(model,u,mode='vector')
-
-#remove trivial components from initial values
-   if not y0.size:
-      y0=numpy.zeros(u.shape)
-   y0p=removeScaled(model,y0,mode='vector')
-
-#overload w/ rhs
-#RHS has a time component
-#typically of shape n by N
-   if rhs.size:
-#to N by n shape
-#      rhs=rhs.transpose()
-#remove trivial components
-      rhsp=removeScaled(model,rhs,mode='columnsOnly')
-#to (N,n,1)
-      rhsp=numpy.expand_dims(rhsp,axis=2) 
-   else:
-      rhsp=numpy.zeros(0)
-   #print(numpy.linalg.matrix_rank(A))
-
-#diagonalize matrix
-   lam,P=numpy.linalg.eig(A)
-   P1=numpy.linalg.inv(P)
-   D=numpy.diag(lam)
-
-#next are all (N,n,n)
-   W=numpy.exp(numpy.einsum('a,bc',t,numpy.diag(lam)))
-   tE=numpy.ones(t.shape)
-   invL=numpy.einsum('a,bc',tE,numpy.diag(1/lam))
-   E=numpy.ones(W.shape)
-
-#(n,n) @ (N,n,n) @ (n,n) @ (n,)
-#y1 is (N,n)
-#sum works for (N,n)+(n,)
-   #print(up)
-   #print(P1 @ up)
-   #print( (W-E)*invL @ P1 @ up)
-   y= y0p + P @ ((W-E)*invL) @ P1 @ up
-
-   if rhsp.size:
-
-#time dependent component
-#t-t'
-#we have exp(lam*(t-t')) as Green's function
-#construct matrix (N,N,n,n) where the last part is diag(lam)
-#W[i,j,k,l]=delta(k,l)*exp(lam_k*(t_i-t'_j))
-#delta(i,j) is Kronecker's delta
-      T=numpy.outer(t,tE)-numpy.outer(tE,t)
-      T1=numpy.einsum('ab,cd',T,numpy.eye(n))
-      F=numpy.einsum('ab,cd',T>=0,numpy.eye(n))
-      L1=numpy.einsum('ab,cd',numpy.ones(T.shape),numpy.diag(lam))
-      W1=numpy.zeros(F.shape,dtype='complex')
-      W1[F>0]=numpy.exp(L1[F>0]*T1[F>0])
-#(a,b,c,e)=(a,b,c,d) @ (b,d,e)
-#(N,N,n,1)=(n,n) @ (N,N,n,n) @ (n,n) @ (N,n,1)
-#Y1 is (N,N,n,1)
-      Y1=P @ W1 @ P1 @ rhsp
-#to apply Green's theorem, integral along axis=1 should be performed
-#(N,n,1)
-      
-      Y2=integrateFull(t,Y1,axis=1)
-#to remove last dimension, squeeze at the end
-#(N,n), same as y1
-      y+=Y2.squeeze()
-     
-#back to original system, add trivial compartments
-   yOut=numpy.zeros((N,n0))
-   
-   for c in updatedMap:
-      i=updatedMap[c]
-      i0=originalMap[c]
-      yOut[:,i0]=y[:,i]
-   
-#here we need full sized y0,u,rhs, so make sure they are not overwritten by reduced versions
-   return getScaled(model,t,yOut,u0,y0,u,rhs)
-
-def getTotal(model,t,u0=0):
-   u=model.u(0)
-   iT=model.lut['total']
-   U=u[iT]
-#*numpy.ones(t.shape)
-   fU=numpy.zeros(t.shape)
-   fU[1:]=u0+integrate(t,U).squeeze()
-   fU[0]=u0
-   return fU
-
-
-def getScaled(model,t,y,u0=0,y0=numpy.zeros(0),u=numpy.zeros(0),rhs=numpy.zeros(0)):
+def setScaled(model,t,y,u0=0,y0=numpy.zeros(0),u=numpy.zeros(0),rhs=numpy.zeros(0)):
 #update y with scaled containers, which were excluded from direct calculation to make M of full rank
 #u has shape (n,)
 #rhs has shape (N,n)
@@ -784,3 +442,35 @@ def getScaled(model,t,y,u0=0,y0=numpy.zeros(0),u=numpy.zeros(0),rhs=numpy.zeros(
    
    return y
 
+
+def setScaledDerivative(model,y,u=numpy.zeros(0),rhs=numpy.zeros(0)):
+#update y with scaled containers, which were excluded from direct calculation to make M of full rank. 
+#Rhs is added to the argument list so that the same routine can be used for Jacobi calculation
+#M is (n,n), the system matrix, M=model.M(0)
+#y is (N,n), the solution to the coupled compartments (ie. excluded excrement)
+#u is (n,), the input part
+#rhs is (N,n), for Jacobi, set it to (grad M) @ y, for k-th component of gradient
+
+#assemble non-homogenous part of equation
+   if not rhs.size:
+#(N,n)
+      rhs=numpy.zeros(y.shape)
+
+   if u.size:
+#rhs is (n,n) so adding (n,) should work, rhs remains (N,n)
+      rhs+=u
+
+#apply couplings and inputs
+#(n,n)
+   M=model.M(0)
+#(N,n)=[(n,n) @ (N,n,1)] + (N,n)
+   fI=( M @ numpy.expand_dims(y,axis=2) ).squeeze() + rhs
+
+#only set values back to y for excrement compartments (scaled)
+   for c in model.scaled:
+      i=model.lut[c]
+      y[:,i]=fI[:,i]
+
+#(N,n)
+   return y
+

setup/setupDayMatrix.json

@@ -0,0 +1,10 @@
+{
+"atol":1e-8,
+"rtol":1e-8,
+"mode":"solveMatrix",
+"method":"separately",
+"nt":201,
+"tmax":1,
+"tUnit":"day",
+"stride":{"length":1,"tUnit":"day"}
+}

setup/setupMinuteMatrix.json

@@ -0,0 +1,10 @@
+{
+"atol":1e-8,
+"rtol":1e-8,
+"mode":"solveMatrix",
+"method":"separately",
+"nt":201,
+"tmax":1,
+"tUnit":"min",
+"stride":{"length":1,"tUnit":"min"}
+}

setup/setupMonthMatrix.json

@@ -0,0 +1,10 @@
+{
+"atol":1e-8,
+"rtol":1e-8,
+"mode":"solveMatrix",
+"method":"separately",
+"nt":201,
+"tmax":1,
+"tUnit":"month",
+"stride":{"length":1,"tUnit":"month"}
+}

setup/setupYearMatrix.json

@@ -0,0 +1,10 @@
+{
+"atol":1e-8,
+"rtol":1e-8,
+"mode":"solveMatrix",
+"method":"separately",
+"nt":201,
+"tmax":2,
+"tUnit":"year",
+"stride":{"length":2,"tUnit":"year"}
+}