Python numpy.trapz() Examples

The following are 30 code examples of numpy.trapz(). 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 numpy , or try the search function .
Example #1
Source File: qrnn.py    From typhon with MIT License 7 votes vote down vote up
def posterior_mean(self, x):
        r"""
        Computes the posterior mean by computing the first moment of the
        estimated posterior CDF.

        Arguments:

            x(np.array): Array of shape `(n, m)` containing `n` inputs for which
                         to predict the posterior mean.
        Returns:

            Array containing the posterior means for the provided inputs.
        """
        y_pred, qs = self.cdf(x)
        mus = y_pred[-1] - np.trapz(qs, x=y_pred)
        return mus 
Example #2
Source File: test_tlcorrection.py    From rampy with GNU General Public License v2.0 6 votes vote down vote up
def test_tlcorrection(self):
        # testing long correction function
        x_for_long = np.array([20.,21.,22.,23.,24.,25.])
        y_for_long = np.array([1.0,1.0,1.0,1.0,1.0,1.0])

        nu0 = 1.0/(514.532)*1e9 #laser wavenumber at 514.532
        nu = 100.0*x_for_long # cm-1 to m-1
        T = 23.0+273.15 # the temperature in K

        x_long,long_res,eselong = rp.tlcorrection(x_for_long, y_for_long,23.0,514.532,correction = 'long',normalisation='area') # using the function
        t0 = nu0**3.0*nu/((nu0-nu)**4)
        t1= 1.0 - np.exp(-h*c*nu/(k*T)) # c in m/s  : t1 dimensionless
        long_calc= y_for_long*t0*t1 # pour les y
        long_calc = long_calc/np.trapz(long_calc,x_for_long) # area normalisation

        np.testing.assert_equal(long_res,long_calc)
        np.testing.assert_equal(x_for_long,x_long)

        x_long,long_res,eselong = rp.tlcorrection(x_for_long, y_for_long,23.0,514.532,correction = 'long',normalisation='no') # using the function
        t0 = nu0**3.0*nu/((nu0-nu)**4)
        t1= 1.0 - np.exp(-h*c*nu/(k*T)) # c in m/s  : t1 dimensionless
        long_calc= y_for_long*t0*t1 # pour les y

        np.testing.assert_equal(long_res,long_calc)
        np.testing.assert_equal(x_for_long,x_long) 
Example #3
Source File: common.py    From typhon with MIT License 6 votes vote down vote up
def integrate_column(y, x=None, axis=0):
    """Integrate array along an arbitrary axis.

    Note:
        This function is just a wrapper for :func:`numpy.trapz`.

    Parameters:
        y (ndarray): Data array.
        x (ndarray): Coordinate array.
        axis (int): Axis to integrate along for multidimensional input.

    Returns:
        float or ndarray: Column integral.

    Examples:
        >>> import numpy as np
        >>> x = np.linspace(0, 1, 5)
        >>> y = np.arange(5)
        >>> integrate_column(y)
        8.0
        >>> integrate_column(y, x)
        2.0
    """
    return np.trapz(y, x, axis=axis) 
Example #4
Source File: test_soderman.py    From dmipy with MIT License 6 votes vote down vote up
def test_RTOP_to_diameter_soderman(samples=1000):
    """This tests if the RTAP of the sphere relates correctly to the diameter
    of the sphere."""
    diameter = 10e-6

    delta = np.tile(1e-10, samples)  # delta towards zero
    Delta = np.tile(1e10, samples)  # Delta towards infinity
    qvals_perp = np.linspace(0, 10e6, samples)
    n_perp = np.tile(np.r_[1., 0., 0.], (samples, 1))
    scheme = acquisition_scheme_from_qvalues(
        qvals_perp, n_perp, delta, Delta)

    soderman = sphere_models.S2SphereStejskalTannerApproximation(
        diameter=diameter)

    E_soderman = soderman(scheme)

    rtop_soderman = 4 * np.pi * np.trapz(
        E_soderman * qvals_perp ** 2, x=qvals_perp
    )

    sphere_volume = 1. / rtop_soderman
    diameter_soderman = 2 * (sphere_volume / ((4. / 3.) * np.pi)) ** (1. / 3.)

    assert_almost_equal(diameter_soderman, diameter, 7) 
