Updates

pythonScripts/convolveLN.py

@@ -121,7 +121,7 @@ def getZArray(n=2000):
     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)
+    fa=complex(-0.1,1)
     
     #internal structures
     qh=numpy.linspace(u0,u1,2*n+1)
@@ -140,7 +140,7 @@ def inverseL(fx,qz,qPhi,fa,u0,h,n):
       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)
+      fI[i]=complex(fRc-fIs,fIc+fRs)
    fI1=fa*fI
    return numpy.imag(fI1)
 

pythonScripts/propagateErrorLN.py

@@ -18,7 +18,76 @@ def getEdges(a):
    e[-1]=2*e[-2]-e[-3]
    return e
 
-def calculateDistribution(fy,muTarget,mu,cv,dydp):
+def deltaDistribution(xTarget):
+   yTarget=numpy.zeros(len(xTarget))
+#fill bin at 0
+   n=numpy.argmin(numpy.abs(xTarget))
+   yTarget[n]=1
+   return yTarget
+
+
+def scale(xTarget,xOrig,yOrig,scale,sign=1):
+   f = scipy.interpolate.CubicSpline(xOrig, yOrig)
+   x0=numpy.min(xOrig)
+   x1=numpy.max(xOrig)
+   yTarget=[f(sign*scale*x) if sign*x>x0 and sign*scale*x<x1 else 0 for x in xTarget]
+   s=numpy.sum(yTarget)
+   if s==0:
+      yTarget=deltaDistribution(xTarget)
+   else:
+      yTarget/=s
+   return yTarget
+
+def cDistScaled(gx,qa,cv,sign=+1,nbins=100,umax=5):
+   sigmaS=numpy.sqrt(numpy.log(1+cv*cv))
+   muS=numpy.log(qa/numpy.sqrt(1+cv*cv))
+   q=1/numpy.sum(qa)
+   muS1=muS+numpy.log(q)
+#temporary space to convolve
+   qx=numpy.linspace(0,umax,nbins)
+   qy=convolveLN.convolveLN(qx,sigmaS,muS1)
+   gy=scale(gx,qx,qy,q,sign)
+   return gy
+
+def shift(xTarget,xOrig,yOrig,shift):
+   f = scipy.interpolate.CubicSpline(xOrig, yOrig)
+   x0=numpy.min(xOrig)
+   x1=numpy.max(xOrig)
+   yTarget=[f(x-shift) if x-shift>x0 and x-shift<x1 else 0 for x in xTarget]
+   s=numpy.sum(yTarget)
+   if s==0:
+        #fill bin at 0
+      n=numpy.argmin(numpy.abs(xTarget))
+      yTarget[n]=1
+   else:
+      yTarget/=s
+   return yTarget
+
+def calculateDistribution(xTarget,muTarget,mu,cv,dydp,nbins=100,umax=5):
+   #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 distribution
+   #xTarget is the set of values x where distribution is evaluated at
+   #Parameters mu, cv and dydp are given as vectors/lists
+   
+   qa=mu*dydp
+   y0=numpy.max(numpy.abs(qa))
+   fx=numpy.linspace(-5*y0,5*y0,201)
+   
+   sign=qa>=0
+   fy1=cDistScaled(fx,qa[sign],cv[sign],nbins=nbins,umax=umax)
+   fy2=cDistScaled(fx,-qa[~sign],cv[~sign],sign=-1,nbins=nbins,umax=umax)
+
+   fy=scipy.signal.convolve(fy1,fy2,'same')
+
+   a=numpy.average(fx,weights=fy)
+   return shift(xTarget,fx,fy,muTarget-a)
+
+
+
+def calculateDistribution1(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) 
@@ -27,7 +96,8 @@ def calculateDistribution(fy,muTarget,mu,cv,dydp):
    #fy is the set of values y where distribution is evaluated at
    #Parameters mu, cv and dydp are given as vectors/lists
    
-
+   print(f'muTarget={muTarget}')
+   print('mu={}'.format(mu))
    #parameters qa should be around 1
    qa=mu*dydp
    print('qa={}'.format(qa))