def test_x2_in_x1_2():
    """
    x2 has a couple of bins, each of which span more than one original bin
    """
    # old size
    m = 10

    # bin edges
    x_old = np.linspace(0., 1., m+1)
    x_new = np.array([0.25, 0.55, 0.75])

    # some arbitrary distribution
    y_old = 1. + np.sin(x_old[:-1]*np.pi) / np.ediff1d(x_old)

    y_old = unp.uarray(y_old, 0.1*y_old*uniform((m,)))

    # rebin
    y_new = rebin.rebin(x_old, y_old, x_new, interp_kind='piecewise_constant')

    # compute answer here to check rebin
    y_new_here = unp.uarray(np.zeros(2), np.zeros(2))
    y_new_here[0] = 0.5 * y_old[2] + y_old[3] + y_old[4] + 0.5 * y_old[5]
    y_new_here[1] = 0.5 * y_old[5] + y_old[6] + 0.5 * y_old[7]

    assert_allclose(unp.nominal_values(y_new),
                   unp.nominal_values(y_new_here))

    # mean or nominal value comparison
    assert_allclose(unp.std_devs(y_new),
                       unp.std_devs(y_new_here))

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 other(g=True):
    freq = unumpy.uarray([100, 500, 1000, 5000, 10000, 50000], np.array([100, 500, 1000, 5000, 10000, 50000]) * 0.01)
    vin = ufloat(1.01, 0.01)
    vout = unumpy.uarray([0.640, 3.02, 5.27, 9.2, 9.6, 6.4], [0.01, 0.01, 0.01, 0.1, 0.1, 0.1])
    fase = unumpy.uarray([92, 108, 123, 166, 178, -125], [1, 2, 1, 1, 1, 1])

    Gs = vout / vin
    dB = 10 * unumpy.log10(Gs)

    if not g:
        return None

    f = plt.figure(figsize=(8, 8))
    f.suptitle("Differenziatore", fontsize=15, y=0.98)
    
    ax = f.add_subplot(211)

    ax.errorbar(x=unumpy.nominal_values(freq),
        y=unumpy.nominal_values(Gs),
        c='black', fmt='o-')

    ax.set_xlabel('Frequenza', fontsize=14)
    ax.set_ylabel('Guadagno G', fontsize=14)

    ax.set_xscale('log')
    #ax.set_ylim((-13, 1))
    #ax.set_yticklabels(('', 2, 4, 6, 8, 10, 12))
    ax.set_xticklabels(('', '100 Hz', u"1 kHz", u"10 kHz", u"100 kHz", u"1 MHz"))
    
    ax.minorticks_on()
    ax.grid(b=True, which='major', color='0.7', linestyle='-', zorder=-5)
    ax.grid(b=True, which='minor', color='0.9', linestyle='-', zorder=-9)
    ax.set_axisbelow(True)
    
    ax2 = f.add_subplot(212)

    ax2.errorbar(x=unumpy.nominal_values(freq),
        y=unumpy.nominal_values(fase),
        c='black', fmt='o-')

    ax2.set_ylabel('Sfasamento [Gradi]', fontsize=14)

    ax2.set_xscale('log')
    #ax2.set_yticklabels(('', 25, 50, 75, 100))
    ax2.set_xticklabels(('', '100 Hz', u"1 kHz", u"10 kHz", u"100 kHz"))
    
    ax2.minorticks_on()
    ax2.grid(b=True, which='major', color='0.7', linestyle='-', zorder=-5)
    ax2.grid(b=True, which='minor', color='0.9', linestyle='-', zorder=-9)
    ax2.set_axisbelow(True)

    ax3 = ax2.twiny()
    ax3.set_xticks((0, 0.333, 0.666, 1))
    ax3.set_xticklabels(('', "1 kHz", u"10 kHz", u"100 kHz"))

    f.subplots_adjust(top=0.93, hspace=0.25, bottom=0.07, right=0.95)

    plt.savefig("../latex/diff.pdf")

    plt.show()

