PBPK_public/pythonScripts/propagateErrorLN.py

import numpy
import scipy.interpolate
import importlib
import convolveLN
importlib.reload(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
   

   #parameters qa should be around 1
   qa=mu*dydp
   print('qa={}'.format(qa))
   sign=qa>=0
   Ax=numpy.sum(qa[sign])
   Ay=numpy.sum(-qa[~sign])
   print('A={},{}'.format(Ax,Ay))
   if Ax>0:
      qa[sign]/=Ax
   if Ay>0:
      qa[~sign]/=Ay
   cvPrime=cv
   print('cvPrime={}'.format(cvPrime))

   #convolve qa multiplied (1,cv) distributions 
   sigmaS=numpy.sqrt(numpy.log(1+cvPrime*cvPrime))
   print('sigmaS={}'.format(sigmaS))
   muS=numpy.zeros(qa.shape)
   muS[sign]=numpy.log(qa[sign]/numpy.sqrt(1+cvPrime[sign]*cvPrime[sign]))
   muS[~sign]=numpy.log(-qa[~sign]/numpy.sqrt(1+cvPrime[~sign]*cvPrime[~sign]))
   print('muS={}'.format(muS))

   #do the convolution
   nbins=100
   umax=5
   fw=numpy.linspace(0,umax,nbins)
   wx=numpy.array([])
   wy=numpy.array([])
   
   if Ax>0:
      wx=convolveLN.convolveLN(fw,sigmaS[sign],muS[sign])
   
   if Ay>0:
      wy=convolveLN.convolveLN(fw,sigmaS[~sign],muS[~sign])

   #for correlate, binning of wx and wy should be equal, but one has scale A1, the other A2
   #how do you deal with that?
   if Ax>0:
      fwx=Ax*fw
   if Ay>0:
      fwy=Ay*fw

   fw=fwx
   if not Ay>Ax:
      if wy.size>0:
      #re-sample wy
         wy = scipy.interpolate.CubicSpline(fwy, wy)(fw)

   if Ay>Ax:
      fw=fwy
      #resample wx
      if wx.size>0:
         wx = scipy.interpolate.CubicSpline(fwx, wx)(fw)
   
   #output is twice the length of w and symmetric around 0
   fw1=numpy.concatenate((numpy.flip(-fw[1:]),fw))
   w=numpy.zeros(fw1.size)

   if wx.size==0:
      w[:len(wy)]=numpy.flip(wy)
   if wy.size==0:
      w[-len(wx):]=wx
   if wx.size>0 and wy.size>0:
      w=scipy.signal.correlate(wx,wy)
      #w is a convolution of input LN distributions with a mean of A


   #now shift it to the target meat muTarget
   B=muTarget-(Ax-Ay)
   print('B={}'.format(B))

   #interpolate on w and fill y
   f = scipy.interpolate.CubicSpline(fw1, w)
     
   y=numpy.zeros(fy.shape)
   
   for i in range(len(y)):
      try:
         y[i]=f(fy[i]-B)
      except ValueError:
         pass
   
   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*=A
   print('qa={}'.format(qa))
   cvPrime=cv
   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={:.2g} std={:.2g} 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]))
   convSample+=muTarget-A   
   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)