def run(data, sweep_rate=20.,
        c_range=(0.4, 0.6), co_range=(0.6, 0.9),
        exe='', graph=True, baseline=False, copy=False):
    # INIT data
    params = {}
    area_CO = None
    area_H = None

    cycle_CO = data.get_scan(1)
    cycle_baseline = data.get_scan(2)

    if data.sweep_rate:
        sr = data.sweep_rate
    else:
        sr = sweep_rate

    # RUN stuff
    if "CO" in exe:
        params["CO"] = CO(cycle_CO, cycle_baseline, c_range, co_range, add_baseline=baseline, copy=copy)
        x, y = params["CO"][0]
        Q_CO = np.trapz(y, x)  # V C / s
        factor_CO = 420e-6 * sr  # C V / s cm2
        area_CO = Q_CO / factor_CO  # cm2
    if "H" in exe:
        params["H"] = H(cycle_CO, cycle_baseline, co_range[0], copy)
        x, y = params["H"]
        Q_H = np.trapz(y, x)  # V C / s
        factor_H = 210e-6 * sr * 1.e-3  # C V / s cm2
        area_H = Q_H / factor_H  # cm2
    if graph:
        plot(cycle_CO, cycle_baseline, paramsCO=params.get("CO"),
             paramsH=params.get("H"), exe=exe, graph=graph)

    return area_CO, area_H

 def __init__(self, fct, makenorm=True, makeunlog=False, fct_log=False,
              npts=1000, *args, **kwargs):
     if not callable(fct):
         raise Exception("Not callable")
     self.fct = fct
     self.args = args
     self.kwargs = kwargs
     # apply exp on a ln-function at callback or not
     self.makeunlog = makeunlog is True
     # renorm at callback or not
     self.makenorm = makenorm is True
     # whether to add or divide the normalization
     self.fct_log = fct_log is not False and makeunlog is False
     self.prior_log = kwargs.get('prior_log', False)  # jeffrey
     if kwargs.get('prior_bounds', None) is None:
         raise Exception("No prior bounds found")
     if self.prior_log:
         self.x = np.logspace(np.log(kwargs.get('prior_bounds')[0]),
                              np.log(kwargs.get('prior_bounds')[1]),
                              npts,
                              base=np.e)
     else:
         self.x = np.linspace(kwargs.get('prior_bounds')[0],
                              kwargs.get('prior_bounds')[1],
                              npts)
     if makenorm is True:
         resfct = self.fct(self.x, *self.args, **self.kwargs)
         if self.fct_log:
             self.normval = np.log(np.trapz(np.exp(resfct), self.x))
         elif self.makeunlog:
             self.normval = np.trapz(np.exp(resfct), self.x)
         else:
             self.normval = np.trapz(resfct, self.x)

    def beam_phase_sharpWindow(self):
        '''
        *Beam phase measured at the main RF frequency and phase. The beam is 
        averaged over a window. The coefficients of sine and cosine components
        determine the 
        beam phase, projected to the range -Pi/2 to 3/2 Pi. Note that this beam
        phase is already w.r.t. the instantaneous RF phase.*
        '''    
        
        # Main RF frequency at the present turn
        omega_rf = self.rf_station.omega_rf[0,self.rf_station.counter[0]]
        phi_rf = self.rf_station.phi_rf[0,self.rf_station.counter[0]]
        
        if self.alpha != 0.0:
            left_boundary = self.profile.bin_centers \
                >= self.time_offset - np.pi/omega_rf
            right_boundary = self.profile.bin_centers \
                <= -1/self.alpha + self.time_offset - 2*np.pi/omega_rf
            indexes = left_boundary * right_boundary
        else:
            indexes = np.ones(self.profile.n_slices, dtype=bool)
        
        # Convolve with window function
        scoeff = np.trapz( np.sin(omega_rf*self.profile.bin_centers[indexes]\
                                   + phi_rf) \
                           *self.profile.n_macroparticles[indexes],
                           dx=self.profile.bin_size )
        ccoeff = np.trapz( np.cos(omega_rf*self.profile.bin_centers[indexes]\
                                   + phi_rf) \
                           *self.profile.n_macroparticles[indexes],
                           dx=self.profile.bin_size )

        # Project beam phase to (pi/2,3pi/2) range
        self.phi_beam = np.arctan(scoeff/ccoeff) + np.pi

