我们从Python开源项目中，提取了以下49个代码示例，用于说明如何使用numpy.trapz()。
def integrate(self, attrname, radio, x_box, y_box): element = radio.labels[radio.active] source_local = getattr(self, attrname+"_"+element+"_source") lower_xlim = float(x_box.value) lower_ylim = float(y_box.value) x = np.array(source_local.data["x"]) y = np.array(source_local.data["y"]) x_change = x[x>lower_xlim]*1e-4 y_change = y[len(y)-len(x_change):] integral = np.trapz(y_change, x = x_change) comparsion = np.sum(y) * x[-1]*1e-4 print(integral, comparsion)
def band_area(spectrum, low_endmember=None, high_endmember=None): """ Compute the area under the curve between two endpoints where the y-value <= 1. """ x = spectrum.index y = spectrum if not low_endmember: low_endmember = x[0] if not high_endmember: high_endmember = x[-1] ny = y[low_endmember:high_endmember] return np.trapz(-ny[ny <= 1.0])
def add_filtered_planck(self, wavelength, trans): vals = np.zeros(self.x.shape, 'float64') nu = clight/(wavelength*1e-8) nu = nu[::-1] for i,logT in enumerate(self.x): T = 10**logT # Black body at this nu, T Bnu = ((2.0 * hcgs * nu**3) / clight**2.0) / \ (np.exp(hcgs * nu / (kboltz * T)) - 1.0) # transmission f = Bnu * trans[::-1] # integrate transmission over nu vals[i] = np.trapz(f,nu) # normalize by total transmission over filter self.y = vals/trans.sum() #/np.trapz(trans[::-1],nu) #self.y = np.clip(np.maximum(vals, self.y), 0.0, 1.0)
def convolve(self, wavelengths, densities): # define short names for the involved wavelength grids wa = wavelengths wb = self._Wavelengths # create a combined wavelength grid, restricted to the overlapping interval w1 = wa[ (wa>=wb[0]) & (wa<=wb[-1]) ] w2 = wb[ (wb>=wa[0]) & (wb<=wa[-1]) ] w = np.unique(np.hstack((w1,w2))) if len(w) < 2: return 0 # log-log interpolate SED and transmission on the combined wavelength grid # (use scipy interpolation function for SED because np.interp does not support broadcasting) F = np.exp(interp1d(np.log(wa), _log(densities), copy=False, bounds_error=False, fill_value=0.)(np.log(w))) T = np.exp(np.interp(np.log(w), np.log(wb), _log(self._Transmission), left=0., right=0.)) # perform the integration if self._PhotonCounter: return np.trapz(x=w, y=w*F*T) / self._IntegratedTransmission else: return np.trapz(x=w, y=F*T) / self._IntegratedTransmission ## This function calculates and returns the integrated value for a given spectral energy distribution over the # filter's wavelength range,
def integrate(self, wavelengths, densities): # define short names for the involved wavelength grids wa = wavelengths wb = self._Wavelengths # create a combined wavelength grid, restricted to the overlapping interval w1 = wa[(wa >= wb[0]) & (wa <= wb[-1])] w2 = wb[(wb >= wa[0]) & (wb <= wa[-1])] w = np.unique(np.hstack((w1, w2))) if len(w) < 2: return 0 # log-log interpolate SED and transmission on the combined wavelength grid # (use scipy interpolation function for SED because np.interp does not support broadcasting) F = np.exp(interp1d(np.log(wa), _log(densities), copy=False, bounds_error=False, fill_value=0.)(np.log(w))) T = np.exp(np.interp(np.log(w), np.log(wb), _log(self._Transmission), left=0., right=0.)) # perform the integration if self._PhotonCounter: return np.trapz(x=w, y=w * F * T) else: return np.trapz(x=w, y=F * T) ## This private helper function returns the natural logarithm for positive values, and a large negative number # (but not infinity) for zero or negative values.
def obj_func_disc_nofit(xx, e_g, e_p, g, rho, Q_g, Q_p, head_g, head_p, optimizing = True): H_T = int(price_duration.Frequency.max()) # total duration (100%) prc_g, prc_p, freq_g, freq_p = [],[],[],[] for i,x in enumerate(price_duration.Frequency): if x < xx: # Power Generation price and duration prc_g.append(price_duration.Price[i]), freq_g.append(x) if H_T - xx < x < H_T: # Pumping price and duration prc_p.append(price_duration.Price[i]), freq_p.append(x) prc_g = np.array(prc_g) # generation price prc_p = np.array(prc_p) # pumping price freq_g = np.array(freq_g) # generation duration freq_p = np.array(freq_p) # pumping duration # Use numerical integration to integrate (Trapezoidal rule) Power_Revenue = np.trapz(prc_g, freq_g, dx=0.1, axis = -1)*e_g*rho*g*Q_g*head_g/(10**6) Pumping_Cost = np.trapz(prc_p, freq_p, dx=0.1, axis = -1)/e_p*rho*g*Q_p*head_p/(10**6) z = Power_Revenue - Pumping_Cost # profit return z if optimizing else -z # fit a curve
def pan_tompkins(self,data_buf): ''' 1) 3rd differential of the data_buffer (512 samples) 2) square of that 3) integrate -> integrate(data_buffer) 4) take the pulse from that ''' order = 3 diff = np.diff(data_buf,3) for i in range(order): diff = np.insert(diff,0,0) square = np.square(diff) window = 64 integral = np.zeros((len(square))) for i in range(len(square)): integral[i] = np.trapz(square[i:i+window]) # print(integral) return integral
def nena(locs, info, callback=None): bin_centers, dnfl_ = next_frame_neighbor_distance_histogram(locs, callback) def func(d, a, s, ac, dc, sc): f = a * (d / s**2) * _np.exp(-0.5 * d**2 / s**2) fc = ac * (d / sc**2) * _np.exp(-0.5 * (d**2 + dc**2) / sc**2) * _iv(0, d * dc / sc) return f + fc pdf_model = _lmfit.Model(func) params = _lmfit.Parameters() area = _np.trapz(dnfl_, bin_centers) median_lp = _np.mean([_np.median(locs.lpx), _np.median(locs.lpy)]) params.add('a', value=area/2, min=0) params.add('s', value=median_lp, min=0) params.add('ac', value=area/2, min=0) params.add('dc', value=2*median_lp, min=0) params.add('sc', value=median_lp, min=0) result = pdf_model.fit(dnfl_, params, d=bin_centers) return result, result.best_values['s']
def test_gamma(): tsample = 0.005 / 1000 with pytest.raises(ValueError): t, g = utils.gamma(0, 0.1, tsample) with pytest.raises(ValueError): t, g = utils.gamma(2, -0.1, tsample) with pytest.raises(ValueError): t, g = utils.gamma(2, 0.1, -tsample) for tau in [0.001, 0.01, 0.1]: for n in [1, 2, 5]: t, g = utils.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))
def mag2diam(ref_mag, ref_band, teff, nbr_pts=10): """ Returns the DIAMETER (not radius) of a black body of a given magnitude (U,B,V,R,I,J,H,K) """ filt = ALL_FILTERS[ref_band] if nbr_pts != 10: lam = _gen_lam_for_filter(filt, nbr_pts) else: lam = filt['span_10pts'] # flux in watts from magnitude res = filt['flux_wm2m']*10**(-0.39809*ref_mag)*filt['delta_wl'] # integrated flux res /= np.trapz(blackbody_spectral_irr(teff, lam), lam) return 2*RAD2MAS*np.sqrt(res)
def checkFacts_pdf_Mu(nu=0, tau=0, mu_grid=None, **kwargs): cPrior = c_Prior(nu=nu, tau=tau) if mu_grid is None: mu_grid = make_mu_grid(**kwargs) phi_grid = mu2phi(mu_grid) logpdf_grid = tau * phi_grid - nu * c_Phi(phi_grid) - cPrior logJacobian_grid = logJacobian_mu(mu_grid) pdf_grid = np.exp(logpdf_grid + logJacobian_grid) Integral = np.trapz(pdf_grid, mu_grid) E_mu_numeric = np.trapz(pdf_grid * mu_grid, mu_grid) E_mu_formula = tau/nu E_phi_numeric = np.trapz(pdf_grid * phi_grid, mu_grid) E_phi_formula = E_phi(nu,tau) print "nu=%7.3f tau=%7.3f" % (nu, tau) print " Integral=% 7.3f should be % 7.3f" % (Integral, 1.0) print " E[mu]=% 7.3f should be % 7.3f" % (E_mu_numeric, E_mu_formula) print " E[phi]=% 7.3f should be % 7.3f" % (E_phi_numeric, E_phi_formula)
def checkFacts_pdf_Phi( nu=0, tau=0, phi_grid=None, **kwargs): pdf_grid, phi_grid = make_pdfgrid_phi(nu, tau, phi_grid, **kwargs) mu_grid = phi2mu(phi_grid, **kwargs) IntegralVal = np.trapz(pdf_grid, phi_grid) E_phi_numeric = np.trapz(pdf_grid * phi_grid, phi_grid) E_phi_formula = E_phi(nu=nu, tau=tau, **kwargs) E_mu_numeric = np.trapz(pdf_grid * mu_grid, phi_grid) E_mu_formula = tau/nu mode_phi_numeric = phi_grid[np.argmax(pdf_grid)] mode_phi_formula = mu2phi(tau/nu, **kwargs) print "nu=%7.3f tau=%7.3f" % (nu, tau) print " Integral=% 7.3f should be % 7.3f" % (IntegralVal, 1.0) print " E[mu]=% 7.3f should be % 7.3f" % (E_mu_numeric, E_mu_formula) print " E[phi]=% 7.3f should be % 7.3f" % (E_phi_numeric, E_phi_formula) print " mode[phi]=% 7.3f should be % 7.3f" % ( mode_phi_numeric, mode_phi_formula)
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
def make_curve(x, y, xlabel, ylabel, filename): """Draw and save ROC or PR curves.""" auc = np.trapz(x, y, dx=0.001) # area under the curve sns.plt.figure() sns.plt.clf() sns.plt.plot(x, y) sns.plt.xlim([0, 1]) sns.plt.ylim([0, 1]) sns.plt.xlabel(xlabel) sns.plt.ylabel(ylabel) sns.plt.title("AUC: {}".format(auc)) sns.plt.savefig(filename) return auc
def plotROC(self): for classifier in self.structures: x, y = [], [] for threshold in np.arange(self.minScore, self.maxScore, 0.1): tp, tn, fp, fn = 0, 0, 0, 0 for realS in self.structures: for score in self.scores[(classifier, realS)]: if score >= threshold: # Positive if classifier == realS: # True tp += 1 else: # False fp += 1 else: # Negative if classifier != realS: # True tn += 1 else: # False fn += 1 tpr = tp / (tp + fn) fpr = fp / (tn + fp) x.append(fpr) y.append(tpr) print("----- Receiver Operating Characteristic Curve -----") print("Classifier:", classifier) x.sort() y.sort() auc = np.trapz(y, x) # Area Under Curve print("AUC:", auc) pyplot.title("Classifier " + classifier) pyplot.xlabel('False Positive Rate') pyplot.ylabel('True Positive Rate') pyplot.plot([0, 1], [0, 1]) pyplot.plot(x, y) pyplot.show()
def test_limbdark(tol = 1e-4): ''' Test the limb darkening normalization by comparing the total flux of the star to what you get with the Stefan-Boltzmann law. ''' # Instantiate the star r = 0.1 teff = 3200 star = Star('A', m = 0.1, r = r, nz = 31, color = 'k', limbdark = [u1], teff = teff) # True luminosity truth = 4 * np.pi * (r * RSUN) ** 2 * SBOLTZ * teff ** 4 # System system = System(star, quiet = True) system.distance = 12.2 # Compute the light curve time = np.arange(0., 1., 10) system.compute(time, lambda1 = 0.01, lambda2 = 1000, R = 3000) # Luminosity bol = np.trapz(system.flux[0], system.A.wavelength * 1e6) lum = 4 * np.pi * bol * (12.2 * PARSEC) ** 2 # Check! assert np.abs(lum - truth) / truth < tol, \ "Incorrect bolometric flux: %.10e != %.10e." % (lum, truth)
def autocorrelation_function(self, r): """compute the real space autocorrelation function for the Gaussian random field model """ # compute the cut-level parameter beta beta = np.sqrt(2) * erfinv(2 * (1-self.frac_volume) - 1) # the covariance of the GRF acf_psi = (np.exp(-r/self.corr_length) * (1 + r / self.corr_length) * np.sinc(2*r/self.repeat_distance)) # integral discretization. henning says: the resolution 1e-2 is ad hoc, test required, # the integrand has a (integrable) singularity for t=1 and acf_psi = 1, so an adaptive # discretization seems preferable -> TODO dt = 1e-2 t = np.arange(0, 1, dt) # the gridded integrand, via change of integration variable # compared to the wp-2 docu, to enable array-based computation t_gridded, acf_psi_gridded = np.meshgrid(t, acf_psi) integrand_gridded = (acf_psi_gridded / np.sqrt(1 - (t_gridded * acf_psi_gridded)**2) * np.exp( - beta**2 / (1 + t_gridded * acf_psi_gridded))) acf = 1.0 / (2 * np.pi) * np.trapz(integrand_gridded, x=t_gridded) return acf # ft not known analytically: deligate # def ft_autocorrelation_function(self, k):
def trapz2(f, x=None, y=None, dx=1.0, dy=1.0): """Double integrate.""" return numpy.trapz(numpy.trapz(f, x=y, dx=dy), x=x, dx=dx)
def evaluate_recall(self, candidate_boxes, ar_thresh=0.5): 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.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.05 thresholds = np.minimum(np.arange(0.2, 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
def compute_ess_transmission(beam_sigma, slits, dispersion): """Compute the transmission as a function of the momentum offset (in %) from a simple analytical model.""" n_steps = 10000 dx = 3.0/n_steps sigma = beam_sigma/2.8 slits_at = slits/dispersion error = np.arange(-1.5, 1.5, dx) slits = np.zeros(n_steps) slits[np.where((error < slits_at) & (error > -slits_at))] = 1.0 beam = np.exp(-(error/sigma)**2/2) return np.roll(np.convolve(slits, beam, mode="same"), -1)/np.trapz(beam)
def bpm_fit(a, p): if a is None: return None x = a['channel']['data'][p][:, 1] y = a['channel']['data'][p][:, 2] ar = np.trapz(y / np.sum(y) * len(y), x) mean = np.mean(x * y / np.sum(y) * len(y)) rms = np.std(x * y / np.sum(y) * len(y)) popt, pcov = curve_fit(gaussian, x, y, p0=[ar, mean, rms]) return [popt, np.sqrt(pcov.diagonal())]
def avg_mole_per_sec(raw, start_point=0, end_point='eq', constV=False): #print type(end_point) T = raw['temperature'] x = raw['axis0'] rr = raw['net_reaction_rate'] V = raw['volume'] n_rxn = rr.shape[1] - 1 i_start = start_point if type(end_point) is int: i_end = end_point if i_end == i_start: mole_per_sec = rr[i_start, :] * float(V[i_start]) return np.transpose(mole_per_sec) else: if end_point is 'ign': T_end = T[0] + 400 elif end_point is 'eq': T_end = T[-1] - 50 i_end = np.argmin(abs(T - T_end)) avg_rr = np.ndarray([n_rxn+1,1]) for id_rxn in range(1, n_rxn + 1): if constV: mole_per_sec = rr[i_start:i_end,id_rxn] * V[0] else: mole_per_sec = np.multiply(rr[i_start:i_end, id_rxn], V[i_start:i_end]) int_rr = np.trapz(np.transpose(mole_per_sec), np.transpose(x[i_start:i_end])) avg_rr[id_rxn] = float(int_rr) / (x[i_end] - x[i_start]) return avg_rr
def WRelease(self,kr,r0,kt,t0,kb,b0,T): """ Do the analytical "release" of guest to standard concentration """ ### SETUP # Example: print WRelease(5.0,24.00,100.0,180.0,100.0,180.0,298.15) # kr, r0: distance restraint (r) force constant, target value # kt, t0: angle restraint (theta) force constant, target value # kb, b0: angle restraint (b) force constant, target value # T: Temperature R = 1.987204118e-3 # kcal/mol-K, a.k.a. boltzman constant beta = 1/(T*R) rlb,rub,rst = [0.0,100.0,0.0001] # r lower bound, upper bound, step size tlb,tub,tst = [0.0,np.pi,0.00005] # theta ", ", " blb,bub,bst = [0.0,np.pi,0.00005] # b ", ", " def fr(val): return (val**2)*np.exp(-beta*kr*(val-r0)**2) def ft(val): return np.sin(val)*np.exp(-beta*kt*(val-np.radians(t0))**2) def fb(val): return np.sin(val)*np.exp(-beta*kb*(val-np.radians(b0))**2) ### Integrate rint,tint,bint = [0.0,0.0,0.0] intrange = np.arange(rlb,rub,rst) rint = np.trapz(fr(intrange),intrange) intrange = np.arange(tlb,tub,tst) tint = np.trapz(ft(intrange),intrange) intrange = np.arange(blb,bub,bst) bint = np.trapz(fb(intrange),intrange) return R*T*np.log(np.pi*(1.0/1660.0)*rint*tint*bint)
def quad(fintegrand, xmin, xmax, n=1e4): spacings = np.logspace(np.log10(xmin), np.log10(xmax), n) integrand_arr = fintegrand(spacings) val = np.trapz(integrand_arr, dx=np.diff(spacings)) return val
def a2t(at, Om0=0.27, Oml0=0.73, h=0.700): integrand = lambda x: 1./(x*sqrt(Oml0+Om0*x**-3.0)) # current_time,err = si.quad(integrand,0.0,at,epsabs=1e-6,epsrel=1e-6) current_time = quad(integrand, 1e-4, at) # spacings = np.logspace(-5,np.log10(at),1e5) # integrand_arr = integrand(spacings) # current_time = np.trapz(integrand_arr,dx=np.diff(spacings)) current_time *= 9.779/h return current_time
def integrate_inf(fcn, error=1e-3, initial_guess=10): """ Integrate a function *fcn* from zero to infinity, stopping when the answer changes by less than *error*. Hopefully someday we can do something better than this! """ xvals = np.logspace(0,np.log10(initial_guess), initial_guess+1)-.9 yvals = fcn(xvals) xdiffs = xvals[1:] - xvals[:-1] # Trapezoid rule, but with different dxes between values, so np.trapz # will not work. areas = (yvals[1:] + yvals[:-1]) * xdiffs / 2.0 area0 = np.sum(areas) # Next guess. next_guess = 10 * initial_guess xvals = np.logspace(0,np.log10(next_guess), 2*initial_guess**2+1)-.99 yvals = fcn(xvals) xdiffs = xvals[1:] - xvals[:-1] # Trapezoid rule. areas = (yvals[1:] + yvals[:-1]) * xdiffs / 2.0 area1 = np.sum(areas) # Now we refine until the error is smaller than *error*. diff = area1 - area0 area_last = area1 one_pow = 3 while diff > error: next_guess *= 10 xvals = np.logspace(0,np.log10(next_guess), one_pow*initial_guess**one_pow+1) - (1 - 0.1**one_pow) yvals = fcn(xvals) xdiffs = xvals[1:] - xvals[:-1] # Trapezoid rule. areas = (yvals[1:] + yvals[:-1]) * xdiffs / 2.0 area_next = np.sum(areas) diff = area_next - area_last area_last = area_next one_pow+=1 return area_last
def trapzint(f, a, b, bins=10000): zbins = np.logspace(np.log10(a + 1), np.log10(b + 1), bins) - 1 return np.trapz(f(zbins[:-1]), x=zbins[:-1], dx=np.diff(zbins))
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(0, nStep): errs[p] = np.shape(np.where(distances > thres[p]))[-1] score = np.trapz(errs) return score
def obj_func_cont(xx, e_g, e_p, g, rho, Q_g, Q_p, head_g, head_p, optimizing = True): H_T = int(price_duration.Frequency.max()) # total duration (100%) x1 = np.arange(0,xx) y1 = f(x1) x2 = np.arange(H_T-xx,H_T) y2 = f(x2) Power_Revenue = np.trapz(y1, x1, dx=0.1, axis = -1)*e_g*rho*g*Q_g*head_g/(10**6) Pumping_Cost = np.trapz(y2, x2, dx=0.1, axis = -1)/e_p*rho*g*Q_p*head_p/(10**6) z = Power_Revenue - Pumping_Cost # profit return -z if optimizing else z # objective function to maximize - discrete
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
def _get_indice(cls, w, flux, blue, red, band=None, unit='ew', degree=1, **kwargs): """ compute spectral index after continuum subtraction Parameters ---------- w: ndarray (nw, ) array of wavelengths in AA flux: ndarray (N, nw) array of flux values for different spectra in the series blue: tuple(2) selection for blue continuum estimate red: tuple(2) selection for red continuum estimate band: tuple(2), optional select region in this band only. default is band = (min(blue), max(red)) unit: str `ew` or `mag` wether equivalent width or magnitude degree: int (default 1) degree of the polynomial fit to the continuum Returns ------- ew: ndarray (N,) equivalent width array """ wi, fi = cls.continuum_normalized_region_around_line(w, flux, blue, red, band=band, degree=degree) if unit in (0, 'ew', 'EW'): return np.trapz(1. - fi, wi, axis=-1) else: m = np.trapz(fi, wi, axis=-1) m = -2.5 * np.log10(m / np.ptp(wi)) return m
def __init__(self, wavelength, transmit, name='', dtype="photon", unit=None): """Constructor""" self.name = name try: # get units from the inputs self._wavelength = wavelength.magnitude self.wavelength_unit = str(wavelength.units) except AttributeError: self._wavelength = wavelength self.wavelength_unit = unit self.transmit = transmit self.norm = trapz(transmit, self._wavelength) self._lT = trapz(self._wavelength * transmit, self._wavelength) self._lpivot = np.sqrt( self._lT / trapz(transmit / self._wavelength, self._wavelength) ) self._cl = self._lT / self.norm self.set_dtype(dtype)
def leff(self): """ Unitwise Effective wavelength leff = int (lamb * T * Vega dlamb) / int(T * Vega dlamb) """ with Vega() as v: s = self.reinterp(v.wavelength) w = s._wavelength leff = np.trapz(w * s.transmit * v.flux.magnitude, w, axis=-1) leff /= np.trapz(s.transmit * v.flux.magnitude, w, axis=-1) if self.wavelength_unit is not None: return leff * unit[self.wavelength_unit] else: return leff
def energy_radiated_approximation_and_farfield2(self,trajectory, gamma, x, y, distance): # N = trajectory.shape[1] N = trajectory.nb_points() if distance == None: # in radian : n_chap = np.array([x, y, 1.0 - 0.5 * (x ** 2 + y ** 2)]) X = np.sqrt(x ** 2 + y ** 2 )#TODO a changer #in meters : else : X = np.sqrt(x ** 2 + y ** 2 + distance ** 2) n_chap = np.array([x, y, distance]) / X E = np.zeros((3,), dtype=np.complex) integrand = np.zeros((3, N), dtype=np.complex) A1 = (n_chap[0] * trajectory.a_x + n_chap[1] * trajectory.a_y + n_chap[2] * trajectory.a_z) A2 = (n_chap[0] * (n_chap[0] - trajectory.v_x) + n_chap[1] * (n_chap[1] - trajectory.v_y) + n_chap[2] * (n_chap[2] - trajectory.v_z)) Alpha2 = np.exp( 0. + 1j * self.photon_frequency * (trajectory.t + X / codata.c - n_chap[0] * trajectory.x - n_chap[1] * trajectory.y - n_chap[2] * trajectory.z)) Alpha1 = (1.0 / (1.0 - n_chap[0] * trajectory.v_x - n_chap[1] * trajectory.v_y - n_chap[2] * trajectory.v_z)) ** 2 integrand[0] += (A1 * (n_chap[0] - trajectory.v_x) - A2 * trajectory.a_x) * Alpha2 * Alpha1 integrand[1] += (A1 * (n_chap[1] - trajectory.v_y) - A2 * trajectory.a_y) * Alpha2 * Alpha1 integrand[2] += (A1 * (n_chap[2] - trajectory.v_z) - A2 * trajectory.a_z) * Alpha2 * Alpha1 for k in range(3): # E[k] = np.trapz(integrand[k], self.trajectory.t) E[k] = np.trapz(integrand[k], trajectory.t) return E
def energy_radiated_farfield2(self, trajectory, gamma, x, y, distance): N = trajectory.nb_points() if distance == None: # in radian : n_chap = np.array([x, y, 1.0 - 0.5 * (x ** 2 + y ** 2)]) X = np.sqrt(x ** 2 + y ** 2 )#TODO a changer #in meters : else : X = np.sqrt(x ** 2 + y ** 2 + distance ** 2) n_chap = np.array([x, y, distance]) / X R= X/codata.c - n_chap[0] * trajectory.x- n_chap[1] * trajectory.y - n_chap[2] * trajectory.z E = np.zeros((3,), dtype=np.complex) integrand = np.zeros((3, N), dtype=np.complex) A1 = (n_chap[0] * trajectory.a_x + n_chap[1] * trajectory.a_y + n_chap[2] * trajectory.a_z) A2 = (n_chap[0] * (n_chap[0] - trajectory.v_x) + n_chap[1] * (n_chap[1] - trajectory.v_y) + n_chap[2] * (n_chap[2] - trajectory.v_z)) Alpha2 = np.exp(0. + 1j * self.photon_frequency * (trajectory.t + R)) Alpha1 = (1.0 / (1.0 - n_chap[0] * trajectory.v_x - n_chap[1] * trajectory.v_y - n_chap[2] * trajectory.v_z)) ** 2 cst=codata.c/(R*gamma**2) integrand[0] += (A1 * (n_chap[0] - trajectory.v_x) - A2 * trajectory.a_x + cst * (n_chap[0] - trajectory.v_x) ) * Alpha2 * Alpha1 integrand[1] += (A1 * (n_chap[1] - trajectory.v_y) - A2 * trajectory.a_y + cst * (n_chap[1] - trajectory.v_y) ) * Alpha2 * Alpha1 integrand[2] += (A1 * (n_chap[2] - trajectory.v_z) - A2 * trajectory.a_z + cst * (n_chap[2] - trajectory.v_z) ) * Alpha2 * Alpha1 for k in range(3): # E[k] = np.trapz(integrand[k], self.trajectory.t) E[k] = np.trapz(integrand[k], trajectory.t) E *= -1j * np.exp(1j * self.photon_frequency / codata.c * (n_chap[0] * x + n_chap[1] * y + n_chap[2] * distance)) return E
def energy_radiated_approx2(self,trajectory, gamma, x, y, distance): N = trajectory.nb_points() n_chap = np.array([x - trajectory.x * codata.c, y - trajectory.y * codata.c, distance - trajectory.z * codata.c]) # R = np.zeros(n_chap.shape[1]) # for i in range(n_chap.shape[1]): # R[i] = np.linalg.norm(n_chap[:, i]) # n_chap[:, i] /= R[i] R = np.sqrt( n_chap[0]**2 + n_chap[1]**2 + n_chap[2]**2 ) n_chap[0,:] /= R n_chap[1,:] /= R n_chap[2,:] /= R E = np.zeros((3,), dtype=np.complex) integrand = np.zeros((3, N), dtype=np.complex) A1 = (n_chap[0] * trajectory.a_x + n_chap[1] * trajectory.a_y + n_chap[2] * trajectory.a_z) A2 = (n_chap[0] * (n_chap[0] - trajectory.v_x) + n_chap[1] * (n_chap[1] - trajectory.v_y) + n_chap[2] * (n_chap[2] - trajectory.v_z)) Alpha2 = np.exp(0. + 1j * self.photon_frequency * (trajectory.t + R / codata.c)) Alpha1 = (1.0 / (1.0 - n_chap[0] * trajectory.v_x - n_chap[1] * trajectory.v_y - n_chap[2] * trajectory.v_z)) ** 2 integrand[0] += ((A1 * (n_chap[0] - trajectory.v_x) - A2 * trajectory.a_x) ) * Alpha2 * Alpha1 integrand[1] += ((A1 * (n_chap[1] - trajectory.v_y) - A2 * trajectory.a_y) ) * Alpha2 * Alpha1 integrand[2] += ((A1 * (n_chap[2] - trajectory.v_z) - A2 * trajectory.a_z) ) * Alpha2 * Alpha1 for k in range(3): E[k] = np.trapz(integrand[k], trajectory.t) #E[k] = integrate.simps(integrand[k], trajectory.t) return E
def integration(self,is_quadrant=0,use_flux_per_mrad2_or_mm2=1): if (self.X is None or self.Y is None): raise Exception(" X and Y must be array for integration") if self.X.shape != self.Y.shape: raise Exception(" X and Y must have the same shape") if len(self.X.shape)==2 : X = np.linspace(self.X.min(), self.X.max(), self.X.shape[0]) Y = np.linspace(self.Y.min(), self.Y.max(), self.Y.shape[1]) res1=integrate.trapz(self.intensity, X) res = integrate.trapz(res1, Y) # res = integrate.simps(integrate.simps(self.intensity, X), Y) # res = self.intensity.sum() * (self.X[1,0] - self.X[0,0]) * (self.Y[0,1] - self.X[0,0]) # regular grid only else : # X and Y are 1d array if len(self.X) == 1: res = self.intensity[0] else: # choix arbitraire XY = np.zeros_like(self.X) for i in range(1, len(self.X)): XY[i] = XY[i-1]+np.sqrt((self.X[i] - self.X[i-1]) * 2 + (self.Y[i] - self.Y[i-1]) ** 2) res = np.trapz(self.intensity, XY) # Note that the value of flux is in phot/s/0.1%bw/mrad2 (or .../mm2) and # our grid is in rad (or m), therefore we must account for this: if use_flux_per_mrad2_or_mm2: res *= 1e6 # in case the calculation is for a quadrant, the integral is four times the calculated value if is_quadrant: res *= 4 return res
def kde(self): """ Calculate kernel density estimator distribution function """ plot = True # Input data x_data = np.linspace(-math.pi, math.pi, 1000) # Kernels, centered at input data points kernels = [] for datapoint in self.data: # Make the basis function as a von mises PDF kernel = self.vonMisesPDF(x_data, mu=datapoint) kernels.append(kernel) # Handle weights if len(self.weights) > 0: kernels = np.asarray(kernels) weighted_kernels = np.multiply(kernels, self.weights[:, None]) else: weighted_kernels = kernels # Normalize pdf vmkde = np.sum(weighted_kernels, axis=0) vmkde = vmkde / np.trapz(vmkde, x=x_data) self.fn = interp1d(x_data, vmkde)
