Python scipy.integrate.simps() Examples

The following are 30 code examples of scipy.integrate.simps(). You can vote up the ones you like or vote down the ones you don't like, and go to the original project or source file by following the links above each example. You may also want to check out all available functions/classes of the module scipy.integrate , or try the search function .
Example #1
Source File: acc_utils.py    From ocelot with GNU General Public License v3.0 8 votes vote down vote up
def slice_bunching(tau, charge, lambda_mod, smooth_sigma=None):
    """
    Function calculates bunching factor for wavelength lambda_mod

    $b(\lambda) = \frac{1}{N_0}\left| \langle e^{- i\frac{ 2 \pi}{\lambda} s} N(s)\rangle \right|$

    :param p_array: ParticleArray
    :param lambda_mod: wavelength
    :param smooth_sigma: smoothing parameter
    :return: bunching factor
    """
    if smooth_sigma is None:
        smooth_sigma = lambda_mod/10.

    B = s_to_cur(tau, sigma=smooth_sigma, q0=charge, v=speed_of_light)

    b = np.abs(simps(B[:, 1] / speed_of_light * np.exp(-1j * 2 * np.pi / lambda_mod * B[:, 0]), B[:, 0])) / charge
    return b 
Example #2
Source File: acc_utils.py    From ocelot with GNU General Public License v3.0 7 votes vote down vote up
def bunching(p_array, lambda_mod, smooth_sigma=None):
    """
    Function calculates bunching factor for wavelength lambda_mod

    $b(\lambda) = \frac{1}{N_0}\left| \langle e^{- i\frac{ 2 \pi}{\lambda} s} N(s)\rangle \right|$

    :param p_array: ParticleArray
    :param lambda_mod: wavelength
    :param smooth_sigma: smoothing parameter
    :return: bunching factor
    """
    if smooth_sigma is None:
        smooth_sigma = min(np.std(p_array.tau())*0.01, lambda_mod*0.1)

    B = s_to_cur(p_array.tau(), sigma=smooth_sigma, q0=np.sum(p_array.q_array), v=speed_of_light)

    b = np.abs(simps(B[:, 1] / speed_of_light * np.exp(-1j * 2 * np.pi / lambda_mod * B[:, 0]), B[:, 0])) / np.sum(
        p_array.q_array)
    return b 
Example #3
Source File: states.py    From strawberryfields with Apache License 2.0 7 votes vote down vote up
def p_quad_values(self, mode, xvec, pvec):

        r"""Calculates the discretized p-quadrature probability distribution of the specified mode.

        Args:
            mode (int): the mode to calculate the p-quadrature probability values of
            xvec (array): array of discretized :math:`x` quadrature values
            pvec (array): array of discretized :math:`p` quadrature values

        Returns:
            array: 1D array of size len(pvec), containing reduced p-quadrature
            probability values for a specified range of x and p.
        """

        W = self.wigner(mode, xvec, pvec)
        y = []
        for i in range(0, len(pvec)):
            res = simps(W[i, : len(xvec)], xvec)
            y.append(res)
        return np.array(y) 
Example #4
Source File: aerodynamic_pdf.py    From AeroPy with MIT License 6 votes vote down vote up
def expected(data, airFrame):
    alpha, V, lift_to_drag = data

    pdf = airFrame.pdf.score_samples(np.vstack([alpha.ravel(), V.ravel()]).T)
    pdf = np.exp(pdf.reshape(lift_to_drag.shape))
    expected_value = 0
    numerator_list = []
    denominator_list = []
    for i in range(len(lift_to_drag)):
        numerator = simps(lift_to_drag[i]*pdf[i], alpha[i])
        denominator = simps(pdf[i], alpha[i])
        numerator_list.append(numerator)
        denominator_list.append(denominator)
    numerator = simps(numerator_list, V[:, 0])
    denominator = simps(denominator_list, V[:, 0])
    expected_value = numerator/denominator
    return(expected_value) 
Example #5
Source File: states.py    From strawberryfields with Apache License 2.0 6 votes vote down vote up
def x_quad_values(self, mode, xvec, pvec):

        r"""Calculates the discretized x-quadrature probability distribution of the specified mode.

        Args:
            mode (int): the mode to calculate the x-quadrature probability values of
            xvec (array): array of discretized :math:`x` quadrature values
            pvec (array): array of discretized :math:`p` quadrature values

        Returns:
            array: 1D array of size len(xvec), containing reduced x-quadrature
            probability values for a specified range of x and p.
        """

        W = self.wigner(mode, xvec, pvec)
        y = []
        for i in range(0, len(xvec)):
            res = simps(W[: len(pvec), i], pvec)
            y.append(res)
        return np.array(y) 
