Python numpy.trapz函数代码示例

def run(data, sweep_rate=20.,
        c_range=(0.4, 0.6), co_range=(0.6, 0.9),
        exe='', graph=True, baseline=False, copy=False):
    # INIT data
    params = {}
    area_CO = None
    area_H = None

    cycle_CO = data.get_scan(1)
    cycle_baseline = data.get_scan(2)

    if data.sweep_rate:
        sr = data.sweep_rate
    else:
        sr = sweep_rate

    # RUN stuff
    if "CO" in exe:
        params["CO"] = CO(cycle_CO, cycle_baseline, c_range, co_range, add_baseline=baseline, copy=copy)
        x, y = params["CO"][0]
        Q_CO = np.trapz(y, x)  # V C / s
        factor_CO = 420e-6 * sr  # C V / s cm2
        area_CO = Q_CO / factor_CO  # cm2
    if "H" in exe:
        params["H"] = H(cycle_CO, cycle_baseline, co_range[0], copy)
        x, y = params["H"]
        Q_H = np.trapz(y, x)  # V C / s
        factor_H = 210e-6 * sr * 1.e-3  # C V / s cm2
        area_H = Q_H / factor_H  # cm2
    if graph:
        plot(cycle_CO, cycle_baseline, paramsCO=params.get("CO"),
             paramsH=params.get("H"), exe=exe, graph=graph)

    return area_CO, area_H

 def __init__(self, fct, makenorm=True, makeunlog=False, fct_log=False,
              npts=1000, *args, **kwargs):
     if not callable(fct):
         raise Exception("Not callable")
     self.fct = fct
     self.args = args
     self.kwargs = kwargs
     # apply exp on a ln-function at callback or not
     self.makeunlog = makeunlog is True
     # renorm at callback or not
     self.makenorm = makenorm is True
     # whether to add or divide the normalization
     self.fct_log = fct_log is not False and makeunlog is False
     self.prior_log = kwargs.get('prior_log', False)  # jeffrey
     if kwargs.get('prior_bounds', None) is None:
         raise Exception("No prior bounds found")
     if self.prior_log:
         self.x = np.logspace(np.log(kwargs.get('prior_bounds')[0]),
                              np.log(kwargs.get('prior_bounds')[1]),
                              npts,
                              base=np.e)
     else:
         self.x = np.linspace(kwargs.get('prior_bounds')[0],
                              kwargs.get('prior_bounds')[1],
                              npts)
     if makenorm is True:
         resfct = self.fct(self.x, *self.args, **self.kwargs)
         if self.fct_log:
             self.normval = np.log(np.trapz(np.exp(resfct), self.x))
         elif self.makeunlog:
             self.normval = np.trapz(np.exp(resfct), self.x)
         else:
             self.normval = np.trapz(resfct, self.x)

    def beam_phase_sharpWindow(self):
        '''
        *Beam phase measured at the main RF frequency and phase. The beam is 
        averaged over a window. The coefficients of sine and cosine components
        determine the 
        beam phase, projected to the range -Pi/2 to 3/2 Pi. Note that this beam
        phase is already w.r.t. the instantaneous RF phase.*
        '''    
        
        # Main RF frequency at the present turn
        omega_rf = self.rf_station.omega_rf[0,self.rf_station.counter[0]]
        phi_rf = self.rf_station.phi_rf[0,self.rf_station.counter[0]]
        
        if self.alpha != 0.0:
            left_boundary = self.profile.bin_centers \
                >= self.time_offset - np.pi/omega_rf
            right_boundary = self.profile.bin_centers \
                <= -1/self.alpha + self.time_offset - 2*np.pi/omega_rf
            indexes = left_boundary * right_boundary
        else:
            indexes = np.ones(self.profile.n_slices, dtype=bool)
        
        # Convolve with window function
        scoeff = np.trapz( np.sin(omega_rf*self.profile.bin_centers[indexes]\
                                   + phi_rf) \
                           *self.profile.n_macroparticles[indexes],
                           dx=self.profile.bin_size )
        ccoeff = np.trapz( np.cos(omega_rf*self.profile.bin_centers[indexes]\
                                   + phi_rf) \
                           *self.profile.n_macroparticles[indexes],
                           dx=self.profile.bin_size )

        # Project beam phase to (pi/2,3pi/2) range
        self.phi_beam = np.arctan(scoeff/ccoeff) + np.pi

