def dist2(x):
  R, ETA = pylab.meshgrid(r[r<fit_rcutoff], eta)
  g = pylab.zeros_like(ETA)
  g = evalg(x, ETA, R)
  gfit = pylab.reshape(g, len(eta)*len(r[r<fit_rcutoff]))
  return gfit - pylab.reshape([g[r<fit_rcutoff] for g in ghs],
                              len(eta)*len(r[r<fit_rcutoff]))

def dist2(x):
  R, GSIGMAS = pylab.meshgrid(r[r<fit_rcutoff], gsigmas)
  g = pylab.zeros_like(GSIGMAS)
  g = evalg(x, GSIGMAS, R)
  gfit = pylab.reshape(g, len(eta)*len(r[r<fit_rcutoff]))
  return gfit - pylab.reshape([g[r<fit_rcutoff] for g in ghs],
                              len(gsigmas)*len(r[r<fit_rcutoff]))

def loadMNISTImages(filename):
  f = open(filename, 'rb')

  # Verify Magic Number
  s = f.read(4)
  magic = int(s.encode('hex'),16)
  assert(magic == 2051)

  # Get Number of Images
  s = f.read(4)
  numImages = int(s.encode('hex'),16)
  s = f.read(4)
  numRows = int(s.encode('hex'),16)
  s = f.read(4)
  numCols = int(s.encode('hex'),16)

  # Get Data
  s = f.read()
  a = frombuffer(s, uint8)

  # Use 'F' to ensure that we read by column
  a = reshape(a, (numCols , numRows, numImages), order='F');
  images = transpose(a, (1, 0, 2))
  f.close()

  # Reshape to #pixels * #examples
  images  = reshape(a, (shape(images)[0] * shape(images)[1], numImages),
          order='F');
  images = double(images)/255
  return images

def getParamCovMat(prefix,dlogpower = 2, theoconstmult = 1.,dlogfilenames = ['dlogpnldloga.dat'],volume=256.**3,startki = 0, endki = 0, veff = [0.]):
    """
    Calculates parameter covariance matrix from the power spectrum covariance matrix and derivative term
    in the prefix directory
    """
    nparams = len(dlogfilenames)

    kpnl = M.load(prefix+'pnl.dat')
    k = kpnl[startki:,0]

    nk = len(k)
    if (endki == 0):
        endki = nk
        
    pnl = M.array(kpnl[startki:,1],M.Float64)
    covarwhole = M.load(prefix+'covar.dat')
    covar = covarwhole[startki:,startki:]
    if len(veff) > 1:
        sqrt_veff = M.sqrt(veff[startki:])
    else:
        sqrt_veff = M.sqrt(volume*M.ones(nk))

    dlogs = M.reshape(M.ones(nparams*nk,M.Float64),(nparams,nk))
    paramFishMat = M.reshape(M.zeros(nparams*nparams*(endki-startki),M.Float64),(nparams,nparams,endki-startki))
    paramCovMat = paramFishMat * 0.

    # Covariance matrices of dlog's
    for param in range(nparams):
        if len(dlogfilenames[param]) > 0:
            dlogs[param,:] = M.load(prefix+dlogfilenames[param])[startki:,1]

    normcovar = M.zeros(M.shape(covar),M.Float64)
    for i in range(nk):
        normcovar[i,:] = covar[i,:]/(pnl*pnl[i])

    M.save(prefix+'normcovar.dat',normcovar)

    f = k[1]/k[0]

    if (volume == -1.):
        volume = (M.pi/k[0])**3

    #theoconst = volume * k[1]**3 * f**(-1.5)/(12.*M.pi**2) #1 not 0 since we're starting at 1
    for ki in range(1,endki-startki):
        for p1 in range(nparams):
            for p2 in range(nparams):
                paramFishMat[p1,p2,ki] = M.sum(M.sum(\
                M.inverse(normcovar[:ki+1,:ki+1]) *
                M.outerproduct(dlogs[p1,:ki+1]*sqrt_veff[:ki+1],\
                               dlogs[p2,:ki+1]*sqrt_veff[:ki+1])))
                
                
        paramCovMat[:,:,ki] = M.inverse(paramFishMat[:,:,ki])

    return k[1:],paramCovMat[:,:,1:]

    def retrieve_result(self, filename):
	"""
	retrieve_many_result(filename)

	Import the result of the dynamics as a class variable.

	Parameters
	----------
	filename : name of the result file.
	"""
	
	if filename[-3:] == "txt":
	    raw_result = numpy.loadtxt("%s/%s"%(self.path_to_output_dir, filename))
	    numpy.save("%s/%s"%(self.path_to_output_dir, filename[:-4]), raw_result)
	else:
	    raw_result = numpy.load("%s/%s"%(self.path_to_output_dir, filename))
	
	self.times = raw_result[:,0]
	self.field = raw_result[:,1]
	
	raw_result = raw_result[:,2:]


	self.psi = zeros([self.vib_basis_size, self.el_basis_size, 
	    raw_result.shape[0]], dtype = complex)
	
	for index in range(raw_result.shape[0]):
	    self.psi[:,:,index] = reshape(raw_result[index, :self.basis_size] + 
		1j * raw_result[index, self.basis_size:], 
		[self.vib_basis_size, self.el_basis_size], order = "F")
	self.psi_final = self.psi[:,:,-1]