Example #6
Source File: expected_value_MonteCarlos.py    From AeroPy with MIT License 6 votes vote down vote up
def expected(data, airFrame):
    alpha, V, lift_to_drag = data

    pdf = airFrame.pdf.score_samples(np.vstack([alpha.ravel(), V.ravel()]).T)
    pdf = np.exp(pdf.reshape(lift_to_drag.shape))
    expected_value = 0
    numerator_list = []
    denominator_list = []
    for i in range(len(lift_to_drag)):
        numerator = simps(lift_to_drag[i]*pdf[i], alpha[i])
        denominator = simps(pdf[i], alpha[i])
        numerator_list.append(numerator)
        denominator_list.append(denominator)
    numerator = simps(numerator_list, V[:, 0])
    denominator = simps(denominator_list, V[:, 0])
    expected_value = numerator/denominator
    return(expected_value) 
Example #7
Source File: expected_value.py    From AeroPy with MIT License 6 votes vote down vote up
def expected(data, airFrame):
    alpha, V, lift_to_drag = data

    pdf = airFrame.pdf.score_samples(np.vstack([alpha.ravel(), V.ravel()]).T)
    pdf = np.exp(pdf.reshape(lift_to_drag.shape))
    expected_value = 0
    numerator_list = []
    denominator_list = []
    for i in range(len(lift_to_drag)):
        numerator = simps(lift_to_drag[i]*pdf[i], alpha[i])
        denominator = simps(pdf[i], alpha[i])
        numerator_list.append(numerator)
        denominator_list.append(denominator)
    numerator = simps(numerator_list, V[:, 0])
    denominator = simps(denominator_list, V[:, 0])
    expected_value = numerator/denominator
    return(expected_value) 
Example #8
Source File: test_polya_gamma_identity.py    From pgmult with MIT License 6 votes vote down vote up
def simps_approx():
    # Compute the left hand side analytically
    loglhs = psi*a - b * np.log1p(np.exp(psi))
    lhs = np.exp(loglhs)

    # Compute the right hand side with quadrature
    from scipy.integrate import simps
    # Lay down a grid of omegas
    # TODO: How should we choose the bins?
    omegas = np.linspace(1e-15, 5, 1000)
    # integrand = lambda om: 2**-b \
    #                        * np.exp((a-b/2.)*psi 0 0.5*om*psi**2) \
    #                        * polya_gamma_density(om, b, 0)
    # y = map(integrand, omegas)
    # rhs = simps(integrand, y)
    logy = -b * np.log(2) + (a-b/2.) * psi -0.5*omegas*psi**2
    logy += log_polya_gamma_density(omegas, b, 0, trunc=21)
    y = np.exp(logy)
    rhs = simps(y, omegas)

    print("Numerical Quadrature")
    print("log LHS: ", loglhs)
    print("log RHS: ", np.log(rhs)) 
Example #9
Source File: test_polya_gamma_identity.py    From pgmult with MIT License 6 votes vote down vote up
def plot_density():
    omegas = np.linspace(1e-16, 5, 1000)
    logpomega = log_polya_gamma_density(omegas, b, 0, trunc=1000)
    pomega = np.exp(logpomega).real

    import matplotlib.pyplot as plt
    plt.ion()
    plt.figure()
    plt.plot(omegas, pomega)

    y = -b * np.log(2) + (a-b/2.) * psi -0.5*omegas*psi**2


    from scipy.integrate import simps
    Z = simps(pomega, omegas)
    print("Z: ", Z) 
Example #10
Source File: thermo.py    From pwtools with BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
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) 
Example #11
Source File: analysis.py    From ChainConsumer with MIT License 5 votes vote down vote up
def _get_smoothed_histogram(self, chain, parameter):
        data = chain.get_data(parameter)
        smooth = chain.config["smooth"]
        if chain.grid:
            bins = get_grid_bins(data)
        else:
            bins = chain.config["bins"]
            bins, smooth = get_smoothed_bins(smooth, bins, data, chain.weights)

        hist, edges = np.histogram(data, bins=bins, density=True, weights=chain.weights)
        if chain.power is not None:
            hist = hist ** chain.power
        edge_centers = 0.5 * (edges[1:] + edges[:-1])
        xs = np.linspace(edge_centers[0], edge_centers[-1], 10000)

        if smooth:
            hist = gaussian_filter(hist, smooth, mode=self.parent._gauss_mode)
        kde = chain.config["kde"]
        if kde:
            kde_xs = np.linspace(edge_centers[0], edge_centers[-1], max(200, int(bins.max())))
            ys = MegKDE(data, chain.weights, factor=kde).evaluate(kde_xs)
            area = simps(ys, x=kde_xs)
            ys = ys / area
            ys = interp1d(kde_xs, ys, kind="linear")(xs)
        else:
            ys = interp1d(edge_centers, hist, kind="linear")(xs)
        cs = ys.cumsum()
        cs /= cs.max()
        return xs, ys, cs 
