Python integrate.cumtrapz方法代码示例

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import cumtrapz [as 别名]
def get_velocity_displacement(time_step, acceleration, units="cm/s/s",
                              velocity=None, displacement=None):
    '''
    Returns the velocity and displacment time series using simple integration
    :param float time_step:
        Time-series time-step (s)
    :param numpy.ndarray acceleration:
        Acceleration time-history
    :returns:
        velocity - Velocity Time series (cm/s)
        displacement - Displacement Time series (cm)
    '''
    acceleration = convert_accel_units(acceleration, units)
    if velocity is None:
        velocity = time_step * cumtrapz(acceleration, initial=0.)
    if displacement is None:
        displacement = time_step * cumtrapz(velocity, initial=0.)
    return velocity, displacement 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import cumtrapz [as 别名]
def K0_fin_anf_opt(self, i, traj, wmin, gamma):
        s = np.zeros(i)
        n = np.zeros((i, 3))
        R = np.zeros(i)
        w = np.zeros(i)
        self.K0_0(i, traj[0], traj[1], traj[2], traj[3], gamma, s, n, R, w)
        j = np.where(w <= wmin)[0]

        if len(j) > 0:
            j = j[-1]
            w = w[j:i]
            s = s[j:i]
        else:
            j = 0
        K = self.K0_1(i, j, R, n, traj[4], traj[5], traj[6], w, gamma)

        if len(K) > 1:
            a = np.append(0.5 * (K[0:-1] + K[1:]) * np.diff(s), 0.5 * K[-1] * s[-1])
            KS = np.cumsum(a[::-1])[::-1]
            # KS = cumsum_inv_jit(a)
            # KS = cumtrapz(K[::-1], -s[::-1], initial=0)[::-1] + 0.5*K[-1]*s[-1]
        else:
            KS = 0.5 * K[-1] * s[-1]
        return w, KS 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import cumtrapz [as 别名]
def test_nd(self):
        x = np.arange(3 * 2 * 4).reshape(3, 2, 4)
        y = x
        y_int = cumtrapz(y, x, initial=0)
        y_expected = np.array([[[0., 0.5, 2., 4.5],
                                [0., 4.5, 10., 16.5]],
                               [[0., 8.5, 18., 28.5],
                                [0., 12.5, 26., 40.5]],
                               [[0., 16.5, 34., 52.5],
                                [0., 20.5, 42., 64.5]]])

        assert_allclose(y_int, y_expected)

        # Try with all axes
        shapes = [(2, 2, 4), (3, 1, 4), (3, 2, 3)]
        for axis, shape in zip([0, 1, 2], shapes):
            y_int = cumtrapz(y, x, initial=3.45, axis=axis)
            assert_equal(y_int.shape, (3, 2, 4))
            y_int = cumtrapz(y, x, initial=None, axis=axis)
            assert_equal(y_int.shape, shape) 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import cumtrapz [as 别名]
def test_x_none(self):
        y = np.linspace(-2, 2, num=5)

        y_int = cumtrapz(y)
        y_expected = [-1.5, -2., -1.5, 0.]
        assert_allclose(y_int, y_expected)

        y_int = cumtrapz(y, initial=1.23)
        y_expected = [1.23, -1.5, -2., -1.5, 0.]
        assert_allclose(y_int, y_expected)

        y_int = cumtrapz(y, dx=3)
        y_expected = [-4.5, -6., -4.5, 0.]
        assert_allclose(y_int, y_expected)

        y_int = cumtrapz(y, dx=3, initial=1.23)
        y_expected = [1.23, -4.5, -6., -4.5, 0.]
        assert_allclose(y_int, y_expected) 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import cumtrapz [as 别名]
def from_file(self, filename="REPORT"):
        """
        Reads values from files and stores it in the `parse_dict` attribute

        Args:
            filename (str): Path to the file that needs to be parsed
        """
        with open(filename, "r") as f:
            lines = f.readlines()
        rel_lines = [lines[i + 2] for i, line in enumerate(lines) if "Blue_moon" in line]
        if len(rel_lines) > 0:
            [lam, _, _, _] = [val for val in np.genfromtxt(rel_lines, usecols=[1, 2, 3, 4]).T]
            rel_lines = [lines[i] for i, line in enumerate(lines) if "cc>" in line]
            cv = np.genfromtxt(rel_lines, usecols=[2])
            fe = cumtrapz(lam, cv)
            self.parse_dict["cv_full"] = cv
            self.parse_dict["derivative"] = lam
            self.parse_dict["cv"] = cv[:-1]
            self.parse_dict["free_energy"] = fe 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import cumtrapz [as 别名]