Example #5
Source File: test_transformations.py    From pybids with MIT License 6 votes vote down vote up
def test_resample(collection):
    coll = collection.clone()

    transform.ToDense(coll, 'parametric gain', output='pg_dense')
    pg = coll.variables['pg_dense']
    old_shape = pg.values.shape
    old_auc = np.trapz(np.abs(pg.values.values.squeeze()), dx=0.1)
    transform.Resample(coll, 'pg_dense', 1)
    pg = coll.variables['pg_dense']
    new_shape = pg.values.shape
    # Spacing (dx) is 10* larger when downsampled fro 10hz to 1hz
    new_auc = np.trapz(np.abs(pg.values.values.squeeze()), dx=1)

    # Shape from 10hz to 1hz
    assert new_shape[0] == old_shape[0] / 10

    # Assert that the auc is more or less the same (not exact, rounding error)
    # Values are around 0.25
    assert np.allclose(old_auc, new_auc, rtol=0.05) 
Example #6
Source File: plot_hawkes.py    From tick with BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
def _normalize_functions(y_values_list, t_values):
    """Normalize list of functions by their integral value

    Parameters
    ----------
    y_values_list : `list` of np.ndarray
        y values of the list of function we want to normalize

    t_values : `np.ndarray`
        t values shared by all functions given with y_values_list

    Returns
    -------
    normalized_y_values_list : `list` of np.ndarray
        Normalized y values of the given list of function

    normalizations : `np.ndarray`
        Normalization factors that have been used
    """
    y_values_list = np.array(y_values_list)
    normalizations = [
        1. / np.trapz(y_values, t_values) for y_values in y_values_list
    ]
    normalized_y_values_list = (y_values_list.T * normalizations).T
    return normalized_y_values_list, normalizations 
Example #7
Source File: test_soderman.py    From dmipy with MIT License 6 votes vote down vote up
def test_RTPP_to_diameter_soderman(samples=1000):
    """This tests if the RTPP of the plane relates correctly to the diameter
    of the plane."""
    diameter = 10e-6

    delta = np.tile(1e-10, samples)  # delta towards zero
    Delta = np.tile(1e10, samples)  # Delta towards infinity
    qvals_perp = np.linspace(0, 10e6, samples)
    n_perp = np.tile(np.r_[1., 0., 0.], (samples, 1))
    scheme = acquisition_scheme_from_qvalues(
        qvals_perp, n_perp, delta, Delta)

    soderman = plane_models.P2PlaneStejskalTannerApproximation(
        diameter=diameter)

    E_soderman = soderman(scheme)

    rtpp_soderman = 2 * np.trapz(
        E_soderman, x=qvals_perp
    )

    diameter_soderman = 1 / rtpp_soderman

    assert_almost_equal(diameter_soderman, diameter, 7) 
Example #8
Source File: test_base.py    From pulse2percept with BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
def test_gamma():
    tsample = 0.005 / 1000

    with pytest.raises(ValueError):
        t, g = gamma(0, 0.1, tsample)
    with pytest.raises(ValueError):
        t, g = gamma(2, -0.1, tsample)
    with pytest.raises(ValueError):
        t, g = gamma(2, 0.1, -tsample)

    for tau in [0.001, 0.01, 0.1]:
        for n in [1, 2, 5]:
            t, g = gamma(n, tau, tsample)
            npt.assert_equal(np.arange(0, t[-1] + tsample / 2.0, tsample), t)
            if n > 1:
                npt.assert_equal(g[0], 0.0)

            # Make sure area under the curve is normalized
            npt.assert_almost_equal(np.trapz(np.abs(g), dx=tsample), 1.0,
                                    decimal=2)

            # Make sure peak sits correctly
            npt.assert_almost_equal(g.argmax() * tsample, tau * (n - 1)) 