def avgField(time, field, tstart=None, std=False):
    """
    This subroutine computes the time-average (and the std) of a time series

    >>> ts = MagicTs(field='misc', iplot=False, all=True)
    >>> nuavg = avgField(ts.time, ts.topnuss, 0.35)
    >>> print(nuavg)

    :param time: time
    :type time: numpy.ndarray
    :param field: the time series of a given field
    :type field: numpy.ndarray
    :param tstart: the starting time of the averaging
    :type tstart: float
    :param std: when set to True, the standard deviation is also calculated
    :type std: bool
    :returns: the time-averaged quantity
    :rtype: float
    """
    if tstart is not None:
        mask = np.where(abs(time - tstart) == min(abs(time - tstart)), 1, 0)
        ind = np.nonzero(mask)[0][0]
    else:  # the whole input array is taken!
        ind = 0
    fac = 1.0 / (time[-1] - time[ind])
    avgField = fac * np.trapz(field[ind:], time[ind:])

    if std:
        stdField = np.sqrt(fac * np.trapz((field[ind:] - avgField) ** 2, time[ind:]))
        return avgField, stdField
    else:
        return avgField

def test_spectra_get_flux_contributions():
    timestepmin = 40
    timestepmax = 80
    dfspectrum = at.spectra.get_spectrum(
        specfilename, timestepmin=timestepmin, timestepmax=timestepmax, fnufilterfunc=None)

    integrated_flux_specout = np.trapz(dfspectrum['f_lambda'], x=dfspectrum['lambda_angstroms'])

    specdata = pd.read_csv(specfilename, delim_whitespace=True)
    arraynu = specdata.loc[:, '0'].values
    arraylambda_angstroms = const.c.to('angstrom/s').value / arraynu

    contribution_list, array_flambda_emission_total = at.spectra.get_flux_contributions(
        modelpath, timestepmin=timestepmin, timestepmax=timestepmax)

    integrated_flux_emission = -np.trapz(array_flambda_emission_total, x=arraylambda_angstroms)

    # total spectrum should be equal to the sum of all emission processes
    print(f'Integrated flux from spec.out:     {integrated_flux_specout}')
    print(f'Integrated flux from emission sum: {integrated_flux_emission}')
    assert math.isclose(integrated_flux_specout, integrated_flux_emission, rel_tol=4e-3)

    # check each bin is not out by a large fraction
    diff = [abs(x - y) for x, y in zip(array_flambda_emission_total, dfspectrum['f_lambda'].values)]
    print(f'Max f_lambda difference {max(diff) / integrated_flux_specout}')
    assert max(diff) / integrated_flux_specout < 2e-3

def test_AR_LD():
    """

    Test the Levinson Durbin estimate of the AR coefficients against the
    expercted PSD

    """
    arsig,_,_ = utils.ar_generator(N=512)
    avg_pwr = (arsig*arsig.conjugate()).real.mean()
    order = 8
    ak, sigma_v = tsa.AR_est_LD(arsig, order)
    w, psd = tsa.AR_psd(ak, sigma_v)

    # the psd is a one-sided power spectral density, which has been
    # multiplied by 2 to preserve the property that
    # 1/2pi int_{-pi}^{pi} Sxx(w) dw = Rxx(0)

    # evaluate this integral numerically from 0 to pi
    dw = np.pi/len(psd)
    avg_pwr_est = np.trapz(psd, dx=dw) / (2*np.pi)
    npt.assert_almost_equal(avg_pwr, avg_pwr_est, decimal=0)

    # Test for providing the autocovariance as an input:
    ak,sigma_v = tsa.AR_est_LD(arsig, order, utils.autocov(arsig))
    w, psd = tsa.AR_psd(ak, sigma_v)
    avg_pwr_est = np.trapz(psd, dx=dw) / (2*np.pi)
    npt.assert_almost_equal(avg_pwr, avg_pwr_est, decimal=0)