def plot_hist(X,Y,title,name):
    # get list of tracks and list of labels
    xs = X.values
    ys = Y.values
    ys = pl.reshape(ys,[ys.shape[0],])
    
    pl.figure(figsize=(15, 6), dpi=100)
    for i in range(many_features):
        if (i==2):
            counts0, bins0 = pl.histogram(xs[ys==0,i],100,range=(0.,0.08))
            counts1, bins1 = pl.histogram(xs[ys==1,i],100,range=(0.,0.08))
        elif (i==5):
            counts0, bins0 = pl.histogram(xs[ys==0,i],100,range=(1,5))
            counts1, bins1 = pl.histogram(xs[ys==1,i],100,range=(1,5))
        elif (i==6):
            counts0, bins0 = pl.histogram(xs[ys==0,i],100,range=(0,15))
            counts1, bins1 = pl.histogram(xs[ys==1,i],100,range=(0,15))
        elif (i==7):
            counts0, bins0 = pl.histogram(xs[ys==0,i],100,range=(-1.5,1.))
            counts1, bins1 = pl.histogram(xs[ys==1,i],100,range=(-1.5,1.))	      
        else:
            counts0, bins0 = pl.histogram(xs[ys==0,i],100)
            counts1, bins1 = pl.histogram(xs[ys==1,i],100)
        pl.hold()
        pl.subplot(2,4,i+1)
        pl.plot(bins0[0:100],counts0,'r',bins1[0:100],counts1,'b')
        pl.title(feature_names[i])
    pl.tight_layout()
    pl.savefig("../out/{0}/{1}".format(WHICH_EXP,name),bbox_inches='tight')

def read_big_data(filename, type=None, format='ascii'):

    if format == 'ascii':
        data = filename.replace('.cfg', '_' + type + '.dat')
    elif format == 'hdf5':
        data = filename.replace('.cfg', '.' + type + '.h5')

    try:
        A = plt.np.load(data.replace('.dat', '.npy'))
    except IOError:
        if format == 'ascii':
            n_blocks = get_first_newline(data)
            A = plt.loadtxt(data).T
            [n_cols, n_turns, n_blocks] = [A.shape[0], A.shape[1] / n_blocks, n_blocks]
            A = plt.reshape(A, (n_cols, n_turns, n_blocks))
            A = plt.rollaxis(A, 1, 3)

            plt.np.save(data.replace('.dat', '.npy'), A)
        elif format == 'hdf5':
            A = h5py.File(data, 'r')
            turns = natsort.natsorted(A.keys(), signed=False)
            cols = ['z0', 'x', 'xp', 'y', 'yp', 'z', 'zp']
            try:
                n_cols, n_particles, n_turns = len(cols), len(A[turns[0]][cols[0]]), len(turns)

                B = plt.zeros((n_cols, n_particles, n_turns))
                for i, d in enumerate(cols):
                    for j, g in enumerate(turns):
                        B[i, :, j] = A[g][d]
            except KeyError:
                B = A['Fields']['Data'][:].T

            A = B

    return A, data

