Adding log-normal convolution

models/humanHG.json

@@ -110,7 +110,8 @@
         "intestineVolume":{"value":0.98,"unit":"kg"},
         "kRBC":{"value":0.60,"unit":"l/min"},
         "kU":{"value":0,"unit":"l/min"},
-        "kH":{"value":2.823e-6,"unit":"l/min"},
+        "commentkH":"assume pc-times larger input flow, 2.823e-6*248.7",
+        "kH":{"value":7.021e-4,"unit":"l/min"},
         "kBR":{"value":4.033e-3,"unit":"l/min"},
         "kI":{"value":4.033e-6,"unit":"l/min"},
         "kB":{"value":4.033e-5,"unit":"l/min"},

pythonScripts/convolveLN.py

@@ -0,0 +1,126 @@
+import numpy
+import os
+import scipy.interpolate
+import filionQuadrature
+
+def ln(mu,cv):
+   #generate log-normal distributed variable with mean mu and coefficient of variation cv
+    m=numpy.log(mu/numpy.sqrt(1+cv*cv))
+    s=numpy.sqrt(numpy.log(1+cv*cv))
+    return ln_s(m,s)
+
+def ln_s(mu_s,sigma_s):
+   #generate r.v. with log-normal distribution
+   #mu_s and sigma_s are parameters of the normal distribution in the exponent
+   #see ln() for targeted mean and coefficient of variation
+    z=numpy.random.normal()
+    return numpy.exp(mu_s+sigma_s*z)
+
+def readCoef(sigma):
+   #read coefficients of gamma distribution set that approximates LN distro
+   #with mu_s=0 and sigma_s=sigma, valid for 0.05<=sigma<=3
+   #from 
+   #Furman E, Hackmann D, Kuznetsov A. 
+   #On log-normal convolutions: An analytical–numerical method with 
+   #applications to economic capital determination.
+   #Insurance: Mathematics and Economics, 90 (2020) 120-134.
+   #return fa and fb which are alpha and beta for n=20 approximations 
+    
+    fdir=os.path.join(os.path.expanduser('~'),'software','src','PBPK',
+      'lnnCoeficients20','cBase.txt')
+    with open(fdir,'r') as f:
+        lines=f.read().splitlines()
+    fm=[q.split(' ') for q in lines]
+    fm=numpy.array([[float(q) for q in x] for x in fm])
+    fn=(fm.shape[1]-1)//2
+    #print(fn)
+    fa=fm[:,1:fn+1]
+    fb=fm[:,fn+1:2*fn+1]
+    fx=fm[:,0]
+    falpha=[]
+    fbeta=[]
+    for i in range(fn):
+        f2a = scipy.interpolate.interp1d(fx, fa[:,i], kind='cubic')
+        falpha.append(f2a(sigma))
+        f2b = scipy.interpolate.interp1d(fx, fb[:,i], kind='cubic')
+        fbeta.append(f2b(sigma))
+    return numpy.array(falpha),numpy.array(fbeta)
+
+def getLTransform(z,fa,fb, mu=0):
+   #evaluate Laplace transform of a gamma-set convolution with 
+   #coefficients alpha,fa and beta,fb
+   #at point z with an optional parameter mu equal to log-normal parameter mu_s, 
+   #see labeling at ln() and ln_s() functions
+    s=1
+    fn=len(fa)
+    for i in range(fn):
+        s*=numpy.power((1+z*numpy.exp(mu)/fb[i]),-fa[i])
+    return s
+
+def getLTransformGrid(sigma,fz, mu=0):
+   #evaluate the Laplace transform of a gamma-set convolution 
+   #that approximates a log-normal distribution with parameters
+   #sigma_s=sigma, mu_s=mu at an array of complex numbers fz
+    fa,fb=readCoef(sigma)
+    fn=len(fa)
+    return numpy.array([getLTransform(z,fa,fb,mu) for z in fz])
+
+def getLTransformPart(z, fa, fb, i):
+   #get a Laplace transform of a single gamma distribution, specified by index
+   #into a set of alpha and beta arrays fa and fb, respectively
+    return numpy.power((1+z/fb[i]),-fa[i])
+
+def getLTransformGridPart(sigma,fz,i):
+   #evaluate the Laplace transform of a single gamma distribution appearing at
+   #index i of alpha/beta pairs that approximate log-normal distribution 
+   #for sigma_s=sigma  
+    fa,fb=readCoef(sigma)
+    fn=len(fa)
+    print('Using alpha={} beta={}'.format(fa[i],fb[i]))
+    return numpy.array([getLTransformPart(z,fa,fb,i) for z in fz])
+
+def addExpA(fz,fq,x):
+    #apply the exp(c*x)*exp(Re(a)*u*x term, see paper for reasoning
+    return fq*numpy.exp(numpy.real(fz)*x)
+
+def convolveLN(fx,sigma,mu):
+   #convolve a set of LN with parameters sigma and mu as numpy arrays and
+   #return the pdf of convolution evaluated at the set of points fx
+    
+    #internal computation parameters
+    #number of grid points 1M suggested in paper, good fit found for few thousand
+    n=2000
+    #integral range, 250 to 500 suggested in paper
+    u0=0
+    u1=300
+    #constant c (to the right of all poles, saw artefacts for c of 2 or above)
+    c=1
+    #constant a (inclined integral path, seems to work best for Re(a) close to 0,
+   #-0.1>a>-0.5 was the range tested in paper, using upper limit
+    fa=numpy.complex(-0.1,1)
+    
+    #internal structures
+    qh=numpy.linspace(u0,u1,2*n+1)
+    h=(u1-u0)/2/n
+    #curve segmentation
+    qz=c+fa/numpy.absolute(fa)*qh
+    #function values at grid points
+    qPhi=numpy.ones(qh.size,dtype=complex)
+    for i in range(len(sigma)):
+        m=mu[i]
+        s=sigma[i]
+        qPhi*=getLTransformGrid(sigma[i],qz,mu[i])
+        
+    
+    #Filion quadratures by parts of complex numbers
+    fI=numpy.zeros(len(fx),dtype=complex)
+    for i in range(len(fx)):
+        x=fx[i]
+        qPhi1=addExpA(qz,qPhi,x)
+        fRc=filionQuadrature.filionQuadrature(numpy.real(qPhi1),x,u0,h,n,'cos')
+        fIc=filionQuadrature.filionQuadrature(numpy.imag(qPhi1),x,u0,h,n,'cos')
+        fRs=filionQuadrature.filionQuadrature(numpy.real(qPhi1),x,u0,h,n,'sin')
+        fIs=filionQuadrature.filionQuadrature(numpy.imag(qPhi1),x,u0,h,n,'sin')
+        fI[i]=numpy.complex(fRc-fIs,fIc+fRs)
+    fI1=fa*fI
+    return numpy.imag(fI1)