def g_1000hz():
    g_1000hz = np.genfromtxt("../dati/g_1000hz.csv", skip_header=1, delimiter=",")
    vin = unumpy.uarray(g_1000hz[:,0] * 0.001, g_1000hz[:,1] * 0.001)
    vout = unumpy.uarray(g_1000hz[:,2], g_1000hz[:,3])

    Gs = vout / vin
    G = sum(Gs) / float(len(Gs))

    print Gs, "\n", G

    def usum(self, width=False):
        """
        Return the sum of the bin contents and uncertainty.

        See sum().
        """
        if width:
            return np.dot(np.diff(self.bins), uarray((self.hist, self.errs)))
        else:
            return np.sum(uarray((self.hist, self.errs)))

def compute_hwz(N_list, ttor, fit, plotname, title, sl=slice(None,None), Uscale=1, p0=None, eq=None, plabels=None, punits=None, Th_erw=None):
    
    N = np.sum(unp.uarray(N_list,np.sqrt(N_list)), axis=0)
    t = np.arange(len(N))*ttor+ttor/2.

    table = pt.PrettyTable()
    table.add_column('t [s]', t.astype(int), align='r')
    if len(N_list) > 1:
        for i in range(len(N_list)):
            table.add_column('N'+str(i+1), N_list[i].astype(int), align='r')
        table.add_column('Summe', N, align='r')
    else:
        table.add_column('N', N, align='r')
    with open("Resources/table_"+plotname+".txt", "w") as text_file:
        text_file.write(table.get_string())


    global N_U
    N_U = N_U0*Uscale*ttor
    popt, pstats = papstats.curve_fit(fit, t[sl], N[sl], p0=p0)

    # Untergrundfehler
    N_U = (N_U0-N_U0.s)*Uscale*ttor
    popt_min, pstats_min = papstats.curve_fit(fit, t[sl], N[sl], p0=p0)
    N_U = (N_U0+N_U0.s)*Uscale*ttor
    popt_max, pstats_max = papstats.curve_fit(fit, t[sl], N[sl], p0=p0)
    N_U = N_U0*Uscale*ttor
    s_U = unp.nominal_values(((np.abs(popt-popt_min)+np.abs(popt-popt_max))/2.))
    s_corrected = np.sqrt(unp.std_devs(popt)**2 + s_U**2)
    popt_corrected = unp.uarray(unp.nominal_values(popt),s_corrected)
    
    # Halbwertszeit
    Th = popt_corrected[::2]*unc.umath.log(2)
    for i in range(len(Th)):
        papstats.print_rdiff(Th[i]/60, Th_erw[i]/60)

    # Plot
    plt.clf()
    plt.title('Diagramm '+plotname+': '+title)
    plt.xlabel('Messzeit $t \, [s]$')
    plt.ylabel('Ereigniszahl $N$')
    xspace = np.linspace(0, t[-1])
    papstats.plot_data(t, N, label='Messpunkte')
    papstats.plot_fit(fit, popt, pstats, xspace, eq=eq, plabels=plabels, punits=punits)
    plt.fill_between(xspace, fit(xspace, *unp.nominal_values(popt_min)), fit(xspace, *unp.nominal_values(popt_max)), color='g', alpha=0.2)
    Nmin = np.amin(unp.nominal_values(N))
    for i in range(len(Th)):
        plt.hlines(popt[1::2][i].n/2.+N_U.n, 0, Th[i].n, lw=2, label='Halbwertszeit $'+papstats.pformat(Th[i], label=r'T_{\frac{1}{2}}'+('^'+str(i+1) if len(Th)>1 else ''), unit='s')+'$')
    handles, labels = plt.gca().get_legend_handles_labels()
    p = plt.Rectangle((0, 0), 1, 1, color='g', alpha=0.2)
    handles.append(p)
    labels.append('Fit im '+r'$1 \sigma$'+'-Bereich von $N_U$:'+''.join(['\n$'+papstats.pformat(s_U[i], label='\Delta '+plabels[i]+'^{U}', unit=punits[i])+'$' for i in range(len(plabels))]))
    plt.legend(handles, labels)
    papstats.savefig_a4(plotname+'.png')

def I_with_err(f_val=1, f_err=0, g_val=g_DEFAULT, g_err=0,
               e=1, e1=1, e2=E_INF):
    """Wrapper for f so the user doesn't have to know about
    the uncertainties module"""
    from uncertainties import unumpy
    f = unumpy.uarray(f_val, f_err)
    g = unumpy.uarray(g_val, g_err)
    _I = power_law_integral_flux(f, g, e, e1, e2)
    I_val = unumpy.nominal_values(_I)
    I_err = unumpy.std_devs(_I)
    return I_val, I_err