def get_band_phot( spectrum , bandpass, z, a=1, zero_point=48.6 ):
    import numpy as np

    c = (3.0e18)      #  Angs / s

  #  These conditionals import text files if necessary

    if isinstance( spectrum , type('string') ): spectrum = np.loadtxt(spectrum)
    if isinstance( bandpass , type('string') ): bandpass = np.loadtxt(bandpass)

    interp_flux = np.interp(bandpass[:,0], spectrum[:,0]*(1.0+z),
        spectrum[:,a] )

    weighted_flux = []
    nu_arr = []
    weights = []

    nu_arr = c / bandpass[:,0]
    weighted_flux = interp_flux * bandpass[:,1] * bandpass[:,0]/nu_arr/nu_arr
    weights = bandpass[:,1] / nu_arr

    numer = np.trapz( weighted_flux , x=nu_arr )
    denom = np.trapz( weights , x=nu_arr)

    #return -2.5*np.log10( numer / denom ) + zero_point   # <--- "MAGNITUDE"
    return -2.5*np.log10(numer/denom/zero_point)# - zero_point   # <--- "MAGNITUDE"

    def solve_nonlinear(self, params, unknowns, resids):
        '''a VERY coarse approximation that takes the speed profile given in the original proposal as given. It just serves as a place holder for now. It's better than nothing, but a real analysis is needed here'''
        t1 = (300 - 0) * 1609 / (9.81 * 0.5) / 3600 # time needed to accelerate
        t2 = (555 - 300) * 1609 / (9.81 * 0.5) / 3600 # time needed to accelerate
        t3 = (params['max_velocity'] - 555) * 1609 / (9.81 * 0.5) / 3600 # time needed to accelerate
        #speed profile data from hyperloop alpha proposal, pg43
        dataStart = np.array([
            [0, 0],
            [t1, 300 * 1.609], # (300[mi/h] * 1609[m/mi]) / (9.81[m/s] * 0.5) / 3600[s/h]
            [167, 300 * 1.609],
            [167 + t2, 555 * 1.609], # (555 - 300[mi/h] * 1609[m/mi]) / (9.81[m/s] * 0.5) / 3600[s/h]
            [435, 555 * 1.609],
            [435 + t3, params['max_velocity']]])
        startUp = np.trapz(dataStart[:, 1] / 3600, dataStart[:, 0]) * 1000 # km covered during start up Los Angeles Grapevine
        dataEnd = np.array([
            [0, params['max_velocity']],
            [t3, 555 * 1.609],
            [t3 + 100, 555 * 1.609],
            [t3 + 100 + t2, 300 * 1.609],
            [t3 + 100 + t2 + 400, 300 * 1.609],
            [t3 + 100 + t2 + 400 + t1, 0]])
        windDown = np.trapz(dataEnd[:, 1] / 3600, dataEnd[:, 0]) * 1000 # km covered during wind down along I-580 to SF      
        len_middle = params['tube_length'] - (startUp + windDown)
        time_middle = middleLength / params['max_velocity']
        unknowns['time_mission'] = time_middle + 435 + t3 + t3 + 100 + t2 + 400 + t1

        unknowns['energy'] = (params['pwr_req'] * unknowns['time_mission'] / 3600.0) * (1 + params['pwr_marg']) # convert to hours

