Python scipy.integrate.simps() Examples

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 

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 

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) 

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) 

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) 

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) 

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) 

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

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) 

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 _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 

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 

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 

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 

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)])))) 

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 

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) 

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

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

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

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] 

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 

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 

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] 

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() 

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) 

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) 

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 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 

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)