PBPK_public/pythonScripts/runSolver.py

import ivp
import solveMatrix
import os
import cModel
import sys
import json
import numpy
import time
import scipy.interpolate
import shutil
import importlib
importlib.reload(solveMatrix)

defaultValues={\
      'method':'LSODA',\
      'atol':1e-4,\
      'rtol':1e-4,\
      't0':0,\
      'tmax':1,\
      'mode':'IVPSimultaneous',\
      'nt':201,
      'tUnit':'day',\
      'network':'labkey-klimt.json',\
      'startFromRef':'None',\
      'project':'Analysis/Run',\
      'schema':'lists',\
      'query':'runs',\
      'view':'WithStrings',
      'exposureScaleRef':1}

def get(setup,par):
   try:
      return setup[par]
   except KeyError:
      pass
   return defaultValues[par]


def parseSetup(setupFile):
   with open(setupFile,'r') as f:
      setup=json.load(f)
   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':
      return 60
   if tUnit=='day':
      return 24*60
   if tUnit=='month':
      return 24*60*30
   if tUnit=='year':
      return 24*60*365
   return 1

def getStrides(setup):
   tmax=float(setup['tmax'])*convertUnit(setup['tUnit'])
   t0=setup['t0']
   print('t0 {} tmax {}'.format(t0,tmax))
   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 getFiles(setup, getAll=False):
   strides=getStrides(setup)
   readable=['t','sol','se','qt','sOut']
   #get last stride
   steps=[strides[-1]]
   if getAll:
      steps=strides
   inFiles=[]
   for s in steps:
      l=s['label']
      inFiles.append({x:'{}{}.txt'.format(x,l) for x in readable})
   return inFiles

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']))
      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
   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:
         #this could be modified for different ordering of lut as well
         sIn[:,model.lutSE[x]]=s0[:,startPoint['lutSE'][x]]
   #copy t0 to setup to help get the strides right
   exposureScale=setup['exposureScaleRef']
   y0*=exposureScale
   print('sIn {}'.format(sIn.shape))
   try:
      for x in model.lutSE:
         sIn[:,model.lutSE[x]]*=exposureScale
   except IndexError:
      #if there are no initial conditions, skip this part
      pass
   setup['t0']=t0 
   strides=getStrides(setup)
   full_t=[t0]
   full_sol=[y0]
   full_sOut=[sIn]
   full_qt=[t0]
   for step in strides:
      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(step['label']))

      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,t1,\
            nt=setup['nt'],\
            atol=setup['atol'],\
            rtol=setup['rtol'],method=setup['method'],\
            t0=t0,y0=y0,sIn=sIn)

      if setup['mode']=='solveMatrix':
         t,sol,se,s1=solveMatrix.solveMatrix(model,t1,\
            nt=setup['nt'],\
            t0=t0,y0=y0,sIn=sIn,method=setup['method'])

      end_time=time.time()
      #interpolate on s1
      #s1 has shape ft x n x m
      #where ft is LSODE driven number of points
      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]
      sIn=sOut[-1]

      #store
      writeable={'t':t,'sol':sol,'se':se,'qt':qt,'sOut':sOut}
      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])

      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:
      f.write(json.dumps(setup))
      #this is 3D, so new routines
   print('Written {}'.format(setupOut))
   #write model and parameter file
   origFiles=[modelFile,parameterFile]
   outFiles=["model.json","parameters.json"]
   for x in zip(origFiles,outFiles):
      xOut=os.path.join(jobDir,x[1])
      shutil.copyfile(x[0],xOut)
      print('Written {}'.format(xOut))


def interpolate3d(qt,t,s1):
   #interpolate on s1
   #s1 has shape ft x n x m
   #where ft is LSODE driven number of points
   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 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)
   for x in lut:
      fcode[lut[x]]=x
   return fcode

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 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

def parseShape(line):
   shape=line.replace('#Shape ','').\
      replace('\n','').\
      replace('(','').\
      replace(')','').\
      split(',')
   return [int(x) for x in shape]

      
def read3D(fName):
   new_data=numpy.loadtxt(fName)
#loadtxt will ignore lines starting with #, so that is where we store contextual information like shape
   #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


