我们从Python开源项目中，提取了以下23个代码示例，用于说明如何使用scipy.integrate.trapz()。
def fastMie_SD(m, wavelength, dp, ndp): # http://pymiescatt.readthedocs.io/en/latest/inverse.html#fastMie_SD dp = coerceDType(dp) ndp = coerceDType(ndp) _length = np.size(dp) Q_sca = np.zeros(_length) Q_abs = np.zeros(_length) Q_back = np.zeros(_length) aSDn = np.pi*((dp/2)**2)*ndp*(1e-6) for i in range(_length): Q_sca[i],Q_abs[i],Q_back[i] = fastMieQ(m,wavelength,dp[i]) Bsca = trapz(Q_sca*aSDn) Babs = trapz(Q_abs*aSDn) Bback = trapz(Q_back*aSDn) return Bsca, Babs, Bback
def density_energy_electron(spectrum, gamma): """ Calculate the energy density of relativistic electrons. Parameters ---------- spectrum : 1D float `~numpy.ndarray` The number density of the electrons w.r.t. Lorentz factors Unit: [cm^-3] gamma : 1D float `~numpy.ndarray` The Lorentz factors of electrons Returns ------- e_re : float The energy density of the relativistic electrons. Unit: [erg cm^-3] """ e_re = integrate.trapz(spectrum*gamma*AU.mec2, gamma) return e_re
def _integrate(self, y, f): """ Integrate `y` along axis=1, i.e. over freq axis for all T. Parameters ---------- y : 2d array (nT, ndos) where nT = len(self.T), ndos = len(self.dos) f : self.f, (len(self.dos),) Returns ------- array (nT,) """ mask = np.isnan(y) if mask.any(): self._printwarn("HarmonicThermo._integrate: warning: " " %i NaNs found in y!" %len(mask)) if self.fixnan: self._printwarn("HarmonicThermo._integrate: warning: " "fixing %i NaNs in y!" %len(mask)) y[mask] = self.nanfill # this call signature works for scipy.integrate,{trapz,simps} return self.integrator(y, x=f, axis=1)
def debye_func(x, nstep=100, zero=1e-8): """Debye function :math:`f(x) = 3 \int_0^1 t^3 / [\exp(t x) - 1] dt` Parameters ---------- x : float or 1d array nstep : int number of points for integration zero : float approximate the 0 in the integral by this (small!) number """ x = np.atleast_1d(x) if x.ndim == 1: x = x[:,None] else: raise StandardError("x is not 1d array") tt = np.linspace(zero, 1.0, nstep) return 3.0 * trapz(tt**3.0 / (np.exp(tt*x) - 1.0), tt, axis=1)
def calculate_auc(self, y, dx): """ calculate the Area Under the Curve using the average of the estimated areas by the two composites the Simpsons's and the trapezoidal rules. """ AUC = (simps(y, dx=dx) + trapz(y, dx=dx)) / 2. return AUC
def trapez_area(sefl, y, dx): """ Calculate the area under a curve using the composite trapezoidal rule. """ return trapz(y, dx=dx)
def Mie_SD(m, wavelength, dp, ndp, interpolate=False, asDict=False): # http://pymiescatt.readthedocs.io/en/latest/forward.html#Mie_SD dp = coerceDType(dp) ndp = coerceDType(ndp) _length = np.size(dp) Q_ext = np.zeros(_length) Q_sca = np.zeros(_length) Q_abs = np.zeros(_length) Q_pr = np.zeros(_length) Q_back = np.zeros(_length) Q_ratio = np.zeros(_length) g = np.zeros(_length) # scaling of 1e-6 to cast in units of inverse megameters - see docs aSDn = np.pi*((dp/2)**2)*ndp*(1e-6) # _logdp = np.log10(dp) for i in range(_length): Q_ext[i], Q_sca[i], Q_abs[i], g[i], Q_pr[i], Q_back[i], Q_ratio[i] = AutoMieQ(m,wavelength,dp[i]) Bext = trapz(Q_ext*aSDn) Bsca = trapz(Q_sca*aSDn) Babs = Bext-Bsca Bback = trapz(Q_back*aSDn) Bratio = trapz(Q_ratio*aSDn) bigG = trapz(g*Q_sca*aSDn)/trapz(Q_sca*aSDn) Bpr = Bext - bigG*Bsca if asDict: return dict(Bext=Bext, Bsca=Bsca, Babs=Babs, G=bigG, Bpr=Bpr, Bback=Bback, Bratio=Bratio) else: return Bext, Bsca, Babs, bigG, Bpr, Bback, Bratio
def SF_SD(m, wavelength, dp, ndp, minAngle=0, maxAngle=180, angularResolution=0.5, space='theta', angleMeasure='radians', normalization=None): # http://pymiescatt.readthedocs.io/en/latest/forward.html#SF_SD _steps = int(1+(maxAngle-minAngle)/angularResolution) ndp = coerceDType(ndp) dp = coerceDType(dp) SL = np.zeros(_steps) SR = np.zeros(_steps) SU = np.zeros(_steps) kwargs = {'minAngle':minAngle, 'maxAngle':maxAngle, 'angularResolution':angularResolution, 'space':space, 'normalization':None} for n,d in zip(ndp,dp): measure,l,r,u = ScatteringFunction(m,wavelength,d,**kwargs) SL += l*n SR += r*n SU += u*n if normalization in ['n','N','number','particles']: _n = trapz(ndp,dp) SL /= _n SR /= _n SU /= _n elif normalization in ['m','M','max','MAX']: SL /= np.max(SL) SR /= np.max(SR) SU /= np.max(SU) elif normalization in ['t','T','total','TOTAL']: SL /= trapz(SL,measure) SR /= trapz(SR,measure) SU /= trapz(SU,measure) return measure,SL,SR,SU
def main(num_samples, burn, lag, w): alpha = 1.0 N = 25 Z = gu.simulate_crp(N, alpha) K = max(Z) + 1 # CRP with gamma prior. log_pdf_fun = lambda alpha : gu.logp_crp_unorm(N, K, alpha) - alpha proposal_fun = lambda : np.random.gamma(1.0, 1.0) D = (0, float('Inf')) samples = su.slice_sample(proposal_fun, log_pdf_fun, D, num_samples=num_samples, burn=burn, lag=lag, w=w) minval = min(samples) maxval = max(samples) xvals = np.linspace(minval, maxval, 100) yvals = np.array([math.exp(log_pdf_fun(x)) for x in xvals]) yvals /= trapz(xvals, yvals) fig, ax = plt.subplots(2,1) ax[0].hist(samples, 50, normed=True) ax[1].hist(samples, 100, normed=True) ax[1].plot(xvals,-yvals, c='red', lw=3, alpha=.8) ax[1].set_xlim(ax[0].get_xlim()) ax[1].set_ylim(ax[0].get_ylim()) plt.show()
def _area(self, y): from scipy.integrate import trapz # Before the integration: # flip the crs curve to x axis # and start at y=0 yy = [abs(p-1) for p in y] deltax = self.wvl[1] - self.wvl[0] area = trapz(yy, dx=deltax) return area
def __init__(self, wavelength, transmit, name='', dtype="photon", unit=None): """Constructor""" self.name = name try: # get units from the inputs self._wavelength = wavelength.magnitude self.wavelength_unit = str(wavelength.units) except AttributeError: self._wavelength = wavelength self.wavelength_unit = unit self.transmit = transmit self.norm = trapz(transmit, self._wavelength) self._lT = trapz(self._wavelength * transmit, self._wavelength) self._lpivot = np.sqrt( self._lT / trapz(transmit / self._wavelength, self._wavelength) ) self._cl = self._lT / self.norm self.set_dtype(dtype)
def leff(self): """ Unitwise Effective wavelength leff = int (lamb * T * Vega dlamb) / int(T * Vega dlamb) """ with Vega() as v: s = self.reinterp(v.wavelength) w = s._wavelength leff = np.trapz(w * s.transmit * v.flux.magnitude, w, axis=-1) leff /= np.trapz(s.transmit * v.flux.magnitude, w, axis=-1) if self.wavelength_unit is not None: return leff * unit[self.wavelength_unit] else: return leff
def get_hist_density(data, nbins=25): y, x = np.histogram(data, nbins) x = x[1:] - 0.5*(x[1] - x[0]) p = y / trapz(y, x) return [p, x]
def auprc(self): recall, precision = self.precision_recall() return trapz(precision, recall)
def auroc(self): fp, tp = self.roc_curve() return trapz(tp, fp)
def integration(self,is_quadrant=0,use_flux_per_mrad2_or_mm2=1): if (self.X is None or self.Y is None): raise Exception(" X and Y must be array for integration") if self.X.shape != self.Y.shape: raise Exception(" X and Y must have the same shape") if len(self.X.shape)==2 : X = np.linspace(self.X.min(), self.X.max(), self.X.shape[0]) Y = np.linspace(self.Y.min(), self.Y.max(), self.Y.shape[1]) res1=integrate.trapz(self.intensity, X) res = integrate.trapz(res1, Y) # res = integrate.simps(integrate.simps(self.intensity, X), Y) # res = self.intensity.sum() * (self.X[1,0] - self.X[0,0]) * (self.Y[0,1] - self.X[0,0]) # regular grid only else : # X and Y are 1d array if len(self.X) == 1: res = self.intensity[0] else: # choix arbitraire XY = np.zeros_like(self.X) for i in range(1, len(self.X)): XY[i] = XY[i-1]+np.sqrt((self.X[i] - self.X[i-1]) * 2 + (self.Y[i] - self.Y[i-1]) ** 2) res = np.trapz(self.intensity, XY) # Note that the value of flux is in phot/s/0.1%bw/mrad2 (or .../mm2) and # our grid is in rad (or m), therefore we must account for this: if use_flux_per_mrad2_or_mm2: res *= 1e6 # in case the calculation is for a quadrant, the integral is four times the calculated value if is_quadrant: res *= 4 return res
def norm_int(y, x, area=1.0, scale=True, func=simps): """Normalize integral area of y(x) to `area`. Parameters ---------- x,y : numpy 1d arrays area : float scale : bool, optional Scale x and y to the same order of magnitude before integration. This may be necessary to avoid numerical trouble if x and y have very different scales. func : callable Function to do integration (like scipy.integrate.{simps,trapz,...} Called as ``func(y,x)``. Default: simps Returns ------- scaled y Notes ----- The argument order y,x might be confusing. x,y would be more natural but we stick to the order used in the scipy.integrate routines. """ if scale: fx = np.abs(x).max() fy = np.abs(y).max() sx = x / fx sy = y / fy else: fx = fy = 1.0 sx, sy = x, y # Area under unscaled y(x). _area = func(sy, sx) * fx * fy return y*area/_area
def ScatteringFunction(m, wavelength, diameter, minAngle=0, maxAngle=180, angularResolution=0.5, space='theta', angleMeasure='radians', normalization=None): # http://pymiescatt.readthedocs.io/en/latest/forward.html#ScatteringFunction x = np.pi*diameter/wavelength _steps = int(1+(maxAngle-minAngle)/angularResolution) # default 361 if angleMeasure in ['radians','RADIANS','rad','RAD']: adjust = np.pi/180 elif angleMeasure in ['gradians','GRADIANS','grad','GRAD']: adjust = 1/200 else: adjust = 1 if space in ['q','qspace','QSPACE','qSpace']: _steps *= 10 if minAngle==0: minAngle = 1e-5 measure = np.logspace(np.log10(minAngle),np.log10(maxAngle),_steps)*np.pi/180 _q = True else: measure = np.linspace(minAngle,maxAngle,_steps)*adjust _q = False if x == 0: return measure,0,0,0 _measure = np.linspace(minAngle,maxAngle,_steps)*np.pi/180 SL = np.zeros(_steps) SR = np.zeros(_steps) SU = np.zeros(_steps) for j in range(_steps): u = np.cos(_measure[j]) S1, S2 = MieS1S2(m,x,u) SL[j] = (np.sum(np.conjugate(S1)*S1)).real SR[j] = (np.sum(np.conjugate(S2)*S2)).real SU[j] = (SR[j]+SL[j])/2 if normalization in ['m','M','max','MAX']: SL /= np.max(SL) SR /= np.max(SR) SU /= np.max(SU) elif normalization in ['t','T','total','TOTAL']: SL /= trapz(SL,measure) SR /= trapz(SR,measure) SU /= trapz(SU,measure) if _q: measure = (4*np.pi/wavelength)*np.sin(measure/2)*(diameter/2) return measure,SL,SR,SU
def getFlux(self, slamb, sflux, axis=-1): """getFlux Integrate the flux within the filter and return the integrated energy If you consider applying the filter to many spectra, you might want to consider extractSEDs. Parameters ---------- slamb: ndarray(dtype=float, ndim=1) spectrum wavelength definition domain sflux: ndarray(dtype=float, ndim=1) associated flux Returns ------- flux: float Energy of the spectrum within the filter """ _wavelength = self._get_filter_in_units_of(slamb) _slamb = _drop_units(slamb) #clean data for inf values by interpolation if True in np.isinf(sflux): indinf = np.where(np.isinf(sflux)) indfin = np.where(np.isfinite(sflux)) sflux[indinf] = np.interp(_slamb[indinf], _slamb[indfin], sflux[indfin]) # reinterpolate transmission onto the same wavelength def as the data ifT = np.interp(_slamb, _wavelength, self.transmit, left=0., right=0.) # if the filter is null on that wavelength range flux is then 0 ind = ifT > 0. if True in ind: try: _sflux = sflux[:, ind] except: _sflux = sflux[ind] # limit integrals to where necessary ind = ifT > 0. if 'photon' in self.dtype: a = np.trapz(_slamb[ind] * ifT[ind] * _sflux, _slamb[ind], axis=axis) b = np.trapz(_slamb[ind] * ifT[ind], _slamb[ind] ) elif 'energy' in self.dtype: a = np.trapz( ifT[ind] * _sflux, _slamb[ind], axis=axis ) b = np.trapz( ifT[ind], _slamb[ind]) if (np.isinf(a).any() | np.isinf(b).any()): print(self.name, "Warn for inf value") return a / b else: return 0.
def filters1( model_nus, model_fluxes, filterdict, z ): """ Projects the model SEDs into the filter curves of each photometric band. ##input: - model_nus: template frequencies [log10(nu)] - model_fluxes: template fluxes [F_nu] - filterdict: dictionary with all band filter curves' information. To change this, add one band and filter curve, etc, look at DICTIONARIES_AGNfitter.py - z: redshift ##output: - bands [log10(nu)] - Filtered fluxes at these bands [F_nu] """ bands, files_dict, lambdas_dict, factors_dict = filterdict filtered_model_Fnus = [] # Costumize model frequencies and fluxes [F_nu] # to same units as filter curves (to wavelengths [angstrom] and F_lambda) model_lambdas = nu2lambda_angstrom(model_nus) * (1+z) model_lambdas = model_lambdas[::-1] model_fluxes_nu = model_fluxes[::-1] model_fluxes_lambda = fluxnu_2_fluxlambda(model_fluxes_nu, model_lambdas) mod2filter_interpol = interp1d(model_lambdas, model_fluxes_lambda, kind = 'nearest', bounds_error=False, fill_value=0.) # For filter curve at each band. # (Vectorised integration was not possible -> different filter-curve-arrays' sizes) for iband in bands: # Read filter curves info for each data point # (wavelengths [angstrom] and factors [non]) lambdas_filter = np.array(lambdas_dict[iband]) factors_filter = np.array(factors_dict[iband]) iband_angst = nu2lambda_angstrom(iband) # Interpolate the model fluxes to #the exact wavelengths of filter curves modelfluxes_at_filterlambdas = mod2filter_interpol(lambdas_filter) # Compute the flux ratios, equivalent to the filtered fluxes: # F = int(model)/int(filter) integral_model = trapz(modelfluxes_at_filterlambdas*factors_filter, x= lambdas_filter) integral_filter = trapz(factors_filter, x= lambdas_filter) filtered_modelF_lambda = (integral_model/integral_filter) # Convert all from lambda, F_lambda to Fnu and nu filtered_modelFnu_atfilter_i = fluxlambda_2_fluxnu(filtered_modelF_lambda, iband_angst) filtered_model_Fnus.append(filtered_modelFnu_atfilter_i) return bands, np.array(filtered_model_Fnus)