pythonScripts/copy_data.py

@@ -0,0 +1,9 @@
+import os
+import sys
+
+def loadData(pars):
+    print(pars)
+
+
+if __name__=="__main__":
+    loadData(sys.argv)

pythonScripts/dataTransfer.py

@@ -0,0 +1,89 @@
+import os
+import sys
+import json
+
+def loadData(parFile):
+    #parameterFile is a json with settings
+
+    with open(parFile,'r') as f:
+        pars=json.load(f)
+
+    ds=openData(pars['data']['filename'],pars['data']['format'])
+    #var names are here
+    #print(ds[0])
+    #open data
+    transferToDestination(ds,pars)
+
+def openData(qFile,qFormat="DTA"):
+    if qFormat=="DTA":
+        import pandas as pd
+        ds=pd.read_stata(qFile)
+        #print(ds.dtypes)
+        rows=[]
+        n=len(ds.columns)
+        rng=range(len(ds.columns))
+        for row in ds.itertuples(index=False):
+            varRow={ds.columns[i]:row[i] for i in rng if not pd.isna(row[i])}
+            rows.append(varRow)
+        return rows
+
+def transferToDestination(rows,pars):
+    dest=pars['destination']
+    if dest['format']=='LabKey':
+        return transferToLabKey(rows,pars)
+
+def transferToLabKey(rows,pars):
+    src=pars['data']
+    dest=pars['destination']
+    fsetup=os.path.join(os.path.expanduser('~'),'.labkey','setup.json')
+    with open(fsetup,'r') as f:
+        setup=json.load(f)
+
+    sys.path.append(setup['paths']['nixWrapper'])
+    import nixWrapper
+    nixWrapper.loadLibrary('labkeyInterface')
+    import labkeyInterface
+    import labkeyDatabaseBrowser
+
+    net=labkeyInterface.labkeyInterface()
+    net.init(dest['server'])
+    db=labkeyDatabaseBrowser.labkeyDB(net)
+
+    for r in rows:
+        idField=dest['participantField']
+        idValue=r[src['idField']]
+        idFilter={'variable':idField,'value':'{}'.format(idValue),'oper':'eq'}
+
+        try: 
+            sField=dest['visitField']
+        except KeyError:
+            sField='SequenceNum'
+        try:
+            seqNum=src['visitField']
+        except KeyError:
+            seqNum="0.0"
+        
+        visitFilter={'variable':sField,'value':seqNum,'oper':'eq'}
+
+        ds=db.selectRows(dest['project'],dest['schema'],dest['query'],
+                [idFilter,visitFilter])
+        if len(ds['rows'])>0:
+            mode='update'
+            outRow=ds['rows'][0]
+        else:
+            mode='insert'
+            outRow={idField:idValue,sField:seqNum}
+        for var in r.keys():
+            outRow[var]=r[var]
+        
+        resp=db.modifyRows(mode,dest['project'],dest['schema'],dest['query'],
+                [outRow])
+        print(resp)
+    print('Done')
+
+        
+
+
+
+if __name__=="__main__":
+    loadData(sys.argv[1])