Example #12
Source File: thermal_properties.py    From DynaPhoPy with MIT License 5 votes vote down vote up
def get_free_energy_correction_dos(temperature, frequency, dos, dos_r):

    def n(temp, freq):
        return pow(np.exp(freq*h_bar/(k_b*temp))-1, -1)

    free_energy_1 = np.nan_to_num([ dos_r[i] * -h_bar/2 * freq*(n(temperature, freq) + 1 / 2.)
                                      for i, freq in enumerate(frequency)])

    free_energy_2 = np.nan_to_num([ dos[i] * -h_bar/2 * freq*(n(temperature, freq) + 1 / 2.)
                                      for i, freq in enumerate(frequency)])

    free_energy_c = free_energy_1 - free_energy_2

    free_energy_c = integrate.simps(free_energy_c, frequency) * N_a / 1000 # KJ/K/mol
    return free_energy_c 
Example #13
Source File: util.py    From recurrent-slds with MIT License 5 votes vote down vote up
def compute_psi_cmoments(alphas):
    K = alphas.shape[0]
    psi = np.linspace(-10,10,1000)

    mu = np.zeros(K-1)
    sigma = np.zeros(K-1)
    for k in range(K-1):
        density = get_density(alphas[k], alphas[k+1:].sum())
        mu[k] = simps(psi*density(psi),psi)
        sigma[k] = simps(psi**2*density(psi),psi) - mu[k]**2

    return mu, sigma 
Example #14
Source File: utils.py    From pgmult with MIT License 5 votes vote down vote up
def compute_psi_cmoments(alphas):
    K = alphas.shape[0]
    psi = np.linspace(-10,10,1000)

    mu = np.zeros(K-1)
    sigma = np.zeros(K-1)
    for k in range(K-1):
        density = get_density(alphas[k], alphas[k+1:].sum())
        mu[k] = simps(psi*density(psi),psi)
        sigma[k] = simps(psi**2*density(psi),psi) - mu[k]**2
        # print '%d: mean=%0.3f var=%0.3f' % (k, mean, s - mean**2)

    return mu, sigma 
Example #15
Source File: utils.py    From pgmult with MIT License 5 votes vote down vote up
def dirichlet_to_psi_meanvar(alphas,psigrid=np.linspace(-10,10,1000)):
    K = alphas.shape[0]

    def meanvar(k):
        density = get_marginal_psi_density(alphas[k], alphas[k+1:].sum())
        mean = simps(psigrid*density(psigrid),psigrid)
        s = simps(psigrid**2*density(psigrid),psigrid)
        return mean, s - mean**2

    return list(map(np.array, list(zip(*[meanvar(k) for k in range(K-1)])))) 
Example #16
Source File: thermal_properties.py    From DynaPhoPy with MIT License 5 votes vote down vote up
def get_entropy2(temperature, frequency, dos):

    def n(temp, freq):
        return pow(np.exp(freq*h_bar/(k_b*temp))-1, -1)

    entropy = np.nan_to_num([dos[i] * k_b * ((n(temperature, freq) + 1) * np.log(n(temperature, freq) + 1)
                                             - n(temperature, freq) * np.log(n(temperature, freq)))
                         for i, freq in enumerate(frequency)])
    entropy = integrate.simps(entropy, frequency) * N_a # J/K/mol
    return entropy 
Example #17
Source File: thermal_properties.py    From DynaPhoPy with MIT License 5 votes vote down vote up
def get_entropy(temperature, frequency, dos):

    def coth(x):
        return  np.cosh(x)/np.sinh(x)

    entropy = np.nan_to_num([dos[i]*(1.0 / (2. * temperature) * h_bar * freq * coth(h_bar * freq / (2 * k_b * temperature))
                                     - k_b * np.log(2 * np.sinh(h_bar * freq / (2 * k_b * temperature))))
                             for i, freq in enumerate(frequency)])
    entropy = integrate.simps(entropy, frequency) * N_a # J/K/mol
    return entropy