def plot_u(cal_file, fm_file, offset_file, description, accidental_offset,
        results_file):
    M = np.genfromtxt(cal_file)
    N = np.genfromtxt(fm_file)
    O = np.genfromtxt(offset_file)

    i_helm = M[:,1] #current applied to helmholtz for calibration measurement
    b_helm = M[:,2] #field applied to helmholtz coil for calibration measurement
    p, cov = np.polyfit(i_helm, b_helm, 1,  cov=True) #fit a line to calibration measurement so that we get a calibration

    i_fm = N[:,1] #current applied to helmmholtz for shielding measurement
    b_fm = unumpy.uarray(N[:,2],0.0005) - accidental_offset #field measured inside of ferromagnet shield


    B_earth = np.polyval(p,0) #We get the Earths magnetic field from i=0 of the Helmholtz calibration
    B_fm_no_i = unumpy.uarray(np.mean(O[:,2]), np.std(O[:,2])) #Get average and error for initial magnetization

    mag = B_fm_no_i - B_earth #initial magnetization is the field inside of the ferromagnet before any field is applied minus the earths magnetic field 

    Bin = b_fm - mag #internal magnetization is the measured internal field minus the initial magnetization. This correction might not be necessary for a soft ferromagnet
    
    Bext = unumpy.uarray(np.polyval(p,i_fm), 0.0005) #external field
    Bext_nom = unumpy.nominal_values(Bext)
    Bext_err = unumpy.std_devs(Bext)

    B = Bext/Bin
    c = a/b
    u=(-2*B + c**2 - 2*unumpy.sqrt(B**2 - B*c**2 - B + c**2) + 1)/(c**2 - 1)

    u_nom = unumpy.nominal_values(u)
    u_err = unumpy.std_devs(u)

    #cakculate uerr with just point to point uncertainties. I define this as just
    #uncertainty from the field measurements
    u_pp=(-2*B + c.nominal_value**2 - 2*unumpy.sqrt(B**2 - B*c.nominal_value**2 - B + c.nominal_value**2) + 1)/(c.nominal_value**2 - 1)
    
    #calculate uerr from just geometry uncertainties
    u_geom=(-2*unumpy.nominal_values(B) + c**2 - 2*unumpy.sqrt(unumpy.nominal_values(B)**2 - unumpy.nominal_values(B)*c**2 - unumpy.nominal_values(B) + c**2) + 1)/(c**2 - 1)

    ##obtain uncertainties from field
    u_err_pp = unumpy.std_devs(u_pp)

    ##obtain uncertainties from geometry
    u_err_geom = unumpy.std_devs(u_geom)

    with open(results_file, "w") as myfile:
        myfile.write('#Bext, sig_Bext, ur, sig_ur, sig_ur_pp, sig_ur_corr\n')
        for j in range(0, len(u_nom)):
            myfile.write('%s\t%s\t%s\t%s\t%s\t%s\n' %(
                Bext_nom[j], Bext_err[j],
                u_nom[j], u_err[j], u_err_pp[j], u_err_geom[j]))
    

    plt.errorbar(Bext_nom, u_nom, u_err, marker = '.', label = description)

def f_with_err(I_val=1, I_err=0, g_val=g_DEFAULT, g_err=0,
               e=1, e1=1, e2=E_INF):
    """Wrapper for f so the user doesn't have to know about
    the uncertainties module"""
    from uncertainties import unumpy
    I = unumpy.uarray(I_val, I_err)
    g = unumpy.uarray(g_val, g_err)
    _f = power_law_flux(I, g, e, e1, e2)
    f_val = unumpy.nominal_values(_f)
    f_err = unumpy.std_devs(_f)
    return f_val, f_err

