def time_fitting(x_fit,y_fit):
    """Fit a linear relation to the x_fit and y_fit parameters

    Returns the actual fit and the parameters of the fit,
    """
    import numpy as np
    x_fit = np.array( x_fit )
    y_fit = np.array( y_fit )

    ###First fit iteration and remove outliers
    POLY_FIT_ORDER = 1

    slope,intercept = scipy.polyfit(x_fit,y_fit,POLY_FIT_ORDER)
    fit = scipy.polyval((slope,intercept),x_fit)
    fit_sigma = fit.std()
    include_index = np.where(np.abs(fit-y_fit) < 1.5*fit_sigma)[0]

    if len(include_index) < 4:
        return None,None,False

    ###Final Fit
    x_fit_clipped = x_fit[include_index]
    y_fit_clipped = y_fit[include_index]

    parameters = scipy.polyfit(x_fit_clipped,y_fit_clipped,POLY_FIT_ORDER)
    fit = scipy.polyval(parameters,x_fit)

    return fit,parameters,True

    def x_do_contact(self, args):
        """
        DEBUG COMMAND to be activated in the future
        """
        xext, ysub, contact = self.find_contact_point2(debug=True)

        contact_plot = self.plots[0]
        contact_plot.add_set(xext, ysub)
        contact_plot.add_set([xext[contact]], [self.plots[0].vectors[0][1][contact]])
        # contact_plot.add_set([first_point[0]],[first_point[1]])
        # contact_plot.add_set([last_point[0]],[last_point[1]])
        contact_plot.styles = [None, None, None, "scatter"]
        self._send_plot([contact_plot])
        return

        index, regr, regr_contact = self.find_contact_point2(debug=True)
        print regr
        print regr_contact
        raw_plot = self.current.curve.default_plots()[0]
        xret = raw_plot.vectors[0][0]
        # nc_line=[(item*regr[0])+regr[1] for item in x_nc]
        nc_line = scipy.polyval(regr, xret)
        c_line = scipy.polyval(regr_contact, xret)

        contact_plot = self.current.curve.default_plots()[0]
        contact_plot.add_set(xret, nc_line)
        contact_plot.add_set(xret, c_line)
        contact_plot.styles = [None, None, None, None]
        # contact_plot.styles.append(None)
        contact_plot.destination = 1
        self._send_plot([contact_plot])

def test_scipy():
    #Sample data creation
    #number of points
    n=50
    t=linspace(-5,5,n)
    #parameters
    a=0.8; b=-4
    x=polyval([a,b],t)
    #add some noise
    xn=x+randn(n)

    #Linear regressison -polyfit - polyfit can be used other orders polys
    (ar,br)=polyfit(t,xn,1)
    xr=polyval([ar,br],t)
    #compute the mean square error
    err=sqrt(sum((xr-xn)**2)/n)

    print('Linear regression using polyfit')
    print('parameters: a=%.2f b=%.2f \nregression: a=%.2f b=%.2f, ms error= %.3f' % (a,b,ar,br,err))

    #matplotlib ploting
    title('Linear Regression Example')
    plot(t,x,'g.--')
    plot(t,xn,'k.')
    plot(t,xr,'r.-')
    legend(['original','plus noise', 'regression'])

    show()

    #Linear regression using stats.linregress
    (a_s,b_s,r,tt,stderr)=stats.linregress(t,xn)
    print('Linear regression using stats.linregress')
    print('parameters: a=%.2f b=%.2f \nregression: a=%.2f b=%.2f, std error= %.3f' % (a,b,a_s,b_s,stderr))