def render_network(A):
    [L, M] = shape(A)
    sz = int(sqrt(L))
    buf = 1
    A = asarray(A)

    if floor(sqrt(M)) ** 2 != M:
        m = int(sqrt(M / 2))
        n = M / m
    else:
        m = int(sqrt(M))
        n = m

    array = -ones([buf + m * (sz + buf), buf + n * (sz + buf)], "d")

    k = 0
    for i in range(m):
        for j in range(n):
            clim = max(abs(A[:, k]))
            x_offset = buf + i * (sz + buf)
            y_offset = buf + j * (sz + buf)
            array[x_offset : x_offset + sz, y_offset : y_offset + sz] = reshape(A[:, k], [sz, sz]) / clim

            k += 1
    return array

def homog2D(xPrime, x):
    """
    
    Compute the 3x3 homography matrix mapping a set of N 2D homogeneous 
    points (3xN) to another set (3xN)

    """

    numPoints = xPrime.shape[1]
    assert numPoints >= 4

    A = None
    for i in range(0, numPoints):
        xiPrime = xPrime[:, i]
        xi = x[:, i]

        Ai_row0 = pl.concatenate((pl.zeros(3), -xiPrime[2] * xi, xiPrime[1] * xi))
        Ai_row1 = pl.concatenate((xiPrime[2] * xi, pl.zeros(3), -xiPrime[0] * xi))
        Ai = pl.row_stack((Ai_row0, Ai_row1))

        if A is None:
            A = Ai
        else:
            A = pl.vstack((A, Ai))

    U, S, V = pl.svd(A)
    V = V.T
    h = V[:, -1]
    H = pl.reshape(h, (3, 3))
    return H

def homog3D(points2d, points3d):
    """
    
    Compute a matrix relating homogeneous 3D points (4xN) to homogeneous
    2D points (3xN)

    Not sure why anyone would do this.  Note that the returned transformation 
    *NOT* an isometry.  But it's here... so deal with it.

    """

    numPoints = points2d.shape[1]
    assert numPoints >= 4

    A = None
    for i in range(0, numPoints):
        xiPrime = points2d[:, i]
        xi = points3d[:, i]

        Ai_row0 = pl.concatenate((pl.zeros(4), -xiPrime[2] * xi, xiPrime[1] * xi))
        Ai_row1 = pl.concatenate((xiPrime[2] * xi, pl.zeros(4), -xiPrime[0] * xi))
        Ai = pl.row_stack((Ai_row0, Ai_row1))

        if A is None:
            A = Ai
        else:
            A = pl.vstack((A, Ai))

    U, S, V = pl.svd(A)
    V = V.T
    h = V[:, -1]
    P = pl.reshape(h, (3, 4))
    return P

def f_l2_wd(x0, S, I, gamma):
    M = shape(S)[0]
    L, batch_size = shape(I)
    A = matrix(reshape(x0,[L,M]))
    E = I - A*S
    f = 0.5*(E.A**2).sum()/batch_size + 0.5*gamma*(A.A**2).sum()
    return f

def g_l2_wd(x0, S, I, gamma):
    M = shape(S)[0]
    L, batch_size = shape(I)
    A = matrix(reshape(x0,[L,M]))
    E = I - A*S
    g = -E*S.T/batch_size + gamma*A
    return g.A1

    def simulate(self, T):

        """Simulates the full neural field model

		Arguments
		----------

		T: ndarray
				simulation time instants
		Returns
		----------
		V: list of matrix
			each matrix is the neural field at a time instant

		Y: list of matrix
			each matrix is the observation vector corrupted with noise at a time instant
		"""

        Y = []
        V = []

        spatial_location_num = (len(self.field_space)) ** 2
        sim_field_space_len = len(self.field_space)

        # initial field
        v0 = self.Sigma_e_c * pb.matrix(np.random.randn(spatial_location_num, 1))
        v_membrane = pb.reshape(v0, (sim_field_space_len, sim_field_space_len))
        firing_rate = self.act_fun.fmax / (1.0 + pb.exp(self.act_fun.varsigma * (self.act_fun.v0 - v_membrane)))

        for t in T[1:]:

            v = self.Sigma_varepsilon_c * pb.matrix(np.random.randn(len(self.obs_locns), 1))
            w = pb.reshape(
                self.Sigma_e_c * pb.matrix(np.random.randn(spatial_location_num, 1)),
                (sim_field_space_len, sim_field_space_len),
            )

            print "simulation at time", t
            g = signal.fftconvolve(self.K, firing_rate, mode="same")
            g *= self.spacestep ** 2
            v_membrane = self.Ts * pb.matrix(g) + self.xi * v_membrane + w
            firing_rate = self.act_fun.fmax / (1.0 + pb.exp(self.act_fun.varsigma * (self.act_fun.v0 - v_membrane)))
            # Observation
            Y.append((self.spacestep ** 2) * (self.C * pb.reshape(v_membrane, (sim_field_space_len ** 2, 1))) + v)
            V.append(v_membrane)

        return V, Y