Example #9
Source File: test_callaghan.py    From dmipy with MIT License 6 votes vote down vote up
def test_RTPP_to_length_callaghan(samples=1000):
    length = 10e-6

    delta = 1e-10  # delta towards zero
    Delta = 1e10  # Delta towards infinity
    qvals_perp = np.linspace(0, 10e6, samples)
    n_perp = np.tile(np.r_[1., 0., 0.], (samples, 1))
    scheme = acquisition_scheme_from_qvalues(qvals_perp, n_perp, delta, Delta)

    plane = plane_models.P3PlaneCallaghanApproximation(diameter=length)
    E_callaghan = plane(scheme)

    rtpp_callaghan = 2 * np.trapz(E_callaghan, x=qvals_perp)

    length_callaghan = 1 / rtpp_callaghan
    assert_almost_equal(length_callaghan, length, 7) 
Example #10
Source File: test_callaghan.py    From dmipy with MIT License 6 votes vote down vote up
def test_RTAP_to_diameter_callaghan(samples=10000):
    mu = [0, 0]
    lambda_par = 1.7
    diameter = 10e-6

    delta = np.tile(1e-10, samples)  # delta towards zero
    Delta = np.tile(1e10, samples)  # Delta towards infinity
    qvals_perp = np.linspace(0, 10e6, samples)
    n_perp = np.tile(np.r_[1., 0., 0.], (samples, 1))
    scheme = acquisition_scheme_from_qvalues(qvals_perp, n_perp, delta, Delta)

    callaghan = cylinder_models.C3CylinderCallaghanApproximation(
        mu=mu, lambda_par=lambda_par, diameter=diameter)

    E_callaghan = callaghan(scheme)

    rtap_callaghan = 2 * np.pi * np.trapz(
        E_callaghan * qvals_perp, x=qvals_perp
    )

    diameter_callaghan = 2 / np.sqrt(np.pi * rtap_callaghan)
    assert_almost_equal(diameter_callaghan, diameter, 7) 
Example #11
Source File: test_model.py    From pyhawkes with MIT License 6 votes vote down vote up
def test_generate_statistics():
    K = 1
    T = 100
    dt = 1.0
    network_hypers = {'c': np.zeros(K, dtype=np.int), 'p': 1.0, 'kappa': 10.0, 'v': 10*5.0}
    true_model = DiscreteTimeNetworkHawkesModelSpikeAndSlab(K=K, dt=dt,
                                                            network_hypers=network_hypers)
    S,R = true_model.generate(T=T)

    E_N = np.trapz(R, dt * np.arange(T), axis=0)
    std_N = np.sqrt(E_N)
    N = S.sum(axis=0)

    assert np.all(N >= E_N - 3*std_N), "N less than 3std below mean"
    assert np.all(N <= E_N + 3*std_N), "N more than 3std above mean"

    print("Expected number of events: ", E_N)
    print("Actual number of events:   ", S.sum(axis=0)) 
Example #12
Source File: genesis4.py    From ocelot with GNU General Public License v3.0 6 votes vote down vote up
def calc_spec(self, zi=None, loc='near', npad=1, estimate_ph_sp_dens=1):

        field = self.rad_field(zi=zi, loc=loc)
        axis = field.ndim - 1

        spec = np.abs(np.fft.fft(field, axis=axis))**2
        spec = np.fft.fftshift(spec, axes=axis)

        scale_ev = h_eV_s * speed_of_light * (np.fft.fftfreq(self.nSlices, d=self.s[1]-self.s[0]) + 1 / self.lambdaref)
        scale_ev = np.fft.fftshift(scale_ev)

        if estimate_ph_sp_dens:
            tt=np.trapz(spec, scale_ev, axis=axis)
            if axis==1:
                tt[tt==0] = np.inf
                spec *= (self.n_photons / tt)[:, np.newaxis]
            else:
                if tt==0:
                    tt = np.inf
                spec *= (self.n_photons[zi] / tt)

        return scale_ev, spec 