def main(FileName):
    
    for trial in kep.iofiles.passFileToTrialList(FileName):
        tr = kep.pipelinepars.keptrial(trial)
        kw = kep.keplc.kw(**tr.kw)
        X = kep.keplc.keplc(tr.kid)
        X.runPipeline(kw)
        lc = X.lcData.lcData
        Y = kep.qats.qatslc(X.lcFinal,X.KID)
        idx = num.where((X.lcFinal['eMask'] == False))[0]
        Y.padLC(flagids=idx)
        Y.addNoise()
        Y.runQATS(f=0.01)
        
        P,ignored = calcPeriodogram(Y.periods, lc['x'], lc['ydt'], lc['yerr']) ### ADDED ###
                
        coeff = num.polyfit(num.log10(Y.periods),num.log10(Y.snrLC),1)
        outy = scipy.polyval(coeff,num.log10(Y.periods))
        normalizedPower = 10**(outy)
                
        coeff2 = num.polyfit(num.log10(Y.periods),num.log10(Y.snrFLAT),1) ### ADDED ###
        outy2 = scipy.polyval(coeff2,num.log10(Y.periods)) ### ADDED ###
        normalizedPower2 = 10**(outy2) ### ADDED ###
                
        fittedr=Y.snrFLAT/normalizedPower2 ### ADDED ###
        fittedb=Y.snrLC/normalizedPower2 ### ADDED ###
        normalizedPower3 = normalizedPower2/normalizedPower2 ### ADDED ###
                
        plt.subplot(211) ### ADDED ###
        plt.title('unflipped,'+' KID = '+X.KID)
        plt.plot(Y.periods,Y.snrFLAT/normalizedPower,'r-') 
        plt.plot(Y.periods,Y.snrLC/normalizedPower,'b-')
        plt.plot(Y.periods,normalizedPower/normalizedPower,'k-')
                #plt.plot(Y.periods,normalizedPower2, 'k-') ### ADDED ###
        plt.setp(plt.gca().set_xscale('log'))
        plt.ylabel('Signal Power')
        plt.subplot(212)
        plt.plot(Y.periods,P,'g-') ### ADDED ###
        plt.setp(plt.gca().set_xscale('log')) ### ADDED ###
        plt.xlabel('Period (days)') ### MOVED ###
        plt.ylabel('F Transform')
        plt.savefig('sn.SG0XX.unflipped.'+X.KID+'.png') ### MOVED ###
        #coeff = num.polyfit(num.log10(Y.periods),num.log10(Y.snrLC),1)
        #outy = scipy.polyval(coeff,num.log10(Y.periods))
        #normalizedPower = 10**(outy)
        
        #plt.plot(Y.periods,Y.snrFLAT,'r-') 
        #plt.plot(Y.periods,Y.snrLC,'b-')
        #plt.plot(Y.periods,normalizedPower,'k-')
        #plt.setp(plt.gca().set_xscale('log'))
        #plt.savefig('sn.'+X.KID+'.png')
        
        dfile = open('signal.'+X.KID+'.data','w')
        print >> dfile,'#',X.KID,'|', Y.periods[num.argmax(Y.snrLC/normalizedPower)],\
                Y.periods[num.argmax(Y.snrLC)], max(Y.snrLC/normalizedPower), max(Y.snrLC)
        for i in range(len(Y.SignalPower)):
            print >> dfile, Y.periods[i],'|',Y.snrLC[i],'|',\
                            Y.snrFLAT[i],'|',normalizedPower[i]
        dfile.close()

