Python uncertainties.ufloat函數代碼示例

def getPerformanceOnMC(
    files,
    histogramForEstimation="QCDStudy/PFIsolation_WithMETCutAndAsymJetCuts_DR03",
    bjetBin="",
    function="expo",
    fitRange=(0.3, 1.6),
    additionFitRanges=[(0.2, 1.6), (0.4, 1.6)],
):
    if bjetBin:
        histogramForEstimation = histogramForEstimation + "_" + bjetBin

    hists = [histogramForEstimation]
    hists = getHistsFromFiles(hists, files)
    hists = addSampleSum(hists)

    histogramForComparison = hists["qcd"][histogramForEstimation]
    histogramForEstimation = hists["allMC"][histogramForEstimation]

    estimate, absoluteError = getQCDEstimateFor(
        histogramForEstimation, function, fitRange=fitRange, additionFitRanges=additionFitRanges
    )

    qcdInSignalRegion, qcdError = getIntegral(histogramForComparison, (0, 0.1))

    N_est = ufloat((estimate, absoluteError))
    N_qcd = ufloat((qcdInSignalRegion, qcdError))

    relativeDeviation = N_est / N_qcd

    result = {}
    result["performance"] = (relativeDeviation.nominal_value, relativeDeviation.std_dev())
    result["estimate"] = (estimate, absoluteError)
    result["qcdInSignalRegion"] = (qcdInSignalRegion, qcdError)

    return result

def get_vhtt_yield(mass):
    log.warning("Get VHTT yield %s", mass)
    mmt_yield = results_file.Get('mmt_mumu_final_count/VH%i' % mass).Integral()
    emt_yield = results_file.Get('emt_emu_final_count/VH%i' % mass).Integral()

    _, (mmt_passed, mmt_total) = get_stat_error('VH%i' % mass, mmt_yield)
    _, (emt_passed, emt_total) = get_stat_error('VH%i' % mass, emt_yield)

    assert(mmt_total == emt_total)

    mmt_passed = ufloat( (mmt_passed, math.sqrt(mmt_passed) ))
    emt_passed = ufloat( (emt_passed, math.sqrt(emt_passed) ))
    total = ufloat( (mmt_total, math.sqrt(mmt_total) ))

    all_passed = mmt_passed + emt_passed

    log.warning("total %s", total)
    log.warning("all %s", all_passed)

    #multply by fraction of ZH_WH_ttH that is WH
    wh_xsec = ROOT.HiggsCSandWidth().HiggsCS(3, mass, 7)
    zh_xsec = ROOT.HiggsCSandWidth().HiggsCS(4, mass, 7)
    tth_xsec = ROOT.HiggsCSandWidth().HiggsCS(5, mass, 7)

    wh_fraction = wh_xsec/(wh_xsec + zh_xsec + tth_xsec)

    log.warning("wh_frac %s", wh_fraction)

    br_bonus = 1.0/(0.1075+0.1057+0.1125)

    log.warning("br_bonus %s", br_bonus)

    return 100*br_bonus*(all_passed/total)/wh_fraction