def avgField(time, field, tstart=None, std=False):
    """
    This subroutine computes the time-average (and the std) of a time series

    >>> ts = MagicTs(field='misc', iplot=False, all=True)
    >>> nuavg = avgField(ts.time, ts.topnuss, 0.35)
    >>> print(nuavg)

    :param time: time
    :type time: numpy.ndarray
    :param field: the time series of a given field
    :type field: numpy.ndarray
    :param tstart: the starting time of the averaging
    :type tstart: float
    :param std: when set to True, the standard deviation is also calculated
    :type std: bool
    :returns: the time-averaged quantity
    :rtype: float
    """
    if tstart is not None:
        mask = np.where(abs(time - tstart) == min(abs(time - tstart)), 1, 0)
        ind = np.nonzero(mask)[0][0]
    else:  # the whole input array is taken!
        ind = 0
    fac = 1.0 / (time[-1] - time[ind])
    avgField = fac * np.trapz(field[ind:], time[ind:])

    if std:
        stdField = np.sqrt(fac * np.trapz((field[ind:] - avgField) ** 2, time[ind:]))
        return avgField, stdField
    else:
        return avgField

def test_spectra_get_flux_contributions():
    timestepmin = 40
    timestepmax = 80
    dfspectrum = at.spectra.get_spectrum(
        specfilename, timestepmin=timestepmin, timestepmax=timestepmax, fnufilterfunc=None)

    integrated_flux_specout = np.trapz(dfspectrum['f_lambda'], x=dfspectrum['lambda_angstroms'])

    specdata = pd.read_csv(specfilename, delim_whitespace=True)
    arraynu = specdata.loc[:, '0'].values
    arraylambda_angstroms = const.c.to('angstrom/s').value / arraynu

    contribution_list, array_flambda_emission_total = at.spectra.get_flux_contributions(
        modelpath, timestepmin=timestepmin, timestepmax=timestepmax)

    integrated_flux_emission = -np.trapz(array_flambda_emission_total, x=arraylambda_angstroms)

    # total spectrum should be equal to the sum of all emission processes
    print(f'Integrated flux from spec.out:     {integrated_flux_specout}')
    print(f'Integrated flux from emission sum: {integrated_flux_emission}')
    assert math.isclose(integrated_flux_specout, integrated_flux_emission, rel_tol=4e-3)

    # check each bin is not out by a large fraction
    diff = [abs(x - y) for x, y in zip(array_flambda_emission_total, dfspectrum['f_lambda'].values)]
    print(f'Max f_lambda difference {max(diff) / integrated_flux_specout}')
    assert max(diff) / integrated_flux_specout < 2e-3

def test_AR_LD():
    """

    Test the Levinson Durbin estimate of the AR coefficients against the
    expercted PSD

    """
    arsig,_,_ = utils.ar_generator(N=512)
    avg_pwr = (arsig*arsig.conjugate()).real.mean()
    order = 8
    ak, sigma_v = tsa.AR_est_LD(arsig, order)
    w, psd = tsa.AR_psd(ak, sigma_v)

    # the psd is a one-sided power spectral density, which has been
    # multiplied by 2 to preserve the property that
    # 1/2pi int_{-pi}^{pi} Sxx(w) dw = Rxx(0)

    # evaluate this integral numerically from 0 to pi
    dw = np.pi/len(psd)
    avg_pwr_est = np.trapz(psd, dx=dw) / (2*np.pi)
    npt.assert_almost_equal(avg_pwr, avg_pwr_est, decimal=0)

    # Test for providing the autocovariance as an input:
    ak,sigma_v = tsa.AR_est_LD(arsig, order, utils.autocov(arsig))
    w, psd = tsa.AR_psd(ak, sigma_v)
    avg_pwr_est = np.trapz(psd, dx=dw) / (2*np.pi)
    npt.assert_almost_equal(avg_pwr, avg_pwr_est, decimal=0)