Example #13
Source File: metrics.py    From rainymotion with MIT License 6 votes vote down vote up
def AUC(x, y):
    '''
    AUC - area under curve

    Note: area under curve wich had been computed by standard trapezial
          method (np.trapz)

    Args:

        x (numpy.ndarray): array of one metric rate (1D)
        y (numpy.ndarray): array of another metric rate (1D)

    Returns:

        float - area under curve

    '''

    return np.trapz(y, x) 
Example #14
Source File: test_model.py    From pyhawkes with MIT License 6 votes vote down vote up
def test_compute_rate():
    K = 1
    T = 100
    dt = 1.0
    network_hypers = {'c': np.zeros(K, dtype=np.int), 'p': 1.0, 'kappa': 10.0, 'v': 10*5.0}
    true_model = DiscreteTimeNetworkHawkesModelSpikeAndSlab(K=K, dt=dt,
                                                            network_hypers=network_hypers)
    S,R = true_model.generate(T=T)

    print("Expected number of events: ", np.trapz(R, dt * np.arange(T), axis=0))
    print("Actual number of events:   ", S.sum(axis=0))

    print("Lambda0:  ", true_model.bias_model.lambda0)
    print("W:        ", true_model.weight_model.W)
    print("")

    R_test = true_model.compute_rate()
    assert np.allclose(R, R_test) 
Example #15
Source File: basis.py    From pyhawkes with MIT License 6 votes vote down vote up
def interpolate_basis(self, basis, dt, dt_max,
                          norm=True):
        # Interpolate basis at the resolution of the data
        L,B = basis.shape
        t_int = np.arange(0.0, dt_max, step=dt)
        t_bas = np.linspace(0.0, dt_max, L)

        ibasis = np.zeros((len(t_int), B))
        for b in np.arange(B):
            ibasis[:,b] = np.interp(t_int, t_bas, basis[:,b])

        # Normalize so that the interpolated basis has volume 1
        if norm:
            # ibasis /= np.trapz(ibasis,t_int,axis=0)
            ibasis /= (dt * np.sum(ibasis, axis=0))

        if not self.allow_instantaneous:
            # Typically, the impulse responses are applied to times
            # (t+1:t+R). That means we need to prepend a row of zeros to make
            # sure the basis remains causal
            ibasis = np.vstack((np.zeros((1,B)), ibasis))

        return ibasis 
Example #16
Source File: eval_utils.py    From siameseFC-pytorch-vot with Apache License 2.0 5 votes vote down vote up
def precisionAuc(positions, groundTruth, radius, nStep):
    thres = np.linspace(0, radius, nStep)
    errs = np.zeros([nStep], dtype=np.float32)

    distances = np.sqrt(np.power(positions[:, 0] - groundTruth[:, 0], 2) +
                        np.power(positions[:, 1] - groundTruth[:, 1], 2))
    distances[np.where(np.isnan(distances))] = []

    for p in range(nStep):
        errs[p] = np.shape(np.where(distances > thres[p]))[-1]
    score = np.trapz(errs)
    return score 
Example #17
Source File: stable_solution.py    From AeroPy with MIT License 5 votes vote down vote up
def work(self):
        u = self.g.x3(self.g.x1_grid)
        if self.load_type == 'concentrated':
            self.W = self.load*u[-1]
        elif self.load_type == 'distributed':
            self.W = np.trapz(self.load*u, self.g.x1_grid) 