def split_dat(data):
    """
        Takes the row in the form r_n e_r_n ext_n e_ext_n and returns two lists
        one with the radii and the other with the extinction values both lists
        are of ufloats

        Parameters
        --------
                data : list
                A list of floats of the radii and errors and extinctions and 
                errors.

        Returns 
        --------
                rad : list
                List of ufloats of the radii.

                ext : list
                List of ufloats of the extinctions.
    """
    rad, ext = [], []

    for i in range(0, len(data) - 1, 4):
        vec = [data[i], data[i+1], data[i+2], data[i+3]]
        
        if vec[0].strip() != '':
            radius = ufloat(vec[0], vec[1])
            extinct = ufloat(vec[2], vec[3])

            rad.append(radius)
            ext.append(extinct)

    return ext, rad

 def fit_single_line_BS(self, x, y, zero_lev, err_continuum, fitting_parameters, bootstrap_iterations = 1000):
     
     #Declare parameters and containers for the fit
     params_dict     = fitting_parameters.valuesdict()
     initial_values  = [params_dict['A0'], params_dict['mu0'], params_dict['sigma0']]
     area_array      = empty(bootstrap_iterations)
     params_array    = empty([3, bootstrap_iterations])
     n_points        = len(x)
                         
     #Perform the fit
     for i in range(bootstrap_iterations):
         y_new               = y + np_normal_dist(0, err_continuum, n_points)
         area_array[i]       = simps(y_new, x) - simps(zero_lev, x)
         best_vals, covar    = curve_fit(gaussian_curveBS, (x, zero_lev), y_new, p0=initial_values, maxfev = 1600)
         params_array[:,i]   = best_vals
                 
     #Compute Bootstrap output
     mean_area, std_area = mean(area_array), std(area_array)
     mean_params_array, stdev_params_array = params_array.mean(1), params_array.std(1)
     
     #Store the data
     self.fit_dict['area_intg'],     self.fit_dict['area_intg_er']   = mean_area, std_area
     self.fit_dict['A0_norm'],       self.fit_dict['A0_norm_er']     = mean_params_array[0], stdev_params_array[0]
     self.fit_dict['mu0_norm'],      self.fit_dict['mu0_norm_er']    = mean_params_array[1], stdev_params_array[1]  
     self.fit_dict['sigma0_norm'],   self.fit_dict['sigma0_norm_er'] = mean_params_array[2], stdev_params_array[2]  
             
     A = ufloat(mean_params_array[0], stdev_params_array[0]) 
     sigma = ufloat(mean_params_array[2], stdev_params_array[2])         
     fwhm0_norm = 2.354820045 * sigma
     areaG0_norm = A * sigma * self.sqrt2pi
   
     self.fit_dict['fwhm0_norm'], self.fit_dict['fwhm0_norm_er'] = fwhm0_norm.nominal_value, fwhm0_norm.std_dev
     self.fit_dict['area_G0_norm'], self.fit_dict['area_G0_norm_er'] = areaG0_norm.nominal_value, areaG0_norm.std_dev
                        
     return

def line_fit(x, y, errors=True):
    """
    Simple helper function that speeds up line fitting. Uses lmfit
    
    Parameters
    ----------
    x, y (float)
        Sample length x, y arrays of data to be fitted

    Returns
    -------
    Returns slope and intercept, with uncertainties (use uncertainties package if availabe)
    Also returns fit object out, can be dropped
    """
    from lmfit.models import LinearModel
    mod = LinearModel()
    par = mod.guess(y, x=x)
    out = mod.fit(y, par, x=x)

    s = out.params['slope']
    i = out.params['intercept']

    if errors:
        try:
            from uncertainties import ufloat
            return ufloat(s.value, s.stderr), ufloat(i.value, i.stderr), out
        except:
            return s.value, s.stderr, i.value, i.stderr, out
    else:
        return s.value, i.value, out

def calculate_latencies(version_dates):
    linux_latencies = latency(version_dates['linux'], OrderedDict(avo.os_to_kernel))
    set_latex_value('linuxMeanUpdateLatency', ufloat(statistics.mean(linux_latencies.values()),statistics.stdev(linux_latencies.values())))
    openssl_latencies = latency(version_dates['openssl'], OrderedDict(avo.os_to_project['openssl']))
    set_latex_value('opensslMeanUpdateLatency', ufloat(statistics.mean(openssl_latencies.values()),statistics.stdev(openssl_latencies.values())))
    bouncycastle_latencies = latency(version_dates['bouncycastle'], OrderedDict(avo.os_to_project['bouncycastle']))
    set_latex_value('bouncycastleMeanUpdateLatency',ufloat(statistics.mean(bouncycastle_latencies.values()),statistics.stdev(bouncycastle_latencies.values())))