def energy(powers, default_timestep=8):
    """Compute the energy associated to a list of measures (in W)
    and associated timestamps (in s).
    """
    energy = {"night_rate": 0, "day_rate": 0, "value": 0}
    if len(powers) == 1:
        if powers[0].night_rate == 1:
            energy["night_rate"] = powers[0].value / 1000 * default_timestep / 3600
        else:
            energy["day_rate"] = powers[0].value / 1000 * default_timestep / 3600
        energy["value"] = energy["day_rate"] + energy["night_rate"]
    else:
        x = []
        day_rate = []
        night_rate = []
        for i in powers:
            x.append(i.timestamp)
            if i.night_rate == 1:
                night_rate.append(i.value)
                day_rate.append(0)
            else:
                day_rate.append(i.value)
                night_rate.append(0)
        energy["night_rate"] = numpy.trapz(night_rate, x) / 1000 / 3600
        energy["day_rate"] = numpy.trapz(day_rate, x) / 1000 / 3600
        energy["value"] = energy["day_rate"] + energy["night_rate"]
    return energy

    def set_angle_width(self, width):
        # Calculates the normalizing coefficient for the distribution. From
        # Bringi & Chandrasekar 2001, pg. 68 (which cite Mardia 1972) for an
        # axial distribution.
        self._angle_width = width
        if width <= 0.:
            return
        t = np.linspace(0, 1, 500)
        self._angle_b = 1. / (2. * np.trapz(
            np.exp(-self._angle_width * t**2), t))

        # We neglect here the factor of 1/2pi, as without it:
        # * The graphs on page 69 can be reproduced correctly
        # * The distribution is normalized (sums to 1)
        # Size 90 below *must* match the number returned from T-matrix
        self._angle_vals = np.linspace(np.pi/2, np.pi, TMATRIX_ANGLES)
        self._angle_weights = self._angle_b * (np.exp(
            -self._angle_width * np.cos(self._angle_vals)**2)
            * np.sin(self._angle_vals))

        # Since the distribution and scattering calculations are symmetric, we
        # just calculate the right half of the distribution and double the
        # weights
        self._angle_weights[1:] *= 2
        norm_factor = np.trapz(self._angle_weights, self._angle_vals)
        self._angle_weights /= norm_factor

def integrate_3D(x, y, z, xlin, ylin, zlin, csd):
    # TODO: NOT YET WORKING AS EXPECTED
    X, Y, Z = np.meshgrid(xlin, ylin, zlin)
    Nz = zlin.shape[0]
    Ny = ylin.shape[0]

    m = np.sqrt((x - xlin)**2 + (y - ylin)**2 + (z - zlin)**2)
    m[m < 0.00001] = 0.00001
    csd = csd / m
    #print 'CSD:'
    #print csd

    J = np.zeros((Ny, Nz))
    for i in xrange(Nz):
        J[:, i] = np.trapz(csd[:, :, i], zlin)

    #print '1st integration'
    #print J

    Ny = ylin.shape[0]
    I = np.zeros(Ny)
    for i in xrange(Ny):
        I[i] = np.trapz(J[:, i], ylin)

    #print '2nd integration'
    #print I

    norm = np.trapz(I, xlin)
    
    #print '3rd integration'
    #print norm
    
    return norm

    def effective_wavelength(self, binned=True, wavelengths=None,
                             mode='efflerg'):
        """Calculate :ref:`effective wavelength <synphot-formula-effwave>`.

        Parameters
        ----------
        binned : bool
            Sample data in native wavelengths if `False`.
            Else, sample binned data (default).

        wavelengths : array-like, `~astropy.units.quantity.Quantity`, or `None`
            Wavelength values for sampling.
            If not a Quantity, assumed to be in Angstrom.
            If `None`, ``self.waveset`` or `binset` is used, depending
            on ``binned``.

        mode : {'efflerg', 'efflphot'}
            Flux is first converted to the unit below before calculation:

                * 'efflerg' - FLAM
                * 'efflphot' - PHOTLAM (deprecated)

        Returns
        -------
        eff_lam : `~astropy.units.quantity.Quantity`
            Observation effective wavelength.

        Raises
        ------
        synphot.exceptions.SynphotError
            Invalid mode.

        """
        mode = mode.lower()
        if mode == 'efflerg':
            flux_unit = units.FLAM
        elif mode == 'efflphot':
            warnings.warn(
                'Usage of EFFLPHOT is deprecated.', AstropyDeprecationWarning)
            flux_unit = units.PHOTLAM
        else:
            raise exceptions.SynphotError(
                'mode must be "efflerg" or "efflphot"')

        if binned:
            x = self._validate_binned_wavelengths(wavelengths).value
            y = self.sample_binned(wavelengths=x, flux_unit=flux_unit).value
        else:
            x = self._validate_wavelengths(wavelengths).value
            y = units.convert_flux(x, self(x), flux_unit).value

        num = np.trapz(y * x ** 2, x=x)
        den = np.trapz(y * x, x=x)

        if den == 0.0:  # pragma: no cover
            eff_lam = 0.0
        else:
            eff_lam = abs(num / den)

        return eff_lam * self._internal_wave_unit

    def get_cracking_history():
        XK = []  # position of the first crack
        sig_c_K = [0.]
        eps_c_K = [0.]

        idx_0 = np.argmin(sig_mu_x)
        XK.append(x[idx_0])
        sig_c_0 = sig_mu_x[idx_0] * Ec / Em
        sig_c_K.append(sig_c_0)
        eps_c_K.append(sig_mu_x[idx_0] / Em)

        while True:
            z_x = get_z_x(x, XK)
            sig_c_k, y_i = get_sig_c_K(z_x, sig_c_K[-1])
            if sig_c_k == sig_cu:
                break
            XK.append(y_i)
            sig_c_K.append(sig_c_k)
            eps_c_K.append(
                np.trapz(eps_f(get_z_x(x, XK), sig_c_k), x) / np.amax(x))  # Eq. (10)

            # save the figure