def go(
    mags=(15, 30.1, 0.5), redshifts=(0.01, 12, 0.01), coeff_file="prior_K_zmax7_coeff.dat", outfile="prior_K_extend.dat"
):

    fp = open(coeff_file)
    lines = fp.readlines()
    fp.close()

    mag_list = np.cast[float](lines[0].split()[2:])
    z0 = np.cast[float](lines[1].split()[1:])
    gamma = np.cast[float](lines[2].split()[1:])

    z_grid = np.arange(redshifts[0], redshifts[1], redshifts[2])
    NZ = z_grid.shape[0]

    mag_grid = np.arange(mags[0], mags[1], mags[2])
    NM = mag_grid.shape[0]

    #### Polynomial extrapolation not reliable
    # p_z0 = scipy.polyfit(mag_list, z0, order)
    # z0_grid = np.maximum(scipy.polyval(p_z0, mag_grid), 0.05)
    # p_gamma = scipy.polyfit(mag_list, gamma, order)
    # gamma_grid = np.maximum(scipy.polyval(p_gamma, mag_grid), 0.05)

    #### Interpolations on the defined grid
    z0_grid = np.interp(mag_grid, mag_list, z0)
    gamma_grid = np.interp(mag_grid, mag_list, gamma)

    #### Linear extrapolations of fit coefficients
    p_z0 = scipy.polyfit(mag_list[:3], z0[:3], 1)
    z0_grid[mag_grid < mag_list[0]] = np.maximum(scipy.polyval(p_z0, mag_grid[mag_grid < mag_list[0]]), 0.05)
    p_z0 = scipy.polyfit(mag_list[-3:], z0[-3:], 1)
    z0_grid[mag_grid > mag_list[-1]] = np.maximum(scipy.polyval(p_z0, mag_grid[mag_grid > mag_list[-1]]), 0.05)

    p_gamma = scipy.polyfit(mag_list[:3], gamma[:3], 1)
    gamma_grid[mag_grid < mag_list[0]] = np.maximum(scipy.polyval(p_gamma, mag_grid[mag_grid < mag_list[0]]), 0.05)
    p_gamma = scipy.polyfit(mag_list[-3:], gamma[-3:], 1)
    gamma_grid[mag_grid > mag_list[-1]] = np.maximum(scipy.polyval(p_gamma, mag_grid[mag_grid > mag_list[-1]]), 0.05)

    out_matrix = np.zeros((NZ, NM + 1))
    out_matrix[:, 0] = z_grid

    for i in range(NM):
        pz = z_grid * np.exp(-(z_grid / z0_grid[i]) ** gamma_grid[i])
        pz /= np.trapz(pz, z_grid)
        plt.plot(z_grid, pz, label=mag_grid[i])
        out_matrix[:, i + 1] = pz

    plt.legend(ncol=3, loc="upper right")

    header = "# z "
    for m in mag_grid:
        header += "%6.1f" % (m)

    fp = open(outfile, "w")
    fp.write(header + "\n")
    np.savetxt(fp, out_matrix, fmt="%6.3e")
    fp.close()

def concat_extrap_ends(x, npts, polyorder=1, lowside=True, highside=True):
    i=numpy.arange(npts, dtype='float64')
    if lowside:
        ans=scipy.polyfit(-1*(i+1.), x[:npts], polyorder)
        x=numpy.concatenate([scipy.polyval(list(ans), i[::-1]), x])
    if highside:
        ans=scipy.polyfit(-1*(i[::-1]-1.), x[-1*npts:], polyorder)
        x=numpy.concatenate([x, scipy.polyval(list(ans), i)])
    return x    

	def mNm_2_raw(self,calib,tq):
		#ticks: calibrated sensor value
		#calib: yaml config calibration map
		if self.ctype=='sea_vertx_14bit':
			val=(tq-calib['cb_bias'])/calib['cb_scale']
			val=float(polyval(calib['cb_inv_torque'],val))
			return max(0,min(val,VERTX_14BIT_MAX))
		if self.ctype=='adc_poly':
			val=(tq-calib['cb_bias'])/calib['cb_scale']
			val=float(polyval(calib['cb_inv_torque'],val))
			return max(0,min(val,M3EC_ADC_TICKS_MAX))

    def calculo(self):
        entrada = self.kwargs["entrada"]
        self.rendimientoCalculado = Dimensionless(self.kwargs["rendimiento"])

        if self.kwargs["Pout"]:
            DeltaP = Pressure(self.kwargs["Pout"]-entrada.P)
        elif self.kwargs["deltaP"]:
            DeltaP = Pressure(self.kwargs["deltaP"])
        elif self.kwargs["Carga"]:
            DeltaP = Pressure(self.kwargs["Carga"]*entrada.Liquido.rho*g)
        else:
            DeltaP = Pressure(0)

        if self.kwargs["usarCurva"]:
            b1 = self.kwargs["diametro"] != self.kwargs["curvaCaracteristica"][0]  # noqa
            b2 = self.kwargs["velocidad"] != self.kwargs["curvaCaracteristica"][1]  # noqa
            if b1 or b2:
                self.curvaActual = self.calcularCurvaActual()
            else:
                self.curvaActual = self.kwargs["curvaCaracteristica"]
            self.Ajustar_Curvas_Caracteristicas()

        if not self.kwargs["usarCurva"]:
            head = Length(DeltaP/g/entrada.Liquido.rho)
            power = Power(head*g*entrada.Liquido.rho*entrada.Q /
                          self.rendimientoCalculado)
            P_freno = Power(power*self.rendimientoCalculado)
        elif not self.kwargs["incognita"]:
            head = Length(polyval(self.CurvaHQ, entrada.Q))
            DeltaP = Pressure(head*g*entrada.Liquido.rho)
            power = Power(entrada.Q*DeltaP)
            P_freno = Power(polyval(self.CurvaPotQ, entrada.Q))
            self.rendimientoCalculado = Dimensionless(power/P_freno)
        else:
            head = Length(self.DeltaP/g/entrada.Liquido.rho)
            poli = [self.CurvaHQ[0], self.CurvaHQ[1], self.CurvaHQ[2]-head]
            Q = roots(poli)[0]
            power = Power(Q*self.DeltaP)
            entrada = entrada.clone(split=Q/entrada.Q)
            P_freno = Power(polyval(self.CurvaPotQ, Q))
            self.rendimientoCalculado = Dimensionless(power/P_freno)

        self.deltaP = DeltaP
        self.headCalculada = head
        self.power = power
        self.P_freno = P_freno
        self.salida = [entrada.clone(P=entrada.P+DeltaP)]
        self.Pin = entrada.P
        self.PoutCalculada = self.salida[0].P
        self.volflow = entrada.Q
        self.cp_cv = entrada.Liquido.cp_cv