def loadSolutionFromRef(setup, loadAll=False):
   
   if setup['startFromRef']=='None':
      return {'t':0,'sol':numpy.array([]),'sOut':numpy.array([]),'x':numpy.array([]),'y':numpy.array([])}

   ref=setup['startFromRef']
   return loadSolutionFromRunRef(setup,ref,loadAll)

def connect(setup):
   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
   import labkeyDatabaseBrowser

   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)
   return db,fb

def getRunRef(setup,key):
   #get solution from key rather than runRef (easier to copy, understand)
   db,fb=connect(setup)
   keyFilter={'variable':'Key','value':'{}'.format(key),'oper':'eq'}
   ds=db.selectRows(setup['project'],setup['schema'],setup['query'],[keyFilter],setup['view'])
   try:
      if len(ds['rows'])==0:
         return 'NONE'
   except KeyError:
      print(ds)
      return 'NONE'
   print('Found {} rows for {}'.format(len(ds['rows']),key))
   return ds['rows'][0]['runRef']

def loadSolutionFromRunRef(setup,ref,loadAll=False):

   db,fb=connect(setup)

   remoteDir=fb.formatPathURL(setup['project'],'/'.join(['jobs',ref]))

   try:
      localDir=setup['localDir']
   except KeyError:
      localDir=os.path.join(os.path.expanduser('~'),'temp','jobDir')
   if not os.path.isdir(localDir):
      os.mkdir(localDir)

   #download setup
   fs=['setup.json','model.json','parameters.json']
   for f in fs:
      localPath=os.path.join(localDir,f)
      remotePath='/'.join([remoteDir,f])
      fb.readFileToFile(remotePath,localPath)
   
   setupFile=os.path.join(localDir,fs[0])

   setupOld=parseSetup(setupFile)
   inFiles=getFiles(setupOld,loadAll)

   #inFiles.extend(parFiles)
   i=0
   n=len(inFiles)
   for fSet in inFiles:
      for x in fSet:
         f=fSet[x]
         localPath=os.path.join(localDir,f)
         remotePath='/'.join([remoteDir,f])
         fb.readFileToFile(remotePath,localPath)
      print('Loaded {}/{} {}'.format(i+1,n,localDir))
      i+=1
   return loadSolutionFromDir(localDir,loadAll)

def merge(a1,a2):
   try: 
      return numpy.concatenate(a1,a2)
   except ValueError:
      return a2

def loadSolutionFromDir(jobDir,loadAll=False):
   setupFile=os.path.join(jobDir,'setup.json')
   setup=parseSetup(setupFile)
   inFiles=getFiles(setup,loadAll)
   return loadSolutionFromFiles(setup,jobDir,inFiles)

def loadSolutionFromFiles(setup,jobDir,inFiles):
   isFirst=True 
   operators={x:numpy.loadtxt for x in inFiles[0]}
   operators['sOut']=read3D
   data={}
   i=0
   n=len(inFiles)
   for fSet in inFiles:
      idata={}
      print('Parsing [{}/{}]'.format(i+1,n))
      for f in fSet:
         path=os.path.join(jobDir,fSet[f])
         idata=operators[f](os.path.join(jobDir,fSet[f]))
         if f=='sOut':
            lut=idata[1]
            lutSE=idata[2]
            idata=idata[0]
         try:
            data[f]=numpy.concatenate((data[f],idata))
         except KeyError:
            data[f]=idata
      i+=1
   data['lut']=lut
   data['lutSE']=lutSE
   data['setup']=setup
   data['model']=os.path.join(jobDir,'model.json')
   data['parameters']=os.path.join(jobDir,'parameters.json')
   return data

def getStartPointFromDir(jobDir):
   data=loadSolutionFromDir(jobDir)
   return startPointObject(data)

def startPointObject(data):
   vrs={'t':'t0','sol':'y0','sOut':'s0'}
   out={vrs[x]:data[x][-1] for x in vrs}
   addVrs=['lut','lutSE']
   for x in addVrs:
      out[x]=data[x]
   return out

def getStartPoint(setup):
   if setup['startFromRef']=='None':
      return {'t0':0,'y0':numpy.array([]),'s0':numpy.array([]),'lut':[],'lutSE':[]}
   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]
    main(parFiles,jobDir)