def combine_complex_df( df1, df2 ):
	'''
	Takes a 2 pandii dataframes of the form:
		A 	| 	B   	 	A 	| 	B 
	  (v,e) | (v,e)		  (v,e) | (v,e) 

	Returns 1 pandas dataframe of the form
			      A   |   B 
	  			(v,e) | (v,e)
	'''
	from uncertainties import ufloat
	l1=df1.columns.tolist()
	l2=df2.columns.tolist()
	if l1 != l2:
		print "Trying to combine two non compatible dataframes"
		print l1
		print l2
		return

	combined_result = {}
	for sample in l1:
		results = []
		for entry1, entry2 in zip(df1[sample], df2[sample]):
			v1 = ufloat(entry1[0], entry1[1])
			v2 = ufloat(entry2[0], entry2[1])
			s = v1 + v2
			results.append( ( s.nominal_value, s.std_dev ) )
		combined_result[sample] = results
	df = dict_to_df(combined_result)
	return df

def calculate_flux(f, age, arar_constants=None):
    """
        #rad40: radiogenic 40Ar
        #k39: 39Ar from potassium
        f: F value rad40Ar/39Ar
        age: age of monitor in years

        solve age equation for J
    """
    # if isinstance(rad40, (list, tuple)):
    # rad40 = ufloat(*rad40)
    # if isinstance(k39, (list, tuple)):
    # k39 = ufloat(*k39)

    if isinstance(f, (list, tuple)):
        f = ufloat(*f)

    if isinstance(age, (list, tuple)):
        age = ufloat(*age)
        #    age = (1 / constants.lambdak) * umath.log(1 + JR)
    try:
        # r = rad40 / k39
        if arar_constants is None:
            arar_constants = ArArConstants()

        j = (umath.exp(age * arar_constants.lambda_k.nominal_value) - 1) / f
        return j.nominal_value, j.std_dev
    except ZeroDivisionError:
        return 1, 0

def small_sep02 ( mx , mxErr, lag=0):
	""" 
	Given a frequency matrix and errors calculates the (scaled) small 
	separation d02 and propagates errors.
	
	Notes
	-----
	The parameter lag is the difference between the radial orders of
	the first modes of l=0 and l=2.  
	"""

	(Nn, Nl) = np.shape(mx)
	d02    = np.zeros((1,Nn))
	d02Err = np.zeros((1,Nn))
	d02.fill(None)      # values that can't be calculated are NaN

	
	for n in range(1-lag,Nn):
		if (mx[n,0] != 0. and mx[n-1+lag,2] != 0.):
		    a = un.ufloat( (mx[n,0], mxErr[n,0]) )
		    b = un.ufloat( (mx[n-1+lag,2], mxErr[n-1+lag,2]) )
		    result = (a-b) / 3.	
		    d02[0,n-1+lag]    = un.nominal_value(result)
		    d02Err[0,n-1+lag] = un.std_dev(result)
		
	return d02, d02Err

def day3_vacuum():
    plt.clf()
    ux = unp.uarray([800, 1000, 500], ux_error)  # V
    ug_crit = unp.uarray([110, 140, 30], ug_error)  # V
    omega = 2 * math.pi * ufloat(48, 1)  # Hz

    ug_crit = _apply_additional_proportional_error(ug_crit, stability_uncertainty_vacuum)
    ux *= ux_correction
    ug_crit *= ug_correction

    stability(ux, ug_crit, omega, label="Particle V1")

    plt.title("p = 300 mbar")
    plt.legend(loc=2)
    plt.savefig("images/stability_vacuum_1.pdf")

    plt.clf()
    ux = unp.uarray([700, 800], ux_error)  # V
    ug_crit = unp.uarray([100, 140], ug_error)  # V
    omega = 2 * math.pi * ufloat(45, 1)  # Hz

    ug_crit = _apply_additional_proportional_error(ug_crit, stability_uncertainty_vacuum)
    ux *= ux_correction
    ug_crit *= ug_correction

    stability(ux, ug_crit, omega, label="Particle V2")

    plt.title("p = 180 mbar")
    plt.legend(loc=2)
    plt.savefig("images/stability_vacuum_2.pdf")