def estimate_rate_func(t, T, N, plot_flag=False, method='central diff'):

    t_est_pts = scipy.linspace(t.min(), t.max(), N+2) 
    interp_func = scipy.interpolate.interp1d(t,T,'linear')
    T_est_pts = interp_func(t_est_pts)
 
    if plot_flag == True:
        pylab.figure()
        pylab.subplot(211)
        pylab.plot(t_est_pts, T_est_pts,'or')

    # Estimate slopes
    slope_pts = scipy.zeros((N,)) 
    T_slope_pts = scipy.zeros((N,))  

    if method == 'local fit':
        for i in range(1,(N+1)):
            mask0 = t > 0.5*(t_est_pts[i-1] + t_est_pts[i])
            mask1 = t < 0.5*(t_est_pts[i+1] + t_est_pts[i])
            mask = scipy.logical_and(mask0, mask1)
            t_slope_est = t[mask]
            T_slope_est = T[mask]
            local_fit = scipy.polyfit(t_slope_est,T_slope_est,2)
            dlocal_fit = scipy.polyder(local_fit)
            slope_pts[i-1] = scipy.polyval(dlocal_fit,t_est_pts[i]) 
            T_slope_pts[i-1] = scipy.polyval(local_fit,t_est_pts[i])
            if plot_flag == True:
                t_slope_fit = scipy.linspace(t_slope_est[0], t_slope_est[-1], 100)
                T_slope_fit = scipy.polyval(local_fit,t_slope_fit)
                pylab.plot(t_slope_fit, T_slope_fit,'g')
    elif method == 'central diff':
        dt = t_est_pts[1] - t_est_pts[0]
        slope_pts = (T_est_pts[2:] - T_est_pts[:-2])/(2.0*dt)
        T_slope_pts = T_est_pts[1:-1]
    else:
        raise ValueError, 'unkown method %s'%(method,)
        
    # Fit line to slope estimates
    fit = scipy.polyfit(T_slope_pts, slope_pts,1)

    if plot_flag == True:
        T_slope_fit = scipy.linspace(T_slope_pts.min(), T_slope_pts.max(),100)
        slope_fit = scipy.polyval(fit,T_slope_fit)
        pylab.subplot(212)
        pylab.plot(T_slope_fit, slope_fit,'r')
        pylab.plot(T_slope_pts, slope_pts,'o')
        pylab.show()

    rate_slope = fit[0]
    rate_offset = fit[1]
    return rate_slope, rate_offset 

def fit_curve():
	global data, running
	if running == True or len(data[0])==0:
		return
	aa = []
	bb = []
	for k in range(len(data[0])):
		if data[1][k] > 1.0:
			aa.append(data[0][k])
			bb.append(data[1][k])
	x  = array(bb)
	y  = array(aa)
	from scipy import polyfit, polyval
	(ar,br)=polyfit(x,y,1)
	print polyval([ar,br],[0])

