Solution reuse infrastructure (to continue on a time curve)

pythonScripts/ivp.py

@@ -59,74 +59,82 @@ def jacobiSEFull(t,S,system):
         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)
+def solveSimultaneous(model,tmax,atol,rtol,method='LSODA',t0=0,y0=None, S1=None):
+    if S1==None:
+       S1=numpy.zeros((model.n,model.m+1))
+    if y0==None:
+       y0=numpy.zeros(model.n)
     #set initial condition
     S1[:,0]=y0
     S1=S1.ravel()
-    sol=scipy.integrate.solve_ivp(dfdySFull,[0, tmax],S1, args=(sys,), jac=jacobiSEFull,
+    sol=scipy.integrate.solve_ivp(dfdySFull,[t0, tmax],S1, args=(model,), 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))
+    sFull=numpy.reshape(numpy.transpose(sol.y),(len(t),model.n,model.m+1))
     s1=sFull[:,:,1:]
     ysol=sFull[:,:,0]
-    se=sys.calculateUncertainty(ysol,s1)
+    se=model.calculateUncertainty(ysol,s1)
     print('Done simultaneous LSODA SE')
     return t,ysol,se,s1
 
-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,
+def solveSequential(model,tmax,atol,rtol,method='LSODA',t0=0,y0=None, S1=None):
+   if y0==None:
+      y0=numpy.zeros(model.n)
+   solIVP=scipy.integrate.solve_ivp(dfdy,[t0, tmax],y0, args=(model,), jac=jacobi,
                                   method=method, atol=atol, rtol=rtol)
-#y is n x nt (odeint nt x n)
+   #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()
+   model.setY(t,sol)
+   if S1==None:
+      S1=numpy.zeros((model.n,model.m))
+   S1=S1.ravel()
     
-   solIVPSE=scipy.integrate.solve_ivp(dfdyS,[0, tmax],S0, args=(sys,), jac=jacobiSE,
+   solIVPSE=scipy.integrate.solve_ivp(dfdyS,[0, tmax],S1, args=(model,), 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):
+   sraw=numpy.reshape(numpy.transpose(solIVPSE.y),(len(solIVPSE.t),model.n,model.m))
+   #interpolate on t
+   s1=numpy.zeros((len(t),model.n,model.m))
+   for i in range(model.n):
+      for j in range(model.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)
+   se=model.calculateUncertainty(sol,s1)
    return t,sol,se,s1
 
-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
+def solveSimultaneousOdeint(model,tmax,nt=201,t0=0,y0=None, S1=None):
+   t = numpy.linspace(t0,tmax, nt)
+   if y0==None:
+      y0=numpy.zeros(model.n)
+   if S1==None:
+      S1=numpy.zeros((model.n,model.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))
+   solSE1=scipy.integrate.odeint(dfdySFull, S1, t, args=(model,),Dfun=jacobiSEFull,tfirst=True)
+   sFull=numpy.reshape(solSE1,(len(t),model.n,model.m+1))
    s1=sFull[:,:,1:]
    sol=sFull[:,:,0]
    se=sys.calculateUncertainty(sol,s1)
    print('Done simultaneous SE')
    return t,sol,se,s1
 
-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)
+def solveSequentialOdeint(model,tmax,nt=201,t0=0,y0=None,S1=None):
+   t = numpy.linspace(t0,tmax, nt)
+   if y0==None:
+      y0=numpy.zeros(sys.n)
+   sol = scipy.integrate.odeint(dfdy, y0=y0, t=t, args=(model,),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)
+   model.setY(t,sol)
+   if S1==None:
+      S1=numpy.zeros((model.n,model.m))
+   S1=S1.ravel()
+   solSE=scipy.integrate.odeint(dfdyS, S1, t, args=(model,),Dfun=jacobiSE,tfirst=True)
+   s1=numpy.reshape(solSE,(len(t),model.n,model.m))
+   se=model.calculateUncertainty(sol,s1)
    print('Done sequential SE')
    return t,sol,se,s1
 

