Python pylab.ones函数代码示例

def validate_once(true_cf = [pl.ones(3)/3.0, pl.ones(3)/3.0], true_std = 0.01*pl.ones(3), std_bias = [1., 1., 1.], save=False, dir='', i=0):
    """
    Generate a set of simulated estimates for the provided true cause fractions; Fit the bad model and 
    the latent simplex model to this simulated data and calculate quality metrics. 
    """ 
    
    # generate simulation data
    X = data.sim_data_for_validation(1000, true_cf, true_std, std_bias)

    # fit bad model, calculate fit metrics 
    bad_model = models.bad_model(X)
    bad_model_metrics = calc_quality_metrics(true_cf, true_std, std_bias, bad_model)
    retrieve_estimates(bad_model, True, 'bad_model', dir, i)
    
    # fit latent simplex model, calculate fit metrics 
    m, latent_simplex = models.fit_latent_simplex(X)
    latent_simplex_metrics = calc_quality_metrics(true_cf, true_std, std_bias, latent_simplex)
    retrieve_estimates(latent_simplex, True, 'latent_simplex', dir, i)
    
    # either write results to disk or return them 
    if save: 
        pl.rec2csv(bad_model_metrics, '%s/metrics_bad_model_%i.csv' % (dir, i)) 
        pl.rec2csv(latent_simplex_metrics, '%s/metrics_latent_simplex_%i.csv' % (dir, i))
    else: 
        return bad_model_metrics, latent_simplex_metrics

    def tempo_search(db, Key, tempo):
        """
        ::

            Static tempo-invariant search
            Returns search results for query resampled over a range of tempos.
        """
        if not db.configCheck():
            print "Failed configCheck in query spec."
            print db.configQuery
            return None
        prop = 1.0 / tempo  # the proportion of original samples required for new tempo
        qconf = db.configQuery.copy()
        X = db.retrieve_datum(Key)
        P = db.retrieve_datum(Key, powers=True)
        X_m = pylab.mat(X.mean(0))
        X_resamp = pylab.array(adb.resample_vector(X - pylab.mat(pylab.ones(X.shape[0])).T * X_m, prop))
        X_resamp += pylab.mat(pylab.ones(X_resamp.shape[0])).T * X_m
        P_resamp = pylab.array(adb.resample_vector(P, prop))
        seqStart = int(pylab.around(qconf["seqStart"] * prop))
        qconf["seqStart"] = seqStart
        seqLength = int(pylab.around(qconf["seqLength"] * prop))
        qconf["seqLength"] = seqLength
        tmpconf = db.configQuery
        db.configQuery = qconf
        res = db.query_data(featData=X_resamp, powerData=P_resamp)
        res_resorted = adb.sort_search_result(res.rawData)
        db.configQuery = tmpconf
        return res_resorted

def plotInit(Plotting, Elements):
	if (Plotting == 2):
		loc = [i.xy for i in Elements]
		x = [i.real for i in loc]
		y = [i.imag for i in loc]
		x = list(sorted(set(x))) 
		x.remove(-10)
		y = list(sorted(set(y)))

		X, Y = pylab.meshgrid(x, y)
		U = pylab.ones(shape(X))
		V = pylab.ones(shape(Y))

		pylab.ion()
		fig, ax = pylab.subplots(1,1)
		graph = ax.quiver(X, Y, U, V)
		pylab.draw()
	else:
		pylab.ion()
		graph, = pylab.plot(1, 'ro', markersize = 2) 
		x = 2
		pylab.axis([-x,x,x,-x])

		graph.set_xdata(0)
		graph.set_ydata(0)
		pylab.draw()

	return graph

def filter2d(x, y, axes=['y'], algos=['2sigma']):
    """
    Perform 2D data filtration by selected exes.
    In:
        x : ndarray, X vector
        y : ndarray, Y vector
        axes : list, axes names which are used to choose filtered values. x, y or any combination
    Out:
        xnew : ndarray, filtered X
        ynew : ndarray, filtered Y
    """
    xnew = pl.array(x, dtype='float')
    ynew = pl.array(y, dtype='float')
    mask_x = pl.ones(len(x), dtype='bool')
    mask_y = pl.ones(len(y), dtype='bool')
    if 'y' in axes:
        mask_y = filter1d(y,algos=algos)        
    if 'x' in axes:
        mask_x = filter1d(x,algos=algos)
    mask = mask_x * mask_y
    xnew *= mask
    ynew *= mask
    
    xnew = pl.ma.masked_equal(xnew,0)
    xnew = pl.ma.compressed(xnew)
    ynew = pl.ma.masked_equal(ynew,0)
    ynew = pl.ma.compressed(ynew)

    assert pl.shape(xnew) == pl.shape(ynew)
    return xnew, ynew