def get_band_phot( spectrum , bandpass, z, a=1, zero_point=48.6 ):
    import numpy as np

    c = (3.0e18)      #  Angs / s

  #  These conditionals import text files if necessary

    if isinstance( spectrum , type('string') ): spectrum = np.loadtxt(spectrum)
    if isinstance( bandpass , type('string') ): bandpass = np.loadtxt(bandpass)

    interp_flux = np.interp(bandpass[:,0], spectrum[:,0]*(1.0+z),
        spectrum[:,a] )

    weighted_flux = []
    nu_arr = []
    weights = []

    nu_arr = c / bandpass[:,0]
    weighted_flux = interp_flux * bandpass[:,1] * bandpass[:,0]/nu_arr/nu_arr
    weights = bandpass[:,1] / nu_arr

    numer = np.trapz( weighted_flux , x=nu_arr )
    denom = np.trapz( weights , x=nu_arr)

    #return -2.5*np.log10( numer / denom ) + zero_point   # <--- "MAGNITUDE"
    return -2.5*np.log10(numer/denom/zero_point)# - zero_point   # <--- "MAGNITUDE"

    def solve_nonlinear(self, params, unknowns, resids):
        '''a VERY coarse approximation that takes the speed profile given in the original proposal as given. It just serves as a place holder for now. It's better than nothing, but a real analysis is needed here'''
        t1 = (300 - 0) * 1609 / (9.81 * 0.5) / 3600 # time needed to accelerate
        t2 = (555 - 300) * 1609 / (9.81 * 0.5) / 3600 # time needed to accelerate
        t3 = (params['max_velocity'] - 555) * 1609 / (9.81 * 0.5) / 3600 # time needed to accelerate
        #speed profile data from hyperloop alpha proposal, pg43
        dataStart = np.array([
            [0, 0],
            [t1, 300 * 1.609], # (300[mi/h] * 1609[m/mi]) / (9.81[m/s] * 0.5) / 3600[s/h]
            [167, 300 * 1.609],
            [167 + t2, 555 * 1.609], # (555 - 300[mi/h] * 1609[m/mi]) / (9.81[m/s] * 0.5) / 3600[s/h]
            [435, 555 * 1.609],
            [435 + t3, params['max_velocity']]])
        startUp = np.trapz(dataStart[:, 1] / 3600, dataStart[:, 0]) * 1000 # km covered during start up Los Angeles Grapevine
        dataEnd = np.array([
            [0, params['max_velocity']],
            [t3, 555 * 1.609],
            [t3 + 100, 555 * 1.609],
            [t3 + 100 + t2, 300 * 1.609],
            [t3 + 100 + t2 + 400, 300 * 1.609],
            [t3 + 100 + t2 + 400 + t1, 0]])
        windDown = np.trapz(dataEnd[:, 1] / 3600, dataEnd[:, 0]) * 1000 # km covered during wind down along I-580 to SF      
        len_middle = params['tube_length'] - (startUp + windDown)
        time_middle = middleLength / params['max_velocity']
        unknowns['time_mission'] = time_middle + 435 + t3 + t3 + 100 + t2 + 400 + t1

        unknowns['energy'] = (params['pwr_req'] * unknowns['time_mission'] / 3600.0) * (1 + params['pwr_marg']) # convert to hours

def energy(powers, default_timestep=8):
    """Compute the energy associated to a list of measures (in W)
    and associated timestamps (in s).
    """
    energy = {"night_rate": 0, "day_rate": 0, "value": 0}
    if len(powers) == 1:
        if powers[0].night_rate == 1:
            energy["night_rate"] = powers[0].value / 1000 * default_timestep / 3600
        else:
            energy["day_rate"] = powers[0].value / 1000 * default_timestep / 3600
        energy["value"] = energy["day_rate"] + energy["night_rate"]
    else:
        x = []
        day_rate = []
        night_rate = []
        for i in powers:
            x.append(i.timestamp)
            if i.night_rate == 1:
                night_rate.append(i.value)
                day_rate.append(0)
            else:
                day_rate.append(i.value)
                night_rate.append(0)
        energy["night_rate"] = numpy.trapz(night_rate, x) / 1000 / 3600
        energy["day_rate"] = numpy.trapz(day_rate, x) / 1000 / 3600
        energy["value"] = energy["day_rate"] + energy["night_rate"]
    return energy

    def set_angle_width(self, width):
        # Calculates the normalizing coefficient for the distribution. From
        # Bringi & Chandrasekar 2001, pg. 68 (which cite Mardia 1972) for an
        # axial distribution.
        self._angle_width = width
        if width <= 0.:
            return
        t = np.linspace(0, 1, 500)
        self._angle_b = 1. / (2. * np.trapz(
            np.exp(-self._angle_width * t**2), t))

        # We neglect here the factor of 1/2pi, as without it:
        # * The graphs on page 69 can be reproduced correctly
        # * The distribution is normalized (sums to 1)
        # Size 90 below *must* match the number returned from T-matrix
        self._angle_vals = np.linspace(np.pi/2, np.pi, TMATRIX_ANGLES)
        self._angle_weights = self._angle_b * (np.exp(
            -self._angle_width * np.cos(self._angle_vals)**2)
            * np.sin(self._angle_vals))

        # Since the distribution and scattering calculations are symmetric, we
        # just calculate the right half of the distribution and double the
        # weights
        self._angle_weights[1:] *= 2
        norm_factor = np.trapz(self._angle_weights, self._angle_vals)
        self._angle_weights /= norm_factor