def safe_div(n, d):
    if n.size == 0:
        return None
    if n.shape != d.shape:
        return None
    if isinstance(n[0, 0], UFloat) is False:
        return None
    z = np.zeros_like(n)
    z.fill(ufloat(np.nan, np.nan))
    it = np.nditer(d, op_flags=["readonly"], flags=["multi_index", "refs_ok"])
    while not it.finished:
        i = it.multi_index[0]
        j = it.multi_index[1]
        if abs(d[i, j].n) == 0.0:
            # z[i,j] = ufloat(0.,0.)
            z[i, j] = ufloat(np.nan, np.nan)
        #            if abs(n[i,j].n) > 0. :
        #                z[i,j] = ufloat(np.nan,np.nan)
        #            else :
        #                err = np.sqrt( n[i,j].s**2. + d[i,j].s**2. )
        #                z[i,j] = ufloat(0.,err)
        # z[i,j] = ufloat(np.nan,np.nan)
        # z[i,j] = ufloat(np.inf,np.inf)
        else:
            z[i, j] = n[i, j] / d[i, j]
        it.iternext()
    return z

def day2_air():
    plt.clf()
    ux = unp.uarray([1000, 860, 800, 760, 1000, 920, 620], ux_error)  # V
    ug_crit_1 = unp.uarray([390, 170, 600, 120, 225, 175, 53], ug_error)  # V
    ug_crit_2 = unp.uarray([390, 170, 600, 120, 160, 145, 67], ug_error)  # V
    ug_crit_3 = unp.uarray([390, 370, 600, 120, 225, 175, 67], ug_error)  # V
    omega = 2 * math.pi * ufloat(28, 1)  # Hz

    ux *= ux_correction
    for i, ug_crit in enumerate((ug_crit_1, ug_crit_2, ug_crit_3), 1):
        ug_crit = _apply_additional_proportional_error(ug_crit, stability_uncertainty_air)
        ug_crit *= ug_correction

        stability(ux, ug_crit, omega, "Particle A%d" % i)

    plt.legend(loc=2)
    plt.savefig("images/stability_air_1.pdf")

    plt.clf()
    ux = unp.uarray([550, 740, 870, 1040], ux_error)
    ug_crit = unp.uarray([39, 84, 133, 147], ug_error)
    omega = 2 * math.pi * ufloat(32, 1)  # Hz

    ug_crit = _apply_additional_proportional_error(ug_crit, stability_uncertainty_air)
    ux *= ux_correction
    ug_crit *= ug_correction

    stability(ux, ug_crit, omega, label="Particle A4")

    plt.legend(loc=2)
    plt.savefig("images/stability_air_2.pdf")