pythonScripts/runSolver.py

@@ -11,10 +11,17 @@ defaultValues={\
       'method':'LSODA',\
       'atol':1e-4,\
       'rtol':1e-4,\
+      't0':0,\
       'tmax':1,\
       'mode':'IVPSimultaneous',\
       'nt':201,
-      'tUnit':'day'}
+      'tUnit':'day',\
+      'network':'labkey-klimt.json',\
+      'startFromRef':'None',\
+      'project':'Analysis/Run',\
+      'schema':'lists',\
+      'query':'runs',\
+      'view':'withStrings'}
 
 def get(setup,par):
    try:
@@ -57,20 +64,23 @@ def main(parFiles,jobDir):
    setup=parseSetup(setupFile)
    scale=getScale(setup)
    tmax=setup['tmax']*scale
+   t0,y0,S1=loadSolutionFromRef(setup)
 
    start_time=time.time()
    if setup['mode']=='SequentialOdeint':
-      t,sol,se,s1=ivp.solveSequentialOdeint(model,tmax,setup['nt'])
+      t,sol,se,s1=ivp.solveSequentialOdeint(model,tmax,setup['nt'],t0,y0,S1)
    if setup['mode']=='SimultaneousOdeint':
-      t,sol,se,s1=ivp.solveSimultaneousOdeint(model,tmax,setup['nt'])
+      t,sol,se,s1=ivp.solveSimultaneousOdeint(model,tmax,setup['nt'],t0,y0,S1)
 
    if setup['mode']=='IVP':
-      t,sol,se,s1=ivp.solveSequential(model,tmax,atol=setup['atol'],
-         rtol=setup['rtol'],method=setup['method'])
+      t,sol,se,s1=ivp.solveSequential(model,tmax,atol=setup['atol'],\
+         rtol=setup['rtol'],method=setup['method'],\
+         t0,y0,S1)
         
    if setup['mode']=='IVPSimultaneous':
-      t,sol,se,s1=ivp.solveSimultaneous(model,tmax,atol=setup['atol'],
-         rtol=setup['rtol'],method=setup['method'])
+      t,sol,se,s1=ivp.solveSimultaneous(model,tmax,atol=setup['atol'],\
+         rtol=setup['rtol'],method=setup['method'],\
+         t0,y0,S1)
 
    end_time=time.time()
    print('Time: {:.3f} s'.format(end_time-start_time))
@@ -130,4 +140,46 @@ if __name__=="__main__":
     jobDir=sys.argv[2]
     main(parFiles,jobDir)
 
+def loadSolutionFromRef(setup):
+   if setup['startFromRef']=='None':
+      return 0,None,None
 
+   nixSuite=os.path.join(os.path.expanduser('~'),'software','src','nixSuite')
+   sys.path.append(os.path.join(nixSuite,'wrapper'))
+   import nixWrapper
+   nixWrapper.loadLibrary('labkeyInterface')
+      
+   import labkeyInterface
+   import labkeyFileBrowser
+
+   net=labkeyInterface.labkeyInterface()
+   fconfig=os.path.join(os.path.expanduser('~'),'.labkey',setup['network'])
+   net.init(fconfig)
+   #debug
+   net.getCSRF()
+
+   #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']
+   #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']
+   #inFiles.extend(parFiles)
+   for f in inFiles:
+       localPath=os.path.join(localDir,f)
+       remotePath='/'.join([remoteDir,f])
+       fb.readFileToFile(remotePath,localPath)
+
+   t=numpy.loadtxt(os.path.join(localDir,'t.txt'))
+   sol=numpy.loadtxt(os.path.join(localDir,'sol.txt'))
+   se=numpy.loadtxt(os.path.join(localDir,'se.txt'))
+   sOut=read3D(os.path.join(localDir,'sOut.txt'))
+   qt=numpy.loadtxt(os.path.join(localDir,'qt.txt'))
+   return t[-1],sol[-1],sOut[-1]