Example #18
Source File: imdb2.py    From TF_Deformable_Net with MIT License 5 votes vote down vote up
def evaluate_recall(self, candidate_boxes, ar_thresh=0.5):
        # Record max overlap value for each gt box
        # Return vector of overlap values
        gt_overlaps = np.zeros(0)
        for i in range(self.num_images):
            gt_inds = np.where(self.roidb[i]['gt_classes'] > 0)[0]
            gt_boxes = self.roidb[i]['boxes'][gt_inds, :]

            boxes = candidate_boxes[i]
            if boxes.shape[0] == 0:
                continue
            overlaps = bbox_overlaps(boxes.astype(np.float),
                                     gt_boxes.astype(np.float))

            # gt_overlaps = np.hstack((gt_overlaps, overlaps.max(axis=0)))
            _gt_overlaps = np.zeros((gt_boxes.shape[0]))
            for j in range(gt_boxes.shape[0]):
                argmax_overlaps = overlaps.argmax(axis=0)
                max_overlaps = overlaps.max(axis=0)
                gt_ind = max_overlaps.argmax()
                gt_ovr = max_overlaps.max()
                assert(gt_ovr >= 0)
                box_ind = argmax_overlaps[gt_ind]
                _gt_overlaps[j] = overlaps[box_ind, gt_ind]
                assert(_gt_overlaps[j] == gt_ovr)
                overlaps[box_ind, :] = -1
                overlaps[:, gt_ind] = -1

            gt_overlaps = np.hstack((gt_overlaps, _gt_overlaps))

        num_pos = gt_overlaps.size
        gt_overlaps = np.sort(gt_overlaps)
        step = 0.001
        thresholds = np.minimum(np.arange(0.5, 1.0 + step, step), 1.0)
        recalls = np.zeros_like(thresholds)
        for i, t in enumerate(thresholds):
            recalls[i] = (gt_overlaps >= t).sum() / float(num_pos)
        ar = 2 * np.trapz(recalls, thresholds)

        return ar, gt_overlaps, recalls, thresholds 
Example #19
Source File: imdb2.py    From Faster-RCNN_TF with MIT License 5 votes vote down vote up
def evaluate_recall(self, candidate_boxes, ar_thresh=0.5):
        # Record max overlap value for each gt box
        # Return vector of overlap values
        gt_overlaps = np.zeros(0)
        for i in xrange(self.num_images):
            gt_inds = np.where(self.roidb[i]['gt_classes'] > 0)[0]
            gt_boxes = self.roidb[i]['boxes'][gt_inds, :]

            boxes = candidate_boxes[i]
            if boxes.shape[0] == 0:
                continue
            overlaps = bbox_overlaps(boxes.astype(np.float),
                                     gt_boxes.astype(np.float))

            # gt_overlaps = np.hstack((gt_overlaps, overlaps.max(axis=0)))
            _gt_overlaps = np.zeros((gt_boxes.shape[0]))
            for j in xrange(gt_boxes.shape[0]):
                argmax_overlaps = overlaps.argmax(axis=0)
                max_overlaps = overlaps.max(axis=0)
                gt_ind = max_overlaps.argmax()
                gt_ovr = max_overlaps.max()
                assert(gt_ovr >= 0)
                box_ind = argmax_overlaps[gt_ind]
                _gt_overlaps[j] = overlaps[box_ind, gt_ind]
                assert(_gt_overlaps[j] == gt_ovr)
                overlaps[box_ind, :] = -1
                overlaps[:, gt_ind] = -1

            gt_overlaps = np.hstack((gt_overlaps, _gt_overlaps))

        num_pos = gt_overlaps.size
        gt_overlaps = np.sort(gt_overlaps)
        step = 0.001
        thresholds = np.minimum(np.arange(0.5, 1.0 + step, step), 1.0)
        recalls = np.zeros_like(thresholds)
        for i, t in enumerate(thresholds):
            recalls[i] = (gt_overlaps >= t).sum() / float(num_pos)
        ar = 2 * np.trapz(recalls, thresholds)

        return ar, gt_overlaps, recalls, thresholds 