def get_1d_loose_efficiencies( int_stat, int_syst, lead_reg, lead_ptrange, systematics=None) :

    eff_stat = {}
    eff_syst = {}

    if int_stat['lead']['real']['loose'].n == 0 :
        eff_stat['eff_R_T_lead'] = ufloat( 1.0, int_stat['lead']['real']['tight'].s/int_stat['lead']['real']['tight'].n )
    else :
        eff_stat['eff_R_T_lead'] = int_stat['lead']['real']['tight'] / (int_stat['lead']['real']['tight']+int_stat['lead']['real']['loose'])
    eff_stat['eff_F_T_lead'] = int_stat['lead']['fake']['tight'] / (int_stat['lead']['fake']['tight']+int_stat['lead']['fake']['loose'])
    eff_stat['eff_R_L_lead'] = int_stat['lead']['real']['loose'] / (int_stat['lead']['real']['tight']+int_stat['lead']['real']['loose'])
    eff_stat['eff_F_L_lead'] = int_stat['lead']['fake']['loose'] / (int_stat['lead']['fake']['tight']+int_stat['lead']['fake']['loose'])

    eff_syst['eff_R_T_lead'] = int_syst['lead']['real']['tight'] / (int_syst['lead']['real']['tight']+int_syst['lead']['real']['loose'])
    eff_syst['eff_F_T_lead'] = int_syst['lead']['fake']['tight'] / (int_syst['lead']['fake']['tight']+int_syst['lead']['fake']['loose'])
    eff_syst['eff_R_L_lead'] = int_syst['lead']['real']['loose'] / (int_syst['lead']['real']['tight']+int_syst['lead']['real']['loose'])
    eff_syst['eff_F_L_lead'] = int_syst['lead']['fake']['loose'] / (int_syst['lead']['fake']['tight']+int_syst['lead']['fake']['loose'])


    # Do systematics
    # the integrals may already have some systematics
    # that are propagated to the eff_*
    # therefore, don't overwrite the existing 
    # systematics, but make a ufloat with a
    # zero value, and non-zero syst
    eff_syst['eff_R_L_lead'] = eff_syst['eff_R_L_lead'] + ufloat( 0.0, math.fabs(eff_syst['eff_R_L_lead'].n)*get_syst_uncertainty( 'RealTemplateNom', lead_reg, lead_ptrange, 'real', 'loose' ), 'Template_lead_real_loose')
    eff_syst['eff_F_L_lead'] = eff_syst['eff_F_L_lead'] + ufloat( 0.0, math.fabs(eff_syst['eff_F_L_lead'].n)*get_syst_uncertainty( 'FakeTemplate%s'%systematics, lead_reg, lead_ptrange, 'fake', 'loose' ), 'Template_lead_fake_loose' )
    eff_syst['eff_R_T_lead'] = eff_syst['eff_R_T_lead'] + ufloat( 0.0, math.fabs(eff_syst['eff_R_T_lead'].n)*get_syst_uncertainty( 'RealTemplateNom', lead_reg, lead_ptrange, 'real', 'loose' ), 'Template_lead_real_loose')
    eff_syst['eff_F_T_lead'] = eff_syst['eff_F_T_lead'] + ufloat( 0.0, math.fabs(eff_syst['eff_F_T_lead'].n)*get_syst_uncertainty( 'FakeTemplate%s'%systematics, lead_reg, lead_ptrange, 'fake', 'loose' ), 'Template_lead_fake_loose' )

    return eff_stat, eff_syst

    def _set_age_values(self, f, include_decay_error=False):
        arc = self.arar_constants
        j = copy(self.j)
        if j is None:
            j = ufloat(1e-4, 1e-7)
        j.tag = 'Position'
        j.std_dev = self.position_jerr or 0
        age = age_equation(j, f, include_decay_error=include_decay_error, arar_constants=arc)
        self.uage_w_position_err = age

        j = self.j
        if j is None:
            j = ufloat(1e-4, 1e-7, tag='J')

        age = age_equation(j, f, include_decay_error=include_decay_error, arar_constants=arc)
        self.uage_w_j_err = age

        j = copy(self.j)
        if j is None:
            j = ufloat(1e-4, 1e-7, tag='J')

        j.std_dev = 0
        age = age_equation(j, f, include_decay_error=include_decay_error, arar_constants=arc)
        self.uage = age

        self.age = nominal_value(age)
        self.age_err = std_dev(age)
        self.age_err_wo_j = std_dev(age)

        for iso in self.itervalues():
            iso.age_error_component = self.get_error_component(iso.name)

 def test_decimal_places(self):
     cases = [(0, '0'), (-1, '-1'), (0.012, '0.01'), (0.016, '0.02')
             , (ufloat(0,0), r'$0.0 \pm 0.0$'), (ufloat(10,0), r'$10.0 \pm 0.0$')
             , (ufloat(0.12, 0.012), r'$0.12 \pm 0.01$'), (ufloat(0.12, 0.016), r'$0.12 \pm 0.02$')
             , (ufloat(-0.12, 0.016), r'$-0.12 \pm 0.02$')]
     for test_case, result in cases:
         self.assertEqual(display_num(test_case, decimal_places=2), result)