def interpMaskedAreas(x, dtfunc):
    idx = num.where(dtfunc == -1)[0]
    # find consecutive points of dtfunc set at -1
    diffs = idx - num.roll(idx, 1)
    setIdx1 = num.where(diffs > 1)[0]
    setIdx2 = num.where(diffs > 1)[0] - 1
    setIdx1 = num.hstack( (num.array([0]), setIdx1) )
    setIdx2 = num.hstack( (setIdx2, num.array([-1])) )
    # each int in setIdx1 is beginning of a set of points
    # each int in setIdx2 is end of a set of points
    for i in range(len(setIdx1)):
        i1 = setIdx1[i]
        i2 = setIdx2[i]
        vals_to_interp = dtfunc[ idx[i1]:idx[i2]+1]
        
        # construct interpolation from points on either side of vals
        # first ydata
        nVals = len(vals_to_interp)
        left = num.array( \
            [ dtfunc[idx[i1]-1-nVals], dtfunc[idx[i1]-1] ])
        right = num.array( \
            [ dtfunc[idx[i2]+1], dtfunc[idx[i2]+1+nVals] ])
        vals_to_construct_interp = num.hstack( (left, right) )
        vtci = vals_to_construct_interp
        
        # now xdata
        leftx = num.array([ x[idx[i1]-1-nVals], x[idx[i1]-1] ])
        rightx = num.array([ x[idx[i2]+1], x[idx[i2]+1+nVals] ])
        x_vtci = num.hstack( (leftx, rightx) )
        # conduct interpolation
        intpCoeffs = scipy.polyfit(x_vtci, vtci, 3)
        dtfunc[ idx[i1]:idx[i2]+1 ] = \
            scipy.polyval(intpCoeffs, x[ idx[i1]:idx[i2]+1 ])
        
    return dtfunc

def construct_linear_regression(x_learn, y_dev_learn, x_test, y_dev_test):

    regression = sp.polyfit(x_learn, y_dev_learn, 1)

    y_exp_learn = sp.polyval(regression, x_learn)
    y_exp_test = sp.polyval(regression, x_test)

    print("For train dataset:\n\ty = a*x + b\n\ta = %f\n\tb = %f" % (regression[0], regression[1]))

    mse_learn = np.sqrt(np.mean((y_exp_learn - y_dev_learn) ** 2))
    mse_test = np.sqrt(np.mean((y_exp_test - y_dev_test) ** 2))
    mse_total = np.sqrt((((mse_learn**2) * learn_set_size)
                         + ((mse_test**2) * (dataset_size - learn_set_size)))
                        / dataset_size)
    print("Train MSE = %f\nTest MSE = %f\nTotal MSE = %f" % (mse_learn, mse_test,mse_total ))
    return regression, mse_learn, mse_test, mse_total, y_exp_learn, y_exp_test

def calculatePolyfit(x, y):
    """ Calculate a linear regression using polyfit instead """

    # FIXME(pica) This doesn't work for large ranges (tens of orders of
    # magnitude). No idea why as the fixRange() function above should make
    # all values greater than one. The trouble mainly seems to happen when
    # log10(min(x)) negative.

    # Check the smallest value in x or y isn't below 2x10-7 otherwise
    # we hit numerical instabilities.
    xFactorFix, xFactor = False, 0
    yFactorFix, yFactor = False, 0
    minX = min(x)
    minY = min(y)
    if minX < 2e-7:
        x, xFactor, xFactorFix = fixRange(x)
    if minY < 2e-7:
        y, yFactor, yFactorFix = fixRange(y)

    (ar, br) = polyfit(x, y, 1)
    yf = polyval([ar, br], x)
    xf = x

    if xFactorFix:
        xf = xf / 10**xFactor
    if yFactorFix:
        yf = yf / 10**yFactor

    return xf, yf