Example #20
Source File: hmetrics.py    From topmodel with MIT License 5 votes vote down vote up
def auc(fprs, tprs):
    xs = np.concatenate(([0], fprs[::-1], [1]))
    ys = np.concatenate(([0], tprs[::-1], [1]))
    return np.trapz(ys, xs) 
Example #21
Source File: test_interaction.py    From pySINDy with MIT License 5 votes vote down vote up
def test_trapz_matrix():
    # Test to make sure matrices give the same answer as ndarrays
    # 2018-04-29: moved here from core.tests.test_function_base.
    x = np.linspace(0, 5)
    y = x * x
    r = np.trapz(y, x)
    mx = np.matrix(x)
    my = np.matrix(y)
    mr = np.trapz(my, mx)
    assert_almost_equal(mr, r) 
Example #22
Source File: test_interaction.py    From predictive-maintenance-using-machine-learning with Apache License 2.0 5 votes vote down vote up
def test_trapz_matrix():
    # Test to make sure matrices give the same answer as ndarrays
    # 2018-04-29: moved here from core.tests.test_function_base.
    x = np.linspace(0, 5)
    y = x * x
    r = np.trapz(y, x)
    mx = np.matrix(x)
    my = np.matrix(y)
    mr = np.trapz(my, mx)
    assert_almost_equal(mr, r) 
Example #23
Source File: contactloss.py    From obman_train with GNU General Public License v3.0 5 votes vote down vote up
def meshiou(gt_dists, pred_dists, threshs=[1, 2, 3, 4, 5, 6, 7, 8, 9, 10]):
    """
    For each thresholds, computes thresh_ious and averages accross batch dim
    """
    all_ious = []
    for thresh in threshs:
        ious = thresh_ious(gt_dists, pred_dists, thresh)
        all_ious.append(ious)
    iou_auc = np.mean(
        np.trapz(torch.stack(all_ious).cpu().numpy(), axis=0, x=threshs)
    )
    batch_ious = torch.stack(all_ious).mean(1)
    return batch_ious, iou_auc 
Example #24
Source File: science_utils.py    From oopt-gnpy with BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
def _fast_generalized_psi(self, carrier_cut, pump_carrier, f_eval, f_cut_resolution):
        """ It computes the generalized psi function similarly to the one used in the GN model
        :return: generalized_psi
        """
        # Fiber parameters
        alpha0 = self.fiber.alpha0(f_eval)
        beta2 = self.fiber.params.beta2
        beta3 = self.fiber.params.beta3
        f_ref_beta = self.fiber.params.ref_frequency
        z = self.stimulated_raman_scattering.z
        frequency_rho = self.stimulated_raman_scattering.frequency
        rho_norm = self.stimulated_raman_scattering.rho * np.exp(np.abs(alpha0) * z / 2)
        if len(frequency_rho) == 1:
            def rho_function(f): return rho_norm[0, :]
        else:
            rho_function = interp1d(frequency_rho, rho_norm, axis=0, fill_value='extrapolate')
        rho_norm_pump = rho_function(pump_carrier.frequency)

        f1_array = np.array([pump_carrier.frequency - (pump_carrier.baud_rate * (1 + pump_carrier.roll_off) / 2),
                             pump_carrier.frequency + (pump_carrier.baud_rate * (1 + pump_carrier.roll_off) / 2)])
        f2_array = np.arange(carrier_cut.frequency,
                             carrier_cut.frequency + (carrier_cut.baud_rate * (1 + carrier_cut.roll_off) / 2),
                             f_cut_resolution)  # Only positive f2 is used since integrand_f2 is symmetric

        integrand_f1 = np.zeros(len(f1_array))
        for f1_index, f1 in enumerate(f1_array):
            delta_beta = 4 * np.pi**2 * (f1 - f_eval) * (f2_array - f_eval) * \
                (beta2 + np.pi * beta3 * (f1 + f2_array - 2 * f_ref_beta))
            integrand_f2 = self._generalized_rho_nli(delta_beta, rho_norm_pump, z, alpha0)
            integrand_f1[f1_index] = 2 * np.trapz(integrand_f2, f2_array)  # 2x since integrand_f2 is symmetric in f2
        generalized_psi = 0.5 * sum(integrand_f1) * pump_carrier.baud_rate
        return generalized_psi 
