PBPK_public/pythonScripts/function.py

import numpy
import functools

def fval(t,W):
   if callable(W):
      return W(t)
   return W

def L(t,t0,W):
   ft0=fval(t,t0)
   fW=fval(t,W)
   x=(t-ft0)/fW
   if x<-18:
      return 0
   if x>18:
      return 1
   return 1/(1+numpy.exp(-x))
def dL(t,t0,W):
   return L(t,t0,W)*(1-L(t,t0,W))

def dLdW(t,t0,W):
   fW=fval(t,W)
   ft0=fval(t,t0)
   return dL(t,t0,W)*(t-ft0)/fW/fW

def dLdt0(t,t0,W):
   fW=fval(t,W)
   return dL(t,t0,W)/fW

def logisticPar(par):
   t0=par['t0']['value']
   W=par['W']['value']
   f=functools.partial(L,t0=t0,W=W)
   
   #D1[x] could be a function !!!
   Dt0=par['t0']['derivatives']
   dfdt0=functools.partial(dLdt0,t0=t0,W=W)
   dt0=generate(dfdt0,Dt0)
   #D2[x] could be a function !!!
   DW=par['W']['derivatives']
   dfdW=functools.partial(dLdW,t0=t0,W=W)
   dW=generate(dfdW,DW)
   #merge
   return Object(f,[dt0,dW])

def deltaW(t,t1,t2):
   ft1=fval(t,t1)
   ft2=fval(t,t2)
   return ft2-ft1

def deltaT(t,t1):
   ft1=fval(t,t1)
   return t-t1

def fq(t,t1,t2):
   return deltaT(t,t1)/deltaW(t,t1,t2)

def f1(t,C,t1,t2,L1,L2):
   fC=fval(t,C)
   return fC*fq(t,t1,t2)*L1(t)*(1-L2(t))

def f2(t,C,L1,L2):
   fC=fval(t,C)
   return fC*L1(t)*L2(t)

def lg(t,C,t1,t2,eps,L1,L2):
   fC=fval(t,C)
   return eps*fC+f1(t,C,t1,t2,L1,L2)+f2(t,C,L1,L2)

def dlgdt1(t,C,t1,t2,L1,L2):
   fC=fval(t,C)
   return -fC*L1(t)*(1-L2(t))*(1+fq(t,t1,t2))/deltaW(t,t1,t2)

def dlgdt2(t,C,t1,t2,L1,L2):
   return -f1(t,C,t1,t2,L1,L2)/deltaW(t,t1,t2)

def dlgdL1(t,C,t1,t2,L1,L2):
   fC=fval(t,C)
   fL2=L2(t)
   return fC*(fq(t,t1,t2)*(1-fL2)+fL2)

def dlgdL2(t,C,t1,t2,L1,L2):
   fC=fval(t,C)
   return fC*L1(t)*(1-fq(t,t1,t2))

def dlgdC(t,C,t1,t2,L1,L2):
   fL2=L2(t)
   return L1(t)*(fq(t,t1,t2)*(1-fL2)+fL2)

def lgFixedSlope(t,alpha,t1,t2,eps,L1,L2):
   fA=fval(t,alpha)
   ft1=fval(t,t1)
   dW=deltaW(t,t1,t2)
   return fA*(dW*eps+L1(t)*((t-ft1)*(1-L2(t))+dW*L2(t)))

def dlgFixedSlopedt1(t,alpha,L1,L2):
   fA=fval(t,alpha)
   return -fA*L1(t)

def dlgFixedSlopedt2(t,alpha,L1,L2):
   fA=fval(t,alpha)
   return -fA*L1(t)*L2(t)

def dlgFixedSlopedL1(t,alpha,t1,t2,L1,L2):
   return lgFixedSlope(t,alpha,t1,t2,L1,L2)/L1(t)

def dlgFixedSlopedL2(t,alpha,t1,t2,L1,L2):
   fA=fval(t,alpha)
   ft2=fval(t,t2)
   return fA*L1(t)*(ft2-t)

def dlgFixedSlopedalpha(t,alpha,t1,t2,L1,L2):
   fA=fval(t,alpha)
   return lgFixedSlope(t,alpha,t1,t2,L1,L2)/fA