def integrate_3D(x, y, z, xlin, ylin, zlin, csd):
    # TODO: NOT YET WORKING AS EXPECTED
    X, Y, Z = np.meshgrid(xlin, ylin, zlin)
    Nz = zlin.shape[0]
    Ny = ylin.shape[0]

    m = np.sqrt((x - xlin)**2 + (y - ylin)**2 + (z - zlin)**2)
    m[m < 0.00001] = 0.00001
    csd = csd / m
    #print 'CSD:'
    #print csd

    J = np.zeros((Ny, Nz))
    for i in xrange(Nz):
        J[:, i] = np.trapz(csd[:, :, i], zlin)

    #print '1st integration'
    #print J

    Ny = ylin.shape[0]
    I = np.zeros(Ny)
    for i in xrange(Ny):
        I[i] = np.trapz(J[:, i], ylin)

    #print '2nd integration'
    #print I

    norm = np.trapz(I, xlin)
    
    #print '3rd integration'
    #print norm
    
    return norm

    def effective_wavelength(self, binned=True, wavelengths=None,
                             mode='efflerg'):
        """Calculate :ref:`effective wavelength <synphot-formula-effwave>`.

        Parameters
        ----------
        binned : bool
            Sample data in native wavelengths if `False`.
            Else, sample binned data (default).

        wavelengths : array-like, `~astropy.units.quantity.Quantity`, or `None`
            Wavelength values for sampling.
            If not a Quantity, assumed to be in Angstrom.
            If `None`, ``self.waveset`` or `binset` is used, depending
            on ``binned``.

        mode : {'efflerg', 'efflphot'}
            Flux is first converted to the unit below before calculation:

                * 'efflerg' - FLAM
                * 'efflphot' - PHOTLAM (deprecated)

        Returns
        -------
        eff_lam : `~astropy.units.quantity.Quantity`
            Observation effective wavelength.

        Raises
        ------
        synphot.exceptions.SynphotError
            Invalid mode.

        """
        mode = mode.lower()
        if mode == 'efflerg':
            flux_unit = units.FLAM
        elif mode == 'efflphot':
            warnings.warn(
                'Usage of EFFLPHOT is deprecated.', AstropyDeprecationWarning)
            flux_unit = units.PHOTLAM
        else:
            raise exceptions.SynphotError(
                'mode must be "efflerg" or "efflphot"')

        if binned:
            x = self._validate_binned_wavelengths(wavelengths).value
            y = self.sample_binned(wavelengths=x, flux_unit=flux_unit).value
        else:
            x = self._validate_wavelengths(wavelengths).value
            y = units.convert_flux(x, self(x), flux_unit).value

        num = np.trapz(y * x ** 2, x=x)
        den = np.trapz(y * x, x=x)

        if den == 0.0:  # pragma: no cover
            eff_lam = 0.0
        else:
            eff_lam = abs(num / den)

        return eff_lam * self._internal_wave_unit

    def get_cracking_history():
        XK = []  # position of the first crack
        sig_c_K = [0.]
        eps_c_K = [0.]

        idx_0 = np.argmin(sig_mu_x)
        XK.append(x[idx_0])
        sig_c_0 = sig_mu_x[idx_0] * Ec / Em
        sig_c_K.append(sig_c_0)
        eps_c_K.append(sig_mu_x[idx_0] / Em)

        while True:
            z_x = get_z_x(x, XK)
            sig_c_k, y_i = get_sig_c_K(z_x, sig_c_K[-1])
            if sig_c_k == sig_cu:
                break
            XK.append(y_i)
            sig_c_K.append(sig_c_k)
            eps_c_K.append(
                np.trapz(eps_f(get_z_x(x, XK), sig_c_k), x) / np.amax(x))  # Eq. (10)

            # save the figure
#             plt.figure()
#             plt.plot(x, sig_m(get_z_x(x, XK), sig_c_k))
#             plt.plot(x, sig_mu_x)
#             plt.savefig("D:\\cracking_history\\" + str(len(sig_c_K)) + ".png")
#             plt.close()

        sig_c_K.append(sig_cu)
        eps_c_K.append(np.trapz(eps_f(get_z_x(x, XK), sig_cu), x) / np.amax(x))
        return sig_c_K, eps_c_K

