def calculateFDunc(self):
        #Calculates the uncertainty of the FFT according to:
        #   - J. M. Fornies-Marquina, J. Letosa, M. Garcia-Garcia, J. M. Artacho, "Error Propagation for the transformation of time domain into frequency domain", IEEE Trans. Magn, Vol. 33, No. 2, March 1997, pp. 1456-1459
        #return asarray _tdData
        #Assumes tha the amplitude of each time sample is statistically independent from the amplitude of the other time
        #samples

        # Calculates uncertainty of the real and imaginary part of the FFT and ther covariance
        unc_E_real = []
        unc_E_imag = []
        cov = []
        for f in self.getfreqs():
            unc_E_real.append(py.sum((py.cos(2*py.pi*f*self._tdData.getTimes())*self._tdData.getUncEX())**2))
            unc_E_imag.append(py.sum((py.sin(2*py.pi*f*self._tdData.getTimes())*self._tdData.getUncEX())**2))
            cov.append(-0.5*sum(py.sin(4*py.pi*f*self._tdData.getTimes())*self._tdData.getUncEX()**2))
        
        unc_E_real = py.sqrt(py.asarray(unc_E_real))
        unc_E_imag = py.sqrt(py.asarray(unc_E_imag))
        cov = py.asarray(cov)
        
        # Calculates the uncertainty of the modulus and phase of the FFT
        unc_E_abs = py.sqrt((self.getFReal()**2*unc_E_real**2+self.getFImag()**2*unc_E_imag**2+2*self.getFReal()*self.getFImag()*cov)/self.getFAbs()**2)
        unc_E_ph = py.sqrt((self.getFImag()**2*unc_E_real**2+self.getFReal()**2*unc_E_imag**2-2*self.getFReal()*self.getFImag()*cov)/self.getFAbs()**4)
        
        t=py.column_stack((self.getfreqs(),unc_E_real,unc_E_imag,unc_E_abs,unc_E_ph))
        return self.getcroppedData(t)  

def projectCl(lvec,P,D,dNdz,z,growthFac=None):
    """
    project C_l's given a power spectrum class P (Camb or
    BBKS) and Distance class D together

    arguments:
    lvec: vector of l values
    P: p.pk p.k contains the power spectrum, e.g. pt.Camb instance
    D: frw.Distance instance
    dNdz,z, growthFac: vectors suitable for trapezoid z-integration

    presently it crashes if z=0.0 is included, start from a small z value
    """
    lvec = M.asarray(lvec)
    dNdz2 = M.asarray(dNdz)**2
    z = M.asarray(z)
    da1 = 1./D.rtc(z)/D.h #comoving Da in h^1Mpc

    dNdz2vc = dNdz2/D.vc(z)/D.h**3 # comovin volume in (h^-1Mpc)^3
    #`use growth factor if given
    if growthFac:
        dNdz2vc = dNdz2vc*(growthFac**2)
    lk = M.log(P.k)
    pk = P.pk
    
##     return M.asarray([utils.trapz(utils.splineResample(pk,lk,
##                      M.log(l*da1))*dNdz2vc,z) for l in lvec])
    return M.asarray([utils.trapz(utils.interpolateLin(pk,lk,
                     M.log(l*da1))*dNdz2vc,z) for l in lvec])

def show_first_and_last_frames(vid,f1=1,f2=2):
	vidlen = vidtools.vid_duration(vid)
	out = os.path.splitext(vid)[0]
	os.system('vid2png.py %s %s 0 1 1' % (vid,out))
	os.system('vid2png.py %s %s %s 1 1' % (vid,out,vidlen-1))
	m1 = pylab.asarray(Image.open(sorted(glob(out+'/*.png'))[0]).convert('L'))
	m2 = pylab.asarray(Image.open(sorted(glob(out+'/*.png'))[-1]).convert('L'))
	pylab.matshow(m1,fignum=f1)
	pylab.matshow(m2,fignum=f2)
	return(m1.mean()-m2.mean())