def show_dialog_len():
    import numpy as np
    from scipy import polyfit,polyval
    import matplotlib.pyplot as plt
    import LetsgoCorpus as lc
    import operator

    corpus = lc.Corpus('D:/Data/training',prep=True)
    
    l = [];avg_cs = []
    for dialog in corpus.dialogs():
        if len(dialog.turns) > 40:
            continue
        l.append(len(dialog.turns))
        avg_cs.append(reduce(operator.add,map(lambda x:x['CS'],dialog.turns))/l[-1])
    (ar,br) = polyfit(l,avg_cs,1)
    csr = polyval([ar,br],l)
    plt.plot(l,avg_cs,'g.',alpha=0.75)
    plt.plot(l,csr,'r-',alpha=0.75,linewidth=2)
    plt.axis([0,50,0,1.0])
    plt.xlabel('Dialog length (turn)')
    plt.ylabel('Confidence score')
    plt.title('Dialog length vs. Confidence score')
    plt.grid(True)
    plt.savefig('img/'+'Dialog length vs Confidence score'+'.png')
#    plt.show()
    plt.clf()

	def raw_2_mA(self,calib,ticks_a,ticks_b):
		#ticks: adc sensor value
		#calib: yaml config calibration map
		if self.ctype=='adc_linear_5V':
			mV_a = 5000.0*(ticks_a-calib['cb_ticks_at_zero_a'])/M3EC_ADC_TICKS_MAX
			mV_b = 5000.0*(ticks_b-calib['cb_ticks_at_zero_b'])/M3EC_ADC_TICKS_MAX
			i_a = 1000.0*mV_a/calib['cb_mV_per_A']
			i_b = 1000.0*mV_b/calib['cb_mV_per_A']
			val=max(abs(i_a), abs(i_b))
			return max(0,val*calib['cb_scale']+calib['cb_bias'])	
		if self.ctype=='adc_poly':
			i_a= float(polyval(calib['cb_current_a'],ticks_a))
			i_b= float(polyval(calib['cb_current_b'],ticks_b))
			val=max(abs(i_a), abs(i_b))
			return val*calib['cb_scale']+calib['cb_bias']
		return 0.0

def main(FileName):
    
    for trial in kep.iofiles.passFileToTrialList(FileName):
        tr = kep.pipelinepars.keptrial(trial)
        kw = kep.keplc.kw(**tr.kw)
        X = kep.keplc.keplc(tr.kid)
        X.runPipeline(kw)
        Y = kep.qats.qatslc(X.lcFinal,X.KID,maske=kw.maske)
        Y.padLC()
        Y.addNoise()
        Y.runQATS(f=0.01)
        
        coeff = num.polyfit(num.log10(Y.periods),num.log10(Y.snrLC),1)
        outy = scipy.polyval(coeff,num.log10(Y.periods))
        normalizedPower = 10**(outy)
        
        plt.plot(Y.periods,Y.snrFLAT,'r-') 
        plt.plot(Y.periods,Y.snrLC,'b-')
        plt.plot(Y.periods,normalizedPower,'k-')
        plt.setp(plt.gca().set_xscale('log'))
        plt.savefig('sn.'+X.KID+'.png')
        
        dfile = open('signal.'+X.KID+'.data','w')
        print >> dfile,'#',X.KID,'|', Y.periods[num.argmax(Y.snrLC/normalizedPower)],\
                Y.periods[num.argmax(Y.snrLC)], max(Y.snrLC/normalizedPower), max(Y.snrLC)
        for i in range(len(Y.SignalPower)):
            print >> dfile, Y.periods[i],'|',Y.snrLC[i],'|',\
                            Y.snrFLAT[i],'|',normalizedPower[i],'|',\
                            Y.tmin[i],'|',Y.tmax[i],'|',Y.q[i]
        dfile.close()

def linearRegression(histogram):
    x=[]
    y=[]
    for i in range(256):
        for k in range(256):
            if histogram.GetBinContent(i,k)!=0:
                x.append(i)
                y.append(k)
    (ar,br)=polyfit(x,y,1)
    xr=polyval([ar,br],x)
    err=sqrt(sum((xr-y)**2)/len(y))
    if abs(ar)>1.0:
        (ar,br)=polyfit(y,x,1)
        xr=polyval([ar,br],y)
        err=sqrt(sum((xr-x)**2)/len(y))
    return err

