def main(data_files1, data_files2,
         err, B=0.5003991, l=2e-2, w=1e-2,
         outfile=None):
    mu_H = lambda V_5, V_1: V_5 * l / ( V_1 * B * w)
    mu_He = lambda V_5, V_1, V_5e, V_1e: pylab.sqrt(\
                ( (l / ( V_1 * B * w)) * V_5e )**2 +\
                ( (V_5 / ( V_1 * B * w)) * 1e-4 )**2 +\
                ( (V_5 * l / ( V_1 * B * w)) * 1e-4 )**2 +\
                ( (V_5 * l / ( V_1**2 * B * w)) * V_1e )**2 +\
                ( (V_5 * l / ( V_1 * B**2 * w)) * 2e-8 )**2 +\
                ( (V_5 * l / ( V_1 * B * w**2)) * 1e-4 )**2)
    x5, x1, V5, V1, N5, N1 = [], [], [], [], [], []
    for df in data_files1:
        databox = spinmob.data.load(df)
        x, V, N = interpolate(databox)
        x5 += x
        V5 += V
        N5 += N
    for df in data_files2:
        databox = spinmob.data.load(df)
        x, V, N = interpolate(databox)
        x1 += x
        V1 += V
        N1 += N
    min_len = min([len(x5), len(x1)])
    xs = pylab.array(x5[:min_len])
    V5, V1 = pylab.array(V5[:min_len]), pylab.array(V1[:min_len])
    N5, N1 = pylab.array(N5[:min_len]), pylab.array(N1[:min_len])
    e5, e1 = err / pylab.sqrt(N5), err / pylab.sqrt(N1)
    ys, es = mu_H(V5, V1), mu_He(V5, V1, e5, e1)
    make_fig(xs, ys, es, outfile)

def fitting(x,y,f,p0,c0,yerr='none',qval='none'):
   if yerr=='none':
      yerr=y*0+y/y
   global func
   func=f
   npar=len(p0)
   func=def_func(npar,c0)
   print 'fitting with funcion: ', func
   print 'no of parameters: ', len(p0)
#   plsq = leastsq(residuals, p0, args=(y,x,yerr), col_deriv=1, maxfev=20000)
   if 'duri' in func:
      plsq= leastsq(residuals_duri_global, p0, args=(y,x,yerr,qval), col_deriv=0, ftol=1e-4, maxfev=2000000)
   else:
      plsq = leastsq(residuals, p0, args=(y,x,yerr), col_deriv=0, maxfev=20000)
   if npar==1:
      final_par = array([plsq[0]])
   else:
      final_par = plsq[0]
   if 'duri' in func:
      yfit=0*y
      for i,q in enumerate(qval):
         q=float(q)
         yfit[i,:]=pyl.array(peval_duri_global(x,final_par,q),typecode='d')
   else:
      yfit=pyl.array(peval(x,final_par),typecode='d')
   return yfit, final_par, func

 def updateParticules(self, X, Y, T, freq=10):
     """ rajouter une ligne dans le groupe de particules"""
     l0 = self.getObjFromType("particules")
     lignes = l0.object
     d = l0.getData()
     if d == None:
         data = []
     else:
         data = d * 1
     t = self.getObjFromType("partTime")
     if t != None:
         txt = t.object * 1
     else:
         txt = []
     p, = pl.plot(pl.array(X), pl.array(Y), "r")
     lignes.append(p)
     obj = GraphicObject("particules", lignes, True, None)
     data.append([X, Y, T])
     obj.setData(data)
     self.addGraphicObject(obj)
     if freq > 0:
         tx, ty, tt = X[0::freq], Y[0::freq], T[0::freq]
         for i in range(len(tx)):
             a = str(tt[i])
             b = a.split(".")
             ln = max(4, len(b[0]))
             txt.append(pl.text(tx[i], ty[i], a[:ln], fontsize="8"))
         obj = GraphicObject("partTime", txt, False, None)
         self.addGraphicObject(obj)
     self.gui_repaint()  # bug matplotlib v2.6 for direct draw!!!
     self.draw()

 def __init__(self, bins = 10, rnge = (0,100), temp_file = 'temp.txt',
                    do_histograms = True, do_coincidences = False):
     """ Thread constructor
     
     :arg temp_file: Name of the temporary file
     :arg do_histogram:
     :arg bins:
     :arg rnge:
     :arg do_coincidences: Generate an additional 
                           histogram with the coincidences t2-t1
     """
     threading.Thread.__init__(self)
     self.bins = bins
     self.rnge = rnge
     self.data = [pl.array([]),pl.array([])]
     self.hist = [generateHist([], self.bins, rnge= self.rnge),
                  generateHist([], self.bins, rnge= self.rnge)]
     self.n_starts = 0
     if temp_file != None:
         self.temp_file = open(temp_file,'w')
     else: 
         self.temp_file = None
     self.do_coincidences = do_coincidences
     self.do_histograms = do_histograms
     if do_coincidences:
         self.rnge_c = (self.rnge[0] - self.rnge[1], self.rnge[1] - self.rnge[0])
         self.hist_c = generateHist([], 2*self.bins, rnge= self.rnge_c)