def objfun_l2_wd(x0, S, I, gamma):
    M = shape(S)[0]
    L, batch_size = shape(I)
    A = matrix(reshape(x0,[L,M]))
    E = I - A*S 
    f = 0.5*(E.A**2).sum()/batch_size + 0.5*gamma*(A.A**2).sum()
    g = -E*S.T/batch_size + gamma*A   
    return (f,g.A1)

def readDatDirectory(key, directory):
    global stats
    #Don't read data in if it's already read
    if not key in DATA["mean"]:
        data = defaultdict(array)

        #Process the dat files
        for datfile in glob.glob(directory + "/*.dat"):
            fileHandle = open(datfile, 'rb')
            keys, dataDict = csvExtractAllCols(fileHandle)
            stats = union(stats, keys)
            for aKey in keys:
                if not aKey in data:
                    data[aKey] = reshape(array(dataDict[aKey]),
                                         (1, len(dataDict[aKey])))
                else:
                    data[aKey] = append(data[aKey],
                                        reshape(array(dataDict[aKey]),
                                                (1, len(dataDict[aKey]))),
                                        axis=0)

        #Process the div files'
        for datfile in glob.glob(directory + "/*.div"):
            fileHandle = open(datfile, 'rb')
            keys, dataDict = csvExtractAllCols(fileHandle)
            stats = union(stats, keys)
            for aKey in keys:
                if not aKey in data:
                    data[aKey] = reshape(array(dataDict[aKey]),
                                         (1, len(dataDict[aKey])))
                else:
                    data[aKey] = append(data[aKey],
                                        reshape(array(dataDict[aKey]),
                                                (1, len(dataDict[aKey]))),
                                        axis=0)

        #Iterate through the stats and calculate mean/standard deviation
        for aKey in stats:
            if aKey in data:
                DATA["mean"][key][aKey] = mean(data[aKey], axis=0)
                DATA["median"][key][aKey] = median(data[aKey], axis=0)
                DATA["std"][key][aKey] = std(data[aKey], axis=0)
                DATA["ste"][key][aKey] = std(data[aKey], axis=0)/ sqrt(len(data[aKey]))
                DATA["min"][key][aKey] = mean(data[aKey], axis=0)-amin(data[aKey], axis=0)
                DATA["max"][key][aKey] = amax(data[aKey], axis=0)-mean(data[aKey], axis=0)
                DATA["actual"][key][aKey] = data[aKey]

def load_data(filename):
    state="end"
    dat=Data()
    valdep=[]
    infile=open(filename)
    for line in infile.readlines():
        tagfnd=re.search("\<([^<>]*)\>",line)
        if tagfnd:
            tag=tagfnd.group(1)
            tag.strip()
            if re.search("Qucs Dataset ",line):
                continue
            if tag[0]=="/":
                state="end"
                print "Number of Dimensions:",len(valdep)
                if len(valdep)>1:
                    shape=[]
                    for i in range(len(valdep),0,-1):
                        shape.append(dat.__dict__[valdep[i-1]].len)
                    val=pylab.array(val)                        
                    val=pylab.reshape(val,shape)
                dat.__dict__[name]=val
            else:
                state="start"
                words=tag.split()
                type=words[0]
                name=words[1].replace(".","_")
                name=name.replace(",","")
                name=name.replace("[","")
                name=name.replace("]","")
                if type=="indep":
                    val=Val("f")
                    val.len=int(words[2])
                else:
                    val=[]
                    valdep=words[2:]
        else:
            if state=="start":
                if "j" in line:
                    print line
                    line=line.replace("j","")
                    line="%sj"%line.strip()
                    try:
                        val.append(complex(line))                
                    except:
                        traceback.print_exc()
                        print line # add nan check
                        print name
                        print len(val)
                else:                    
                    val.append(float(line))                
            else:
                print "Parser Error:",line

    return dat      

	def simulate(self,T):

		"""Simulates the full neural field model

		Arguments
		----------

		T: ndarray
				simulation time instants
		Returns
		----------
		V: list of matrix
			each matrix is the neural field at a time instant

		Y: list of matrix
			each matrix is the observation vector corrupted with noise at a time instant
		"""

		Y=[]
		V=[]  

		spatial_location_num=(len(self.simulation_space_x_y))**2
		sim_field_space_len=len(self.simulation_space_x_y) 

		#initial field
		v0=pb.dot(self.Sigma_e_c,np.random.randn(spatial_location_num,1))
		v_membrane=pb.reshape(v0,(sim_field_space_len,sim_field_space_len))
        
		for t in T[1:]:
            
			v = pb.dot(self.Sigma_varepsilon_c,np.random.randn(len(self.obs_locns),1))
			w = pb.reshape(pb.dot(self.Sigma_e_c,np.random.randn(spatial_location_num,1)),(sim_field_space_len,sim_field_space_len))
			#print "simulation at time",t
			g=signal.fftconvolve(self.K,v_membrane,mode='same') 
			g*=(self.spacestep**2)
			v_membrane=g+w
            
			#Observation
			Y.append((self.spacestep**2)*(pb.dot(self.C,pb.reshape(v_membrane,(sim_field_space_len**2,1))))+v)
			V.append(v_membrane)

		return V,Y

