Making output reloadable; figuring out variable ordering

pythonScripts/ivp.py

@@ -59,10 +59,11 @@ 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(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:
+def solveSimultaneous(model,tmax,atol,rtol,method='LSODA',t0=0,y0=numpy.array([]), Sin=numpy.array([])):
+    S1=numpy.zeros((model.n,model.m+1))
+    if Sin.size!=0:
+       S1[:,1:]=Sin
+    if y0.size==0:
        y0=numpy.zeros(model.n)
     #set initial condition
     S1[:,0]=y0
@@ -77,8 +78,8 @@ def solveSimultaneous(model,tmax,atol,rtol,method='LSODA',t0=0,y0=None, S1=None)
     print('Done simultaneous LSODA SE')
     return t,ysol,se,s1
 
-def solveSequential(model,tmax,atol,rtol,method='LSODA',t0=0,y0=None, S1=None):
-   if y0==None:
+def solveSequential(model,tmax,atol,rtol,method='LSODA',t0=0,y0=numpy.array([]), S1=numpy.array([])):
+   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)
@@ -87,7 +88,7 @@ def solveSequential(model,tmax,atol,rtol,method='LSODA',t0=0,y0=None, S1=None):
    t=solIVP.t
    print('shape (y) {}'.format(sol.shape))
    model.setY(t,sol)
-   if S1==None:
+   if S1.size==0:
       S1=numpy.zeros((model.n,model.m))
    S1=S1.ravel()
     
@@ -104,12 +105,14 @@ def solveSequential(model,tmax,atol,rtol,method='LSODA',t0=0,y0=None, S1=None):
    se=model.calculateUncertainty(sol,s1)
    return t,sol,se,s1
 
-def solveSimultaneousOdeint(model,tmax,nt=201,t0=0,y0=None, S1=None):
+def solveSimultaneousOdeint(model,tmax,nt=201,t0=0,y0=numpy.array([]), Sin=numpy.array([])):
    t = numpy.linspace(t0,tmax, nt)
-   if y0==None:
+   if y0.size==0:
       y0=numpy.zeros(model.n)
-   if S1==None:
-      S1=numpy.zeros((model.n,model.m+1))
+   
+   S1=numpy.zeros((model.n,model.m+1))
+   if Sin.size!=0:
+      S1[:,1:]=Sin
    #set initial condition
    S1[:,0]=y0
    S1=S1.ravel()
@@ -121,15 +124,15 @@ def solveSimultaneousOdeint(model,tmax,nt=201,t0=0,y0=None, S1=None):
    print('Done simultaneous SE')
    return t,sol,se,s1
 
-def solveSequentialOdeint(model,tmax,nt=201,t0=0,y0=None,S1=None):
+def solveSequentialOdeint(model,tmax,nt=201,t0=0,y0=numpy.array([]),S1=numpy.array([])):
    t = numpy.linspace(t0,tmax, nt)
-   if y0==None:
+   if y0.size==0:
       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))
 
    model.setY(t,sol)
-   if S1==None:
+   if S1.size==0:
       S1=numpy.zeros((model.n,model.m))
    S1=S1.ravel()
    solSE=scipy.integrate.odeint(dfdyS, S1, t, args=(model,),Dfun=jacobiSE,tfirst=True)

pythonScripts/runSolver.py

@@ -64,23 +64,30 @@ def main(parFiles,jobDir):
    setup=parseSetup(setupFile)
    scale=getScale(setup)
    tmax=setup['tmax']*scale
-   t0,y0,S1=loadSolutionFromRef(setup)
+   t0,y0,S1,lut,lutSE=getStartPoint(setup)
+   #we should reorder S1 to match model.lutSE
+   #S1 is nxm, so we should reorder columns
+   Sin=numpy.zeros(S1.shape)
+   if S1.size>0:
+      for x in model.lutSE:
+         #this could be modified for different ordering of lut as well
+         Sin[:,model.lutSE[x]]=S1[:,lutSE[x]]
 
    start_time=time.time()
    if setup['mode']=='SequentialOdeint':
-      t,sol,se,s1=ivp.solveSequentialOdeint(model,tmax,setup['nt'],t0,y0,S1)
+      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,S1)
+      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,atol=setup['atol'],\
          rtol=setup['rtol'],method=setup['method'],\
-         t0,y0,S1)
+         t0=t0,y0=y0,Sin=Sin)
         
    if setup['mode']=='IVPSimultaneous':
       t,sol,se,s1=ivp.solveSimultaneous(model,tmax,atol=setup['atol'],\
          rtol=setup['rtol'],method=setup['method'],\
-         t0,y0,S1)
+         t0=t0,y0=y0,Sin=Sin)
 
    end_time=time.time()
    print('Time: {:.3f} s'.format(end_time-start_time))