def plot_color_vs_mass_vs_icd():
    galaxies=mk_galaxy_struc()
    # Add the figures

    # Mass vs color plot I-H
    f1 = pyl.figure(1,figsize=(6,4))
    f1s1 = f1.add_subplot(212)


    color1 = []
    mass1 = []
    icd1 = []

    for i in range(len(galaxies)):
        if galaxies[i].ston_I >30.0:
            if -0.05 < galaxies[i].ICD_IH and galaxies[i].ICD_IH < 0.25:
                mass1.append(galaxies[i].Mass)
                color1.append(galaxies[i].Imag-galaxies[i].Hmag)
                icd1.append(galaxies[i].ICD_IH)
            else:
                mass1.append(galaxies[i].Mass)
                color1.append(galaxies[i].Imag-galaxies[i].Hmag)
                icd1.append(0.25)

    # Sort the arrays by ICD
    mass1 = pyl.asarray(mass1)
    color1 = pyl.asarray(color1)
    icd1 = pyl.asarray(icd1)

    IH_array = pyl.column_stack((mass1,color1,icd1))

    IH_array = colsort(IH_array,3)

    sc1 = f1s1.scatter(IH_array[:,0], IH_array[:,1], c=IH_array[:,2], s=50,
    cmap='spectral')


    ############
    # FIGURE 1 #
    ############

    bar = pyl.colorbar(sc1)
    bar.set_label(r"$\xi[I,H]$")

    f1s1.set_xscale('log')
    f1s1.set_xlim(3e7,1e12)
    f1s1.set_ylim(0.0,4.0)
    f1s1.set_xlabel(r"Mass $[M_{\odot}]$")
    f1s1.set_ylabel("$(I-H)_{Observed}$")

#    pyl.subplots_adjust(left=0.15, bottom=0.15, right=.75)
#    pyl.savefig('color_vs_mass_vs_icd_IH.eps',bbox='tight')
    return f1s1

    def __init__( self, t, x, y, A=[1., 0.85], a=[0.25, 0.85] ):
        '''
        
        Initializing generalized thalamo-cortical loop. Full
        functionality is only obtained in the subclasses, like DOG,
        full_eDOG, etc.
        
        Parameters
        ----------
        
        t : array
            1D Time vector
        x : np.array
            1D vector for x-axis sampling points
        y : np.array
            1D vector for y-axis sampling points

        Keyword arguments
        -----------------
        
        A : sequence (default A = [1., 0.85])
            Amplitudes for DOG receptive field for center and surround, 
            respectively
        a : sequence (default a = [0.25, 0.85])
            Width parameter for DOG receptive field for center and surround,
            respectively

        Usage
        -----
            Look at subclasses for example usage            
        '''
        # Set parameteres as attributes
        self.name = 'pyDOG Toolbox'
        self.t = t
        self.A = A
        self.a = a
        self.x = x
        self.y = y
        # Find sampling rates and sampling freqs and 
        self.nu_xs = 1./(x[1]-x[0])
        self.nu_ys = 1./(y[1]-y[0])
        self.fs = 1./(t[1]-t[0])
        self.f = fft.fftfreq(pl.asarray(t).size, t[1]-t[0])
        # fftshift spatial frequency,
        self.nu_x = fft.fftfreq(pl.asarray(x).size, x[1]-x[0])
        self.nu_y = fft.fftfreq(pl.asarray(y).size, y[1]-y[0])
        # Make meshgrids, may come in handy
        self._xx, self._yy = pl.meshgrid(self.x, self.y)
        self._nu_xx, self._nu_yy = pl.meshgrid(self.nu_x, self.nu_y)
        # r is needed for all circular rfs
        self.r = pl.sqrt(self._xx**2 + self._yy**2)
        self.k = 2 * pl.pi * pl.sqrt(self._nu_xx**2 + self._nu_yy**2)

