def standing(self, version, plot):
        dx = 0.1
        dy = 0.1
        self.dx = dx
        self.dy = dy
        b = 0.0
        V = 0.0
        Lx = 5
        Ly = 5
        T = 6

        test = "plug"
        A = 1
        B = 1
        mx = 2.0
        my = 2.0
        kx = mx * np.pi / Lx
        ky = my * np.pi / Ly

        I = np.vectorize(lambda x, y: A * np.cos(kx * x) * np.cos(ky * y))
        q = np.vectorize(lambda x, y: 2)
        f = np.vectorize(lambda x, y, t: 0)
        c = np.sqrt(q(1, 1))

        w = np.sqrt(q(1, 1) * (kx ** 2 + ky ** 2))
        ue = np.vectorize(lambda x, y, t: A * np.cos(kx * x) * np.cos(ky * y) * np.cos(w * t))

        m = 6
        h = 1.0
        err = np.zeros(m + 1)
        R, T, V = convergence_rates(version, solver, ue, m, h, Lx, Ly, T, b, V, I, q, f, plot=False)
        for i in range(len(R)):
            print "For dt = %.4f, the convergence rate is %.4f" % (T[i], R[i])

def make_gaussian_conv_plots(m, nvals):
    """
    Compute mean of means, std of means for many values of n
    Plot against expected mean, std of mean as function of n
    """
    mean_of_means = np.vectorize(lambda x: np.mean(mean_unif_samp(m,x)))
    std_of_means = np.vectorize(lambda x: np.std(mean_unif_samp(m,x)))
    
    mns = mean_of_means(nvals)
    stds = std_of_means(nvals)
    
    mu = 0.5
    std_expect = std_of_mean_unif(m)
    
    plt.plot(nvals,mns, 'ko')
    plt.axhline(mu)
    plt.xscale('log')
    plt.xlabel('Number of %d-sample sums' % (m))
    plt.ylabel('Mean of means')
    plt.show()
    
    plt.plot(nvals,stds, 'ko')
    plt.axhline(std_expect)
    plt.xscale('log')
    plt.xlabel('Number of %d-sample sums' % (m))
    plt.ylabel('Standard deviations of means')
    plt.show()

    def fire(self, x):
        sigmoid = lambda x: 1. / (1. + np.exp(-x))

        z = np.vectorize(sigmoid)(self.hidden_weight.dot(np.r_[np.array([1]), x]))
        y = np.vectorize(sigmoid)(self.output_weight.dot(np.r_[np.array([1]), z]))

        return (z, y)

