runSolver.py 6.5 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232
  1. import ivp
  2. import os
  3. import cModel
  4. import sys
  5. import json
  6. import numpy
  7. import time
  8. import scipy.interpolate
  9. defaultValues={\
  10. 'method':'LSODA',\
  11. 'atol':1e-4,\
  12. 'rtol':1e-4,\
  13. 't0':0,\
  14. 'tmax':1,\
  15. 'mode':'IVPSimultaneous',\
  16. 'nt':201,
  17. 'tUnit':'day',\
  18. 'network':'labkey-klimt.json',\
  19. 'startFromRef':'None',\
  20. 'project':'Analysis/Run',\
  21. 'schema':'lists',\
  22. 'query':'runs',\
  23. 'view':'withStrings'}
  24. def get(setup,par):
  25. try:
  26. return setup[par]
  27. except KeyError:
  28. pass
  29. return defaultValues[par]
  30. def parseSetup(setupFile):
  31. with open(setupFile,'r') as f:
  32. setup=json.load(f)
  33. out={}
  34. for p in defaultValues:
  35. out[p]=get(setup,p)
  36. return out
  37. def getScale(setup):
  38. tUnit=setup['tUnit']
  39. if tUnit=='min':
  40. return 1
  41. if tUnit=='hour':
  42. return 60
  43. if tUnit=='day':
  44. return 24*60
  45. if tUnit=='month':
  46. return 24*60*30
  47. if tUnit=='year':
  48. return 24*60*30*365
  49. return 1
  50. def main(parFiles,jobDir):
  51. model=cModel.model()
  52. setupFile=parFiles[0]
  53. modelFile=parFiles[1]
  54. parameterFile=parFiles[2]
  55. model.parse(modelFile,parameterFile)
  56. setup=parseSetup(setupFile)
  57. scale=getScale(setup)
  58. tmax=setup['tmax']*scale
  59. print('Using {}'.format(jobDir))
  60. t0,y0,S1,lut,lutSE=getStartPoint(setup)
  61. #we should reorder S1 to match model.lutSE
  62. #S1 is nxm, so we should reorder columns
  63. Sin=numpy.zeros(S1.shape)
  64. if S1.size>0:
  65. for x in model.lutSE:
  66. #this could be modified for different ordering of lut as well
  67. Sin[:,model.lutSE[x]]=S1[:,lutSE[x]]
  68. start_time=time.time()
  69. if setup['mode']=='SequentialOdeint':
  70. t,sol,se,s1=ivp.solveSequentialOdeint(model,tmax,setup['nt'],t0,y0,Sin)
  71. if setup['mode']=='SimultaneousOdeint':
  72. t,sol,se,s1=ivp.solveSimultaneousOdeint(model,tmax,setup['nt'],t0,y0,Sin)
  73. if setup['mode']=='IVP':
  74. t,sol,se,s1=ivp.solveSequential(model,tmax,atol=setup['atol'],\
  75. rtol=setup['rtol'],method=setup['method'],\
  76. t0=t0,y0=y0,Sin=Sin)
  77. if setup['mode']=='IVPSimultaneous':
  78. t,sol,se,s1=ivp.solveSimultaneous(model,tmax,atol=setup['atol'],\
  79. rtol=setup['rtol'],method=setup['method'],\
  80. t0=t0,y0=y0,Sin=Sin)
  81. end_time=time.time()
  82. print('Time: {:.3f} s'.format(end_time-start_time))
  83. #store
  84. #interpolate on s1
  85. #s1 has shape ft x n x m
  86. #where ft is LSODE driven number of points
  87. qt,sOut=interpolate(setup,model,t,s1,t0,tmax)
  88. numpy.savetxt(os.path.join(jobDir,'t.txt'),t)
  89. numpy.savetxt(os.path.join(jobDir,'sol.txt'),sol)
  90. numpy.savetxt(os.path.join(jobDir,'se.txt'),se)
  91. numpy.savetxt(os.path.join(jobDir,'qt.txt'),qt)
  92. #this is 3D, so new routines
  93. write3D(os.path.join(jobDir,'sOut.txt'),sOut,model.lut,model.lutSE)
  94. def interpolate(setup,model,t,s1,t0,tmax):
  95. #interpolate on s1
  96. #s1 has shape ft x n x m
  97. #where ft is LSODE driven number of points
  98. sOut=numpy.zeros((setup['nt'],model.n,model.m))
  99. qt=numpy.linspace(t0,tmax,setup['nt'])
  100. for i in range(model.n):
  101. for j in range(model.m):
  102. y=s1[:,i,j]
  103. f = scipy.interpolate.interp1d(t, y, kind='cubic')
  104. sOut[:,i,j]=f(qt)
  105. return qt,sOut
  106. def invertLUT(lut):
  107. fcode=[None]*len(lut)
  108. for x in lut:
  109. fcode[lut[x]]=x
  110. return fcode
  111. def write3D(fName,a,lut,lutSE):
  112. with open(fName,'w') as f:
  113. f.write('#Shape {}\n'.format(a.shape))
  114. f.write('#Lut {}\n'.format(invertLUT(lut)))
  115. f.write('#Lut (SE) {}\n'.format(invertLUT(lutSE)))
  116. i=0
  117. for data_slice in a:
  118. f.write('#\t Slice {}\n'.format(i))
  119. numpy.savetxt(f,data_slice)
  120. i+=1
  121. def parseLUT(line,isSE=False):
  122. SE=''
  123. if isSE:
  124. SE='(SE)'
  125. header='#Lut {}'.format(SE)
  126. names=line.replace(header,'').\
  127. replace('[','').\
  128. replace(']','').\
  129. replace(' ','').\
  130. replace('\n','').\
  131. replace("'",'').\
  132. split(',')
  133. lut={names[i]:i for i in range(len(names))}
  134. return lut
  135. def parseShape(line):
  136. shape=line.replace('#Shape ','').\
  137. replace('\n','').\
  138. replace('(','').\
  139. replace(')','').\
  140. split(',')
  141. return [int(x) for x in shape]
  142. def read3D(fName):
  143. new_data=numpy.loadtxt(fName)
  144. #read and parse first three lines to get shape,lut and lutSE
  145. nLines=3
  146. with open(fName,'r') as f:
  147. lines=[f.readline() for i in range(nLines)]
  148. shape=parseShape(lines[0])
  149. lut=parseLUT(lines[1])
  150. lutSE=parseLUT(lines[2],isSE=True)
  151. return new_data.reshape(shape),lut,lutSE
  152. def loadSolutionFromRef(setup):
  153. if setup['startFromRef']=='None':
  154. return 0,numpy.array([]),numpy.array([]),numpy.array([]),numpy.array([])
  155. nixSuite=os.path.join(os.path.expanduser('~'),'software','src','nixSuite')
  156. sys.path.append(os.path.join(nixSuite,'wrapper'))
  157. import nixWrapper
  158. nixWrapper.loadLibrary('labkeyInterface')
  159. import labkeyInterface
  160. import labkeyFileBrowser
  161. net=labkeyInterface.labkeyInterface()
  162. fconfig=os.path.join(os.path.expanduser('~'),'.labkey',setup['network'])
  163. net.init(fconfig)
  164. #debug
  165. net.getCSRF()
  166. #db=labkeyDatabaseBrowser.labkeyDB(net)
  167. fb=labkeyFileBrowser.labkeyFileBrowser(net)
  168. #refFilter={'variable':'runRef','value':setup['startFromRef'],'oper':'eq'}
  169. #ds=db.selectRows(setup['project'],setup['schema'],setup['query'],[refFilter],setup['view']'WithStrings')
  170. #ordering requires the latest to be first
  171. #ref=ds['rows'][0]['runRef']
  172. #parFileString=ds['rows'][0]['parameterFileString']
  173. #parFiles=parFileString.split(';')
  174. ref=setup['startFromRef']
  175. remoteDir=fb.formatPathURL(setup['project'],'/'.join(['jobs',ref]))
  176. localDir=os.path.join(os.path.expanduser('~'),'temp')
  177. inFiles=['t.txt','sol.txt','se.txt','qt.txt','sOut.txt']
  178. #inFiles.extend(parFiles)
  179. for f in inFiles:
  180. localPath=os.path.join(localDir,f)
  181. remotePath='/'.join([remoteDir,f])
  182. fb.readFileToFile(remotePath,localPath)
  183. t=numpy.loadtxt(os.path.join(localDir,'t.txt'))
  184. sol=numpy.loadtxt(os.path.join(localDir,'sol.txt'))
  185. se=numpy.loadtxt(os.path.join(localDir,'se.txt'))
  186. sOut,lut,lutSE=read3D(os.path.join(localDir,'sOut.txt'))
  187. qt=numpy.loadtxt(os.path.join(localDir,'qt.txt'))
  188. return t,sol,se,sOut,qt,lut,lutSE
  189. def getStartPoint(setup):
  190. if setup['startFromRef']=='None':
  191. return 0,numpy.array([]),numpy.array([]),[],[]
  192. t,sol,se,sOut,qt,lut,lutSE=loadSolutionFromRef(setup)
  193. return t[-1],sol[-1],sOut[-1],lut,lutSE
  194. if __name__=="__main__":
  195. parFiles=sys.argv[1].split(';')
  196. jobDir=sys.argv[2]
  197. main(parFiles,jobDir)