# Alternative way to calculate entropy (not used) 
Example #18
Source File: __init__.py    From DynaPhoPy with MIT License 5 votes vote down vote up
def plot_power_spectrum_wave_vector(self):
        plt.suptitle('Projection onto wave vector')
        plt.plot(self.get_frequency_range(), self.get_power_spectrum_wave_vector(), 'r-')
        plt.xlabel('Frequency [THz]')
        plt.ylabel('eV * ps')
        plt.axhline(y=0, color='k', ls='dashed')
        plt.show()
        total_integral = integrate.simps(self.get_power_spectrum_wave_vector(), x=self.get_frequency_range())
        print("Total Area (Kinetic energy <K>): {0} eV".format(total_integral)) 
Example #19
Source File: __init__.py    From DynaPhoPy with MIT License 5 votes vote down vote up
def write_power_spectrum_full(self, file_name):
        reading.write_curve_to_file(self.get_frequency_range(),
                                    self.get_power_spectrum_full()[None].T,
                                    file_name)
        total_integral = integrate.simps(self.get_power_spectrum_full(), x=self.get_frequency_range())
        print("Total Area (Kinetic energy <K>): {0} eV".format(total_integral)) 
Example #20
Source File: __init__.py    From DynaPhoPy with MIT License 5 votes vote down vote up
def write_power_spectrum_wave_vector(self, file_name):
        reading.write_curve_to_file(self.get_frequency_range(),
                                    self.get_power_spectrum_wave_vector()[None].T,
                                    file_name)
        total_integral = integrate.simps(self.get_power_spectrum_wave_vector(), x=self.get_frequency_range())
        print("Total Area (Kinetic energy <K>): {0} eV".format(total_integral)) 
Example #21
Source File: __init__.py    From DynaPhoPy with MIT License 5 votes vote down vote up
def get_thermal_properties(self, force_constants=None):

        temperature = self.get_temperature()

        phonopy_dos = pho_interface.obtain_phonopy_dos(self.dynamic.structure,
                                                       mesh=self.parameters.mesh_phonopy,
                                                       freq_min=0.01,
                                                       freq_max=self.get_frequency_range()[-1],
                                                       force_constants=force_constants,
                                                       projected_on_atom=self.parameters.project_on_atom,
                                                       NAC=self.parameters.use_NAC)

        free_energy = thm.get_free_energy(temperature, phonopy_dos[0], phonopy_dos[1])
        entropy = thm.get_entropy(temperature, phonopy_dos[0], phonopy_dos[1])
        c_v = thm.get_cv(temperature, phonopy_dos[0], phonopy_dos[1])
        integration = integrate.simps(phonopy_dos[1], x=phonopy_dos[0]) / (
            self.dynamic.structure.get_number_of_atoms() *
            self.dynamic.structure.get_number_of_dimensions())
        total_energy = thm.get_total_energy(temperature, phonopy_dos[0], phonopy_dos[1])

        if force_constants is not None:
            # Free energy correction
            phonopy_dos_h = pho_interface.obtain_phonopy_dos(self.dynamic.structure,
                                                             mesh=self.parameters.mesh_phonopy,
                                                             freq_min=0.01,
                                                             freq_max=self.get_frequency_range()[-1],
                                                             projected_on_atom=self.parameters.project_on_atom,
                                                             NAC=self.parameters.use_NAC)

            free_energy += thm.get_free_energy_correction_dos(temperature, phonopy_dos_h[0], phonopy_dos_h[1],
                                                              phonopy_dos[1])
            total_energy += thm.get_free_energy_correction_dos(temperature, phonopy_dos_h[0], phonopy_dos_h[1],
                                                               phonopy_dos[1])

            # correction = thm.get_free_energy_correction_dos(temperature, phonopy_dos_h[0], phonopy_dos_h[1], phonopy_dos[1])
            # print('Free energy/total energy correction: {0:12.4f} KJ/mol'.format(correction))

        return [free_energy, entropy, c_v, total_energy, integration] 
Example #22
Source File: thermal_properties.py    From DynaPhoPy with MIT License 5 votes vote down vote up
def get_free_energy(temperature, frequency, dos):

    free_energy = np.nan_to_num([dos[i] * k_b * temperature * np.log(2 * np.sinh(h_bar * freq / (2 * k_b * temperature)))
                                 for i, freq in enumerate(frequency)])

    free_energy[0] = 0
    free_energy = integrate.simps(free_energy, frequency) * N_a / 1000  # KJ/K/mol
    return free_energy 
Example #23
Source File: thermal_properties.py    From DynaPhoPy with MIT License 5 votes vote down vote up
def get_free_energy_correction_shift(temperature, frequency, dos, shift):

    def n(temp, freq):
        return pow(np.exp(freq*h_bar/(k_b*temp))-1, -1)

    free_energy_c = np.nan_to_num([dos[i] * -h_bar/2 *shift*(n(temperature, freq) + 1 / 2.)
                                   for i, freq in enumerate(frequency)])

    free_energy_c = integrate.simps(free_energy_c, frequency) * N_a / 1000 # KJ/K/mol
    return free_energy_c 