def solve_linalg(k, T, F0, F1, f):
    N, h = len(X), L/len(X)
    I = np.eye(N)
    S,Y = np.meshgrid(X,X)
    abs_func = np.vectorize(apply_abs)
    F0, F1 = partial(F0, k), partial(F1, k)
    G0 = lambda i,j: abs_func(i, j, F0, F0) - b*h
    G1 = lambda i,j: abs_func(i, j, F1, lambda x: -F1(x)) - b*h*h*(i+j-.5)
    A = weight_matrix(
        lambda i,j: (j-i)*G0(-i,j) - G1(-i,j)/h,
        lambda i,j: G1(-i,j)/h - (j-i-1)*G0(-i,j)
    )
    B = weight_matrix(
        lambda i,j: (j+i)*G0(i,j) - G1(i,j)/h,
        lambda i,j: G1(i,j)/h - (j+i-1)*G0(i,j)
    )
    #splot(X, A*T(np.abs(S-Y)/k)/k)
    #splot(X, B*T((S+Y)/k)/k)
    #py.show()
    phi = solve(a*I - A*T(np.abs(S-Y)/k)/k + B*T((S+Y)/k)/k, f(X))
    p_xy = -(k*(T2(0)-T2(1./k)) + np.trapz((T1((L-X)/k) - T1((L+X)/k))*phi, X))*2/a
    Phi = np.outer(phi,np.ones(N))
    Q = np.trapz(T2((L-X)/k)-T2((L+X)/k) + np.trapz((T1(np.abs(S-Y)/k) - T1((S+Y)/k))*Phi, X)/k, X)/2/a
    #splot(X, K(*XX))
    #py.plot(X, phi)
    w = np.vectorize(lambda x: 0 if x==0 else x*np.log(x)/a)
    ww = lambda x: k*x*x*(2*np.log(x)-1)/4/a
    #print >> sys.stderr, k, np.trapz(phi, X), np.trapz(phi - w((L-X)/k), X) + ww(L/k)
    print >> sys.stderr, k, p_xy, np.trapz(phi, X)/2, Q
    #np.savetxt(sys.stdout, np.transpose((X, phi)), fmt='%1.4e')
    return k, p_xy, np.trapz(phi, X)/2, Q

    def __call__(self, mass, taper=False, integral_form=False, **kwargs):
        """ Unclear if integral_form is right..."""
        if taper:

            def num_func(x, mass_):
                tf = (1-(mass_/x)**(1-self.j))**0.5
                return self.imf(x)*(1./x)**(1-self.j) * (2/((1+self.Rmdot**2*x**1.5)**0.5+1)) * tf

            def integrate(lolim, mass_):
                integral = scipy.integrate.quad(num_func, lolim, self.mmax, args=(mass_,), **kwargs)[0]
                return integral

            numerator = np.vectorize(integrate)(np.where(self.mmin < mass, mass, self.mmin), mass)

        else:
            def num_func(x):
                return self.imf(x)*(1./x)**(1-self.j) * (2/((1+self.Rmdot**2*x**1.5)**0.5+1))

            def integrate(lolim):
                integral = scipy.integrate.quad(num_func, lolim, self.mmax, **kwargs)[0]
                return integral

            numerator = np.vectorize(integrate)(np.where(self.mmin < mass, mass, self.mmin))

        result = (1-self.j) * mass**(1-self.j) * numerator / self.denominator
        if integral_form:
            warnings.warn("The 'integral form' of the Chabrier PMF is not correctly normalized; "
                          "it is just PMF(m) * m")
            return result * self.normfactor * mass
            raise ValueError("Integral version not yet computed")
        else:
            return result * self.normfactor

    def get_dir_tot_ele(x, y, z, slope, aspect):
        sun_gap = get_sun_gap_ele(x, y)
        ele = z

        def get_incidence_ele(zeni, azi, slope_in, aspect_in):
            inc = np.arccos(np.cos(np.deg2rad(zeni)) * np.cos(np.deg2rad(slope_in))
                            + np.sin(np.deg2rad(zeni)) * np.sin(np.deg2rad(slope_in))
                            * np.cos(azi - np.deg2rad(aspect_in)))
            return inc

        f_incidence = np.vectorize(get_incidence_ele, excluded=['slope_in', 'aspect_in'])
        incidence = f_incidence(zeni_c, azi_c, slope, aspect)

        def get_trans_ele(zeni, elevation, transmittance):
            m = np.exp(-0.000118 * elevation - 1.638 * np.power(10, -9) * elevation * elevation) / \
                (np.cos(np.deg2rad(zeni)) + 0.50572 * np.power((96.07995 - zeni), -1.6364))
            return np.power(transmittance, m)

        f_trans = np.vectorize(get_trans_ele, excluded=['elevation', 'transmittance'])
        trans = f_trans(zeni_c, ele, beta)

        def get_dir_ele(time_p_ele, sun_gap_ele, incidence_ele, trans_ele):
            # incidence = get_incidence_angle(zeni_c_ele, azi_c_ele, slope, aspect)
            dir_ele = 1.367 * trans_ele * time_p_ele * sun_gap_ele * np.cos(incidence_ele)
            if dir_ele < 0.0:
                return 0.0
            return dir_ele

        f_dir = np.vectorize(get_dir_ele)
        dir_array = f_dir(time_p, sun_gap, incidence, trans)
        dir_tot = np.sum(dir_array)
        return dir_tot

    def __call__(self, luminosity, taper=False, integral_form=False, **kwargs):
        """ Unclear if integral_form is right..."""
        if taper:

            def num_func(x, luminosity_):
                tf = (1-(luminosity_/x)**(1-self.j))**0.5
                return self.imf(x)*x**(self.j-self.jf-1) * tf

            def integrate(lolim, luminosity_):
                integral = scipy.integrate.quad(num_func, lolim, self.mmax, args=(luminosity_,), **kwargs)[0]
                return integral

            numerator = np.vectorize(integrate)(np.where(self.mmin < luminosity, luminosity, self.mmin), luminosity)

        else:
            def num_func(x):
                return self.imf(x)*x**(self.j-self.jf-1)

            def integrate(lolim):
                integral = scipy.integrate.quad(num_func, lolim, self.mmax, **kwargs)[0]
                return integral

            numerator = np.vectorize(integrate)(np.where(self.mmin < luminosity, luminosity, self.mmin))

        result = (1-self.j) * luminosity**(1-self.j) * numerator / self.denominator
        if integral_form:
            warnings.warn("The 'integral form' of the Chabrier PMF is not correctly normalized; "
                          "it is just PMF(m) * m")
            return result * self.normfactor * luminosity
            raise ValueError("Integral version not yet computed")
        else:
            return result * self.normfactor