def evaluate_permeability( fname_data="samplemeasurement.csv", Bin="B1", Bout="B2", fname_do="fm618_do_cryo.txt", fname_th="fm618_th_cryo.txt" ):

    print ("Evaluating permeability for: ", fname_data, Bin, Bout, fname_do, fname_th )

    # prepare result dataframe
    result = pd.DataFrame(columns = ["Bout","Bout_sdev","mu","mu_err_pp","mu_err_geom"])
    
    # get data
    data = pd.read_csv(fname_data)

    print(data.head(10))

    # set uncertainty for field reading manually to 0 if it does not exist
    if Bin + '_sdev' not in data:
        print ("Use uncertainty of 0.005 for " , Bin+'_sdev')
        data[Bin + '_sdev'] = 0.005 #Gaussmeter precision 0.01 mT above 30 mT)

    # set uncertainty for field reading manually to 0 if it does not exist
    if Bout + '_sdev' not in data:
        print ("Use uncertainty of 0.005 for " , Bout+'_sdev')
        data[Bout + '_sdev'] = 0.005

    # copy Bin and Bout to results
    result['Bout_c'] = unumpy.uarray( abs( data[Bout].values ), data[Bout + '_sdev'].values)
    result['Bin_c'] = unumpy.uarray( abs( data[Bin].values ), data[Bin + '_sdev'].values)
    
    # get inner and outer diameter
    diam_out = mean_from_file('diameter_files/'+fname_do)
    thickness = mean_from_file('diameter_files/'+fname_th)
    diam_in = diam_out - 2*thickness

    # check ideal permeability
    calc_mu_cloak(diam_in, diam_out)

    # calculate permeability and store values and uncertainties in arrays
    (mu_c, mu_err, mu_err_pp, mu_err_geom) = calc_mu(Bin=result['Bin_c'], Bout=result['Bout_c'], radius_inner=diam_in, radius_outer=diam_out)

    result["mu"] = mu_c
    result["mu_err"] = mu_err
    result["mu_err_pp"] = mu_err_pp
    result["mu_err_geom"] = mu_err_geom

    # same for Bout
    result["Bout"] = unumpy.nominal_values(result["Bout_c"])
    result["Bout_sdev"] = unumpy.std_devs(result["Bout_c"])

    #print(result.head(10))

    # drop 'mu_c' and 'Bout_c'
    result.drop('Bin_c', axis=1, inplace=True)
    result.drop('Bout_c', axis=1, inplace=True)

    return (result)

def calcGasMass(DMpc, FHI, FCO, FCOerr):
	# Calculate HI gas mass
	FHI = unumpy.uarray(FHI,(FHI*0.2))
	MHI = (2.356E5)*((DMpc)**2.)*(FHI)
	# Calculate H2 gas mass
	FCO = unumpy.uarray(FCO,FCOerr)
	MH2 = 7845*FCO*((DMpc)**2.)
	where_are_nans = unumpy.isnan(MH2)
	MH2[where_are_nans] = 0
	# Total gas mass
	Mgas = MHI + MH2
	return MHI, MH2, Mgas

def power_law_I_with_err(
    f_val=1, f_err=0, g_val=g_DEFAULT, g_err=0, e=1, e1=1, e2=E_INF
):
    """Evaluate power-law flux and propagate errors."""
    from uncertainties import unumpy

    f = unumpy.uarray(f_val, f_err)
    g = unumpy.uarray(g_val, g_err)
    _I = power_law_integral_flux(f, g, e, e1, e2)
    I_val = unumpy.nominal_values(_I)
    I_err = unumpy.std_devs(_I)
    return I_val, I_err

def power_law_f_with_err(
    I_val=1, I_err=0, g_val=g_DEFAULT, g_err=0, e=1, e1=1, e2=E_INF
):
    """Evaluate power-law ``dnde`` and propagate errors."""
    from uncertainties import unumpy

    I = unumpy.uarray(I_val, I_err)
    g = unumpy.uarray(g_val, g_err)
    _f = power_law_flux(I, g, e, e1, e2)
    f_val = unumpy.nominal_values(_f)
    f_err = unumpy.std_devs(_f)
    return f_val, f_err