Example #25
Source File: background.py    From FRETBursts with GNU General Public License v2.0 5 votes vote down vote up
def _compute_error(residuals, x_residuals, error_metrics):
    assert error_metrics in ['KS', 'CM', None]
    error = None
    if error_metrics == 'KS':
        error = np.abs(residuals).max()*100
    elif error_metrics == 'CM':
        error = np.trapz(residuals**2, x=x_residuals)
    return error 
Example #26
Source File: burst_plot.py    From FRETBursts with GNU General Public License v2.0 5 votes vote down vote up
def _fitted_E_plot(d, i=0, F=1, no_E=False, ax=None, show_model=True,
                   verbose=False, two_gauss_model=False, lw=2.5, color='k',
                   alpha=0.5, fillcolor=None):
    """Plot a fitted model overlay on a FRET histogram."""
    if ax is None:
        ax2 = gca()
    else:
        ax2 = plt.twinx(ax=ax)
        ax2.grid(False)

    if d.fit_E_curve and show_model:
        x = r_[-0.2:1.21:0.002]
        y = d.fit_E_model(x, d.fit_E_res[i, :])
        scale = F*d.fit_E_model_F[i]
        if two_gauss_model:
            assert d.fit_E_res.shape[1] > 2
            if d.fit_E_res.shape[1] == 5:
                m1, s1, m2, s2, a1 = d.fit_E_res[i, :]
                a2 = (1-a1)
            elif d.fit_E_res.shape[1] == 6:
                m1, s1, a1, m2, s2, a2 = d.fit_E_res[i, :]
            y1 = a1*normpdf(x, m1, s1)
            y2 = a2*normpdf(x, m2, s2)
            ax2.plot(x, scale*y1, ls='--', lw=lw, alpha=alpha, color=color)
            ax2.plot(x, scale*y2, ls='--', lw=lw, alpha=alpha, color=color)
        if fillcolor is None:
            ax2.plot(x, scale*y, lw=lw, alpha=alpha, color=color)
        else:
            ax2.fill_between(x, scale*y, lw=lw, alpha=alpha, edgecolor=color,
                             facecolor=fillcolor, zorder=10)
        if verbose:
            print('Fit Integral:', np.trapz(scale*y, x))

    ax2.axvline(d.E_fit[i], lw=3, color=red, ls='--', alpha=0.6)
    xtext = 0.6 if d.E_fit[i] < 0.6 else 0.2
    if d.nch > 1 and not no_E:
        ax2.text(xtext, 0.81, "CH%d: $E_{fit} = %.3f$" % (i+1, d.E_fit[i]),
                 transform=gca().transAxes, fontsize=16,
                 bbox=dict(boxstyle='round', facecolor='#dedede', alpha=0.5)) 
Example #27
Source File: distributions.py    From dmipy with MIT License 5 votes vote down vote up
def __call__(self, **kwargs):
        r"""
        Parameters
        ----------
        diameter : float or array, shape (N)
            cylinder (axon) diameter in meters.

        Returns
        -------
        Pgamma : float or array, shape (N)
            probability of cylinder diameter for given alpha and beta.
        """
        alpha = kwargs.get('alpha', self.alpha)
        beta = kwargs.get('beta', self.beta)

        gamma_dist = stats.gamma(alpha, scale=beta)
        start_point = interpolate.bisplev(alpha, beta, self.start_interpolator)
        end_point = interpolate.bisplev(alpha, beta, self.end_interpolator)
        start_point = max(start_point, 1e-8)
        radii = np.linspace(start_point, end_point, self.Nsteps)
        normalization = self.norm_func(radii)
        radii_pdf = gamma_dist.pdf(radii)
        radii_pdf_area = radii_pdf * normalization
        radii_pdf_normalized = (
            radii_pdf_area /
            np.trapz(x=radii, y=radii_pdf_area)
        )
        return radii, radii_pdf_normalized 