def velocity_field(psi): #takes a symbolic function and returns two lambda functions
#to evaluate the derivatives in both x and y.
   global w
   if velocity_components:
      u = lambdify((x,y), eval(x_velocity), modules='numpy')
      v = lambdify((x,y), eval(y_velocity), modules='numpy')
   else:
      if is_complex_potential:
         print "Complex potential, w(z) given"
         #define u, v symbolically as the imaginary part of the derivatives
         u = lambdify((x, y), sympy.im(psi.diff(y)), modules='numpy')
         v = lambdify((x, y), -sympy.im(psi.diff(x)), modules='numpy')
      else:
         #define u,v as the derivatives 
         print "Stream function, psi given"
         u = sympy.lambdify((x, y), psi.diff(y), 'numpy')
         v = sympy.lambdify((x, y), -psi.diff(x), 'numpy')
   if (branch_cuts): # If it's indicated that there are branch cuts in the mapping,
                      # then we need to return vectorized numpy functions to evaluate
                      # everything numerically, instead of symbolically 
                      # This of course results in a SIGNIFICANT time increase
                      #   (I don't know how to handle more than the primitive root
                      #   (symbolically in Sympy
      return np.vectorize(u),np.vectorize(v)
   else:
       # If there are no branch cuts, then return the symbolic lambda functions (MUCH faster)
      return u,v

    def df_two_degree_counts(self, df, col_i, col_j, operation, file_path = None):
        # if two columns are the same just return one_degree_count
        if col_i == col_j:
            return self.df_one_degree_counts(df, col_i, file_path)

        if operation == 'per':
            task = 'both'
        elif operation == 'num':
            task = 'two'
        else:
            print 'unknown operation'
            return
                
        self._df_get_count_tables(df, task, [col_i, col_j], file_path)

        i_table = one_count_table[col_i]
        ij_table = two_count_table[(col_i, col_j)]

        col_i_data = df[col_i].values
        col_j_data = df[col_j].values
        if operation == 'per':  # 'per': percentage of (elem_i, elem_j) in all (elem_i, col_j)  
            vfunc = np.vectorize(lambda x,y: float(ij_table[x][y])/i_table[x])
            col_new = vfunc(col_i_data, col_j_data)
        elif operation == 'num':    # 'num': number of different kinds of (elem_i, col_j) 
            vfunc = np.vectorize(lambda x: ij_table[x]['unique'])
            col_new = vfunc(col_i_data)

        return col_new

 def readConfig(self):
     configFile = open(self.CNSconfig, 'r')
     self.CNSconfigInfo = np.append(np.vectorize(lambda x: parseConfigFindPath(x,configFile))(['root_folder', 'pathPython',
                              'checkValidity','conservedFastaPath','pickleSkip','pickleName','fasta2phylip','PhyML',
                               'bootstrap','treeFile','treeOut','ratioCopy','outputTreeImages']),
                                    np.vectorize(lambda x: parseConfigFindList(x,configFile))(['masterListSpecies','intragenus','intergenus','subgenome']))
     configFile.close()

