Adding striding for intermediate solution storage

pythonScripts/cModel.py

@@ -137,16 +137,8 @@ class model:
             else:
                #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
+      #build SE part
       self.buildSE()
 
    def buildSE(self):

pythonScripts/ivp.py

@@ -59,7 +59,7 @@ def dfdySFull(t,S,system):
     return numpy.ravel(mOut)
 
 def jacobiSEFull(t,S,system):
-    return system.jacobi(t) 
+   return system.jacobiFull(t) 
 
 def solveSimultaneous(model,tmax,atol,rtol,method='LSODA',nt=201,\
       t0=0,y0=numpy.array([]), sIn=numpy.array([])):

pythonScripts/runSolver.py

@@ -37,10 +37,17 @@ def parseSetup(setupFile):
    out={}
    for p in defaultValues:
       out[p]=get(setup,p)
+   #q=[x if not x in out for x in setup]
+   q=[x for x in setup if not x in out]
+   for x in q:
+      out[x]=setup[x]
    return out
 
 def getScale(setup):
    tUnit=setup['tUnit']
+   return convertUnit(tUnit)
+
+def convertUnit(tUnit):
    if tUnit=='min':
       return 1
    if tUnit=='hour':
@@ -53,7 +60,31 @@ def getScale(setup):
       return 24*60*30*365
    return 1
 
-
+def getStrides(setup):
+   tmax=setup['tmax']*convertUnit(setup['tUnit'])
+   t0=setup['t0']
+   try:
+      stride=setup['stride']
+   except KeyError:
+      return [{'label':'','length':tmax-t0}]
+
+   #cover tmax with steps of length stride['length']
+   strides=[]
+   strideLength=stride['length']*convertUnit(stride['tUnit'])
+   i=0
+   while t0<tmax:
+      strides.append({'length':strideLength,'label':'_step{:03d}'.format(i)})
+      t0+=strideLength
+      i+=1
+   return strides
+
+def getLastStepFiles(setup):
+   strides=getStrides(setup)
+   readable=['t','sol','se','qt','sOut']
+   #get last stride
+   step=strides[-1]
+   l=step['label']
+   return {x:'{}{}.txt'.format(x,l) for x in readable}
 
 def main(parFiles,jobDir,startDir='NONE'):
    model=cModel.model()
@@ -79,41 +110,58 @@ def main(parFiles,jobDir,startDir='NONE'):
       for x in model.lutSE:
          #this could be modified for different ordering of lut as well
          sIn[:,model.lutSE[x]]=s0[:,lutSE[x]]
-
-   start_time=time.time()
-   if setup['mode']=='SequentialOdeint':
-      t,sol,se,s1=ivp.solveSequentialOdeint(model,tmax,setup['nt'],t0,y0,sIn)
-   if setup['mode']=='SimultaneousOdeint':
-      t,sol,se,s1=ivp.solveSimultaneousOdeint(model,tmax,setup['nt'],t0,y0,sIn)
-
-   if setup['mode']=='IVP':
-      t,sol,se,s1=ivp.solveSequential(model,tmax,\
-         nt=setup['nt'],\
-         atol=setup['atol'],\
-         rtol=setup['rtol'],method=setup['method'],\
-         t0=t0,y0=y0,sIn=sIn)
+   #copy t0 to setup to help get the strides right
+   setup['t0']=t0 
+   strides=getStrides(setup)
+   for step in strides:
+      t1=t0+step['length']
+      start_time=time.time()
+      if setup['mode']=='SequentialOdeint':
+         t,sol,se,s1=ivp.solveSequentialOdeint(model,t1,setup['nt'],t0,y0,sIn)
+      if setup['mode']=='SimultaneousOdeint':
+         t,sol,se,s1=ivp.solveSimultaneousOdeint(model,t1,setup['nt'],t0,y0,sIn)
+
+      if setup['mode']=='IVP':
+         t,sol,se,s1=ivp.solveSequential(model,t1,\
+            nt=setup['nt'],\
+            atol=setup['atol'],\
+            rtol=setup['rtol'],method=setup['method'],\
+            t0=t0,y0=y0,sIn=sIn)
         