def example():

    from pylab import rand, ones, concatenate
    import matplotlib.pyplot as plt
    # EXAMPLE data code from:
    # http://matplotlib.sourceforge.net/pyplots/boxplot_demo.py
    # fake up some data
    spread= rand(50) * 100
    center = ones(25) * 50
    flier_high = rand(10) * 100 + 100
    flier_low = rand(10) * -100
    data =concatenate((spread, center, flier_high, flier_low), 0)

    # fake up some more data
    spread= rand(50) * 100
    center = ones(25) * 40
    flier_high = rand(10) * 100 + 100
    flier_low = rand(10) * -100
    d2 = concatenate( (spread, center, flier_high, flier_low), 0 )
    data.shape = (-1, 1)
    d2.shape = (-1, 1)
    #data = [data, d2, d2[::2,0]]
    data = [data, d2]

    fig = plt.figure()
    ax = fig.add_subplot(1,1,1)
    ax.set_xlim(0,4)
    percentile_box_plot(ax, data, [2,3])
    plt.show()

def jetWoGn(reverse=False):
    """
    jetWoGn(reverse=False)
       - returning a colormap similar to cm.jet, but without green.
         if reverse=True, the map starts with red instead of blue.
    """
    m=18 # magic number, which works fine
    m0=pylab.floor(m*0.0)
    m1=pylab.floor(m*0.2)
    m2=pylab.floor(m*0.2)
    m3=pylab.floor(m/2)-m2-m1

    b_ = pylab.hstack( (0.4*pylab.arange(m1)/(m1-1.)+0.6, pylab.ones((m2+m3,)) ) )
    g_ = pylab.hstack( (pylab.zeros((m1,)),pylab.arange(m2)/(m2-1.),pylab.ones((m3,))) )
    r_ = pylab.hstack( (pylab.zeros((m1,)),pylab.zeros((m2,)),pylab.arange(m3)/(m3-1.)))

    r = pylab.hstack((r_,pylab.flipud(b_)))
    g = pylab.hstack((g_,pylab.flipud(g_)))
    b = pylab.hstack((b_,pylab.flipud(r_)))

    if reverse:
        r = pylab.flipud(r)
        g = pylab.flipud(g)
        b = pylab.flipud(b)

    ra = pylab.linspace(0.0,1.0,m)

    cdict = {'red': zip(ra,r,r),
            'green': zip(ra,g,g),
            'blue': zip(ra,b,b)}

    return LinearSegmentedColormap('new_RdBl',cdict,256)