def train_data(N=300):
    sig=3
    X=np.array(np.random.uniform(low=0,high=4*np.math.pi,size=[N,1]))
    f_sin=np.vectorize(sp.sin,otypes=[np.float])
    f_cos=np.vectorize(sp.cos,otypes=[np.float])
    Y=30*f_sin(X)+30*f_cos(2*X+4)+sig*np.array(sp.randn(N,1))
    return [X,Y]

def find_ephemeris_lookup_date(tai_beg, tai_end, obs_md_table):
    '''
    Want to find the 15-minute increment (0, 15, 30, 45) that is sandwiched between the two
    passed datetimes.  However, spans can be less than 15 minutes, so also need handle this
    case; here, will round to whichever increment has the smallest delta between tai_beg
    and tai_end.  I've not seen it, but I have to assume that spans can also be longer than
    15 minutes.
    '''
    vectfunc = np.vectorize(tai_str_to_datetime)
    tai_end_dt = vectfunc(tai_end)
    tai_beg_dt = vectfunc(tai_beg)

    vectfunc = np.vectorize(get_ephemeris_block_in_interval)

    mask = (tai_end_dt - tai_beg_dt) <= ephemeris_max_block_size
    ret = np.zeros((len(tai_end_dt),), dtype=dt.datetime)
    ret[mask] = vectfunc(tai_beg_dt[mask], tai_end_dt[mask])

    def _lookup_str_format(dtval):
        if isinstance(dtval, dt.datetime):
            return dtval.strftime("%Y-%b-%d %H:%M")
        return ""

    vectfunc = np.vectorize(_lookup_str_format)
    ret = vectfunc(ret)

    return ret[mask], obs_md_table[mask]

    def cohort_results(self, group_fields):
        df_grouped = self.dataframe.groupby(group_fields)

        if self.coalation_type == 'best':
            df = df_grouped.max()
            df['pob'] = np.vectorize(percentage_marks)(df['pMark'], df['qMark'])
            df = df.drop(['pMark', 'qMark', 'aID', 'qNum'], axis=1)
            df['Colour'] = np.vectorize(colourise)(df['pob'])
            df.rename(columns={'pob': 'Cohort Max percentage'}, inplace=True)
            if group_fields == ['qTopic']:
                df = df.drop(['qModule'], axis=1)
            elif group_fields == ['qModule']:
                df = df.drop(['qTopic'], axis=1)
        elif self.coalation_type == 'mean':
            df = df_grouped.mean()
            df['pob'] = np.vectorize(percentage_marks)(df['pMark'], df['qMark'])
            df = df.drop(['pMark', 'qMark', 'aID', 'qNum'], axis=1)
            # round off ALL values in the dataframe to 2 decimal places
            df = np.round(df, 2)
            df['Colour'] = np.vectorize(colourise)(df['pob'])
            df.rename(columns={'pob': 'Cohort Mean percentage'}, inplace=True)
        else:
            QtGui.QMessageBox.question(
                self,
                'Uh oh...',
                ("Something went wrong."
                 "Please close the analysis screen and try again."))
        return df

    def class_results(self, tSet, group_fields):
        filtered = self.dataframe[self.dataframe['teaching_set'].isin([tSet])]
        cols = ['qModule', 'qTopic', 'qMark', 'pMark']
        df = filtered[cols]
        df_grouped = df.groupby(group_fields)

        if self.coalation_type == 'best':
            df = df_grouped.max()
            df['pob'] = np.vectorize(percentage_marks)(df['pMark'], df['qMark'])
            df['Colour'] = np.vectorize(colourise)(df['pob'])
            df.rename(columns={'pob': 'Class Max percentage'}, inplace=True)
            df = df.drop(['pMark', 'qMark'], axis=1)
            if group_fields == ['qTopic']:
                df = df.drop(['qModule'], axis=1)
            elif group_fields == ['qModule']:
                df = df.drop(['qTopic'], axis=1)
        elif self.coalation_type == 'mean':
            df = df_grouped.mean()
            df['pob'] = np.vectorize(percentage_marks)(df['pMark'], df['qMark'])
            # Round off ALL values in the dataframe to 2 decimal places
            # otherwise the mean function will give ridiculous accuracy!
            df = np.round(df, 2)
            df['Colour'] = np.vectorize(colourise)(df['pob'])
            df = df.drop(['pMark', 'qMark'], axis=1)
            df.rename(columns={'pob': 'Class Mean percentage'}, inplace=True)
        else:
            QtGui.QMessageBox.question(
                self,
                'Uh oh...',
                ("Something went wrong."
                 "Please close the analysis screen and try again."))
        return df