def plotstuff(cell, electrode):

    figure = mlab.figure(size=(800,800))
    
    l_list = []
    for sec in neuron.h.allsec():
        idx = cell.get_idx_section(sec.name())
        j = 0
        for seg in sec:
            i = idx[j]
            x = pl.array([cell.xstart[i],cell.xend[i]])
            y = pl.array([cell.ystart[i],cell.yend[i]])
            z = pl.array([cell.zstart[i],cell.zend[i]])
            s = pl.array([seg.v, seg.v])
            
            l = mlab.plot3d(x, y, z, s, colormap = 'Spectral',
                            tube_radius=cell.diam[i],
                            representation='surface', vmin=-70, vmax=10)
            l_list.append(l)
            print j
            j += 1
    
    t0 = time()
    ipdb.set_trace()    
    #ms = l_list[0].mlab_source
    while time()-t0 < 10:
        for l in l_list:
            ms = l.mlab_source
            s = pl.rand()*80 -70
            scalars = pl.array([s, s])
            ms.set(scalars = scalars)

def histogram(arguments):

    data = list(map(float, sys.stdin.readlines()))
    data_min = min(data)
    data_avg = pylab.average(pylab.array(data))
    data_max = max(data)
    data_std = pylab.std(pylab.array(data))

    data = filter(
        lambda n: data_avg + arguments.n * data_std > (n**2)**0.5, data)

    pyplot.hist(list(data), bins=arguments.bins)
    pyplot.suptitle(arguments.suptitle)

    if arguments.title is None:
        pyplot.title('min|avg|max|std = {0:0.2f}|{1:0.2f}|{2:0.2f}|{3:0.2f}'
            .format(data_min, data_avg, data_max, data_std))
    else:
        pyplot.title(arguments.title)

    pyplot.xlabel(arguments.xlabel)
    pyplot.ylabel(arguments.ylabel)
    pyplot.grid()

    pyplot.savefig(path(arguments))

