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