def dDCR_moments(SED1, SED2, bandpass):
    zenith_angle = np.pi/4.0 * galsim.radians
    R500 = galsim.dcr.get_refraction(500, zenith_angle)

    # analytic first moment differences
    R = lambda w:(galsim.dcr.get_refraction(w, zenith_angle) - R500) / galsim.arcsec
    x1 = np.union1d(bandpass.wave_list, SED1.wave_list)
    x1 = x1[(x1 >= bandpass.blue_limit) & (x1 <= bandpass.red_limit)]
    x2 = np.union1d(bandpass.wave_list, SED2.wave_list)
    x2 = x2[(x2 >= bandpass.blue_limit) & (x2 <= bandpass.red_limit)]
    numR1 = np.trapz(R(x1) * bandpass(x1) * SED1(x1), x1)
    numR2 = np.trapz(R(x2) * bandpass(x2) * SED2(x2), x2)
    den1 = SED1.calculateFlux(bandpass)
    den2 = SED2.calculateFlux(bandpass)

    R1 = numR1/den1
    R2 = numR2/den2
    dR_analytic = R1 - R2

    # analytic second moment differences
    V1_kernel = lambda w:(R(w) - R1)**2
    V2_kernel = lambda w:(R(w) - R2)**2
    numV1 = np.trapz(V1_kernel(x1) * bandpass(x1) * SED1(x1), x1)
    numV2 = np.trapz(V2_kernel(x2) * bandpass(x2) * SED2(x2), x2)
    V1 = numV1/den1
    V2 = numV2/den2
    dV_analytic = V1 - V2

    return dR_analytic, dV_analytic, len(x2)

	def solve( self, t, D, VdotO20 ):
		self.Time_Grid = t
		self.Demand = D
		self.VdotO2aerobic = self.solveAerobic( VdotO20 )
		self.VdotO2anaerobic = self.solveAnaerobic()
		self.VO2aerobic = trapz( self.VdotO2aerobic, x=self.Time_Grid )
		self.VO2anaerobic = trapz( self.VdotO2anaerobic, x=self.Time_Grid )

def synchroints(sdip, xtwissdip, disperdip, rho, E0):
    # SOME SYNCHRTRON CONSTANTS
    hbar = hb/qe
    Cq = 55*hbar*cl/32/sqrt(3)/E0            # Chao: 3.8319e-13 m
    Ca = cl*re/3/E0**3                       # Chao: 2.113e-24 m²/(eV)³/s
    Cy = 4*pi/cl*Ca                          # Chao: 8.846e-32 m/(eV)³

    # H FUNCTION: H(s) = beta*D'^2 + 2*alpha*D*D' + gamma*D^2
    Hsx = (xtwissdip[1, 1, :]*disperdip[0, :]**2
           - 2*xtwissdip[0, 1, :]*disperdip[0, :]*disperdip[1, :]
           + xtwissdip[0, 0, :]*disperdip[1, :]**2)

    # SYNCHROTRIN INTEGRALS
    SYNIN1 = trapz(disperdip[0, :], sdip)/rho
    SYNIN2 = 2*pi/rho
    SYNIN3 = 2*pi/rho**2
    SYNIN4x = SYNIN1/rho/rho
    SYNIN4y = 0
    SYNIN5x = trapz(Hsx, sdip)/rho/rho/rho

    # DAMPING PARTITION NUMBERS
    Jx = 1 - SYNIN4x/SYNIN2
    Jy = 1 - SYNIN4y/SYNIN2
    Js = 2 + (SYNIN4x+SYNIN4y)/SYNIN2
    return SYNIN1, SYNIN2, SYNIN3, SYNIN4x, SYNIN4y, SYNIN5x, Jx, Jy, Js, Cq, Ca