def getAntLocations(configuration='D'):
    """Return location information for each antenna in array.

    Arguments:
        Takes configuration as argument (choose from 'A', 'C', 'D').
    Returns:
        Tuple of arrays containing antenna diameters (m), names,
        and x, y, and z locations (m).
    """

    # Raw antenna locations in C configuration
    vx = [41.1100006,  134.110001,   268.309998,  439.410004,  644.210022,
          880.309998,  1147.10999,  1442.41003,  1765.41003,  -36.7900009,
          -121.690002,  -244.789993, -401.190002, -588.48999,  -804.690002,
          -1048.48999, -1318.48999, -1613.98999,  -4.38999987,-11.29,  
          -22.7900009, -37.6899986, -55.3899994, -75.8899994, -99.0899963, 
          -124.690002, -152.690002]
    vy = [3.51999998, -39.8300018,  -102.480003, -182.149994, -277.589996,
         -387.839996, -512.119995, -649.76001,  -800.450012, -2.58999991,
         -59.9099998,  -142.889999, -248.410004, -374.690002, -520.599976,
         -685,        -867.099976, -1066.42004,   77.1500015,  156.910004, 
         287.980011,  457.429993,  660.409973,  894.700012,  1158.82996, 
         1451.43005,  1771.48999]
    vz = [0.25,       -0.439999998, -1.46000004, -3.77999997, -5.9000001,
         -7.28999996, -8.48999977, -10.5,       -9.56000042,      0.25, 
         -0.699999988, -1.79999995, -3.28999996, -4.78999996, -6.48999977,
         -9.17000008, -12.5299997, -15.3699999,  1.25999999,   2.42000008, 
         4.23000002,  6.65999985,  9.5,         12.7700005,  16.6800003, 
         21.2299995,  26.3299999]
    d = [25.0,       25.0,         25.0,         25.0,       25.0,
         25.0,       25.0,         25.0,         25.0,       25.0,
         25.0,       25.0,         25.0,         25.0,       25.0,
         25.0,       25.0,         25.0,         25.0,       25.0,
         25.0,       25.0,         25.0,         25.0,       25.0,
         25.0,       25.0]

    # Get scaling factor for antenna locations based on configuration
    if(configuration=='D'):
        scale = 3.0
    else:
        if(configuration=='A'):
            scale = 1.0/9.0
        else:
            scale = 1.0
        print 'Using VLA C-array coords'

    # Move antennas into desired configuration
    nn = len(vx)
    x = (vx - (sum(pl.array(vx))/(nn)))/scale
    y = (vy - (sum(pl.array(vy))/(nn)))/scale
    z = (vz - (sum(pl.array(vz))/(nn)))/scale

    # Label the antenna
    an=[]
    for i in range(0,nn):
        an.append("VLA"+str(i))
    return d,an,x,y,z

    def __init__(self, contact_area_percent=50.0):

        ######################################
        # begin: parameters to be specified

        self.contact_area_percent = contact_area_percent

        # resistor that is in series with the taxel (Ohms)
        self.r1 = 47.0 

        # total voltage across the taxel and r1, which are in serise (Volts)
        self.vtot = 5.0 

        # the maximum resistance of the taxel when no pressure is applied (Ohms)
        self.rtax_max = 50.0 

        # the minimum force that will be applied to the taxel (Newtons)
        self.fz_min = 0.0 

        # the maximum force that will be applied to the taxel (Newtons)
        self.fz_max = 45.0 

        # the number of bits for the analog to digital conversion 
        self.adc_bits = 10 

        # the pressure sensitive area of the taxel (meters^2)
        self.taxel_area = 0.04 * 0.04 

        # pressure that results in minimum resistance after which
        # further pressure does not result in a reduction in the
        # signal, since the sensor is saturated (Pascals = N/m^2)
        self.pressure_max = self.fz_max/(0.4 * self.taxel_area) 

        # hack to specify the minimum resistance of the taxel, which
        # is associated with the maximum pressure. for now, it's
        # specified as a percentage of the maximum resistance, which
        # is associated with 0 applied pressure (no contact)
        self.r_min_percent_of_r_no_contact = 0.001 #

        # end
        ######################################

        self.r_no_contact = self.taxel_area * self.rtax_max 
        self.r_min = self.r_no_contact * (self.r_min_percent_of_r_no_contact/100.0)
        self.fz_array = pl.arange(self.fz_min, self.fz_max, 0.001) # N
        self.adc_range = pow(2.0, self.adc_bits)
        self.volts_per_adc_unit = self.vtot/self.adc_range # V 
        self.contact_area = self.taxel_area * (self.contact_area_percent/100.0) # m^2
        self.no_contact_area = self.taxel_area - self.contact_area # m^2
        self.pressure_array = pl.array([f/self.contact_area for f in self.fz_array]) # Pascals = N/m^2
        self.rtax_array = pl.array([self.rtax(f) for f in self.pressure_array])
        self.vdigi_array = pl.array([self.output_voltage(r) for r in self.rtax_array])
        self.vdigi_max = self.output_voltage(self.rtax_max)
        self.adc_bias = self.vdigi_max/self.volts_per_adc_unit
        self.adc_array = self.vdigi_array/self.volts_per_adc_unit
        self.adc_plot = self.adc_bias - self.adc_array

