cModel.py 23 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761
  1. import numpy
  2. import json
  3. import os
  4. import scipy.interpolate
  5. #for partial function specializations
  6. import functools
  7. import function
  8. import importlib
  9. importlib.reload(function)
  10. class model:
  11. def __init__(self):
  12. self.compartments={}
  13. self.seJ={}
  14. self.scaled=[]
  15. def add_source(self,compartmentName,formula):
  16. self.compartments[compartmentName]['source']=formula
  17. def add_compartment(self,compartmentName):
  18. self.compartments[compartmentName]={}
  19. self.compartments[compartmentName]['targets']={}
  20. self.compartments[compartmentName]['sensTargets']={}
  21. def getTimeUnit(self):
  22. try:
  23. return self.mod['timeUnit']
  24. except KeyError:
  25. return 's'
  26. def bind(self,src,target,qName,pcName):
  27. #establish a flow from source compartment to the target
  28. #the source equation (where we subtract the current)
  29. #in fact, this is the diagonal element
  30. #get volume names
  31. srcVName=self.getVolumePar(src)
  32. #generate coupling object (w/derivatives)
  33. pSrc=self.couplingObject(-1,qName,pcName,srcVName)
  34. #this includes derivatives and value!
  35. self.addValueObject(src,src,pSrc)
  36. #special target which is not part of calculation
  37. if target=='dump':
  38. return
  39. #the target equation (where we add the current)
  40. #get volume names
  41. targetVName=self.getVolumePar(target)
  42. #generate coupling object
  43. pTarget=self.couplingObject(1,qName,pcName,targetVName)
  44. #equation is for target compartment, but binding for source
  45. self.addValueObject(target,src,pTarget)
  46. def addValueObject(self,targetName,srcName,cObject):
  47. #always binds equation id and a variable
  48. targetList=self.compartments[targetName]['targets']
  49. addValue(targetList,srcName,cObject["value"])
  50. der=cObject["derivatives"]
  51. for d in der:
  52. targetSE=self.getSEJ_comp(d,targetName)
  53. addValue(targetSE,srcName,der[d])
  54. def couplingObject(self,sign, qParName, pcParName, vParName):
  55. qPar=self.get(qParName)
  56. pcPar=self.get(pcParName)
  57. vPar=self.get(vParName)
  58. q=qPar['value']
  59. pc=pcPar['value']
  60. v=vPar['value']
  61. DPC=pcPar['derivatives']
  62. DQ=qPar['derivatives']
  63. DV=vPar['derivatives']
  64. if any(['function' in qPar,'function' in pcPar, 'function' in vPar]):
  65. fq=function.to(q)
  66. fpc=function.to(pc)
  67. fv=function.to(v)
  68. f=lambda t,q=fq,pc=fpc,v=fv,s=sign:s*q(t)/v(t)/pc(t)
  69. dfdPC=lambda t,f=f,pc=fpc:-f(t)/pc(t)
  70. dPC=function.generate(dfdPC,DPC)
  71. dfdQ=lambda t,f=f,q=fq: f(t)/q(t)
  72. dQ=function.generate(dfdQ,DQ)
  73. dfdV=lambda t,f=f,v=fv: -f(t)/v(t)
  74. dV=function.generate(dfdV,DV)
  75. return function.Object(f,[dPC,dQ,dV])
  76. else:
  77. f=sign*q/v/pc
  78. return derivedObject(sign*q/v/pc,\
  79. [{'df':-f/pc,'D':DPC},\
  80. {'df':sign/v/pc,'D':DQ},\
  81. {'df':-f/v,'D':DV}])
  82. #derivatives is the combination of the above
  83. def getVolumePar(self,compartment):
  84. #returnis volume name, if found and useVolume is directed,
  85. #or a standard parameter one
  86. try:
  87. return self.mod["volumes"][compartment]
  88. #parV=self.mod["parameters"][parVName]
  89. except KeyError:
  90. pass
  91. return "one"
  92. def build(self):
  93. comps=self.compartments
  94. self.n=len(comps)
  95. self.fu=[lambda t:0]*self.n
  96. self.lut={c:i for (i,c) in zip(range(self.n),comps.keys())}
  97. self.dM={}
  98. self.fM=numpy.zeros((self.n,self.n))
  99. for c in comps:
  100. comp=comps[c]
  101. if 'source' in comp:
  102. self.fu[self.lut[c]]=parseFunction(comp['source'])
  103. for t in comp['targets']:
  104. arr=comp['targets'][t]
  105. if function.contains(arr):
  106. try:
  107. self.dM[self.lut[c]][self.lut[t]]=\
  108. function.sumArray(arr)
  109. except (KeyError,TypeError):
  110. self.dM[self.lut[c]]={}
  111. self.dM[self.lut[c]][self.lut[t]]=\
  112. function.sumArray(arr)
  113. else:
  114. #just set once
  115. self.fM[self.lut[c],self.lut[t]]=sum(arr)
  116. #build SE part
  117. self.buildSE()
  118. def buildSE(self):
  119. #check which parameterst to include
  120. parList=[]
  121. pars=self.parSetup['parameters']
  122. for parName in self.seJ:
  123. #print(par)
  124. par=pars[parName]
  125. usePar=calculateDerivative(par)
  126. #print('[{}]: {}'.format(usePar,par))
  127. if not usePar:
  128. continue
  129. parList.append(parName)
  130. #print(parList)
  131. self.m=len(parList)
  132. self.lutSE={c:i for (i,c) in zip(range(self.m),parList)}
  133. w=self.getWeights(self.lutSE)
  134. w=numpy.sqrt(w)
  135. self.qSS={}
  136. self.SS=numpy.zeros((self.m,self.n,self.n))
  137. for parName in parList:
  138. sources=self.seJ[parName]
  139. for compartment in sources:
  140. targets=sources[compartment]
  141. for t in targets:
  142. k=self.lutSE[parName]
  143. i=self.lut[compartment]
  144. j=self.lut[t]
  145. #print('[{} {} {}] {}'.format(parName,compartment,t,targets[t]))
  146. arr=targets[t]
  147. if not function.contains(arr):
  148. self.SS[k,i,j]=w[k]*sum(arr)
  149. continue
  150. ft=function.sumArray(arr,w[k])
  151. try:
  152. self.qSS[k][i][j]=ft
  153. except (KeyError,TypeError):
  154. try:
  155. self.qSS[k][i]={}
  156. self.qSS[k][i][j]=ft
  157. except (KeyError,TypeError):
  158. self.qSS[k]={}
  159. self.qSS[k][i]={}
  160. self.qSS[k][i][j]=ft
  161. #use fM to build static part of fJ
  162. N=self.n*(self.m+1)
  163. self.fJ=numpy.zeros((N,N))
  164. for i in range(self.m+1):
  165. self.fJ[i*self.n:(i+1)*self.n,i*self.n:(i+1)*self.n]=self.fM
  166. def inspect(self):
  167. comps=self.compartments
  168. pars=self.parSetup['parameters']
  169. #pars=self.mod['parameters']
  170. tu=self.getTimeUnit()
  171. print('Time unit: {}'.format(tu))
  172. print('Compartments')
  173. for c in comps:
  174. print('{}/{}:'.format(c,self.lut[c]))
  175. comp=comps[c]
  176. if 'source' in comp:
  177. print('\tsource\n\t\t{}'.format(comp['source']))
  178. print('\ttargets')
  179. for t in comp['targets']:
  180. print('\t\t{}[{},{}]: {}'.format(t,self.lut[c],self.lut[t],\
  181. comp['targets'][t]))
  182. print('Flows')
  183. for f in self.flows:
  184. fName=self.flows[f]
  185. fParName=self.mod['flows'][fName]
  186. fPar=pars[fParName]
  187. print('\t{}[{}]:{} [{}]'.format(f,fName,fParName,self.get(fParName)))
  188. print('Volumes')
  189. for v in self.mod['volumes']:
  190. vParName=self.mod['volumes'][v]
  191. vPar=pars[vParName]
  192. print('\t{}:{} [{}]'.format(v,vParName,self.get(vParName)))
  193. print('Partition coefficients')
  194. for pc in self.mod['partitionCoefficients']:
  195. pcParName=self.mod['partitionCoefficients'][pc]
  196. pcPar=pars[pcParName]
  197. print('\t{}:{} [{}]'.format(pc,pcParName,self.get(pcParName)))
  198. def inspectSE(self):
  199. print('SE parameters')
  200. for p in self.seJ:
  201. print(p)
  202. sources=self.seJ[p]
  203. for compartment in sources:
  204. targets=sources[compartment]
  205. for t in targets:
  206. print('\t SE bind {}/{}:{}'.format(compartment,t,targets[t]))
  207. def parse(self,setupFile,parameterFile):
  208. with open(setupFile,'r') as f:
  209. self.mod=json.load(f)
  210. with open(parameterFile,'r') as f:
  211. self.parSetup=json.load(f)
  212. self.mod['compartments'].append('total')
  213. for m in self.mod['compartments']:
  214. self.add_compartment(m)
  215. for m in self.mod['scaled']:
  216. self.scaled.append(m)
  217. self.add_default_parameters()
  218. #standard parameters such as one,zero etc.
  219. for s in self.mod['sources']:
  220. #src=mod['sources'][s]
  221. self.add_source(s,self.mod['sources'][s])
  222. self.flows={}
  223. #pars=self.mod['parameters']
  224. pars=self.parSetup['parameters']
  225. for f in self.mod['flows']:
  226. #skip comments
  227. if f.find(':')<0:
  228. continue
  229. comps=f.split(':')
  230. c0=splitVector(comps[0])
  231. c1=splitVector(comps[1])
  232. for x in c0:
  233. for y in c1:
  234. pairName='{}:{}'.format(x,y)
  235. self.flows[pairName]=f
  236. for b in self.mod['bindings']['diffusion']:
  237. #whether to scale transfer constants to organ volume
  238. #default is true, but changing here will assume no scaling
  239. comps=b.split('->')
  240. try:
  241. pcParName=self.mod['partitionCoefficients'][b]
  242. except KeyError:
  243. pcParName="one"
  244. kParName=self.mod['bindings']['diffusion'][b]
  245. #operate with names to allow for value/function/derived infrastructure
  246. self.bind(comps[0],comps[1],kParName,pcParName)
  247. for q in self.mod['bindings']['flow']:
  248. comps=q.split('->')
  249. srcs=splitVector(comps[0])
  250. tgts=splitVector(comps[1])
  251. for cs in srcs:
  252. for ct in tgts:
  253. #get partition coefficient
  254. try:
  255. pcParName=self.mod['partitionCoefficients'][cs]
  256. except KeyError:
  257. pcParName="one"
  258. #get flow (direction could be reversed)
  259. try:
  260. qName=self.flows['{}:{}'.format(cs,ct)]
  261. except KeyError:
  262. qName=self.flows['{}:{}'.format(ct,cs)]
  263. flowParName=self.mod['flows'][qName]
  264. #flowPar=pars[flowParName]
  265. self.bind(cs,ct,flowParName,pcParName)
  266. self.build()
  267. def add_default_parameters(self):
  268. pars=self.parSetup['parameters']
  269. pars['one']={'value':1}
  270. pars['zero']={'value':0}
  271. def M(self,t,y=numpy.array([])):
  272. for i in self.dM:
  273. for j in self.dM[i]:
  274. self.fM[i,j]=self.dM[i][j](t)
  275. #create an array and fill it with outputs of function at t
  276. if (y.size==0):
  277. return self.fM
  278. self.set_scaledM(t,y)
  279. return self.fM
  280. def set_scaledM(self,t,y):
  281. #prevent zero division
  282. eps=1e-8
  283. for c in self.scaled:
  284. i=self.lut[c]
  285. it=self.lut['total']
  286. try:
  287. k=numpy.copy(self.originalK[i])
  288. except AttributeError:
  289. k=numpy.copy(self.fM[i,:])
  290. self.originalK={}
  291. self.originalK[i]=k
  292. #make another copy
  293. k=numpy.copy(self.originalK[i])
  294. except KeyError:
  295. k=numpy.copy(self.fM[i,:])
  296. self.originalK[i]=k
  297. #make another copy
  298. k=numpy.copy(self.originalK[i])
  299. k[i]=k[i]-self.u(t)[it]
  300. #scale all inputs by total input mass
  301. for j in range(self.n):
  302. self.fM[i,j]=k[j]/(y[it]+eps)
  303. def u(self,t):
  304. ub=[f(t) for f in self.fu]
  305. ub[self.lut['total']]=sum(ub)
  306. return numpy.array(ub)
  307. def jacobiFull(self,t):
  308. #update jacobi created during build phase with time dependent values
  309. for i in self.dM:
  310. for j in self.dM[i]:
  311. for k in range(system.m+1):
  312. self.fJ[k*system.n+i,k*system.n+j]=self.dM[i][j](t)
  313. return self.fJ
  314. def fSS(self,t,y=numpy.array([])):
  315. for k in self.qSS:
  316. for i in self.qSS[k]:
  317. for j in self.qSS[k][i]:
  318. #print('[{},{},{}] {}'.format(k,i,j,self.qSS[k][i][j]))
  319. self.SS[k,i,j]=(self.qSS[k][i][j])(t)
  320. if y.size==0:
  321. return self.SS
  322. self.set_scaledSS(t,y)
  323. return self.SS
  324. def set_scaledSS(self,t,y):
  325. #prevent zero division
  326. eps=1e-8
  327. for c in self.scaled:
  328. it=self.lut['total']
  329. i=self.lut[c]
  330. try:
  331. dkdp=numpy.copy(self.originalSS[i])
  332. except AttributeError:
  333. dkdp=numpy.copy(self.SS[:,i,:])
  334. self.originalSS={}
  335. self.originalSS[i]=dkdp
  336. dkdp=numpy.copy(self.originalSS[i])
  337. except KeyError:
  338. dkdp=numpy.copy(self.SS[:,i,:])
  339. self.originalSS[i]=dkdp
  340. dkdp=numpy.copy(self.originalSS[i])
  341. self.SS[:,i,:]=dkdp/(y[it]+eps)
  342. #should add error on u!
  343. def fSY(self,y,t):
  344. #M number of sensitivity parameters
  345. #N number of equations
  346. #fSS is MxNxN
  347. #assume a tabulated solution y(t) at t spaced intervals
  348. qS=self.fSS(t,y).dot(y)
  349. #qS is MxN
  350. #but NxM is expected, so do a transpose
  351. #for simultaneous calculation, a Nx(M+1) matrix is expected
  352. tS=numpy.zeros((self.n,self.m+1))
  353. #columns from 2..M+1 are the partial derivatives
  354. tS[:,1:]=numpy.transpose(qS)
  355. #first column is the original function
  356. tS[:,0]=self.u(t)
  357. return tS
  358. def fS(self,t):
  359. #M number of sensitivity parameters
  360. #N number of equations
  361. #fSS is MxNxN
  362. #assume a tabulated solution y(t) at t spaced intervals
  363. qS=self.fSS(t).dot(self.getY(t))
  364. return numpy.transpose(qS)
  365. def getSEJ(self,parName):
  366. #find the sensitivity (SE) derivative of Jacobi with
  367. #respect to parameter
  368. try:
  369. return self.seJ[parName]
  370. except KeyError:
  371. self.seJ[parName]={}
  372. return self.seJ[parName]
  373. def getSEJ_comp(self,parName,compartmentName):
  374. #find equation dictating concentration in compartmentName
  375. #for jacobi-parameter derivative
  376. seJ=self.getSEJ(parName)
  377. try:
  378. return seJ[compartmentName]
  379. except KeyError:
  380. seJ[compartmentName]={}
  381. return seJ[compartmentName]
  382. def setY(self,t,y):
  383. self.tck=[None]*self.n
  384. for i in range(self.n):
  385. self.tck[i] = scipy.interpolate.splrep(t, y[:,i], s=0)
  386. def getY(self,t):
  387. fY=numpy.zeros(self.n)
  388. for i in range(self.n):
  389. fY[i]=scipy.interpolate.splev(t, self.tck[i], der=0)
  390. return fY
  391. def getWeight(self,parName):
  392. pars=self.parSetup['parameters']
  393. par=pars[parName]
  394. #self.get parses the units
  395. v=self.get(parName)["value"]
  396. if par['dist']=='lognormal':
  397. #this is sigma^2_lnx
  398. sln2=numpy.log(par["cv"]*par["cv"]+1)
  399. #have to multiplied by value to get the derivative
  400. #with respect to lnx
  401. return sln2*v*v
  402. else:
  403. #for Gaussian, cv is sigma/value; get sigma by value multiplication
  404. return par["cv"]*par["cv"]*v*v
  405. def getMax(lutSE):
  406. fm=-1
  407. for x in lutSE:
  408. if int(lutSE[x])>fm:
  409. fm=lutSE[x]
  410. return fm
  411. def getWeights(self,lutSE):
  412. #pars=self.parSetup['parameters']
  413. wts=numpy.zeros((model.getMax(lutSE)+1))
  414. for parName in lutSE:
  415. j=lutSE[parName]
  416. wts[j]=self.getWeight(parName)
  417. return wts
  418. def getDerivatives(self,se,i):
  419. #return latest point derivatives
  420. fse=se[-1][i]
  421. #fse is an m-vector
  422. return fse*fse
  423. def calculateUncertainty(self,se):
  424. s2out=numpy.zeros(se.shape[1:])
  425. se2=numpy.multiply(se,se)
  426. #w=self.getWeights(self.lutSE)
  427. w=numpy.ones((self.m))
  428. return numpy.sqrt(numpy.dot(se2,w))
  429. def get(self,parName):
  430. pars=self.parSetup['parameters']
  431. par=pars[parName]
  432. par['name']=parName
  433. if "value" in par:
  434. return self.getValue(par)
  435. if "function" in par:
  436. return self.getFunction(par)
  437. if "derived" in par:
  438. return self.getDerived(par)
  439. print('Paramter {} not found!'.format(parName))
  440. def getValue(self,par):
  441. v=par["value"]
  442. parName=par['name']
  443. #convert to seconds
  444. try:
  445. parUnits=par['unit'].split('/')
  446. except (KeyError,IndexError):
  447. #no unit given
  448. return valueObject(v,parName)
  449. timeUnit=self.getTimeUnit()
  450. try:
  451. if parUnits[1]==timeUnit:
  452. return valueObject(v,parName)
  453. except IndexError:
  454. #no / in unit name
  455. return valueObject(v,parName)
  456. if parUnits[1]=='min' and timeUnit=='s':
  457. return valueObject(v/60,parName)
  458. if parUnits[1]=='s' and timeUnit=='min':
  459. return valueObject(60*v,parName)
  460. if parUnits[1]=='day' and timeUnit=='min':
  461. return valueObject(v/24/60,parName)
  462. if parUnits[1]=='hour' and timeUnit=='min':
  463. return valueObject(v/60,parName)
  464. #no idea what to do
  465. return valueObject(v,parName)
  466. def getFunction(self,par):
  467. fName=par['function']
  468. #print('[{}]: getFunction({})'.format(par['name'],par['function']))
  469. df=self.parSetup['functions'][fName]
  470. skip=['type']
  471. par1={x:self.get(df[x]) for x in df if x not in skip}
  472. if df['type']=='linearGrowth':
  473. #print(par1)
  474. return function.linearGrowth(par1)
  475. if df['type']=='linearGrowthFixedSlope':
  476. return function.linearGrowthFixedSlope(par1)
  477. print('Function {} not found!'.format(fName))
  478. def getDerived(self,par):
  479. dName=par['derived']
  480. d=self.parSetup['derivedParameters'][dName]
  481. #print('Derived [{}]: type {}'.format(dName,d['type']))
  482. if d['type']=='product':
  483. #print('{}*{}'.format(d['a'],d['b']))
  484. pA=self.get(d['a'])
  485. a=pA['value']
  486. DA=pA['derivatives']
  487. pB=self.get(d['b'])
  488. b=pB['value']
  489. DB=pB['derivatives']
  490. #even more generic -> df/dp=[df/dA*dA/dp+df/dB*dB/dp]
  491. if any(['function' in pA,'function' in pB]):
  492. fa=function.to(a)
  493. fb=function.to(b)
  494. f=lambda t,a=fa,b=fb:a(t)*b(t)
  495. dfdA=lambda t,b=fb: b(t)
  496. dfdB=lambda t,a=fa: a(t)
  497. dA=function.generate(dfdA,DA)
  498. dB=function.generate(dfdB,DB)
  499. return function.Object(f,[dA,dB])
  500. else:
  501. return derivedObject(a*b,[{'df':b,'D':DA},{'df':a,'D':DB}])
  502. if d['type']=='power':
  503. #print('{}^{}'.format(d['a'],d['n']))
  504. pA=self.get(d['a'])
  505. a=pA['value']
  506. DA=pA['derivatives']
  507. pN=self.get(d['n'])
  508. n=pN['value']
  509. DN=pN['derivatives']
  510. if any(['function' in pA,'function' in pN]):
  511. fa=function.to(a)
  512. fn=function.to(n)
  513. f=lambda t,a=fa,n=fn:numpy.power(a(t),n(t))
  514. dfdA=lambda t,n=fn,f=f,a=fa:n(t)*f(t)/a(t)
  515. dfdN=lambda t,a=fa,f=f:numpy.log(a(t))*f(t)
  516. dA=function.generate(dfdA,DA)
  517. dN=function.generate(dfdN,DN)
  518. return function.Object(f,[dA,dN])
  519. else:
  520. f=numpy.power(a,n)
  521. return derivedObject(f,[{'df':n*f/a,'D':DA},{'df':f*numpy.log(a),'D':DN}])
  522. if d['type']=='ratio':
  523. #print('{}/{}'.format(d['a'],d['b']))
  524. pA=self.get(d['a'])
  525. a=pA['value']
  526. DA=pA['derivatives']
  527. pB=self.get(d['b'])
  528. b=pB['value']
  529. DB=pB['derivatives']
  530. #even more generic -> df/dp=[df/dA*dA/dp+df/dB*dB/dp]
  531. if any(['function' in pA,'function' in pB]):
  532. fa=function.to(a)
  533. fb=function.to(b)
  534. f=lambda t,a=fa,b=fb,:a(t)/b(t)
  535. dfdA=lambda t,f=f,a=fa: f(t)/a(t)
  536. dfdB=lambda t,f=f,b=fb: -f(t)/b(t)
  537. dA=function.generate(dfdA,DA)
  538. dB=function.generate(dfdB,DB)
  539. return function.Object(f,[dA,dB])
  540. else:
  541. return derivedObject(a/b,[{'df':1/b,'D':DA},{'df':-a/b/b,'D':DB}])
  542. if d['type']=='sum':
  543. #print('{}+{}'.format(d['a'],d['b']))
  544. pA=self.get(d['a'])
  545. a=pA['value']
  546. DA=pA['derivatives']
  547. pB=self.get(d['b'])
  548. b=pB['value']
  549. DB=pB['derivatives']
  550. #even more generic -> df/dp=[df/dA*dA/dp+df/dB*dB/dp]
  551. if any(['function' in pA,'function' in pB]):
  552. fa=function.to(a)
  553. fb=function.to(b)
  554. f=lambda t,a=fa,b=fb,:a(t)+b(t)
  555. dfdA=lambda t: 1
  556. dfdB=lambda t: 1
  557. dA=function.generate(dfdA,DA)
  558. dB=function.generate(dfdB,DB)
  559. return function.Object(f,[dA,dB])
  560. else:
  561. return derivedObject(a+b,[{'df':1,'D':DA},{'df':1,'D':DB}])
  562. def calculateDerivative(par):
  563. #add derivatives if dist(short for distribution) is specified
  564. return "dist" in par
  565. def sumValues(dArray,x):
  566. s=0
  567. for a in dArray:
  568. try:
  569. s+=a[x]
  570. except KeyError:
  571. continue
  572. return s
  573. def valueObject(v,parName):
  574. #convert everything to functions
  575. d0={parName:1}
  576. return {'value':v,'derivatives':{parName:1}}
  577. def derivedObject(f,ders):
  578. o={'value':f}
  579. DD=[]
  580. for d in ders:
  581. df=d['df']
  582. D=d['D']
  583. DD.append({x:df*D[x] for x in D})
  584. allKeys=[]
  585. for x in DD:
  586. allKeys.extend(x.keys())
  587. allKeys=list(set(allKeys))
  588. o['derivatives']={x:sumValues(DD,x) for x in allKeys}
  589. return o
  590. def splitVector(v):
  591. if v.find('(')<0:
  592. return [v]
  593. return v[1:-1].split(',')
  594. def parseFunction(formula):
  595. if formula['name']=='exponential':
  596. c0=formula['constant']
  597. k=formula['k']
  598. return lambda t,c=c0,k=k:c*numpy.exp(k*t)
  599. if formula['name']=='constant':
  600. c0=formula['value']
  601. return lambda t,c0=c0:c0
  602. if formula['name']=='Heavyside':
  603. t0=formula['limit']
  604. v=formula['value']
  605. return lambda t,v=v,t0=t0:v if t<t0 else 0
  606. return lambda t:1
  607. def addValue(qdict,compName,v):
  608. #add function to compName of dictionary qdict,
  609. #check if compName exists and handle the potential error
  610. #lambda functions can't be summed directly, so qdict is a list
  611. #that will be merged at matrix generation time
  612. try:
  613. qdict[compName].append(v)
  614. except KeyError:
  615. qdict[compName]=[v]
  616. #also add derivatives
  617. #
  618. # for d in dTarget:
  619. # ctarget=self.getSEJ_comp(d,target)
  620. # addValue(ctarget,target,dTarget[d])
  621. def get(timeUnit,par):
  622. v=par["value"]
  623. #convert to seconds
  624. try:
  625. parUnits=par['unit'].split('/')
  626. except (KeyError,IndexError):
  627. #no unit given
  628. return v
  629. try:
  630. if parUnits[1]==timeUnit:
  631. return v
  632. except IndexError:
  633. #no / in unit name
  634. return v
  635. if parUnits[1]=='min' and timeUnit=='s':
  636. return v/60
  637. if parUnits[1]=='s' and timeUnit=='min':
  638. return 60*v
  639. #no idea what to do
  640. return v