def learn(rate,desire,inp,weight):
    """
    desire ... one number;
    input ... one input vector
    wini ... initial weights
    """
    inp = asarray(inp)
    weight = asarray(weight)
    
    if neuron(weight,inp) == desire:
        dw = zeros(len(inp))
    else:
        dw = rate* desire * inp
    return dw

 def _bringToCommonTimeAxis(self,tdDatas):
     #What can happen: 
     #a) datalengthes are not equal due to missing datapoints 
     #   => no equal frequency bins
     #b) time positions might not be equal
         
     miss_points_max=10    
     #check for missing datapoints, allowing 
     #not more than miss_points_max points to miss 
     all_lengthes=[]
     for thisdata in tdDatas:
         all_lengthes.append(len(thisdata[:,0]))
     
     #do it always, just print a warning, if miss_points_max is exceeded
     if min(all_lengthes)!=max(all_lengthes):
         print("Datalength of suceeding measurements not consistent, try to fix")
         if max(all_lengthes)-min(all_lengthes)>miss_points_max:
             print("Warning: Data seems to be corrupted. \n" +\
             "The length of acquired data of repeated measurements differs by \n" + \
                 str(max(all_lengthes)-min(all_lengthes)) + ' datapoints')
     
     #interpolation does no harm, even if everything is consistent (no interpolation in this case)
     commonMIN=max([thistdData[:,0].min() for thistdData in tdDatas])
     commonMAX=min([thistdData[:,0].max() for thistdData in tdDatas])
     commonLENGTH=min([thistdData[:,0].shape[0] for thistdData in tdDatas])
     
     #interpolate the data
     for i in range(self.numberOfDataSets):
         tdDatas[i]=self.getInterData(tdDatas[i],commonLENGTH,commonMIN,commonMAX)
     
     return py.asarray(tdDatas)

    def load_map(self, data, isReversed):
        
        # sort
        data=list(data)
        data.sort( key=lambda x: x[1])
        data.sort( key=lambda x: x[2])
        data=py.asarray(data)
        
        # find nfocus
        firstfocus=data[0][0]
        for row in data[1:]:
            self._nfocus+=1
            if row[0] == firstfocus:
                break
        
        # extract lum data    
        for row in data:
            if isReversed:
                lum = row[3:][::-1]
            else:
                lum = row[3:]
            self._specList.append(Spectrum(wavelen=self._wavelen,lum=lum))

        # split specList into points, ie. [[z1, z1, z3],[z1,z2,z3]]
        self._specList=[self._specList[i:i+self._nfocus] for i in range(0,len(self._specList),self._nfocus)]
        self._z = (data[:,0][0:self._nfocus])
        
        assert len(self._z) == len(self._specList[0]), "len of focuses must match specList sublist length"        

 def callback(self, verts):
     if len(self.Nxy) > 1:
         for j, series in enumerate(self.y):
             # ind = pl.np.nonzero(points_inside_poly(zip(self.x, series), verts))[0]
             ind = pl.np.nonzero(mpl_path(verts).contains_points(zip(self.x, series)))[0]
             for i in range(self.Nxy[1]):
                 if i in ind:
                     self.badpoints.append([i, j, self.y[j, i]])
                     self.axes.collections[j]._offsets[i,1] = pl.nan
                     self.y[j,i] = pl.nan
     else:
         cleanedpoints = self.axes.collections[0]._paths[0].vertices.tolist()
         # ind = pl.np.nonzero(points_inside_poly(cleanedpoints, verts))[0]
         ind = pl.np.nonzero(mpl_path(verts).contains_points(cleanedpoints))[0]
         removedcount = 0
         for i in range(len(cleanedpoints)):
             if i in ind:
                 self.badpoints.append([i, self.x[i], self.y[i]])
                 out = cleanedpoints.pop(i - removedcount)
                 self.axes.collections[0]._paths[0].vertices = pl.asarray(cleanedpoints)
                 original_indx = pl.find(self.x == out[0])
                 self.y[original_indx] = pl.nan
                 removedcount += 1
     self.canvas.draw_idle()
     self.canvas.widgetlock.release(self.lasso)
     del self.lasso