pythonScripts/filionQuadrature.py

@@ -0,0 +1,77 @@
+import numpy
+
+def Taylor(x,coef):
+    q=1
+    s=0
+    for c in coef:
+        s+=c*q
+        q*=x
+    return s
+
+def filionAlpha(th,th3,sth,cth):
+    if abs(th)<0.1:
+        return Taylor(th,[0,0,0,2/45,0,-2/315,0,2/4725])
+    
+    th2=th*th
+    return (th2+sth*cth*th-2*sth*sth)/th3
+
+def filionBeta(th,th3,sth,cth):
+    if abs(th)<0.1:
+        return Taylor(th,[2/3,0,2/15,0,-4/105,0,2/567,0,-4/22275])
+    
+    return 2*((1+cth*cth)*th-2*sth*cth)/th3
+
+def filionGamma(th,th3,sth,cth):
+    if abs(th)<0.1:
+        return Taylor(th,[4/3,0,-2/15,0,1/210,0,-1/11340])
+    
+    return 4*(sth-th*cth)/th3
+
+def cEven(qs,n):
+
+   idx=numpy.arange(0,2*n+1,2)
+   return numpy.sum(qs[idx])
+
+def qFunc(qf,t,x0,h,n,func):
+
+   qx=numpy.linspace(x0,x0+2*h*n,2*n+1)
+   qs=[q*func(t*x) for x,q in zip(qx,qf)]
+   qs[0]*=0.5
+   qs[len(qs)-1]*=0.5
+   return numpy.array(qs)
+        
+def cOdd(qs,n):
+   idx=numpy.arange(1,2*n+1,2)
+   return numpy.sum(qs[idx])
+
+def filionQuadrature(qf,t,x0,h,n,func,debug=0):
+   #calculate integral f(x) func(tx) in [x0,x0+2*n*h] 
+   #where func is either cos or sine 
+    cFunc=numpy.cos
+    iFunc=numpy.sin
+    sign=1
+    if func=='sin':
+        cFunc=numpy.sin
+        iFunc=numpy.cos
+        sign=-1
+    th=t*h
+    sth=numpy.sin(th)
+    cth=numpy.cos(th)
+    th3=th*th*th
+    x2n=x0+2*n*h
+    A=filionAlpha(th,th3,sth,cth)
+    B=filionBeta(th,th3,sth,cth)
+    C=filionGamma(th,th3,sth,cth)
+    qs=qFunc(qf,t,x0,h,n,cFunc)
+    CE=cEven(qs,n)
+    CO=cOdd(qs,n)
+    D=sign*(qf[2*n]*iFunc(t*x2n)-qf[0]*iFunc(t*x0))
+    r=h*(A*D+B*CE+C*CO)
+    if debug:
+      print('fillon r={:.2f} A={:.2f} D={:.2f} B={:.2f} CE={:.2f} C={:.2f} CO={:.2f}'.format(r,A,D,B,CE,C,CO))
+    return r
+    
+def tabulateFunction(f,x0,n,h):
+   qx=[x0+i*h for i in range(2*n+1)]
+   qy=[f(x) for x in qx]
+   return qx,qy