def makeMSFrame(dirname,msname,ra0,dec0,nchan):
  msstokes='RR LL';
  feedtype='perfect R L';


  ## Directory for the MS
  if(not os.path.exists(dirname)):
    cmd = 'mkdir ' + dirname;
    os.system(cmd);

  vx = [41.1100006,  -34.110001,  -268.309998,  439.410004,  -444.210022]
  vy = [3.51999998, 129.8300018,  +102.480003, -182.149994, -277.589996]
  vz = [0.25,       -0.439999998, -1.46000004, -3.77999997, -5.9000001]
  d = [25.0,       25.0,         25.0,         25.0,       25.0]
  an = ['VLA1','VLA2','VLA3','VLA4','VLA5'];
  nn = len(vx)*2.0;
  x = 0.5*(vx - (sum(pl.array(vx))/(nn)));
  y = 0.5*(vy - (sum(pl.array(vy))/(nn)));
  z = 0.5*(vz - (sum(pl.array(vz))/(nn)));

  ####  This call will get locations for all 27 vla antennas.
  #d, an, x, y, z = getAntLocations()


  obspos = me.observatory('EVLA');
  #obspos = me.position('ITRF', '-0.0m', '0.0m', '3553971.510m');

  ## Make MS Frame
  sm.open(ms=msname);
  sm.setconfig(telescopename='EVLA',x=x.tolist(),y=y.tolist(),z=z.tolist(),dishdiameter=d,
               mount=['alt-az'], antname=an,
               coordsystem='local',referencelocation=obspos);
  sm.setspwindow(spwname="CBand",freq="6.0GHz",deltafreq='500MHz',
                 freqresolution='2MHz',nchannels=nchan,stokes=msstokes);
  sm.setfeed(mode=feedtype,pol=['']);
  sm.setfield( sourcename="fake",sourcedirection=me.direction(rf='J2000',v0=ra0,v1=dec0) );
  sm.setlimits(shadowlimit=0.01, elevationlimit='10deg');
  sm.setauto(autocorrwt=0.0);
  sm.settimes(integrationtime='1800s', usehourangle=True,
                       referencetime=me.epoch('UTC','2013/05/10/00:00:00'));
  # Every 30 minutes, from -3h to +3h
  ostep = 0.5
  for loop in pl.arange(-3.0,+3.0,ostep):
    starttime = loop
    stoptime = starttime + ostep
    print starttime, stoptime
    for ch in range(0,nchan):
        sm.observe(sourcename="fake",spwname='CBand',
                   starttime=str(starttime)+'h', stoptime=str(stoptime)+'h');
  sm.close();

  listobs(vis=msname)

  return d

 def get_ind_under_point(self, event):
     "get the index of the vertex under point if within epsilon tolerance"
     x, y = self.lx, self.ly
     d = sqrt((pl.array(x) - event.xdata) ** 2 + (pl.array(y) - event.ydata) ** 2)
     indseq = nonzero(equal(d, amin(d)))
     ind = indseq[0]
     if len(ind) > 1:
         ind = ind[0]
     if d[ind] >= self.epsilon:
         ind = None
     return ind

def fit_power(Ts, Rs, Rerr, pguess):
    'Fit a power function to the data'
    Rs, Rerr = pylab.array(Rs), pylab.array(Rerr)
    model, ps = 'a * (x-x0)**(3/2)', 'a,x0'
    # Make intelligent guesses for the parameters
    a, x0 = pguess
    # Create a spinmob fitter
    fitter = spinmob.data.fitter(model, ps)
    fitter.set_data(Ts, Rs, Rerr)
    fitter.set(a=a, x0=x0) # Set guesses
    # Fit.
    fitter.fit()
    return fitter

def fit_exp(Ts, Rs, Rerr, eguess):
    'Fit an exponential function to the data'
    Rs, Rerr = pylab.array(Rs), pylab.array(Rerr)
    model, ps = 'a * exp(b / x)', 'a,b'
    # Make intelligent guesses for the parameters
    a, b = eguess
    b = pylab.log(Rs[0]/Rs[1]) / ( Rs[-1] - Rs[0] )
    # Create a spinmob fitter
    fitter = spinmob.data.fitter(model, ps)
    fitter.set_data(Ts, Rs, Rerr)
    fitter.set(a=a, b=b) # Set guesses
    # Fit.
    fitter.fit()
    return fitter

def buildCityPoints(fName, scaling):
    cityNames, featureList = readCityData(fName, scaling)
    points = []
    for i in range(len(cityNames)):
        point = City(cityNames[i], pylab.array(featureList[i]))
        points.append(point)
    return points