def mk_grid(llx, ulx, nx, lly, uly, ny):
    # Get the Galaxy info
    #galaxies = mk_galaxy_struc()
    galaxies = pickle.load(open('galaxies.pickle','rb'))
    galaxies = filter(lambda galaxy: galaxy.ston_I > 30., galaxies)
    galaxies = pyl.asarray(filter(lambda galaxy: galaxy.ICD_IH < 0.5, galaxies))

    # Make the low mass grid first
    x = [galaxy.Mass for galaxy in galaxies]
    y = [galaxy.ICD_IH *100 for galaxy in galaxies]
    bins_x =pyl.linspace(llx, ulx, nx)
    bins_y = pyl.linspace(uly, lly, ny)

    grid = []

    for i in range(bins_x.size-1):
        xmin = bins_x[i]
        xmax = bins_x[i+1]
        for j in  range(bins_y.size-1):
            ymax = bins_y[j]
            ymin = bins_y[j+1]

            cond=[cond1 and cond2 and cond3 and cond4 for cond1, cond2, cond3,
                cond4 in zip(x>=xmin, x<xmax, y>=ymin, y<ymax)]

            grid.append(galaxies.compress(cond))

    return grid

def cosmic_rejection(x, y, n):
    '''Try to reject cosmic from a spectrum
    '''
    bla = True
    blabla = False

    if n == 0:
        print " Warning: sigma for rejection = 0. Take 0.01."
        n = 0.1

    msk = abs(y-y.mean()) > n * y.std(ddof=1)

    xrej = x[msk]
    yrej = y[msk]

    if bla:
        print " Rejection of points outside", n, "sigma around the mean."
        print " Number of rejected points:", xrej.size, '/', x.size

    if blabla:
        print " Rejected points:"
        print xrej
        print yrej

    msk = pl.asarray([not i for i in msk])

    return msk, xrej, yrej

def sigma2fromPk(c, r):
    """
    calculate sigma^2 from pk

    this function can be called with vectors or scalars, but always returns a vector
    """
    r = M.asarray(r)
    return 9.0 / r ** 2 * N.trapz(c.k * c.pk * sf.j1(M.outer(r, c.k)) ** 2, M.log(c.k)) / 2.0 / M.pi ** 2

def xi2fromCambPk(c, r):
    """
    calculate 2pt corr. function from Camb instance (with its
    associated k,pk)

    this function can be called with vectors
    """
    r = M.asarray(r)
    return N.trapz(c.k ** 3 * c.pk * sf.j0(M.outer(r, c.k)), M.log(c.k)) / 2.0 / M.pi ** 2

 def matrix_show(arr, nx):
     try: # import MatPlotLib if available
         from pylab import imshow, show, cm, asarray
     except ImportError:
         return "MatPlotLib library (http://matplotlib.sourceforge.net) cannot be imported."
     else:
         matrix = [arr[i:i+nx] for i in xrange(0, len(arr), nx)][::-1]
         matrix = asarray(matrix, dtype="UInt8")
         imshow(matrix, cmap=cm.gray, interpolation="nearest")
         show()
         return "MatPlotLib successiful structure plot show."

def trapz(y, x=None, ax=-1, method=pl.add.reduce):
    """trapz(y,x=None,ax=-1) integrates y along the given dimension of
    the data array using the trapezoidal rule.

    you can call it with method=pl.add.accumulate to yield partial sums
    """
    y = pl.asarray(y)
    if x is None:
        d = 1.0
    else:
        d = pl.diff(x,axis=ax)
    y = pl.asarray(y)
    nd = len(y.shape)
    s1 = nd*[slice(None)]
    s2 = nd*[slice(None)]
    s1[ax] = slice(1,None)
    s2[ax] = slice(None,-1)

    ans = method(0.5* d * (y[s1]+y[s2]),ax)
    return ans