def plot_u(cal_file, fm_file,  description, accidental_offset,
        results_file):
    M = np.genfromtxt(cal_file) #turn calibration file into a matrix
    N = np.genfromtxt(fm_file) #turn fm_scan file into a matrix


    i_helm = M[:,1] #current applied to helmholtz for calibration measurement
    b_helm = M[:,2] #field applied to helmholtz coil for calibration measurement
    p, cov = np.polyfit(i_helm, b_helm, 1,  cov=True) #fit a line to calibration measurement so that we get a calibration

    
    i_fm = N[:,1] #current applied to helmmholtz for shielding measurement
    Bin = unumpy.uarray(N[:,2],0.0005) - accidental_offset #field measured inside of ferromagnet shield


    Bin_nom = unumpy.nominal_values(Bin) 
    Bin_err = unumpy.std_devs(Bin)

    Bext = unumpy.uarray(np.polyval(p,i_fm), 0.0005) #external field
    Bext_nom = unumpy.nominal_values(Bext)
    Bext_err = unumpy.std_devs(Bext)

    B = Bin/Bext #ratio of internal to external field

    #calculate permeability
    u = (B*c**2 + B -2 -2*unumpy.sqrt(B**2*c**2 - B*c**2 - B + 1))/(B*c**2-B) 
    print(u)
    u_nom = unumpy.nominal_values(u)
    u_err = unumpy.std_devs(u)

    #cakculate uerr with just point to point uncertainties. I define this as just
    #uncertainty from the field measurements
    u_pp = (B*c.n**2 + B -2 -2*unumpy.sqrt(B**2*c.n**2 - B*c.n**2 - B +
        1))/(B*c.n**2-B)

    u_err_pp = unumpy.std_devs(u_pp)

    #calculate uerr from just geometry uncertainties
    u_geom = (unumpy.nominal_values(B)*c**2 + unumpy.nominal_values(B) -2 -2*unumpy.sqrt(unumpy.nominal_values(B)**2*c**2 - unumpy.nominal_values(B)*c**2 - unumpy.nominal_values(B) +
        1))/(unumpy.nominal_values(B)*c**2-unumpy.nominal_values(B))

    u_err_geom = unumpy.std_devs(u_geom)


    #write results onto a text file
    with open(results_file, "w") as myfile:
        myfile.write('#Bext, sig_Bext, Bi, sig_Bi, ur, sig_ur, sig_ur_pp, sig_ur_corr\n')
        for j in range(0, len(u_nom)):
            myfile.write('%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\n' %(
                Bext_nom[j], Bext_err[j], Bin_nom[j], Bin_err[j],
                u_nom[j], u_err[j], u_err_pp[j], u_err_geom[j]))
    
    plt.errorbar(Bext_nom, u_nom, u_err, marker = '.', label = description)

	def __init__( self, x, y, text, saveAffix, execludeLast = False):
		from ROOT import TPaveText
		self.x = x
		self.y = y
		self.text = text
		self.saveAffix = saveAffix
		error_ug = 15
		error_ux = 10 * sqrt(2)
		x = uarray(( array( x ), [ error_ug ]*len( x )))
		y = uarray(( array( y ), [ error_ux]*len( y )))
		x = x**2

		self.reg = linearRegression(x, y, execludeLast )
		self.reg.func.SetParNames('a','b')
		self.reg.draw(";U^{2}_{i} [V^{2}];U_{g} [V]" )
		self.reg.canvas.cd()
		label = TPaveText(0.1, 0.95, .86, 1, "NDC")
		label.AddText( text )
		label.SetFillStyle(0)
		label.SetBorderSize(0)
		label.Draw()
		ROOT.kPrint = 0
		'''
		ROOT.kInfo = 1000
		ROOT.kWarning = 2000;
		ROOT.kError = 3000;
		ROOT.Break = 4000;
		ROOT.kSysError = 5000;
		ROOT.kFatal = 6000;
		'''
		self.reg.canvas.SaveAs('linReg%s.pdf'%saveAffix)

		# calculate q_m
		# variable definitions:
		from math import pi
		K = 8
		w =  2. * pi * 30
		e_w = 2. * pi
		r =   0.0305 / 2
		e_r = 0.0002 / 2
		#print( '"{}  +/- {}",'.format(self.reg.func.GetParameter(0), self.reg.func.GetParError(0) ) )
		b = self.reg.func.GetParameter(1)
		if b == 0:
			return
		e_b = self.reg.func.GetParError(1)
		q_m = -2. * w**2 * r**2 * b / ( 3 * K ) * 1e6
		stat = abs( 1.* q_m * e_b / b )
		sys = abs (2. * q_m * sqrt( (e_w/w)**2 + (e_r/r)**2 ) )
		#print('{}: q/m = {:.4e} ± {:.2e} (stat) ± {:.2e} (sys) ± {:.2e} (gesamt) μC/kg\n'.format(saveAffix, q_m, stat, sys, sqrt(stat**2 + sys**2)))
		druck = { 'Luft1': 1000, 'Luft2': 1000, '375bar': 375, '400bar': 400, '425bar': 425, '425bar2': 425 }
		print('{} & {:.3g} ± {:.2g} (stat) ± {:.2g} (sys) {:2g}\\\\'.format(druck[saveAffix], -q_m, stat, sys, sqrt( stat**2 + sys**2) ))