Example #24
Source File: _analogsignalarray.py    From nelpy with MIT License 5 votes vote down vote up
def _pdf(self, bins=None, n_samples=None):
        """Return the probability distribution function for each signal."""
        from scipy import integrate

        if bins is None:
            bins = 100

        if n_samples is None:
            n_samples = 100

        if self.n_signals > 1:
            raise NotImplementedError('multiple signals not supported yet!')

        # fx, bins = np.histogram(self.data.squeeze(), bins=bins, normed=True)
        fx, bins = np.histogram(self.data.squeeze(), bins=bins)
        bin_centers = (bins + (bins[1]-bins[0])/2)[:-1]

        Ifx = integrate.simps(fx, bin_centers)

        pdf = type(self)(abscissa_vals=bin_centers,
                                data=fx/Ifx,
                                fs=1/(bin_centers[1]-bin_centers[0]),
                                support=type(self.support)(self.data.min(), self.data.max())
                               ).simplify(n_samples=n_samples)

        return pdf

        # data = []
        # for signal in self.data:
        #     fx, bins = np.histogram(signal, bins=bins)
        #     bin_centers = (bins + (bins[1]-bins[0])/2)[:-1] 
Example #25
Source File: tests.py    From DeepAlignmentNetwork with MIT License 5 votes vote down vote up
def AUCError(errors, failureThreshold, step=0.0001, showCurve=False):
    nErrors = len(errors)
    xAxis = list(np.arange(0., failureThreshold + step, step))

    ced =  [float(np.count_nonzero([errors <= x])) / nErrors for x in xAxis]

    AUC = simps(ced, x=xAxis) / failureThreshold
    failureRate = 1. - ced[-1]

    print "AUC @ {0}: {1}".format(failureThreshold, AUC)
    print "Failure rate: {0}".format(failureRate)

    if showCurve:
        plt.plot(xAxis, ced)
        plt.show() 
Example #26
Source File: test_distributions.py    From GraphicDesignPatternByPython with MIT License 5 votes vote down vote up
def test_pdf_unity_area(self):
        from scipy.integrate import simps
        # PDF should integrate to one
        p = stats.genexpon.pdf(numpy.arange(0, 10, 0.01), 0.5, 0.5, 2.0)
        assert_almost_equal(simps(p, dx=0.01), 1, 1) 
Example #27
Source File: test_norm_int.py    From pwtools with BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
def test_norm_int():
    # simps(y, x) == 2.0
    x = np.linspace(0, np.pi, 100)
    y = np.sin(x)

    for scale in [True, False]:
        yy = norm_int(y, x, area=10.0, scale=scale)
        assert np.allclose(simps(yy,x), 10.0) 
Example #28
Source File: num.py    From pwtools with BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
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 
Example #29
Source File: electro.py    From VASPy with MIT License 5 votes vote down vote up
def get_dband_center(self, d_cols):
        """
        Get d-band center of the DosX object.

        Parameters:
        -----------
        d_cols: The column number range for d orbitals, int or tuple of int.

        Examples:
        ---------
        # The 5 - 9 columns are state density for d orbitals.
        >>> dos.get_dband_center(d_cols=(5, 10))
        """
        d_cols = (d_cols, d_cols+1) if type(d_cols) is int else d_cols

        # 合并d轨道DOS
        start, end = d_cols
        yd = np.sum(self.data[:, start:end], axis=1)

        #获取feimi能级索引
        for idx, E in enumerate(self.data[:, 0]):
            if E >= 0:
                nfermi = idx
                break
        E = self.data[: nfermi+1, 0]  # negative inf to Fermi
        dos = yd[: nfermi+1]          # y values from negative inf to Fermi
        # Use Simpson integration to get d-electron number
        nelectro = simps(dos, E)
        # Get total energy of dband
        tot_E = simps(E*dos, E)
        dband_center = tot_E/nelectro
        self.dband_center = dband_center

        return dband_center 
Example #30
Source File: test_distributions.py    From Computable with MIT License 5 votes vote down vote up
def test_pdf_unity_area(self):
        from scipy.integrate import simps
        # PDF should integrate to one
        assert_almost_equal(simps(stats.genexpon.pdf(numpy.arange(0,10,0.01),
                                                     0.5, 0.5, 2.0),
                                  dx=0.01), 1, 1)