def log_inv(X): # inverts a 3x3 matrix given by the logscale values
    if (X.shape[0] != X.shape[1]):
        raise Exception("X is not a square matrix and cannot be inverted")
    
    if (X.shape[0] == 1):
        return matrix((-X[0,0]))
    
    ldet = log_det(X)
    if (ldet == nan):
        raise Exception("The determinant of X is 0, cannot calculate the inverse")
     
    if (X.shape[0] == 2): # X is a 2x2 matrix
        I = (-log_det(X)) * ones((2,2))
        I[0,0] += X[1,1]
        I[0,1] += X[0,1] + complex(0, pi)
        I[1,0] += X[1,0] + complex(0, pi)
        I[1,1] += X[0,0]
        return I
    
    if (X.shape[0] == 3): # X is a 3x3 matrix
        I = (-log_det(X)) * ones((3,3))
        I[0,0] += log_subt_exp(X[1,1]+X[2,2], X[1,2]+X[2,1])
        I[0,1] += log_subt_exp(X[0,2]+X[2,1], X[0,1]+X[2,2])
        I[0,2] += log_subt_exp(X[0,1]+X[1,2], X[0,2]+X[1,1])
        I[1,0] += log_subt_exp(X[2,0]+X[1,2], X[1,0]+X[2,2])
        I[1,1] += log_subt_exp(X[0,0]+X[2,2], X[0,2]+X[2,0])
        I[1,2] += log_subt_exp(X[0,2]+X[1,0], X[0,0]+X[1,2])
        I[2,0] += log_subt_exp(X[1,0]+X[2,1], X[2,0]+X[1,1])
        I[2,1] += log_subt_exp(X[2,0]+X[0,1], X[0,0]+X[2,1])
        I[2,2] += log_subt_exp(X[0,0]+X[1,1], X[0,1]+X[1,0])
        return I
    
    raise Exception("log_inv is only implemented for matrices of size < 4")

    def sample(self, model, evidence):
        z = evidence['z']
        T, surfaces, sigma_g, sigma_h = [evidence[var] for var in ['T', 'surfaces', 'sigma_g', 'sigma_h']]
        mu_h, phi, sigma_z_g, sigma_z_h = [model.known_params[var] for var in ['mu_h', 'phi', 'sigma_z_g', 'sigma_z_h']]
        prior_mu_g, prior_cov_g = [model.hyper_params[var] for var in ['prior_mu_g', 'prior_cov_g']]
        prior_mu_h, prior_cov_h = [model.hyper_params[var] for var in ['prior_mu_h', 'prior_cov_h']]
        n = len(g)

        y = ma.asarray(ones((n, 2))*nan)
        if sum(T==1) > 0:
            y[T==1, 0] = z[T==1]
        if sum(T==2) > 0:
            y[T==2, 1] = z[T==2]
        y[isnan(y)] = ma.masked

        kalman = self._kalman
        kalman.initial_state_mean=[prior_mu_g[0], prior_mu_h[0]]
        kalman.initial_state_covariance=diag([prior_cov_g[0,0], prior_cov_h[0,0]])
        kalman.transition_matrices=[[1, 0], [0, phi]]
        kalman.transition_offsets =ones((n, 2))*[0, mu_h*(1-phi)]
        kalman.transition_covariance=[[sigma_g**2, 0], [0, sigma_h**2]]
        kalman.observation_matrices=[[1, 0], [1, 1]]
        kalman.observation_covariance=[[sigma_z_g**2, 0], [0, sigma_z_h**2]]
        sampled_surfaces = forward_filter_backward_sample(kalman, y)

        return sampled_surfaces

