python类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 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 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 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)