def getDataFromModel(ModelName):
	fin=open(ModelName,"r")
	fin.readline()
	fin.readline()
	# meanShape=[]
	# pcaMatrix=None
	cnt=0
	alignedSet=[]

	for line in fin.readlines():
		temp=line.strip().split(":")
		label=int(temp[0])
		data=temp[1].split(" ")
		# print data
		if label==1:
			pcaMatrix=np.array(np.vectorize(float)(data))
		elif label==2:
			meanShape=np.array(np.vectorize(float)(data))
		else:
			alignedSet.append(np.vectorize(float)(data))
	szMean=meanShape.size
	szPca=pcaMatrix.size
	pcaMatrix.reshape(szPca)
	pcaMatrix.shape=(szMean,szPca/szMean)
	meanShape.reshape(1,szMean)
	meanShape.shape=(1,szMean)
	# print meanShape.shape,pcaMatrix.shape
	# print pcaMatrix
			
	return pcaMatrix,meanShape,alignedSet

    def fset_two_degree_counts(self, myfset, col_i, col_j, operation, file_path = None):
        # if two columns are the same just return one_degree_count
        if col_i == col_j:
            return self.fset_one_degree_counts(myfset, col_i, file_path)

        if operation == 'per':
            task = 'both'
        if operation == 'num':
            task = 'two'

        self._fset_get_count_tables(myfset, task, [col_i, col_j], file_path)

        col_i_name = myfset.fname_list[col_i]
        col_j_name = myfset.fname_list[col_j]
        col_i_ind = myfset.find_list[col_i]
        col_j_ind = myfset.find_list[col_j]


        i_table = self.one_count_table[col_i_name]
        ij_table = self.two_count_table[(col_i_name, col_j_name)]

        col_i_data_train = myfset.Xtrain[:, col_i_ind]
        col_i_data_test = myfset.Xtest[:, col_i_ind]
        col_j_data_train = myfset.Xtrain[:, col_j_ind]
        col_j_data_test = myfset.Xtest[:, col_j_ind]
        if operation == 'per':  # 'per': percentage of (elem_i, elem_j) in all (elem_i, col_j)  
            vfunc = np.vectorize(lambda x,y: float(ij_table[x][y])/i_table[x])
            col_new_train = vfunc(col_i_data_train, col_j_data_train)
            col_new_test = vfunc(col_i_data_test, col_j_data_test)
        elif operation == 'num':    # 'num': number of different kinds of (elem_i, col_j) 
            vfunc = np.vectorize(lambda x: ij_table[x]['unique'])
            col_new_train = vfunc(col_i_data_train)
            col_new_test = vfunc(col_i_data_test)

        return col_new_train, col_new_test