def plotRocCurves(file_legend):
	pylab.clf()
	pylab.figure(1)
	pylab.xlabel('1 - Specificity', fontsize=12)
	pylab.ylabel('Sensitivity', fontsize=12)
	pylab.title("Need for Referral")
	pylab.grid(True, which='both')
	pylab.xticks([i/10.0 for i in range(1,11)])
	pylab.yticks([i/10.0 for i in range(0,11)])
	pylab.tick_params(axis="both", labelsize=15)

	for file, legend in file_legend:
		points = open(file,"rb").readlines()
		x = [float(p.split()[0]) for p in points]
		y = [float(p.split()[1]) for p in points]
		dev = [float(p.split()[2]) for p in points]
		x = [0.0] + x
		y = [0.0] + y
		dev = [0.0] + dev
	
		auc = np.trapz(y, x) * 100
		aucDev = np.trapz(dev, x) * 100

		pylab.grid()
		pylab.errorbar(x, y, yerr = dev, fmt='-')
		pylab.plot(x, y, '-', linewidth = 1.5, label = legend + u" (AUC = {0:0.1f}% \xb1 {1:0.1f}%)".format(auc,aucDev))

	pylab.legend(loc = 4, borderaxespad=0.4, prop={'size':12})
	pylab.savefig("referral/referral-curves.pdf", format='pdf')

def qini(y_ordered, percent_targeted = np.linspace(0,1,10), y_rand = None):
    '''
    :param y_ordered: the target variable ordered by uplift model score
    :param percent_targeted:
    :param y_rand: if we want to test a different model, place those ordered outcomes here. If none,
                    plot random outcome
    :return:
    '''

    if y_rand == None:
        y_rand = y_ordered

    nsamp = len(y_ordered)

    lift = []
    random = []
    for pt in percent_targeted:

        n_targ = np.round(nsamp*pt)

        yl = y_ordered[:n_targ].sum()
        yr = y_rand[np.random.randint(0,nsamp,n_targ)].sum()

        lift.append(yl)
        random.append(yr)

    lift = np.array(lift)
    random = np.array(random)

    q = np.trapz(lift, x = percent_targeted) - np.trapz(random, x = percent_targeted)

    return q

  def massLossRate(self, dat, phToSHow = 1):

     #inner boundary
     mdot_a = nm.zeros(dat.nz)
     mdot_a[:] = 2.*nm.pi* dat.Rsc*dat.x[dat.js]* dat.Mx[:, phToSHow,dat.js]*dat.Usc*dat.Dsc                                         
     mdotTot_a =-nm.trapz(mdot_a, dat.z*dat.Rsc)         

     #  outer x boundary
     mout1 = nm.zeros(dat.nz)           
     mout1[:] =  2.*nm.pi* dat.Rsc*dat.x[dat.je]* dat.Mx[:, phToSHow,dat.js]*dat.Usc*self.Dsc                                       
     mout1[:]=[0. if x<0. else x for x in mout1 ]  
     
     mdotTot1 = nm.trapz(mout1, dat.z*dat.Rsc)
     mout2 = nm.zeros(dat.nx)           
     
     #  upper z boundary
     mout2[:] = 2.*nm.pi* dat.Rsc*dat.x[:]*self.Mz[dat.ie, phToSHow,:]*self.Usc*self.Dsc                                    
     mdotTotUp = nm.trapz(mout2, dat.x*dat.Rsc)
     
     mout2[:] = 0.
     #  lower z boundary
     mout2[:] = -2.*nm.pi* dat.Rsc*dat.x[:]*self.Mz[dat.i_s, phToSHow,:]*self.Usc*self.Dsc                           
     mdotTotBot = nm.trapz(mout2, dat.x*dat.Rsc)
          
     return(mdotTot1, mdotTotUp, mdotTotBot, mdotTot_a )

 def test_marginal_likelihood(self):
     """
     Test that the maximum likelihood estimates for the marginal likelihood
     pdfs are correct.
     """
     m = np.linspace(2, 8, 31)
     r = np.linspace(0, 2.0, 21)
     M, R = np.meshgrid(m, r)
     pos = np.empty(M.shape + (2,))
     pos[:, :, 0] = R
     pos[:, :, 1] = M
     mean = np.zeros((2))
     cov = np.zeros((2, 2))
     self.g.process(self.fbdata, 0.5, 0)
     self.g.get_mean_cov(mean, cov)
     rv = stats.multivariate_normal(mean, cov)
     p = rv.pdf(pos)
     mp = np.trapz(p, x=r, axis=0)
     rp = np.trapz(p, x=m, axis=1)
     mp_normed = mp / np.trapz(mp, m)
     rp_normed = rp / np.trapz(rp, r)
     mhat = m[np.argmax(mp_normed)]
     rhat = r[np.argmax(rp_normed)]
     np.testing.assert_almost_equal(mhat, 5.8, decimal=1)
     np.testing.assert_almost_equal(rhat, 1.6, decimal=1)