Python integrate.trapz方法代码示例

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import trapz [as 别名]
def l2(x, funca, funcb):
    """
    Computes the L2 norm for one profile(assumed to be linear between points).
    
    :type x: array of float
    :param x: depths :math:`[m]`.
    :type funca: list of arrays of float
    :param funca: array calculated values
    :type funcb: list of arrays of float
    :param funcb: array of values (observed or calculated) to compare to
    
    :returns: L2 norm
    :rtype: array of floats
    """
    diff = np.array(funca - funcb)
    diff = diff * diff
    
    return integrate.trapz(diff, x) 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import trapz [as 别名]
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) 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import trapz [as 别名]
def debye_func(x, nstep=100, zero=1e-8):
    r"""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 Exception("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) 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import trapz [as 别名]
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,dp)
  Babs = trapz(Q_abs*aSDn,dp)
  Bback = trapz(Q_back*aSDn,dp)

  return Bsca, Babs, Bback 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import trapz [as 别名]
def numeric_integation(func, n_samples=10 ** 5, bound_lower=-10**3, bound_upper=10**3):
  """ Numeric integration over one dimension using the trapezoidal rule

     Args:
       func: function to integrate over - must take numpy arrays of shape (n_samples,) as first argument
             and return a numpy array of shape (n_samples,)
       n_samples: (int) number of samples

     Returns:
       approximated integral - numpy array of shape (ndim_out,)
    """
  # proposal distribution
  y_samples = np.squeeze(np.linspace(bound_lower, bound_upper, num=n_samples))
  values = func(y_samples)
  integral = integrate.trapz(values, y_samples)
  return integral 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import trapz [as 别名]
def __init__(self, wavelength, transmit, name='', dtype="photon",
                 unit=None):
        """Constructor"""
        self.name = name
        self.set_dtype(dtype)
        try:   # get units from the inputs
            self._wavelength = wavelength.value
            unit = str(wavelength.unit)
        except AttributeError:
            self._wavelength = wavelength
        self.set_wavelength_unit(unit)
        self.transmit = np.clip(transmit, 0., np.nanmax(transmit))
        self.norm = trapz(self.transmit, self._wavelength)
        self._lT = trapz(self._wavelength * self.transmit, self._wavelength)
        self._lpivot = self._calculate_lpivot()
        if self.norm > 0:
            self._cl = self._lT / self.norm
        else:
            self._cl = 0. 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import trapz [as 别名]
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
            if s.transmit.max() > 0:
                leff = np.trapz(w * s.transmit * v.flux.value, w, axis=-1)
                leff /= np.trapz(s.transmit * v.flux.value, w, axis=-1)
            else:
                leff = float('nan')
        if s.wavelength_unit is not None:
            leff = leff * Unit(s.wavelength_unit)
            if self.wavelength_unit is not None:
                return leff.to(self.wavelength_unit)
            return leff
        else:
            return leff 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import trapz [as 别名]
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() 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import trapz [as 别名]
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 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import trapz [as 别名]
def Mie_SD(m, wavelength, dp, ndp, nMedium=1.0, interpolate=False, asDict=False):
#  http://pymiescatt.readthedocs.io/en/latest/forward.html#Mie_SD
  nMedium = nMedium.real
  m /= nMedium
  wavelength /= nMedium
  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],nMedium)

  Bext = trapz(Q_ext*aSDn,dp)
  Bsca = trapz(Q_sca*aSDn,dp)
  Babs = Bext-Bsca
  Bback = trapz(Q_back*aSDn,dp)
  Bratio = trapz(Q_ratio*aSDn,dp)
  bigG = trapz(g*Q_sca*aSDn,dp)/trapz(Q_sca*aSDn,dp)
  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 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import trapz [as 别名]
def SF_SD(m, wavelength, dp, ndp, nMedium=1.0, minAngle=0, maxAngle=180, angularResolution=0.5, space='theta', angleMeasure='radians', normalization=None):
#  http://pymiescatt.readthedocs.io/en/latest/forward.html#SF_SD
  nMedium = nMedium.real
  m /= nMedium
  wavelength /= nMedium
  
  _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 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import trapz [as 别名]
def _calculate_lpivot(self):
        if self.transmit.max() <= 0:
            return 0.
        if 'photon' in self.dtype:
            lpivot2 = self._lT / trapz(self.transmit / self._wavelength,
                                       self._wavelength)
        else:
            lpivot2 = self.norm / trapz(self.transmit / self._wavelength ** 2,
                                        self._wavelength)
        return np.sqrt(lpivot2) 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import trapz [as 别名]
def _get_indice(cls, w, flux, blue, red, band=None, unit='ew', degree=1,
                    **kwargs):
        """ compute spectral index after continuum subtraction

        Parameters
        ----------
        w: ndarray (nw, )
            array of wavelengths in AA
        flux: ndarray (N, nw)
            array of flux values for different spectra in the series
        blue: tuple(2)
            selection for blue continuum estimate
        red: tuple(2)
            selection for red continuum estimate
        band: tuple(2), optional
            select region in this band only.
            default is band = (min(blue), max(red))
        unit: str
            `ew` or `mag` wether equivalent width or magnitude
        degree: int (default 1)
            degree of the polynomial fit to the continuum

        Returns
        -------
        ew: ndarray (N,)
            equivalent width array
        """
        wi, fi = cls.continuum_normalized_region_around_line(w, flux, blue,
                                                             red, band=band,
                                                             degree=degree)
        if unit in (0, 'ew', 'EW'):
            return np.trapz(1. - fi, wi, axis=-1)
        else:
            m = np.trapz(fi, wi, axis=-1)
            m = -2.5 * np.log10(m / np.ptp(wi))
            return m 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import trapz [as 别名]
def __init__(self, wavelength, transmit, name='', dtype="photon", unit=None):
        """Constructor"""
        self.name       = name
        self.set_dtype(dtype)
        try:   # get units from the inputs
            self._wavelength = wavelength.magnitude
            unit = str(wavelength.units)
        except AttributeError:
            self._wavelength = wavelength
        self.set_wavelength_unit(unit)
        self.transmit   = np.clip(transmit, 0., np.nanmax(transmit))
        self.norm       = trapz(self.transmit, self._wavelength)
        self._lT        = trapz(self._wavelength * self.transmit, self._wavelength)
        self._lpivot    = self._calculate_lpivot()
        self._cl        = self._lT / self.norm 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import trapz [as 别名]
def _calculate_lpivot(self):
        if 'photon' in self.dtype:
            lpivot2 = self._lT / trapz(self.transmit / self._wavelength, self._wavelength)
        else:
            lpivot2 = self.norm / trapz(self.transmit / self._wavelength ** 2, self._wavelength)
        return np.sqrt(lpivot2)