Example #28
Source File: test_soderman.py    From dmipy with MIT License 5 votes vote down vote up
def test_RTAP_to_diameter_soderman(samples=1000):
    """This tests if the RTAP of the cylinder relates correctly to the diameter
    of the cylinder."""
    mu = [0, 0]
    lambda_par = 1.7
    diameter = 10e-6

    delta = np.tile(1e-10, samples)  # delta towards zero
    Delta = np.tile(1e10, samples)  # Delta towards infinity
    qvals_perp = np.linspace(0, 10e6, samples)
    n_perp = np.tile(np.r_[1., 0., 0.], (samples, 1))
    scheme = acquisition_scheme_from_qvalues(
        qvals_perp, n_perp, delta, Delta)

    soderman = cylinder_models.C2CylinderStejskalTannerApproximation(
        mu=mu, lambda_par=lambda_par, diameter=diameter)

    E_soderman = soderman(scheme)

    rtap_soderman = 2 * np.pi * np.trapz(
        E_soderman * qvals_perp, x=qvals_perp
    )

    diameter_soderman = 2 / np.sqrt(np.pi * rtap_soderman)

    assert_almost_equal(diameter_soderman, diameter, 7) 
Example #29
Source File: vectors.py    From modred with BSD 2-Clause "Simplified" License 5 votes vote down vote up
def inner_product(self, vec1, vec2):
        """Computes inner product."""
        IP = vec1 * vec2
        for grid in reversed(self.grids):
            if not isinstance(grid, np.ndarray):
                raise TypeError(
                    'Each grid must be a numpy array, not a %s' % type(grid))
            IP = np.trapz(IP, x=grid)
        return IP 
Example #30
Source File: test_metrics.py    From edm2016 with Apache License 2.0 5 votes vote down vote up
def test_auc_helper(self):
        """ Test the math meat of the AUC metric computation. """
        num_responses = 50
        num_trials = 100
        for trial in range(num_trials):
            # Create random response correctnesses
            correct_prob = 0.1 + 0.8 * np.random.rand()
            corrects = np.zeros(num_responses, dtype=bool)
            # Make sure there's at least 1 correct and 1 incorrect
            while np.sum(corrects) in (0, num_responses):
                corrects = np.random.rand(num_responses) < correct_prob
                incorrects = np.logical_not(corrects)

            # Create some random response probabilities
            rps = np.random.rand(num_responses)
            num_correct = float(np.sum(corrects))
            num_incorrect = float(np.sum(incorrects))

            # Compute AUC the slow way by iterating through thresholds
            tprs = np.zeros(len(rps) + 2)
            fprs = np.zeros(len(tprs))
            for i, threshold in enumerate(np.r_[-1., np.sort(rps), 2.]):
                tprs[i] = np.sum(np.logical_and(corrects, rps > threshold)) / num_correct
                fprs[i] = np.sum(np.logical_and(incorrects, rps > threshold)) / num_incorrect
            expected_auc = np.trapz(tprs[::-1], fprs[::-1])
            actual_auc = undertest.Metrics.auc_helper(corrects, rps)
            self.assertAlmostEqual(expected_auc, actual_auc)

            # Now compute some edge cases (all rps are 0 or all rps are 1)
            self.assertAlmostEqual(0.5, undertest.Metrics.auc_helper(corrects,
                                                                     np.zeros(num_responses)))
            self.assertAlmostEqual(0.5, undertest.Metrics.auc_helper(corrects,
                                                                     np.ones(num_responses)))

            # Now construct a case where a perfect threshold is possible
            sorted_rps = np.sort(rps)
            corrects = rps > sorted_rps[np.random.randint(1, num_responses-1)]
            self.assertAlmostEqual(1.0, undertest.Metrics.auc_helper(corrects, rps))