pythonScripts/propagateErrorLN.py

@@ -0,0 +1,68 @@
+import numpy
+import convolveLN
+
+def generateSample(mu,cv):
+   return numpy.array([convolveLN.ln(mu,cv) for i in range(10000)])
+
+def getEdges(a):
+   #convert array of bin centers a to array of bin edges with one more element than a
+   a1=numpy.zeros(len(a)+1)
+   a2=numpy.zeros(a1.shape)
+   a1[1:]=a
+   a2[:-1]=a
+   e=0.5*(a1+a2)
+   e[0]=2*e[1]-e[2]
+   e[-1]=2*e[-2]-e[-3]
+   return e
+
+def calculateDistribution(fy,muTarget,mu,cv,dydp):
+   #calculate distribution of random variable y with a functional dependence 
+   #on random parameters p, each with a log-normal 
+   #distribution with a mean mu(p) and coefficients of variation cv(p) 
+   #and dydp giving partial derivatives of y on  p
+   #muTarget is the target mean of y distribution
+   #fy is the set of values y where distribution is evaluated at
+   #Parameters mu, cv and dydp are given as vectors/lists
+    
+   qa=mu*dydp
+   A=numpy.sum(qa)
+   qa*=muTarget/A
+   cvPrime=cv*A/muTarget
+   sigmaS=numpy.sqrt(numpy.log(1+cvPrime*cvPrime))
+   muS=numpy.log(qa/numpy.sqrt(1+cvPrime*cvPrime))
+   y=convolveLN.convolveLN(fy,sigmaS,muS)
+   return y/numpy.sum(y)
+
+def generateDistribution(fx,muTarget,mu,cv,dydp, dbg=0):
+   #sum randomly generated variables to approximate distribution of 
+   #variable y with a functional dependence y(p) on random parameters p, 
+   #each with a log-normal (LN) distribution 
+   #with a mean mu(p) and coefficient of variation cv(p)
+   #target mean of y is muTarget
+   #distribution is evaluated at points fy, 
+   #dydp are derivatives of y with respect to parameters p
+   
+   qa=mu*dydp
+   A=numpy.sum(qa)
+   qa*=muTarget/A
+   cvPrime=cv*A/muTarget
+   samples=[a*generateSample(1,c) for a,c in zip(qa,cvPrime)]
+   mean=[numpy.mean(s) for s in samples]
+   std=[numpy.std(s) for s in samples]
+   convSample=numpy.zeros(len(samples[0]))
+   msg='[{}] mean={:.2f} std={:.2f} cv={:.2f}'.
+   for i in range(len(cv)):
+      convSample+=samples[i]
+      if dbg:
+         print(msg.format(i,mean[i],std[i],std[i]/mean[i]))
+      
+   convM=numpy.mean(convSample)
+   convS=numpy.std(convSample)
+   if dbg:
+      code=x
+      print(msg.format(code,convM,convS,convS/convM))
+   h,be=numpy.histogram(convSample,bins=getEdges(fx))
+   #h should be len(fx)
+   return h/numpy.sum(h)
+    
+