def efficiencies(input, ht_val, alphat_binning, mhtmet_binning):
    vals = np.zeros((len(mhtmet_binning), len(alphat_binning)))
    errs = np.zeros_like(vals)
    diff = unumpy.uarray(vals, errs)
    cumu = unumpy.uarray(vals, errs)
    for imht in range(len(mhtmet_binning)):
        mhtmet_val = mhtmet_binning[imht]
        for ialphat in range(len(alphat_binning)):
            alphat_val = alphat_binning[ialphat]
            diff[imht][ialphat] = efficiency(input, ht_val, alphat_val, mhtmet_val, True)
            cumu[imht][ialphat] = efficiency(input, ht_val, alphat_val, mhtmet_val, False)
    #            if ialphat == len(alphat_binning)-1 :
    #                diff[imht][ialphat] = cumu[imht][ialphat]
    return diff, cumu

    def uintegrate(self, x1, x2, width=False):
        """
        Return the integral of the bin contents from `x1` to `x2` and
        uncertainty.

        See integrate().
        """
        i1, i2 = self.findbin([x1,x2])

        if width:
            return np.dot(np.diff(self.bins[i1:i2+2]),
                          uarray((self.hist[i1:i2+1], self.errs[i1:i2+1])))
        else:
            return np.sum(uarray((self.hist[i1:i2+1], self.errs[i1:i2+1])))

def uarray(x,errx):
	"""
With the new releases of the uncertainties and astropy.io.ascii (0.2.3, the
replacement of asciitable), if I try to create an uncertainties array with
the column of a table imported with ascii I run into trouble. For instance, 
if I use the sequence of commands below:

>>> import astropy.io.ascii as asciitable
>>> raff= asciitable.read('data/rafferty06.dat')
>>> m,errm=raff['mass'],raff['errm']
>>> mass=unumpy.uarray(m,errm)
>>> x=0.2*mass

I get the error message: 

>>> TypeError: unsupported operand type(s) for *: 'float' and 'Column'

which I can only assume is due to the new way ascii handles tables.

I created this method to use as a replacement for unumpy.uarray that handles
the tables created with astropy.io.ascii.

Usage is the same as uncertainties.unumpy.uarray.

:type x,errx: arrays created with astropy.io.ascii.
:returns: uncertainties array.
	"""
	import uncertainties.unumpy as unumpy

	x=numpy.array(x)
	errx=numpy.array(errx)

	return unumpy.uarray(x,errx)

def fit_func(func, x, y, yerr, xlong,ylong):
    start = time.time()
    p, cov = curve_fit(func, x, y, p0 = initial_guess(func,x,y), sigma = yerr,
            maxfev = 1000000)
#    print(p)
    if func == func_power:
        pi = p
        p, cov = curve_fit(func, x, y, p0 = initial_guess(func, x+p[2], y),
                sigma = yerr, maxfev = 1000000)
        while (np.any((pi/p-1)>0.001)):
            pi=p
            p, cov = curve_fit(func, x, y, p0 = initial_guess(func, x+p[2], y),
                    sigma = yerr, maxfev = 1000000)
    perr = np.sqrt(np.diag(cov))

    chi2red, pval = calculate_stats(x, y, yerr, p, perr, func)

    yfit = func(xlong,*p)
    resid = ylong-yfit
    p = unumpy.uarray(p,perr)
    extrapolate_val = func_err_prop(func,extrapolate_times,p)
    print(extrapolate_val)

    end = time.time()
    elapsed = end-start

    return (p, yfit, resid, extrapolate_val, chi2red, pval, elapsed) 

def mean(values, axis=0):
    """Returns mean values and their mean errors of a given array. Return value will be a unp.uarray
    Args:
            values:     (list)  Array containing numbers whose mean is desired.
            axis:       (int)   Axis along which the means are computed. The default is to compute the mean of the flattened array.
    """
    return unp.uarray((np.mean(noms(values), axis=axis), scipy.stats.sem(noms(values), axis=axis)))