def set_filling(h,filling=0.5,nk=10,extrae=0.,delta=1e-1):
    """Set the filling of a Hamiltonian"""
    if h.has_eh: raise
    fill = filling + extrae/h.intra.shape[0] # filling
    n = h.intra.shape[0]
    use_kpm = False
    if n>algebra.maxsize: # use the KPM method
        use_kpm = True
        print("Using KPM in set_filling")
    if use_kpm: # use KPM
        es,ds = h.get_dos(energies=np.linspace(-5.0,5.0,1000),
                use_kpm=True,delta=delta,nk=nk,random=False)
        from scipy.integrate import cumtrapz
        di = cumtrapz(ds,es)
        ei = (es[0:len(es)-1] + es[1:len(es)])/2.
        di /= di[len(di)-1] # normalize
        from scipy.interpolate import interp1d
        f = interp1d(di,ei) # interpolating function
        efermi = f(fill) # get the fermi energy
    else: # dense Hamiltonian, use ED
        es = eigenvalues(h,nk=nk)
        efermi = get_fermi_energy(es,fill)
    h.shift_fermi(-efermi) # shift the fermi energy 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import cumtrapz [as 别名]
def chern_density(h,nk=10,operator=None,delta=0.02,dk=0.02,
        es=np.linspace(-1.0,1.0,40)):
  """Compute the Chern density as a function of the energy"""
  ks = klist.kmesh(h.dimensionality,nk=nk)
  cs = np.zeros(es.shape[0]) # initialize
  f = h.get_gk_gen(delta=delta,canonical_phase=True) # green function generator
  tr = timing.Testimator("CHERN DENSITY",maxite=len(ks))
  from . import parallel
  def fp(k): # compute berry curvatures
    if parallel.cores==1: tr.iterate()
    else: print(k)
#    k = np.random.random(3)
    fgreen = berry_green_generator(f,k=k,dk=dk,operator=operator,full=False) 
    return np.array([fgreen(e).real for e in es])
  out = parallel.pcall(fp,ks) # compute everything
  for o in out: cs += o # add contributions
  cs = cs/(len(ks)*np.pi*2) # normalize
  from scipy.integrate import cumtrapz
  csi = cumtrapz(cs,x=es,initial=0) # integrate
  np.savetxt("CHERN_DENSITY.OUT",np.matrix([es,cs]).T)
  np.savetxt("CHERN_DENSITY_INTEGRATED.OUT",np.matrix([es,csi]).T) 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import cumtrapz [as 别名]
def get_peak_measures(time_step, acceleration, get_vel=False, 
    get_disp=False):
    """
    Returns the peak measures from acceleration, velocity and displacement
    time-series
    :param float time_step:
        Time step of acceleration time series in s
    :param numpy.ndarray acceleration:
        Acceleration time series
    :param bool get_vel:
        Choose to return (and therefore calculate) velocity (True) or otherwise
        (false)
    :returns:
        * pga - Peak Ground Acceleration
        * pgv - Peak Ground Velocity
        * pgd - Peak Ground Displacement
        * velocity - Velocity Time Series
        * dispalcement - Displacement Time series
    """
    pga = np.max(np.fabs(acceleration))
    velocity = None
    displacement = None
    # If displacement is not required then do not integrate to get
    # displacement time series
    if get_disp:
        get_vel = True
    if get_vel:
        velocity = time_step * cumtrapz(acceleration, initial=0.)
        pgv = np.max(np.fabs(velocity))
    else:
        pgv = None
    if get_disp:
        displacement = time_step * cumtrapz(velocity, initial=0.)
        pgd = np.max(np.fabs(displacement))
    else:
        pgd = None
    return pga, pgv, pgd, velocity, displacement 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import cumtrapz [as 别名]
def get_husid(acceleration, time_step):
    """
    Returns the Husid vector, defined as \int{acceleration ** 2.}
    :param numpy.ndarray acceleration:
        Vector of acceleration values
    :param float time_step:
        Time-step of record (s)
    """
    time_vector = get_time_vector(time_step, len(acceleration))
    husid = np.hstack([0., cumtrapz(acceleration ** 2., time_vector)])
    return husid, time_vector 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import cumtrapz [as 别名]
def s(self, n=0):
        self.check(n)

        xp2 = self.xp(n) ** 2
        yp2 = self.yp(n) ** 2
        # zp = np.sqrt(1. / (1. + xp2 + yp2))
        s = cumtrapz(np.sqrt(1. + xp2 + yp2), self.z(n), initial=0)
        return s 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import cumtrapz [as 别名]
def invert_cdf(y, x):
    """
    Invert cumulative distribution function of the probability distribution

    Example:
    --------
    # analytical formula for the beam distribution
    f = lambda x: A * np.exp(-(x - mu) ** 2 / (2. * sigma ** 2))

    # we are interesting in range from -30 to 30 e.g. [um]
    x = np.linspace(-30, 30, num=100)

    # Inverted cumulative distribution function
    i_cdf = invert_cdf(y=f(x), x=x)

    # get beam distribution (200 000 coordinates)
    tau = i_cdf(np.random.rand(200000))


    :param y: array, [y0, y1, y2, ... yn] yi = y(xi)
    :param x: array, [x0, x1, x2, ... xn] xi
    :return: function
    """

    cum_int = integrate.cumtrapz(y, x, initial=0)
    #print("invert", np.max(cum_int), np.min(cum_int))
    cdf = cum_int / (np.max(cum_int) - np.min(cum_int))
    inv_cdf = interpolate.interp1d(cdf, x)
    return inv_cdf 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import cumtrapz [as 别名]