def session_fits(session_id, fit_length=50,
        output_filename='results/check_fit.dat',
        write=True):
    dbs = database.DBSession()
    session = dbs.query(database.Session).get(session_id)
    runs = session.experiments[0].runs

    data = []
    for run in runs:
        parameter = run.parameters['release_rate']
        fitness = run.get_objective('pi_fit')
        if fitness is not None:
            data.append((parameter, fitness))

    fit_sorted_data = sorted(data, key=operator.itemgetter(1))
    px, py = zip(*fit_sorted_data[:fit_length])
    coeffs, R, n, svs, rcond = scipy.polyfit(px, py, 2, full=True)

    x, y = zip(*sorted(data))

    peak = float(-coeffs[1] / (2 * coeffs[0]))
    fit = float(R/fit_length)

    if write:
        p = scipy.polyval(coeffs, x)
        header = '# Parabola peak at %s\n# Parabola R^2/n = %s\n'
        rows = zip(x, y, p)
        _small_writer(output_filename, rows, ['release_rate', 'pi_fit', 'parabola'],
                header=header % (peak, fit))

    return session.parameters.get('release_cooperativity'), peak, fit

def plot_rate_vs_mass(array_vs,ks):
  #array_v=array_vs
  #k_growth = ks
  colors = ['#800000','#c87942']#,'#008080','#82bcd3']
  #col='#82bcd3'
  f1 = p.figure()
  for array_v, k_growth,col in zip(array_vs,ks,colors):
    zahl = max([len(i) for i in array_v['R'].flat if i]) # longest list of X values
    nr_cells=array_v.size
    tmp2=zeros((nr_cells,zahl),dtype=float)
    tmp3=zeros((nr_cells,zahl),dtype=float)

    for k1, vector in enumerate(array_v.flat):
      if vector['R']:
	  for k2, r, am, ad, v in zip(range(zahl-len(vector['R']),zahl),vector['R'],vector['B_Am'],vector['B_Ad'],vector['V']):
	    if vector['S_entry'] != 0:
	      tmp2[k1,k2] = log2(r + am + ad) # cell mass
	      tmp3[k1,k2] = log2(r * (am + ad) * k_growth / v) # metabolic rate
    
    tmp2 = array([i for i in tmp2 if i[-1]!=0])
    tmp3 = array([i for i in tmp3 if i[-1]!=0])
    p.plot(tmp2[:,-1],tmp3[:,-1],'o',color=col)
    
    (a_s,b_s,r,tt,stderr)=stats.linregress(tmp2[:,-1],tmp3[:,-1])
    p.text(min(tmp2[:,-1])+1,min(tmp3[:,-1]),'slope = %s'%(round(a_s,2)))
    p.plot(linspace(min(tmp2[:,-1]),max(tmp2[:,-1]),2),polyval([a_s,b_s],linspace(min(tmp2[:,-1]),max(tmp2[:,-1]),2)),'-',color='#25597c',lw=3)

  p.xlabel('Cell mass')
  p.ylabel('Metabolic rate')
  f1.savefig('rate_vs_mass_0.02_and_0.03_SG2_68_G2_cells.pdf')

	def raw_2_deg(self,calib,qei_on,qei_period,qei_rollover, qs_to_qj = 1.0):
		#qei: raw sensor value
		#calib: yaml config calibration map
			
		if self.ctype=="vertx_14bit":
			val = qs_to_qj*360.0*(qei_on)/VERTX_14BIT_MAX
			return val*calib['cb_scale']+calib['cb_bias']
		
		if self.ctype=="ma3_12bit":
			if qei_period==0.0:
				return 0.0
			val = ((qei_on*MA3_12BIT_PERIOD_US)/qei_period)-1 
			val= qs_to_qj*360.0*(val/MA3_12BIT_MAX_PULSE_US)
			return val*calib['cb_scale']+calib['cb_bias']
		
		if self.ctype=="ma3_10bit":
			if qei_period==0.0:
				return 0.0
			val = ((qei_on*MA3_10BIT_PERIOD_US)/qei_period)-1
			val= qs_to_qj*360.0*(val/MA3_10BIT_MAX_PULSE_US)
			return val*calib['cb_scale']+calib['cb_bias']
		
		if self.ctype=="ma3_12bit_poly":
			if qei_period==0.0:
				return 0.0
			val = qs_to_qj*(((qei_on*MA3_12BIT_PERIOD_US)/qei_period)-1)
			val=float(polyval(calib['cb_theta'],val))
			return val*calib['cb_scale']+calib['cb_bias']