我们从Python开源项目中，提取了以下50个代码示例，用于说明如何使用scipy.integrate.quad()。
def _B_0_function(self, z): """ calculate B_0(z) function defined in: Gould A. 1994 ApJ 421L, 71 "Proper motions of MACHOs http://adsabs.harvard.edu/abs/1994ApJ...421L..71G Yoo J. et al. 2004 ApJ 603, 139 "OGLE-2003-BLG-262: Finite-Source Effects from a Point-Mass Lens" http://adsabs.harvard.edu/abs/2004ApJ...603..139Y """ out = 4. * z / np.pi function = lambda x: (1.-value**2*np.sin(x)**2)**.5 for (i, value) in enumerate(z): if value < 1.: out[i] *= ellipe(value*value) else: out[i] *= integrate.quad(function, 0., np.arcsin(1./value))[0] return out
def _bb_intensity(self, Teff, photon_weighted=False): """ Computes mean passband intensity using blackbody atmosphere: I_pb^E = \int_\lambda B(\lambda) P(\lambda) d\lambda / \int_\lambda P(\lambda) d\lambda I_pb^P = \int_\lambda \lambda B(\lambda) P(\lambda) d\lambda / \int_\lambda \lambda P(\lambda) d\lambda Superscripts E and P stand for energy and photon, respectively. @Teff: effective temperature in K @photon_weighted: photon/energy switch Returns: mean passband intensity using blackbody atmosphere. """ if photon_weighted: pb = lambda w: w*self._planck(w, Teff)*self.ptf(w) return integrate.quad(pb, self.wl[0], self.wl[-1])[0]/self.ptf_photon_area else: pb = lambda w: self._planck(w, Teff)*self.ptf(w) return integrate.quad(pb, self.wl[0], self.wl[-1])[0]/self.ptf_area
def _overpressureatscaledtime(r, y, alt, t): sgr = r / y**(1.0/3) shob = alt / y**(1.0/3) x_m = _scaledmachstemformationrange(y, alt) v = _slantrangescalingfactor(r, y, alt) r1 = _scale_slant_range(r, y, alt) / v ta_air = _scaledfreeairblastwavetoa(r1) * v dp = _scaledopposphasedur(r, y, alt) return _opatscaledtime(r, y, alt, sgr, shob, x_m, ta_air, dp, t) # In lieu of the 20-point Gauss-Legendre quadrature used in the original # BLAST.EXE, this fuction uses scipy.integrate.quad to call the venerable FORTRAN # library QUADPACK and perform a Gauss-Kronod quadrature. This appears # to be more accurate than the BLAST.EXE help file claims for the original # approach, which is not surprising as it uses an adaptive algorithm that # attempts to reduce error to within a particular tolerance. # IPTOTAL
def construct_table(): """:rtype: list""" solutions = [] MAX, _ = quad(integrand, 0, HALF_P) for i in [0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1]: power = MAX*i results = [] for tol in tolerance: if len(results) > 0: break for start in linspace(0, HALF_P, 1000): res, _ = quad(integrand, start, HALF_P) #print res if isclose(res, power, abs_tol=tol): results.append(start) mean = sum(results)/len(results) solutions.append(mean) conv = interp1d([min(solutions), max(solutions)], [DIM_MIN, DIM_MAX]) mapped_sol = conv(solutions) return [int(i) for i in mapped_sol] # = DIM_TABLE
def ShahCondensation_Average(x_min,x_max,Ref,G,D,p,TsatL,TsatV): # ******************************** # Necessary Properties # Calculated outside the quadrature integration for speed # ******************************** mu_f = PropsSI('V', 'T', TsatL, 'Q', 0, Ref) #kg/m-s cp_f = PropsSI('C', 'T', TsatL, 'Q', 0, Ref)#*1000 #J/kg-K k_f = PropsSI('L', 'T', TsatV, 'Q', 0, Ref)#*1000 #W/m-K Pr_f = cp_f * mu_f / k_f #[-] Pstar = p / PropsSI(Ref,'pcrit') h_L = 0.023 * (G*D/mu_f)**(0.8) * Pr_f**(0.4) * k_f / D #[W/m^2-K] def ShahCondensation(x,Ref,G,D,p): return h_L * ((1 - x)**(0.8) + (3.8 * x**(0.76) * (1 - x)**(0.04)) / (Pstar**(0.38)) ) if not x_min==x_max: #A proper range is given return quad(ShahCondensation,x_min,x_max,args=(Ref,G,D,p))[0]/(x_max-x_min) else: #A single value is given return ShahCondensation(x_min,Ref,G,D,p)
def integrate_function(function, interval): """ integrates the given function over given interval :param function: :param interval: :return: """ result = 0 err = 0 for area in interval: # res = integrate.quad(function, area[0], area[1]) res = complex_quadrature(function, area[0], area[1]) result += res[0] err += res[1] return np.real_if_close(result), err
def complex_quadrature(func, a, b, **kwargs): """ wraps the scipy qaudpack routines to handle complex valued functions :param func: callable :param a: lower limit :param b: upper limit :param kwargs: kwargs for func :return: """ def real_func(x): return np.real(func(x)) def imag_func(x): return np.imag(func(x)) real_integral = integrate.quad(real_func, a, b, **kwargs) imag_integral = integrate.quad(imag_func, a, b, **kwargs) return real_integral[0] + 1j * imag_integral[0], real_integral[1] + imag_integral[1]
def cdf_dlf(x, A, m1, a1, m2, a2, start=-26): ''' Cumulative Schechter function. Second LF is set to be 2*A of first LF. @param x: magnitude @param A: Scale factor @param m1: Knee of distribution 1 @param a1: Faint-end turnover of first lf @param m2: Knee of distribution 2 @param a2: Faint-end turnover of second lf @param start: Brightest magnitude @return Probability that galaxy has a magnitude greater than x ''' def integrate(in_x): return quad(dlf, start,in_x,args=(A,m1,a1,m2,a2))[0] if np.isscalar(x): x = np.array([x]) return np.fromiter(map(integrate,x),np.float,count=len(x))
def Ez( z, r, t) : """ Get the 2d Ez field Parameters ---------- z : 1darray t, r : float """ Nz = len(z) Nr = len(r) window_zmax = z.max() ez = np.zeros((Nz, Nr)) for iz in range(Nz) : for ir in range(Nr) : ez[iz, ir] = quad( kernel_Ez, z[iz]-c*t, window_zmax-c*t, args = ( z[iz]-c*t, r[ir] ), limit=30 )[0] return( ez )
def Er( z, r, t) : """ Get the 2d Ez field Parameters ---------- z : 1darray t, r : float """ Nz = len(z) Nr = len(r) window_zmax = z.max() er = np.zeros((Nz, Nr)) for iz in range(Nz) : for ir in range(Nr) : er[iz, ir] = quad( kernel_Er, z[iz]-c*t, window_zmax-c*t, args = ( z[iz]-c*t, r[ir] ), limit=200 )[0] return( er ) # --------------------------- # Comparison plots # ---------------------------
def length(self, t0=0, t1=1, error=LENGTH_ERROR, min_depth=LENGTH_MIN_DEPTH): """Calculate the length of the path up to a certain position""" if t0 == 0 and t1 == 1: if self._length_info['bpoints'] == self.bpoints() \ and self._length_info['error'] >= error \ and self._length_info['min_depth'] >= min_depth: return self._length_info['length'] # using scipy.integrate.quad is quick if _quad_available: s = quad(lambda tau: abs(self.derivative(tau)), t0, t1, epsabs=error, limit=1000)[0] else: s = segment_length(self, t0, t1, self.point(t0), self.point(t1), error, min_depth, 0) if t0 == 0 and t1 == 1: self._length_info['length'] = s self._length_info['bpoints'] = self.bpoints() self._length_info['error'] = error self._length_info['min_depth'] = min_depth return self._length_info['length'] else: return s
def quad_fourier(filtered_enb): """ Integrate fft function for fixed period filtered_enb - info about entropy blocks (size, entropy, shift) """ x = np.linspace(-pi_range, pi_range) y = 0.0 period_length = 10000.0 if len(filtered_enb) == 1 and filtered_enb[0][2] == 8.0: return 0, 0 for i in range(len(filtered_enb)): m = 1 if i % 2 != 0: m = -1 block_length = abs(filtered_enb[i][1] - filtered_enb[i][0]) if block_length == 0: continue T = block_length / period_length y += m * filtered_enb[i][2] * \ np.sin(x * ((2 * np.pi) / T) + block_length / period_length) yf = fft(y) return integrate.quad(lambda a: yf[a], -pi_range, pi_range)
def VCACPFF(engine,app): ''' This method calculates the chemical potential or filling factor. ''' engine.rundependences(app.name) engine.cache.pop('pt_kmesh',None) kmesh,nk=app.BZ.mesh('k'),app.BZ.rank('k') fx=lambda omega,mu: (np.trace(engine.mgf_kmesh(omega=mu+1j*omega,kmesh=kmesh),axis1=1,axis2=2)-engine.ncopt/(1j*omega-app.p)).sum().real if app.task=='CP': gx=lambda mu: quad(fx,0,np.float(np.inf),args=mu)[0]/nk/engine.ncopt/np.pi-app.cf mu=broyden2(gx,app.options.pop('x0',0.0),**app.options) engine.log<<'mu(error): %s(%s)\n'%(mu,gx(mu)) if app.returndata: return mu else: filling=quad(fx,0,np.float(np.inf),args=app.cf)[0]/nk/engine.ncopt/np.pi engine.log<<'Filling factor: %s\n'%filling if app.returndata: return filling
def VCAOP(engine,app): ''' This method calculates the order parameters. ''' engine.rundependences(app.name) engine.cache.pop('pt_kmesh',None) cgf,kmesh,nk=engine.apps[engine.preloads[0]],app.BZ.mesh('k'),app.BZ.rank('k') ops,ms={},np.zeros((len(app.terms),engine.ncopt,engine.ncopt),dtype=np.complex128) for i,term in enumerate(app.terms): order=deepcopy(term) order.value=1.0 for opt in HP.Generator(engine.lattice.bonds,engine.config,table=engine.config.table(mask=engine.mask),terms=[order],half=True).operators.itervalues(): ms[i,opt.seqs[0],opt.seqs[1]]+=opt.value ms[i,:,:]+=ms[i,:,:].T.conjugate() fx=lambda omega,m: (np.trace(engine.mgf_kmesh(omega=app.mu+1j*omega,kmesh=kmesh).dot(m),axis1=1,axis2=2)-np.trace(m)/(1j*omega-app.p)).sum().real for term,m,dtype in zip(app.terms,ms,app.dtypes): ops[term.id]=quad(fx,0,np.float(np.inf),args=m)[0]/nk/engine.ncopt*2/np.pi if dtype in (np.float32,np.float64): ops[term.id]=ops[term.id].real engine.log<<HP.Sheet(corner='Order',rows=['Value'],cols=ops.keys(),contents=np.array(ops.values()).reshape((1,-1)))<<'\n' if app.returndata: return ops
def __init__(self, name, K, L, spaceChargeOn, multipart, twiss, beamdata, nbrOfSplits): LinearElement.__init__(self, "quad " + name) #self.name = name self.K = K self.L = L gamma = gammaFromBeta(beamdata[0]) self.M = self.createMatrixM(K, L, beamdata[0], gamma) # M should be a 6x6 matrix self.T = self.createMatrixT(self.M) # M should be a 9x9 matrix # disunite matrices self.n = nbrOfSplits self.Lsp = self.L/self.n self.Msp = self.createMatrixM(self.K, self.Lsp, beamdata[0], gamma) self.Tsp = self.createMatrixT(self.Msp) # space charge class self.spaceChargeOn = spaceChargeOn if self.spaceChargeOn == 1: #self.sc = SpaceCharge('quad_sc', self.Lsp, multipart, twiss, beamdata) # OLD self.sc = SpaceCharge('quad_sc', self.Lsp, beamdata) # NEW elif self.spaceChargeOn == 2: self.sc = SpaceChargeEllipticalIntegral('quad_sc', self.Lsp, beamdata)
def timeTransitFactor(self, beta): # E_z0_of_s z1z2 = -2*self.halfnbrofoscillations*self.L/3 # /3 since A = 0 and B =/= 0 #print "z1z2: " + str(z1z2) z4z5 = 2*self.halfnbrofoscillations*self.L/3 # /3 since A = 0 and B =/= 0 #print "z4z5: " + str(z4z5) ## Integral # -inf to z1|z2 I1 = quad(lambda s: self.amplitudeB*exp(((s+z1z2)/self.sigma)**self.p)*cos(2*constants.pi/(beta*self.rf_lambda)*s - self.phi_s),-inf,z1z2) #print "I1: " + str(I1) # z2 to z4||z5 I2 = quad(lambda s: (self.amplitudeA*cos(constants.pi*s/self.L)+self.amplitudeB*cos(3*constants.pi*s/self.L))*cos(2*constants.pi/(beta*self.rf_lambda)*s - self.phi_s),z1z2,z4z5) #print "I2: " + str(I2) # z5 to inf I3 = quad(lambda s: self.amplitudeB*exp((-(s-z4z5)/self.sigma)**self.p)*cos(2*constants.pi/(beta*self.rf_lambda)*s - self.phi_s),z4z5,inf) #print "I3: " + str(I3) # sum up res = I1[0]+I2[0]+I3[0] res = res/self.E_0 return res
def __call__(self, sat_pos, args=(), **kwds): """ Return the definite integral of *self.fun*(pos, *args*[0], ..., *args*[-1]) for the line of site from *stn_pos* to *sat_pos* (a :class:`PyPosition`) where pos is a :class:`PyPosition` on the line of site (and note the integration bounds on h defined in __init__). The remaining *kwds* are passed to the quadrature routine (:py:func:`quad`). """ diff = NP.array(sat_pos.xyz) - NP.array(self.stn_pos.xyz) S_los = NP.linalg.norm(diff) / 1e3 def pos(s): """ Return the ECEF vector a distance *s* along the line-of-site (in [km]). """ return PyPosition(*(NP.array(self.stn_pos.xyz) + (s / S_los) * diff)) # determine integration bounds # distance along of line of site at which the geodetic height # is self.height1 s1 = minimize_scalar(lambda l: (pos(l).height / 1e3 - self.height1)**2, bounds=[0, S_los], method='Bounded').x # distance along of line of site at which the geodetic height # is self.height2 s2 = minimize_scalar(lambda l: (pos(l).height / 1e3 - self.height2)**2, bounds=[0, S_los], method='Bounded').x def wrapper(s, *args): return self.fun(pos(s), *args) return quad(wrapper, s1, s2, args=args, **kwds)[0]
def _bjelkeformel_P(mast, j, fh): """Beregner deformasjoner i kontakttrådhøyde grunnet en punklast. Dersom lasten angriper under kontakttrådhøyde: :math:`\\delta = \\frac{P*x^2}{6EI}(3fh-x)` Dersom lasten angriper over kontakttrådhøyde: :math:`\\delta = \\frac{P*fh^2}{6EI}(3x-fh)` :param Mast mast: Aktuell mast som beregnes :param Kraft j: Last som skal påføres ``mast`` :param float fh: Kontakttrådhøyde :math:`[m]` :return: Matrise med forskyvningsbidrag :math:`[mm]` :rtype: :class:`numpy.array` """ D = numpy.zeros((5, 8, 3)) if j.e[0] < 0: E = mast.E delta_topp = mast.h + j.e[0] delta_topp = 0 if delta_topp < 0 else delta_topp L = (mast.h - delta_topp) * 1000 delta_y = integrate.quad(mast.Iy_int_P, 0, L, args=(delta_topp,)) delta_z = integrate.quad(mast.Iz_int_P, 0, L, args=(delta_topp,)) I_y = L ** 3 / (3 * delta_y[0]) I_z = L ** 3 / (3 * delta_z[0]) f_y = j.f[1] f_z = j.f[2] x = -j.e[0] * 1000 fh *= 1000 if fh > x: D[j.type[1], j.type[0], 1] = (f_z * x**2) / (6 * E * I_y) * (3 * fh - x) D[j.type[1], j.type[0], 0] = (f_y * x**2) / (6 * E * I_z) * (3 * fh - x) else: D[j.type[1], j.type[0], 1] = (f_z * fh ** 2) / (6 * E * I_y) * (3 * x - fh) D[j.type[1], j.type[0], 0] = (f_y * fh ** 2) / (6 * E * I_z) * (3 * x - fh) return D
def _bjelkeformel_q(mast, j, fh): """Beregner deformasjoner i kontakttrådhøyde grunnet en fordelt last. Funksjonen beregner horisontale forskyvninger basert på følgende bjelkeformel: :math:`\\delta = \\frac{q*fh^2}{24EI}(fh^2+6h^2-4h*fh)` Lasten antas å være jevnet fordelt over hele mastens høyde :math:`h`, med resultant i høyde :math:`h/2` :param Mast mast: Aktuell mast som beregnes :param Kraft j: Last som skal påføres ``mast`` :param float fh: Kontakttrådhøyde :math:`[m]` :return: Matrise med forskyvningsbidrag :math:`[mm]` :rtype: :class:`numpy.array` """ D = numpy.zeros((5, 8, 3)) if j.b > 0 and j.e[0] < 0: E = mast.E delta_topp = mast.h - j.b L = (mast.h - delta_topp) * 1000 delta_y = integrate.quad(mast.Iy_int_q, 0, L, args=(delta_topp,)) delta_z = integrate.quad(mast.Iz_int_q, 0, L, args=(delta_topp,)) I_y = L ** 4 / (4 * delta_y[0]) I_z = L ** 4 / (4 * delta_z[0]) q_y = j.q[1] / 1000 q_z = j.q[2] / 1000 b = j.b * 1000 fh *= 1000 D[j.type[1], j.type[0], 1] = ((q_z * fh ** 2) / (24 * E * I_y)) * (fh ** 2 + 6 * b ** 2 - 4 * b * fh) D[j.type[1], j.type[0], 0] = ((q_y * fh ** 2) / (24 * E * I_z)) * (fh ** 2 + 6 * b ** 2 - 4 * b * fh) return D
def expected_value_integral(func, bounds, parameters): return quad(func, *bounds, parameters)[0] / (bounds[1] - bounds[0])
def expected(t): h = lambda x:g(x,t) return integrate.quad(h,0,np.inf)[0]
def variance(t): h = lambda x:g2(x,t) return integrate.quad(h,0,np.inf)[0]
def _interp_sync_kernel(xmin=1e-3, xmax=10.0, xsample=256): """ Sample the synchrotron kernel function at the specified X positions and make an interpolation, to optimize the speed when invoked to calculate the synchrotron emissivity. WARNING ------- Do NOT simply bound the synchrotron kernel within the specified [xmin, xmax] range, since it decreases as a power law of index 1/3 at the left end, and decreases exponentially at the right end. Bounding it with interpolation will cause the synchrotron emissivity been *overestimated* on the higher frequencies. Parameters ---------- xmin, xmax : float, optional The lower and upper cuts for the kernel function. Default: [1e-3, 10.0] xsample : int, optional Number of samples within [xmin, xmax] used to do interpolation. Returns ------- F_interp : function The interpolated kernel function ``F(x)``. """ xx = np.logspace(np.log10(xmin), np.log10(xmax), num=xsample) Fxx = [xp * integrate.quad(lambda t: scipy.special.kv(5/3, t), a=xp, b=np.inf)[0] for xp in xx] F_interp = interpolate.interp1d(xx, Fxx, kind="quadratic", bounds_error=True, assume_sorted=True) return F_interp
def growth_factor(self, z): """ Growth factor at redshift z. References: Ref.[randall2002],Eq.(A7) """ x0 = (2 * self.Ode0 / self.Om0) ** (1/3) x = x0 / (1 + z) coef = np.sqrt(x**3 + 2) / (x**1.5) integral = integrate.quad(lambda y: y**1.5 * (y**3+2)**(-1.5), a=0, b=x, epsabs=1e-5, epsrel=1e-5)[0] D = coef * integral return D
def _overpressuretotalimpulse(r, y, alt): sgr = r / y**(1.0/3) shob = alt / y**(1.0/3) x_m = _scaledmachstemformationrange(y, alt) v = _slantrangescalingfactor(r, y, alt) r1 = _scale_slant_range(r, y, alt) / v ta_air = _scaledfreeairblastwavetoa(r1) * v dp = _scaledopposphasedur(r, y, alt) t_p = 13 * (ta_air + dp) / 14 scaled_impulse, _ = quad(lambda t: _opatscaledtime(r, y, alt, sgr, shob, x_m, ta_air, dp, t), ta_air, ta_air + dp) return scaled_impulse * y**(1.0/3)
def quad(self, dt, t): return integrate.quad(self.iuh, max(t-1.+dt, 0.), t+dt)[0]
def update_coefs(self): """(Re)calculate the MA coefficients based on the instantaneous unit hydrograph.""" coefs = [] sum_coefs = 0. for t in itertools.count(0., 1.): coef = integrate.quad(self.quad, 0., 1., args=(t,))[0] sum_coefs += coef if (sum_coefs > .9) and (coef < self.smallest_coeff): coefs = numpy.array(coefs) coefs /= numpy.sum(coefs) self.coefs = coefs break else: coefs.append(coef)
def MinOptionOnSpdCall(F1, F2, dv1, dv2, rho, K1, K2, T): ' min(max(F1-K1),max(F2-K2)) assuming F1 F2 are spread of two assets' v1 = dv1 * numpy.sqrt(T) v2 = dv2 * numpy.sqrt(T) def int_func1(x): return scipy.stats.norm.cdf(((F1-K1)-(F2-K2) + (v1 * rho - v2) * x)/(v1 * numpy.sqrt(1-rho**2))) \ * (v2 * x + F2- K2) * scipy.stats.norm.pdf(x) def int_func2(x): return scipy.stats.norm.cdf(((F2-K2)-(F1-K1) + (v2 * rho - v1) * x)/(v2 * numpy.sqrt(1-rho**2))) \ * (v1 * x + F1- K1) * scipy.stats.norm.pdf(x) res1 = quad(int_func1, (K2-F2)/v2, numpy.inf) res2 = quad(int_func2, (K1-F1)/v1, numpy.inf) return res1[0] + res2[0]
def min_on_call(F1, F2, dv1, dv2, rho, K1, K2, T): v1 = dv1 * np.sqrt(T) v2 = dv2 * np.sqrt(T) def int_func1(x): return scipy.stats.norm.cdf(((F1-K1)-(F2-K2) + (v1 * rho - v2) * x)/(v1 * np.sqrt(1-rho**2))) \ * (v2 * x + F2- K2) * scipy.stats.norm.pdf(x) def int_func2(x): return scipy.stats.norm.cdf(((F2-K2)-(F1-K1) + (v2 * rho - v1) * x)/(v2 * np.sqrt(1-rho**2))) \ * (v1 * x + F1- K1) * scipy.stats.norm.pdf(x) res1 = quad(int_func1, (K2-F2)/v2, np.inf) res2 = quad(int_func2, (K1-F1)/v1, np.inf) return res1[0] + res2[0]
def __init__(self, mcmc_samples, bin_edges): """ Histogram Inference a la Dan Foreman-Mackey Parameters: =========== - mcmc_samples: numpy array of shape (Nobs, Nsamples) MCMC samples for the thing you want to histogram - bin_edges: numpy.ndarray array The edges of the histogram bins to use. """ self.mcmc_samples = mcmc_samples self.bin_edges = bin_edges self.bin_centers = (self.bin_edges[:-1] + self.bin_edges[1:]) / 2 self.bin_widths = np.diff(self.bin_edges) self.Nbins = self.bin_widths.size self.Nobs = self.mcmc_samples.shape[0] # Find which bin each q falls in self.bin_idx = np.digitize(self.mcmc_samples, self.bin_edges) - 1 # Determine the censoring function for each bin (used in the integral) self.censor_integrals = np.array([quad(func=self.censoring_fcn, a=left, b=right)[0] for (left, right) in zip(self.bin_edges[:-1], self.bin_edges[1:])]) # Set values needed for multinest fitting self.n_params = self.Nbins self.param_names = [r'$\theta_{}$'.format(i) for i in range(self.Nbins)]
def integrate_gauss(x1, x2, amp, mean, sigma): """ Integrate a gaussian between the points x1 and x2 """ gauss = lambda x, A, mu, sig: A*np.exp(-(x-mu)**2 / (2.0*sig**2)) if x1 < -1e6: x1 = -np.inf if x2 > 1e6: x2 = np.inf result = quad(gauss, x1, x2, args=(amp, mean, sigma)) return result[0]
def _sanity_test(self): # Marginal of x integrates to one. print quad(lambda x: np.exp(self.logpdf_x(x)), self.D[0], self.D[1]) # Marginal of y integrates to one. print quad(lambda y: np.exp(self.logpdf_y(y)), -1 ,1) # Joint of x,y integrates to one; quadrature will fail for small noise. print dblquad( lambda y,x: np.exp(self.logpdf_xy(x,y)), self.D[0], self.D[1], lambda x: -1, lambda x: 1)
def NPV(rho,g,eff,H,q_process,delta_t,p,Q,a,b): integrand = lambda t: np.exp(-r*t)*power(rho,g,eff,H,q_process)*delta_t*p PV = quad(integrand, 0, np.inf) NPVal = PV[0] - cost(a, b, Q) return NPVal # check the function, it does not look right!!!
def EV_flow(mean,std,Q): EV = quad(lambda x: lognorm_pdf(x,mean,std)*x, 0, Q) EV_2 = (1 - lognorm_cdf(Q,mean,std))*Q EV_total = EV[0] + EV_2 return EV_total
def NPV(rho,g,eff,H,q_process,delta_t,p,Q,a,b): integrand = lambda t: np.exp(-r*t)*power(rho,g,eff,H,q_process)*delta_t*p PV = quad(integrand, 0, np.inf) NPVal = PV[0] - cost(a, b, Q) return NPVal
def Bertsch_MC_Average(x_min,x_max,Ref,G,Dh,q_flux,L,TsatL,TsatV): ''' Returns the average heat transfer coefficient between qualities of x_min and x_max. for Bertsch two-phase evaporation in mico-channel HX ''' if not x_min==x_max: #A proper range is given return quad(Bertsch_MC,x_min,x_max,args=(Ref,G,Dh,q_flux,L))[0]/(x_max-x_min) else: #A single value is given return Bertsch_MC(x_min,Ref,G,Dh,q_flux,L,TsatL,TsatV)
def consumer_surp(self): "Compute consumer surplus" # == Compute area under inverse demand function == # integrand = lambda x: (self.ad/self.bd) - (1/self.bd)* x area, error = quad(integrand, 0, self.quantity()) return area - self.price() * self.quantity()
def producer_surp(self): "Compute producer surplus" # == Compute area above inverse supply curve, excluding tax == # integrand = lambda x: -(self.az/self.bz) + (1/self.bz) * x area, error = quad(integrand, 0, self.quantity()) return (self.price() - self.tax) * self.quantity() - area
def __init__(self, ap_paramList, ra1,dec1, ra2, dec2, background_density, z): ''' @param ap_paramList[seed]: Seed for random number generator @param ra1: Left right ascension @param dec1: Bottom declination @param ra2: Right right ascension @param dec2: Top declination @param background_density: galaxy background density in galaxies/square degree @param z: Redshift of galaxy cluster ''' self.ra1 = ra1 self.dec1 = dec1 self.ra2 = ra2 self.dec2 = dec2 self.background_density = background_density self.z = z self.__H0 = 70 self.__Omega_m = 0.3 self.__R0 = 2 self.__Rs = 0.15 / 0.7 self.__Rcore = 0.1 / 0.7 self.__norm = 1/quad(nfw,0,self.__R0,args=(1,self.__Rs,self.__Rcore))[0] super(CatalogGenerator, self).__init__(ap_paramList)
def nfw_cumulative(self,R): ''' Cumulative radial NFW distribution @param R: Radius @return float: Probability of being within R ''' R0 = self.__R0 norm = self.__norm Rs = self.__Rs Rcore = self.__Rcore def integrate(z): return quad(nfw,0,z,args=(norm,Rs,Rcore))[0] if np.isscalar(R): R = np.array([float(R)]) else: R = R.astype(np.float) res = np.piecewise(R, [R < 0.0, np.logical_and(R >= 0.0, R < R0), R >= R0], [lambda z: z, lambda z: np.fromiter(map(integrate,z),np.float), lambda z: z/R0]) if len(res) == 1: return res[0] else: return res
def get_Kstar_max(self): """ value of K*(dimensionless vortex constant) for which weight flow in maximum See also: NASA TN D-4421 eq. 23-25 """ max_val = 0 while(1): try: f = lambda Mstar: (1 - (max_val/self.Mlstar)**2 * Mstar**2) ** (1/(self.gamma - 1)) integrate.quad(f,self.Mlstar,self.Mustar) except: max_val -= 0.1 break max_val += 0.1 def func(Kstar_max): a1 = (1 - Kstar_max**2)**(1/(self.gamma -1)) a2 = (1 - (Kstar_max**2) * (self.Mustar/self.Mlstar)**2) a3 = (1/(self.gamma - 1)) # This if is to avoid the bug in python3 it self if(a2<0): a2 = ((Kstar_max**2) * (self.Mustar/self.Mlstar)**2 - 1) a = a1 + a2**a3 else : a = a1 - a2 ** a3 f = lambda Mstar: (1 - (Kstar_max/self.Mlstar)**2 * Mstar**2) ** (1/(self.gamma - 1)) b = integrate.quad(f,self.Mlstar,self.Mustar) return a - b[0] Kstar_max = optimize.brentq(func,0.1,max_val) return Kstar_max
def get_C(self): """ get reduction in maximum weight flow due to two-dimensional flow See also: NASA TN D-4421 eq. 34b """ Kstar_max = self.get_Kstar_max() a = 1 - np.sqrt((self.gamma + 1)/(self.gamma - 1)) * ((self.gamma + 1) / 2)**(1 / (self.gamma - 1)) * (self.Mustar / (self.Mustar - self.Mlstar)) f = lambda Mstar: Kstar_max * (1 - (Kstar_max / self.Mlstar)**2 * Mstar**2)**(1 / (self.gamma - 1)) b = integrate.quad(f, self.Mlstar, self.Mustar) C = a * b[0] return C
def get_Q(self): """ get vortex flow parameter Q See also: NASA TN D-4421 eq. 34a """ f = lambda Mstar: ((self.gamma + 1) / 2 - (self.gamma - 1) / 2 * Mstar**2)**(1 / (self.gamma - 1)) a = integrate.quad(f, self.Mlstar, self.Mustar) Q = (self.Mlstar * self.Mustar) / (self.Mustar - self.Mlstar) * a[0] return Q
def maximal_age(z): z = np.double(z) #Cosmological Constants O_m = 0.266 O_r = 0. O_k= 0. O_L = 1. - O_m H_0 = 74.3 #km/s/Mpc H_sec = H_0 / 3.0857e19 secondsinyear = 31556926 ageoftheuniverse = 13.798e9 # Equation for the time elapsed since z and now a = 1/(1+z) E = O_m * (1+z)**3 + O_r *(1+z)**4 + O_k *(1+z) + O_L integrand = lambda z : 1 / (1+z) / sqrt( O_m * (1+z)**3 + O_r *(1+z)**4 + O_k *(1+z) + O_L ) #Integration z_obs = z z_cmb = 1089 #As Beta (not cmb). But 1089 (cmb) would be the exagerated maximun possible redshift for the birth z_now = 0 integral, error = quad( integrand , z_obs, z_cmb) # #t = ageoftheuniverse - (integral * (1 / H_sec) / secondsinyear) t = (integral * (1 / H_sec)) / secondsinyear return t
def z2Dlum(z): """ Calculate luminosity distance from redshift. """ #Cosmo Constants O_m = 0.266 O_r = 0. O_k= 0. O_L = 1. - O_m H_0 = 70. #km/s/Mpc H_sec = H_0 / 3.0857e19 # equation a = 1/(1+z) E = O_m * (1+z)**3 + O_r *(1+z)**4 + O_k *(1+z) + O_L integrand = lambda z : 1 / sqrt(O_m * (1+z)**3 + O_r *(1+z)**4 + O_k *(1+z) + O_L) #integration z_obs = z z_now = 0 c_cm = 2.997e10 integral = quad( integrand , z_now, z_obs) dlum_cm = (1+z)*c_cm/ H_sec * integral[0] dlum_Mpc = dlum_cm/3.08567758e24 return dlum_cm