def splitfit(Ts, Rs, es, a, b, c, d, pguess, eguess, outfile=None):
    ## Split the data in two parts
    x1, x2, y1, y2, e1, e2 = [], [], [], [], [], []
    for T, R, pe, ee in zip(Ts, Rs, es[0], es[1]):
        if a < T < b:
            x1.append(T)
            y1.append(abs(R))
            e1.append(pe)
        elif c < T < d:
            x2.append(T)
            y2.append(abs(R))
            e2.append(ee)
    ## Fit one part with the exponential
    fit1 = fit_power(x1, y1, e1, pguess)
    ## Fit one part with the polynomial
    fit2 = fit_exp(x2, y2, e2, eguess)
    res1 = fit1.results[0]
    res2 = fit2.results[0]
    fct1 = lambda x: res1[0] * (x - res1[1])
    fct2 = lambda x: res2[0] * pylab.exp(res2[1]/x)
    Rs = pylab.array(Rs)
    Ges = [[[e/R**2] for R, e in zip(Rs, es[i])] for i in range(2)]
    pylab.clf()
    make_fig(Ts, Rs, Ges, a, b, c, d, fct1, fct2, outfile)
    return fit1, fit2

def RealtimePloter(arg):
	global values
	#current x axis (cxa) is the length of values from max length to length-100
	cxa = range(len(values)-100,len(values),1)
	Stock1[0].set_data(cxa,pylab.array(values[-100:]))
	ax.axis([min(cxa),max(cxa),high,low])
	manager.canvas.draw()

def scaleFeatures(vals):
    """Assumes vals is a sequence of numbers"""
    result = pylab.array(vals)
    mean = sum(result) / float(len(result))
    result -= mean
    sd = stdDev(result)
    result /= sd
    return result

def testfunction(data):
    # N-D Gaussian or N-D Runge Function
    N, sd = data.shape
    f = ones((N,1))
    for i in range(sd):
        #f = f*array([exp(-15*(data[:,i]-0.5)**2)]).T
        f = f*array([1./(1+(5*data[:,i])**2)]).T
        
    return f

def pruning_plot(costs, figname, max_cost=None):
    if max_cost == None:
        max_cost = max(costs)
    assert max_cost in costs
    costs = PP.array(costs)
    costs /= float(max_cost)
    assert 1 in costs
    PP.plot(range(len(costs)), costs)
    PP.xlabel('time steps')
    PP.ylabel('proportion of edges in use')
    PP.savefig(figname + '.png', format='png')

def main(data_files, a, b, c, d, pguess, eguess, perr=1, eerr=1, outfile=None):
    xs, ys, es, Ns = [], [], [], []
    for df in data_files:
        databox = spinmob.data.load(df)
        x, y, e, N = databox[:4]
        xs += list(x)
        ys += [abs(i) for i in y]
        es += list(e)
        Ns += list(N)
    xs, ys, es, Ns = pylab.array(xs), pylab.array(ys), pylab.array(es), pylab.array(Ns)
    pes, ees = perr / pylab.sqrt(Ns), eerr / pylab.sqrt(Ns)
    if 'current' in databox.hkeys:
        current = databox.h('current')
    else:
        current = 0.001 # A
        databox.h(current=current)
    Rs = ys / current
    Rerrs = pes / current, ees / current
    fits = splitfit(xs, Rs, Rerrs, a, b, c, d, pguess, eguess, outfile)
    return fits

def main(data_files, a, b, c, d, pguess, eguess, perr=1, eerr=1,
         I=0.001, B=0.5003991, sample_thickness=1e-3,
         outfile=None):
    R_H = lambda V_H: V_H*sample_thickness / (I*B) # Vm/AT
    R_He = lambda V_H, V_He: pylab.sqrt(\
                ( sample_thickness/(I*B) * V_He )**2 + \
                ( V_H / (I*B) * 1e-4 )**2 + \
                ( V_H * sample_thickness / (B**2 * I) * 2e-8 )**2 )
    xs, ys, es, Ns = [], [], [], []
    for df in data_files:
        databox = spinmob.data.load(df)
        x, y, e, N = databox[:4]
        xs += list(x)
        ys += [R_H(i) for i in y]
        es += [R_He(i, j) for i, j in zip(y, e)]
        Ns += list(N)
    xs, ys, es, Ns = pylab.array(xs), pylab.array(ys), pylab.array(es), pylab.array(Ns)
    pes, ees = perr / pylab.sqrt(Ns), eerr / pylab.sqrt(Ns)
    fits = splitfit(xs, ys, (pes, ees), a, b, c, d, pguess, eguess, outfile)
    return fits