import ivp import os import cModel import sys import json import numpy import time import scipy.interpolate 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'} 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) return out def getScale(setup): tUnit=setup['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*30*365 return 1 def main(parFiles,jobDir): model=cModel.model() setupFile=parFiles[0] modelFile=parFiles[1] parameterFile=parFiles[2] model.parse(modelFile,parameterFile) 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'],t0,y0,S1) if setup['mode']=='SimultaneousOdeint': 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'],\ t0,y0,S1) if setup['mode']=='IVPSimultaneous': 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)) #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,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) def interpolate(setup,model,t,s1,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']) for i in range(model.n): for j in range(model.m): y=s1[:,i,j] f = scipy.interpolate.interp1d(t, y, kind='cubic') sOut[:,i,j]=f(qt) return qt,sOut def write3D(fName,a): with open(fName,'w') as f: f.write('#Shape {}\n'.format(a.shape)) 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() shape=firstLine.replace('#Shape ','').\ replace('\n','').\ replace('(','').\ replace(')','').\ split(',') shape=[int(x) for x in shape] return new_data.reshape(shape) if __name__=="__main__": parFiles=sys.argv[1].split(';') 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]