@@ -88,21 +95,21 @@ def main(parFiles,jobDir):
    #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,tmax)
+   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)
+   write3D(os.path.join(jobDir,'sOut.txt'),sOut,model.lut,model.lutSE)
 
-def interpolate(setup,model,t,s1,tmax):
+def interpolate(setup,model,t,s1,t0,tmax):
    #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(0,tmax,setup['nt'])
+   qt=numpy.linspace(t0,tmax,setup['nt'])
    for i in range(model.n):
       for j in range(model.m):
          y=s1[:,i,j]
@@ -110,30 +117,60 @@ def interpolate(setup,model,t,s1,tmax):
          sOut[:,i,j]=f(qt)
    return qt,sOut
 
+def invertLUT(lut):
+   fcode=[None]*len(lut)
+   for x in lut:
+      fcode[lut[x]]=x
+   return fcode
 
-
-def write3D(fName,a):
+def write3D(fName,a,lut,lutSE):
    with open(fName,'w') as f:
       f.write('#Shape {}\n'.format(a.shape))
+      f.write('#Lut {}\n'.format(invertLUT(lut)))
+      f.write('#Lut (SE) {}\n'.format(invertLUT(lutSE)))
+
       i=0
       for data_slice in a:
          f.write('#\t Slice {}\n'.format(i))
          numpy.savetxt(f,data_slice)
          i+=1
       
-def read3D(fName):
-   new_data=numpy.loadtxt(fName)
-   #read and parse first row to get shape
-   with open(fName,'r') as f:
-      firstLine=f.readline()
+def parseLUT(line,isSE=False):
+   SE=''
+   if isSE:
+      SE='(SE)'
+   header='#Lut {}'.format(SE)
+   names=line.replace(header,'').\
+      replace('[','').\
+      replace(']','').\
+      replace(' ','').\
+      replace('\n','').\
+      replace("'",'').\
+      split(',')
+   lut={names[i]:i for i in range(len(names))}
+   return lut
 
-   shape=firstLine.replace('#Shape ','').\
+def parseShape(line):
+   shape=line.replace('#Shape ','').\
       replace('\n','').\
       replace('(','').\
       replace(')','').\
       split(',')
-   shape=[int(x) for x in shape]
-   return new_data.reshape(shape)
+   return [int(x) for x in shape]
+
+      
+def read3D(fName):
+   new_data=numpy.loadtxt(fName)
+   #read and parse first three lines to get shape,lut and lutSE
+   nLines=3
+   with open(fName,'r') as f:
+      lines=[f.readline() for i in range(nLines)]
+
+   shape=parseShape(lines[0])
+   lut=parseLUT(lines[1])
+   lutSE=parseLUT(lines[2],isSE=True)
+
+   return new_data.reshape(shape),lut,lutSE
 
 if __name__=="__main__":
     parFiles=sys.argv[1].split(';')
@@ -141,8 +178,9 @@ if __name__=="__main__":
     main(parFiles,jobDir)
 
 def loadSolutionFromRef(setup):
+   
    if setup['startFromRef']=='None':
-      return 0,None,None
+      return 0,numpy.array([]),numpy.array([]),numpy.array([]),numpy.array([])
 
    nixSuite=os.path.join(os.path.expanduser('~'),'software','src','nixSuite')
    sys.path.append(os.path.join(nixSuite,'wrapper'))
@@ -180,6 +218,13 @@ def loadSolutionFromRef(setup):
    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'))
+   sOut,lut,lutSE=read3D(os.path.join(localDir,'sOut.txt'))
    qt=numpy.loadtxt(os.path.join(localDir,'qt.txt'))
-   return t[-1],sol[-1],sOut[-1]
+   return t,sol,se,sOut,qt,lut,lutSE
+
+
+def getStartPoint(setup):
+   if setup['startFromRef']=='None':
+      return 0,numpy.array([]),numpy.array([]),[],[]
+   t,sol,se,sOut,qt,lut,lutSE=loadSolutionFromRef(setup)
+   return t[-1],sol[-1],sOut[-1],lut,lutSE

setup/setup.json

@@ -0,0 +1,7 @@
+{"method":"LSODA",
+"atol":1e-4,
+"rtol":1e-4,
+"mode":"IVPSimultaneous",
+"nt":201,
+"tmax":1,
+"tUnit":"day"}

setup/setupFast.json

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

setup/setupFastLoad.json

@@ -0,0 +1,9 @@
+{"method":"LSODA",
+"atol":1e-4,
+"rtol":1e-4,
+"mode":"IVPSimultaneous",
+"nt":201,
+"tmax":1.05,
+"tUnit":"day",
+"startFromRef":"1668337037"
+}