def linearGrowth(par):
   #make sure we never return zero or negative
   C=par['C']['value']
   t1=par['t1']['value']
   t2=par['t2']['value']
   eps=1e-8

   pL1=logisticPar({'t0':par['t1'],'W':par['W1']})
   pL2=logisticPar({'t0':par['t2'],'W':par['W2']})
   L1=pL1['value']
   L2=pL2['value']
   f=functools.partial(lg,C=C,t1=t1,t2=t2,eps=eps,L1=L1,L2=L2)

   Dt1=par['t1']['derivatives']
   dfdt1=functools.partial(dlgdt1,C=C,t1=t1,t2=t2,L1=L1,L2=L2)
   dt1=generate(dfdt1,Dt1)
   
   Dt2=par['t2']['derivatives']
   dfdt2=functools.partial(dlgdt2,C=C,t1=t1,t2=t2,L1=L1,L2=L2)
   dt2=generate(dfdt2,Dt2)

   DL1=pL1['derivatives']
   dfdL1=functools.partial(dlgdL1,C=C,t1=t1,t2=t2,L1=L1,L2=L2)
   dL1=generate(dfdL1,DL1)

   DL2=pL2['derivatives']
   dfdL2=functools.partial(dlgdL2,C=C,t1=t1,t2=t2,L1=L1,L2=L2)
   dL2=generate(dfdL2,DL2)

   DC=par['C']['derivatives']
   dfdC=functools.partial(dlgdC,C=C,t1=t1,t2=t2,L1=L1,L2=L2)
   dC=generate(dfdC,DC)

   return Object(f,[dt1,dt2,dL1,dL2,dC])

def linearGrowthFixedSlope(par):
   #make sure we never return zero or negative
   alpha=par['alpha']['value']
   t1=par['t1']['value']
   t2=par['t2']['value']
   eps=1e-8

   pL1=logisticPar({'t0':par['t1'],'W':par['W1']})
   pL2=logisticPar({'t0':par['t2'],'W':par['W2']})
   L1=pL1['value']
   L2=pL2['value']
   f=functools.partial(lgFixedSlope,alpha=alpha,t1=t1,t2=t2,eps=eps,L1=L1,L2=L2)

   Dt1=par['t1']['derivatives']
   dfdt1=functools.partial(dlgFixedSlopedt1,alpha=alpha,L1=L1,L2=L2)
   dt1=generate(dfdt2,Dt1)
   
   Dt2=par['t2']['derivatives']
   dfdt2=functools.partial(dlgFixedSlopedt2,alpha=alpha,L1=L1,L2=L2)
   dt2=generate(dfdt2,Dt2)

   DL1=pL1['derivatives']
   dfdL1=functools.partial(dlgFixedSlopedL1,alpha=alpha,t1=t1,t2=t2,L1=L1,L2=L2)
   dL1=generate(dfdL1,DL1)

   DL2=pL2['derivatives']
   dfdL2=functools.partial(dlgFixedSlopedL2,alpha=alpha,t1=t1,t2=t2,L1=L1,L2=L2)
   dL2=generate(dfdL2,DL2)

   DA=par['alpha']['derivatives']
   dfdA=functools.partial(dlgFixedSlopedalpha,alpha=alpha,t1=t1,t2=t2,L1=L1,L2=L2)
   dA=generate(dfdA,DA)

   return Object(f,[dt1,dt2,dL1,dL2,dA])



def Object(v,dArray):
   o={'value':v,'function':True}
   allKeys=[]
   for x in dArray:
      allKeys.extend(x.keys()) 
   allKeys=list(set(allKeys))
   o['derivatives']={x:sumDict(dArray,x) for x in allKeys}
   return o

def sumDict(dArray,x):
   #print('sumFunctions for {}'.format(x))
   fs=[]
   for a in dArray:
      #print('\tChecking {}'.format(a))
      try:
         #check if a[x] is a function with a dummy parameter
         fs.append(a[x])
         #print('\tAdding [{}]'.format(v))
      except KeyError:
         #print('\t{} not in table'.format(x))
         continue
   return sumArray(fs)

def sumArray(fs):
   return lambda t,fs=fs: sum([to(f)(t) for f in fs])

def to(a):
   if callable(a):
      return a
   return lambda t,a=a:a
 
def contains(arr):
   for a in arr:
      if callable(a):
         return True
   return False

def multiply(t,a,b):
   #abstract actual type of a and b
   fa=fval(t,a)
   fb=fval(t,b)
   return fa*fb

def generate(df,D):
   return {x:functools.partial(multiply,a=D[x],b=df) for x in D}