-   if setup['mode']=='IVPSimultaneous':
-      t,sol,se,s1=ivp.solveSimultaneous(model,tmax,\
-         nt=setup['nt'],\
-         atol=setup['atol'],\
-         rtol=setup['rtol'],method=setup['method'],\
-         t0=t0,y0=y0,sIn=sIn)
-
-   end_time=time.time()
-   print('Time: {:.3f} s'.format(end_time-start_time))
-   #store
-   #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,tmax)
-
-   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)
-   numpy.savetxt(os.path.join(jobDir,'qt.txt'),qt)
-   #this is 3D, so new routines
-   write3D(os.path.join(jobDir,'sOut.txt'),sOut,model.lut,model.lutSE)
+      if setup['mode']=='IVPSimultaneous':
+         t,sol,se,s1=ivp.solveSimultaneous(model,t1,\
+            nt=setup['nt'],\
+            atol=setup['atol'],\
+            rtol=setup['rtol'],method=setup['method'],\
+            t0=t0,y0=y0,sIn=sIn)
+
+      end_time=time.time()
+      #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)
+      print('Time: {:.3f} s'.format(end_time-start_time))
+      y0=sol[-1]
+      t0=t[-1]
+      sIn=sOut[-1]
+
+      #store
+      writeable={'t':t,'sol':sol,'se':se,'qt':qt,'sOut':sOut}
+      l=step['label']
+      writeable={os.path.join(jobDir,'{}{}.txt'.format(x,l)):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))
+   print('Completed run')
+   #update setup with new content (particularly t0)
+   setupOut=os.path.join(jobDir,'setup.json')
+   with open(setupOut,'w+') as f:
+      f.write(json.dumps(setup))
+      #this is 3D, so new routines
 
 def interpolate(setup,model,t,s1,t0,tmax):
    #interpolate on s1
@@ -196,6 +244,7 @@ def loadSolutionFromRef(setup):
       
    import labkeyInterface
    import labkeyFileBrowser
+   #import labkeyDatabaseBrowser
 
    net=labkeyInterface.labkeyInterface()
    fconfig=os.path.join(os.path.expanduser('~'),'.labkey',setup['network'])
@@ -206,18 +255,30 @@ def loadSolutionFromRef(setup):
    #db=labkeyDatabaseBrowser.labkeyDB(net)
    fb=labkeyFileBrowser.labkeyFileBrowser(net)
 
-   #refFilter={'variable':'runRef','value':setup['startFromRef'],'oper':'eq'}
-   #ds=db.selectRows(setup['project'],setup['schema'],setup['query'],[refFilter],setup['view']'WithStrings')
-   #ordering requires the latest to be first
-   #ref=ds['rows'][0]['runRef']
+   ref=setup['startFromRef']
+   #refFilter={'variable':'runRef','value':ref,'oper':'eq'}
+   #ds=db.selectRows(setup['project'],setup['schema'],setup['query'],[refFilter],setup['view'])
+
    #parFileString=ds['rows'][0]['parameterFileString']
    #parFiles=parFileString.split(';')
-   ref=setup['startFromRef']
    remoteDir=fb.formatPathURL(setup['project'],'/'.join(['jobs',ref]))
-   localDir=os.path.join(os.path.expanduser('~'),'temp')
-   inFiles=['t.txt','sol.txt','se.txt','qt.txt','sOut.txt']
+   localDir=os.path.join(os.path.expanduser('~'),'temp','jobDir')
+   if not os.path.isdir(localDir):
+      os.mkdir(localDir)
+
+   #download setup
+   f='setup.json'
+   localPath=os.path.join(localDir,f)
+   remotePath='/'.join([remoteDir,f])
+   fb.readFileToFile(remotePath,localPath)
+   setupFile=os.path.join(localDir,f)
+
+   setup=parseSetup(setupFile)
+   inFiles=getLastStepFiles(setup)
+
    #inFiles.extend(parFiles)
-   for f in inFiles:
+   for x in inFiles:
+       f=inFiles[x]
        localPath=os.path.join(localDir,f)
        remotePath='/'.join([remoteDir,f])
        fb.readFileToFile(remotePath,localPath)
@@ -225,11 +286,15 @@ def loadSolutionFromRef(setup):
    return loadSolutionFromDir(localDir)
 
 def loadSolutionFromDir(jobDir):
-   t=numpy.loadtxt(os.path.join(jobDir,'t.txt'))
-   sol=numpy.loadtxt(os.path.join(jobDir,'sol.txt'))
-   se=numpy.loadtxt(os.path.join(jobDir,'se.txt'))
-   sOut,lut,lutSE=read3D(os.path.join(jobDir,'sOut.txt'))
-   qt=numpy.loadtxt(os.path.join(jobDir,'qt.txt'))
+   setupFile=os.path.join(jobDir,'setup.json')
+   setup=parseSetup(setupFile)
+   inFiles=getLastStepFiles(setup)
+
+   t=numpy.loadtxt(os.path.join(jobDir,inFiles['t']))
+   sol=numpy.loadtxt(os.path.join(jobDir,inFiles['sol']))
+   se=numpy.loadtxt(os.path.join(jobDir,inFiles['se']))
+   sOut,lut,lutSE=read3D(os.path.join(jobDir,inFiles['sOut']))
+   qt=numpy.loadtxt(os.path.join(jobDir,inFiles['qt']))
    return t,sol,se,sOut,qt,lut,lutSE
 
 def getStartPointFromDir(jobDir):

setup/setupFast.json

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

setup/setupFastStride.json

@@ -0,0 +1,9 @@
+{"method":"LSODA",
+"atol":1e-8,
+"rtol":1e-8,
+"mode":"IVPSimultaneous",
+"nt":201,
+"tmax":1,
+"tUnit":"hour",
+"stride":{"length":0.5,"tUnit":"hour"}
+}