def test_1d(self):
        x = np.linspace(-2, 2, num=5)
        y = x
        y_int = cumtrapz(y, x, initial=0)
        y_expected = [0., -1.5, -2., -1.5, 0.]
        assert_allclose(y_int, y_expected)

        y_int = cumtrapz(y, x, initial=None)
        assert_allclose(y_int, y_expected[1:]) 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import cumtrapz [as 别名]
def compute_cdf_percentiles(fit, cdf_sigmas=CDF_SIGMAS):
    """
    Compute tabulated percentiles of the PDF
    """
    from scipy.interpolate import Akima1DInterpolator
    from scipy.integrate import cumtrapz
    import scipy.stats

    if cdf_sigmas is None:
        cdf_sigmas = CDF_SIGMAS

    cdf_y = scipy.stats.norm.cdf(cdf_sigmas)

    if len(fit['zgrid']) == 1:
        return np.ones_like(cdf_y)*fit['zgrid'][0], cdf_y

    spl = Akima1DInterpolator(fit['zgrid'], np.log(fit['pdf']), axis=1)
    zrfine = [fit['zgrid'].min(), fit['zgrid'].max()]
    zfine = utils.log_zgrid(zr=zrfine, dz=0.0001)
    ok = np.isfinite(spl(zfine))
    pz_fine = np.exp(spl(zfine))
    pz_fine[~ok] = 0

    cdf_fine = cumtrapz(pz_fine, x=zfine)
    cdf_x = np.interp(cdf_y, cdf_fine/cdf_fine[-1], zfine[1:])

    return cdf_x, cdf_y 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import cumtrapz [as 别名]
def get_traces(self):

        fn = self.config.get_output_filename(self.tempdir)
        data = num.loadtxt(fn, skiprows=1, dtype=num.float)
        nsamples, ntraces = data.shape
        deltat = (data[-1, 0] - data[0, 0]) / (nsamples - 1)
        toffset = data[0, 0]

        tred = self.config.time_reduction
        rec = self.config.receiver
        tmin = rec.tstart + toffset + deltat + tred

        traces = []
        for itrace, comp in enumerate(self.config.components):
            # qseis2d gives velocity-integrate to displacement
            # integration removes one sample, add it again in front
            displ = cumtrapz(num.concatenate(
                (num.zeros(1), data[:, itrace + 1])), dx=deltat)

            tr = trace.Trace(
                '', '%04i' % itrace, '', comp,
                tmin=tmin, deltat=deltat, ydata=displ,
                meta=dict(distance=rec.distance,
                          azimuth=0.0))

            traces.append(tr)

        return traces 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import cumtrapz [as 别名]
def _int_spontaneous_raman(self, z_array, raman_matrix, alphap_fiber, freq_array,
                               cr_raman_matrix, freq_diff, ase_bc, bn_array, temperature):
        spontaneous_raman_scattering = OptimizeResult()

        simulation = Simulation.get_simulation()
        sim_params = simulation.sim_params

        dx = sim_params.raman_params.space_resolution
        h = ph.value('Planck constant')
        kb = ph.value('Boltzmann constant')

        power_ase = np.nan * np.ones(raman_matrix.shape)
        int_pump = cumtrapz(raman_matrix, z_array, dx=dx, axis=1, initial=0)

        for f_ind, f_ase in enumerate(freq_array):
            cr_raman = cr_raman_matrix[f_ind, :]
            vibrational_loss = f_ase / freq_array[:f_ind]
            eta = 1 / (np.exp((h * freq_diff[f_ind, f_ind + 1:]) / (kb * temperature)) - 1)

            int_fiber_loss = -alphap_fiber[f_ind] * z_array
            int_raman_loss = np.sum((cr_raman[:f_ind] * vibrational_loss * int_pump[:f_ind, :].transpose()).transpose(),
                                    axis=0)
            int_raman_gain = np.sum((cr_raman[f_ind + 1:] * int_pump[f_ind + 1:, :].transpose()).transpose(), axis=0)

            int_gain_loss = int_fiber_loss + int_raman_gain + int_raman_loss

            new_ase = np.sum((cr_raman[f_ind + 1:] * (1 + eta) * raman_matrix[f_ind + 1:, :].transpose()).transpose()
                             * h * f_ase * bn_array[f_ind], axis=0)

            bc_evolution = ase_bc[f_ind] * np.exp(int_gain_loss)
            ase_evolution = np.exp(int_gain_loss) * cumtrapz(new_ase *
                                                             np.exp(-int_gain_loss), z_array, dx=dx, initial=0)

            power_ase[f_ind, :] = bc_evolution + ase_evolution

        spontaneous_raman_scattering.x = 2 * power_ase
        return spontaneous_raman_scattering