def filter_xaod_to_numpy(files, max_events=None):
    """Processes some files by converting to numpy and applying filtering"""
    # Branch name remapping for convenience
    branch_dict = {
        'CaloCalTopoClustersAuxDyn.calEta' : 'clusEta',
        'CaloCalTopoClustersAuxDyn.calPhi' : 'clusPhi',
        'CaloCalTopoClustersAuxDyn.calE' : 'clusE',
        'CaloCalTopoClustersAuxDyn.EM_PROBABILITY' : 'clusEM',
        'AntiKt10LCTopoTrimmedPtFrac5SmallR20JetsAux.pt' : 'fatJetPt',
        'AntiKt10LCTopoTrimmedPtFrac5SmallR20JetsAux.eta' : 'fatJetEta',
        'AntiKt10LCTopoTrimmedPtFrac5SmallR20JetsAux.phi' : 'fatJetPhi',
        'AntiKt10LCTopoTrimmedPtFrac5SmallR20JetsAux.m' : 'fatJetM',
        'EventInfoAuxDyn.mcChannelNumber' : 'dsid',
        'EventInfoAuxDyn.mcEventWeights' : 'genWeight',
        'InDetTrackParticlesAuxDyn.theta' : 'trackTheta',
        'InDetTrackParticlesAuxDyn.phi' : 'trackPhi',
    }
    # Convert the data to numpy
    print('Now processing:', files)
    tree = get_tree(files, branch_dict, tree_name='CollectionTree',
                    max_events=max_events)
    if tree is None:
        return None
    # Apply physics
    results = process_events(tree)
    skimTree = results['tree']

    # Get the track coordinates
    vtan = np.vectorize(np.tan, otypes=[np.ndarray])
    vlog = np.vectorize(np.log, otypes=[np.ndarray])
    trackTheta = skimTree['trackTheta']
    results['trackEta'] = -vlog(vtan(trackTheta / 2))
    results['trackPhi'] = skimTree['trackPhi']

    return results

def RESDIFF(x,Q,f0,aleak,ph1,da,ang1,Igain,Qgain,Ioff,Qoff):
#       Q = p[0]          ;  Q
#       f0 = p[1]         ;  resonance frequency
#       aleak = p[2]      ;  amplitude of leakage
#       ph1 = p[3]        ;  phase shift of leakage
#       da = p[4]         ;  variation of carrier amplitude
#       ang1 = p[5]       ;  Rotation angle of data
#       Igain = p[6]      ;  Gain of I channel
#       Qgain = p[7]      ;  Gain of Q channel
#       Ioff = p[8]       ;  Offset of I channel
#       Qoff = p[9]       ;  Offset of Q channel

    l = len(x)
    dx = (x - f0) / f0

    # resonance dip function
    s21a = (np.vectorize(complex)(0,2.0*Q*dx)) / (complex(1,0) + np.vectorize(complex)(0,2.0*Q*dx))
    s21a = s21a - complex(.5,0)
    s21b = np.vectorize(complex)(da*dx,0) + s21a + aleak*np.vectorize(complex)(1.0-np.cos(dx*ph1),-np.sin(dx*ph1))

    # scale and rotate
    Ix1 = s21b.real*Igain
    Qx1 = s21b.imag*Qgain
    nI1 = Ix1*np.cos(ang1) + Qx1*np.sin(ang1)
    nQ1 = -Ix1*np.sin(ang1) + Qx1*np.cos(ang1)

    #scale and offset
    nI1 = nI1 + Ioff
    nQ1 = nQ1 + Qoff

    s21 = np.zeros(l*2)
    s21[:l] = nI1
    s21[l:] = nQ1

    return s21

 def unpack(self):
     """
     creates the vectorized functions
     intention: call this method after loading this object from serialization
     """
     self.fnc = vectorize(self.fnc_nv)
     self.deriv_fnc = vectorize(self.deriv_fnc_nv)

def mat_backprop(weights, x, t):
  gamma = 0.8
  vec_sigmoid = np.vectorize(sigmoid)
  vfunc = np.vectorize(quad_error)

  outputs = [x]

  o = x
  for mat in weights:
    extended_o = extend(o)
    o = vec_sigmoid(extended_o.T.dot(mat)).T
    outputs.append(o)

  e = o - t
  print e
  diags = []

  for oi in outputs[1:]:
    diags.append(np.diagflat(vfunc(oi)))

  deltas = [0]*len(weights)
  deltas[-1] = diags[-1].dot(e)
  for idx in range(len(weights)-1)[::-1]:
    deltas[idx] = diags[idx].dot(weights[idx+1][:-1]).dot(deltas[idx+1])

  wnew = []
  for i in xrange(len(weights)):
    z = deltas[i].dot(extend(outputs[i]).T)
    wnew.append(weights[i] - gamma*z.T)

  return wnew