def run_on_cluster(dir='../data', true_cf = [pl.ones(3)/3.0, pl.ones(3)/3.0], true_std = 0.01*pl.ones(3), std_bias=[1.,1.,1.], reps=5, tag=''):
    """
    Runs validate_once multiple times (as specified by reps) for the given true_cf and 
    true_std. Combines the output and cleans up the temp files. This accomplished in 
    parallel on the cluster. This function requires that the files cluster_shell.sh 
    (which allows for submission of a job for each iteration), cluster_validate.py (which
    runs validate_once for each iteration), and cluster_validate_combine.py (which 
    runs combine_output all exist. The tag argument allows for adding a string to the job 
    names so that this function can be run multiple times simultaneously and not have 
    conflicts between jobs with the same name. 
    """

    T, J = pl.array(true_cf).shape  
    if os.path.exists(dir) == False: os.mkdir(dir)

    # write true_cf and true_std to file
    data.rec2csv_2d(pl.array(true_cf), '%s/truth_cf.csv' % (dir))
    data.rec2csv_2d(pl.array(true_std), '%s/truth_std.csv' % (dir))
    data.rec2csv_2d(pl.array([std_bias]), '%s/truth_bias.csv' % (dir))
    
    # submit all individual jobs to retrieve true_cf and true_std and run validate_once
    all_names = [] 
    for i in range(reps): 
        name = 'cc%s_%i' % (tag, i)
        call = 'qsub -cwd -N %s cluster_shell.sh cluster_validate.py %i "%s"' % (name, i, dir)
        subprocess.call(call, shell=True)
        all_names.append(name)
    
    # submit job to run combine_output and clean_up 
    hold_string = '-hold_jid %s ' % ','.join(all_names)
    call = 'qsub -cwd %s -N cc%s_comb cluster_shell.sh cluster_validate_combine.py %i "%s"' % (hold_string, tag, reps, dir)
    subprocess.call(call, shell=True)  

 def __init__(self):
   self.ai = ones(NN.ni)
   self.ah = ones(NN.nh)
   self.ao = ones(NN.no)
   self.wi = zeros((NN.ni, NN.nh))
   self.wo = zeros((NN.nh, NN.no))
   randomizeMatrix(self.wi, -0.2, 0.2)
   randomizeMatrix(self.wo, -2.0, 2.0)

 def allocate(self,n):
     """Allocate space for the internal state variables.
     `n` is the maximum sequence length that can be processed."""
     ni,ns,na = self.dims
     vars = "cix ci gix gi gox go gfx gf"
     vars += " state output"
     for v in vars.split():
         setattr(self,v,nan*ones((n,ns)))
     self.source = nan*ones((n,na))

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 getDR(self):
     #this function should return the dynamic range
     #this should be the noiselevel of the fft
     noiselevel=py.sqrt(py.mean(abs(py.fft(self._tdData.getAllPrecNoise()[0]))**2))
     #apply a moving average filter on log
     window_size=5
     window=py.ones(int(window_size))/float(window_size)
     hlog=py.convolve(20*py.log10(self.getFAbs()), window, 'valid')
     one=py.ones((2,))
     hlog=py.concatenate((hlog[0]*one,hlog,hlog[-1]*one))
     return hlog-20*py.log10(noiselevel)         

 def X_obs(pi=pi, sigma=sigma, value=X):
     logp = mc.normal_like(pl.array(value).ravel(), 
                           (pl.ones([N,J*T])*pl.array(pi).ravel()).ravel(), 
                           (pl.ones([N,J*T])*pl.array(sigma).ravel()).ravel()**-2)
     return logp
     
     logp = pl.zeros(N)
     for n in range(N):
         logp[n] = mc.normal_like(pl.array(value[n]).ravel(),
                                  pl.array(pi+beta).ravel(),
                                  pl.array(sigma).ravel()**-2)
     return mc.flib.logsum(logp - pl.log(N))

    def __init__(self, r_floop=0.5, z_floop=0.0,
                 i_p_coil_filename='hitpops.05.txt',
                 tris_filename='hitpops.05.t3d'):

        self.r_floop = r_floop
        self.z_floop = z_floop

        # read equilibrium file
        i_p_coils = P.loadtxt(i_p_coil_filename, delimiter=',', dtype=fdtype)
        self.i_p_coils = i_p_coils

        r_p_coils_full = i_p_coils[:, 0]
        z_p_coils_full = i_p_coils[:, 1]
        # ??? what is this scale factor, something to do with mu_0 ???
        beta = i_p_coils[:, 3] * 6.28e7
        i_p_coils_full = i_p_coils[:, 2]

        self.r_p_coils_full = r_p_coils_full
        self.z_p_coils_full = z_p_coils_full
        self.beta = beta
        self.i_p_coils_full = i_p_coils_full

        # choose subset where current is not zero

        sub = P.where(i_p_coils_full != 0.0)

        r_p_coils = r_p_coils_full[sub]
        z_p_coils = z_p_coils_full[sub]
        i_p_coils = i_p_coils_full[sub]
        n_p_coils = len(r_p_coils)

        self.r_p_coils = r_p_coils
        self.z_p_coils = z_p_coils
        self.i_p_coils = i_p_coils
        self.n_p_coils = n_p_coils

        r_p_widths = P.ones(n_p_coils, dtype=fdtype) * 0.05
        z_p_widths = 1.0 * r_p_widths
        n_r_p_filaments = P.ones(n_p_coils, dtype=idtype)
        n_z_p_filaments = 1 * n_r_p_filaments

        self.r_p_widths = r_p_widths
        self.z_p_widths = z_p_widths
        self.n_r_p_filaments = n_r_p_filaments
        self.n_z_p_filaments = n_z_p_filaments

        # read in triangle unstructured mesh information
        rzt, tris, pt = t3dinp(tris_filename)

        self.rzt = rzt
        self.tris = tris
        self.pt = pt