def choose_patches(IMAGES, L, batch_size=1000):
    sz = int(sqrt(L))
    imsz = shape(IMAGES)[0]
    num_images = shape(IMAGES)[2]
    BUFF = 4

    X = matrix(zeros([L,batch_size],'d'))
    for i in range(batch_size):
        j = int(floor(num_images * rand()))
        r = sz/2+BUFF+int(floor((imsz-sz-2*BUFF)*rand()))
        c = sz/2+BUFF+int(floor((imsz-sz-2*BUFF)*rand()))
        X[:,i] = reshape(IMAGES[r-sz/2:r+sz/2,c-sz/2:c+sz/2,j],[L,1])
    return X

def getDissim(data, atype, vbose=0, minRank=0, maxRank=50, NPOZ=50):

    ks = data.keys()

    matr = pylab.ones(len(ks)**2)
    matr = pylab.reshape(matr, (len(ks), len(ks)))
    scs = []
    names = []

    for ik, k_con in enumerate(ks):
        name = ik
        if not k_con in names:
            names.append(k_con)
        for jk, k_pl in enumerate(ks):

            ss1 = computeSimSc(data, k_con, k_pl, vbose=vbose, minRank=minRank, maxRank=maxRank, NPOZ=NPOZ)
            ss2 = computeSimSc(data, k_pl, k_con, vbose=vbose, minRank=minRank, maxRank=maxRank, NPOZ=NPOZ)

            
            if atype=='abs':
                sc1 = sum(ss1)
                sc2 = sum(ss2)
            elif atype=='rms':
                sc1 = pylab.sqrt(pylab.sum(ss1**2))
                sc2 = pylab.sqrt(pylab.sum(ss2**2))
            elif atype=='met':
                sc1 = sum(pylab.logical_and(ss1!=0, True)) 
                sc2 = sum(pylab.logical_and(ss2!=0, True))

            if vbose>=1:
                print 'score for ', k_con, k_pl, ss1, sc1, ss2, sc2

            oldsc = sc1 + sc2
            oldsc *= 0.5


            l1 = len(data[k_con])
            l2 = len(data[k_pl])
            iscale = min(l1, l2)
            nsc = oldsc/(1.0*iscale)
            if vbose>=1:
                print k_con, k_pl, 'oldsc', oldsc, l1, l2, iscale, 'nsc', nsc

            
            matr[ik][jk] = nsc
            if jk<=ik:
                continue

            print nsc, 'xx', ik, k_con, jk, k_pl
            scs.append(nsc)
    return names, pylab.array(scs), matr

def plotit(x, y, z, title):
    plb.figure()
    try:
        X = list()
        Y = list()
        [X.append(i) for i in x if i not in X]
        [Y.append(i) for i in y if i not in Y]
        Z = plb.reshape(z, (len(Y), len(X)))
        plb.contourf(X, Y, Z, 10)
    except:
        plb.scatter(x, y, c=z)
    plb.title(title)
    plb.colorbar()
    plb.show()