#             plt.figure()
#             plt.plot(x, sig_m(get_z_x(x, XK), sig_c_k))
#             plt.plot(x, sig_mu_x)
#             plt.savefig("D:\\cracking_history\\" + str(len(sig_c_K)) + ".png")
#             plt.close()

        sig_c_K.append(sig_cu)
        eps_c_K.append(np.trapz(eps_f(get_z_x(x, XK), sig_cu), x) / np.amax(x))
        return sig_c_K, eps_c_K

def dDCR_moments(SED1, SED2, bandpass):
    zenith_angle = np.pi/4.0 * galsim.radians
    R500 = galsim.dcr.get_refraction(500, zenith_angle)

    # analytic first moment differences
    R = lambda w:(galsim.dcr.get_refraction(w, zenith_angle) - R500) / galsim.arcsec
    x1 = np.union1d(bandpass.wave_list, SED1.wave_list)
    x1 = x1[(x1 >= bandpass.blue_limit) & (x1 <= bandpass.red_limit)]
    x2 = np.union1d(bandpass.wave_list, SED2.wave_list)
    x2 = x2[(x2 >= bandpass.blue_limit) & (x2 <= bandpass.red_limit)]
    numR1 = np.trapz(R(x1) * bandpass(x1) * SED1(x1), x1)
    numR2 = np.trapz(R(x2) * bandpass(x2) * SED2(x2), x2)
    den1 = SED1.calculateFlux(bandpass)
    den2 = SED2.calculateFlux(bandpass)

    R1 = numR1/den1
    R2 = numR2/den2
    dR_analytic = R1 - R2

    # analytic second moment differences
    V1_kernel = lambda w:(R(w) - R1)**2
    V2_kernel = lambda w:(R(w) - R2)**2
    numV1 = np.trapz(V1_kernel(x1) * bandpass(x1) * SED1(x1), x1)
    numV2 = np.trapz(V2_kernel(x2) * bandpass(x2) * SED2(x2), x2)
    V1 = numV1/den1
    V2 = numV2/den2
    dV_analytic = V1 - V2

    return dR_analytic, dV_analytic, len(x2)

	def solve( self, t, D, VdotO20 ):
		self.Time_Grid = t
		self.Demand = D
		self.VdotO2aerobic = self.solveAerobic( VdotO20 )
		self.VdotO2anaerobic = self.solveAnaerobic()
		self.VO2aerobic = trapz( self.VdotO2aerobic, x=self.Time_Grid )
		self.VO2anaerobic = trapz( self.VdotO2anaerobic, x=self.Time_Grid )