def astausgleich(ab2org, mn2org, rhoaorg):
    """shifts the branches of a dc sounding to generate a matching curve."""
    ab2 = P.asarray(ab2org)
    mn2 = P.asarray(mn2org)
    rhoa = P.asarray(rhoaorg)
    um = P.unique(mn2)
    for i in range(len(um) - 1):
        r0, r1 = [], []
        ac = P.intersect1d(ab2[mn2 == um[i]], ab2[mn2 == um[i + 1]])
        for a in ac:
            r0.append(rhoa[(ab2 == a) * (mn2 == um[i])][0])
            r1.append(rhoa[(ab2 == a) * (mn2 == um[i + 1])][0])

        if len(r0) > 0:
            fak = P.mean(P.array(r0) / P.array(r1))
            print(fak)
            if P.isfinite(fak) and fak > 0.:
                rhoa[mn2 == um[i + 1]] *= fak

    return rhoa  # formerly pg as vector

 def calcunc(self,tdDatas):
     #not used anymore, older version, should we remove it???
      #tdDatas is a np array of tdData measurements
     if tdDatas.shape[0]==1:
         repeatability=py.zeros((len(tdDatas[0,:,0]),2))
     else:
         repeatability=py.std(py.asarray(tdDatas[:,:,1:3]),axis=0, ddof = 1)/py.sqrt(self.numberOfDataSets)
     #this line is wrong
     elnNoise=tdDatas[0,:,3:]
     uncarray = py.sqrt(repeatability**2 + elnNoise**2)
     
     return uncarray

def xi2fromPk(k, pk, r):
    """
    calculate 2pt corr. function from k, p(k).
    It doesn't seem to work particularly well, though.

    this function can be called with vectors
    """
    r = M.asarray(r)

    print len(k), len(pk), len(r)

    return N.trapz(k ** 3 * pk * sf.j0(M.outer(r, k)), M.log(k)) / 2.0 / M.pi ** 2

def wPk(c, r, w):
    """
    convolve Pk with an arbitrary window function w(kr)
    
    int k^2 P(k)w(kr) dk/2/pi^2

    e.g., the previous two functions can be realized as
    wPk(c,r,j0)
    wPk(c,r,lambda x: (3*pt.sf.j1(x)/x)**2))
    """
    r = M.asarray(r)
    return N.trapz(c.k ** 3 * c.pk * w(M.outer(r, c.k)), M.log(c.k)) / 2.0 / M.pi ** 2

    def calculaten(self,H,l):
        #this calculates n for a given transferfunction H and a given thickness l
    
        #calculate the initial values
        inits=py.asarray(self.calculateinits(H,l))
        
        res=[] #the array of refractive index
        vals=[] # the array of error function values
        bnds=((1,None),(0,None)) #not used at the moment, with SLSQP bnds can be introduced
        nums=len(H[:,0])
        #minimize the deviation of H_theory to H_measured for each frequency
        for i in range(nums):
#            t=minimize(self.error_func,[inits[0,i],inits[1,i]],args=(H[i,:2],l), method='SLSQP',\
#            bounds=bnds, options={'ftol':1e-9,'maxiter':2000, 'disp': False})
            t=minimize(self.error_func,[inits[0,i],inits[1,i]],args=(H[i,:3],l), method='Nelder-Mead',
            options={'xtol': 1e-6,'disp':False})
            res.append(t.x[0]-1j*t.x[1])
            vals.append(t.fun)
        n=py.asarray(res)
        #self.n is a 5xlengthf array, frequency,n_real,n_imag,n_smoothed_real,n_smoothed_imag
        self.n=py.column_stack((H[:,0],n,n))
        return n