Python pylab.zeros函数代码示例

def gadget_merge_ics( outfile, filename1, filename2, offset1, offset2, voffset1=[0.,0.,0.], voffset2=[0.,0.,0.] ):
	snap1 = gadget_readsnapname( filename1 )
	snap2 = gadget_readsnapname( filename2 )

	for i in range(3):
		snap1.pos[:,i] += offset1[i]
		snap2.pos[:,i] += offset2[i]
		
	for i in range(3):
		snap1.vel[:,i] += voffset1[i]
		snap2.vel[:,i] += voffset2[i]

	npart = snap1.npart + snap2.npart
	data = {}
	data[ 'count' ] = npart
	data[ 'pos' ] = pylab.zeros( [npart, 3] )
	data[ 'pos' ][ 0:snap1.npart, : ] = snap1.pos
	data[ 'pos' ][ snap1.npart:npart, : ] = snap2.pos
	data[ 'vel' ] = pylab.zeros( [npart, 3] )
	data[ 'vel' ][ 0:snap1.npart, : ] = snap1.vel
	data[ 'vel' ][ snap1.npart:npart, : ] = snap2.vel
	data[ 'mass' ] = pylab.zeros( npart )
	data[ 'mass' ][ 0:snap1.npart ] = snap1.data["mass"]
	data[ 'mass' ][ snap1.npart:npart ] = snap2.data["mass"]
	data[ 'u' ] = pylab.zeros( npart )
	data[ 'u' ][ 0:snap1.npart ] = snap1.data["u"]
	data[ 'u' ][ snap1.npart:npart ] = snap2.data["u"]
	nxnuc = pylab.shape( snap1.data["xnuc"] )[1]
	data[ 'xnuc' ] = pylab.zeros( [npart, nxnuc] )
	data[ 'xnuc' ][ 0:snap1.npart, : ] = snap1.data["xnuc"]
	data[ 'xnuc' ][ snap1.npart:npart, : ] = snap2.data["xnuc"]

	gadget_write_ics( outfile, data, transpose=False )
	return

def findcurve(psi1,psi2,n=3,nn_fit=4,nn_out=100):
    '''
    Function to find the elastica curve for start and end orientations
    psi1 and psi2. It finds the best curve across all directions from start
    and end, i.e. the direction independent elastica curve.
    
    Inputs
    ------------
    psi1,psi2: start and end orientations.
    n:     degree of estimation polynomial.
    nn:    number of points on the curve.
             - nn_fit: for fittin purposes
             - nn_out: for the output
    
    Outputs
    ------------
    Returns a tuple (s,psi). 
    s:   points on the curve.
    psi: curvature of the curve as a function of s.
    E:   curvature energy of the curve
    '''
    # 
    
    # define the starting conditions
    a0 = pl.zeros(n+1) 
    
    # Set a high energy: 
    E_best = 10000  
    
    # and predfine output curve
    s       = pl.linspace(0,1,nn_out) # points on the curve
    psi_out = pl.zeros(nn_out)        # curvature at points in curve
    
    
    # across all the start and end directions find the curve with the lowest energy    
    for dpsi1 in (-pl.pi,0,pl.pi):
        for dpsi2 in (-pl.pi,0,pl.pi):
            # For the starting variables,
            # the first two polygon variables can be estimated from the Sharon paper derivation
            # For different starting variables the solution can be hard to find            
            a0[-2] = 4*(   pl.arcsin(- (pl.sin(psi1+dpsi1)+ pl.sin(psi2+dpsi2))/4)    -(psi1+dpsi1+psi2+dpsi2)/2       )
            a0[-1] = 2*a0[-2]/pl.cos( (psi1+dpsi1+psi2+dpsi2)/2 + a0[-2]/4  )               
            
            # find the best variables to minimize the elastica energy
            fit = fsolve(errors,a0,args=(psi1+dpsi1,psi2+dpsi2,nn_fit))
    
            # find the curve and its derivative for the fitted variables
            a    = fit[:-1]
            psi  = Psi(a,s,psi1+dpsi1,psi2+dpsi2)
            dpsi = dPsi(a,s,psi1+dpsi1,psi2+dpsi2)
    
            # find the energy of this curve
            E = sum(dpsi**2)*s[1]
            
            # check against the lowest energy
            if E_best > E:
                E_best = E
                psi_out[:] = pl.copy(psi)    
    
    return (s,psi_out,E_best)

    def Global_Stiffness(self):
        '''
        Generates Global Stiffness Matrix for the plane structure
        '''
        elem = self.element;
        B = py.zeros((6,6))
        for i in range (0,py.size(elem,0)): 
            #for each element find the stifness matrix
            K = py.zeros((self.n_nodes*2,self.n_nodes*2))            
            el = elem[i]
            
            #nodes formatted for input            
            [node1, node2, node3] = el;
            node1x = 2*(node1-1);node2x = 2*(node2-1);node3x = 2*(node3-1);
            node1y = 2*(node1-1)+1;node2y = 2*(node2-1)+1;node3y = 2*(node3-1)+1;
            
            #Area, Strain Matrix and E Matrix multiplied to get element stiffness            
            [J,B] = self.B(el)
            local_k =0.5*abs(J)*py.dot(py.transpose(B),py.dot(self.E_matrix,B))
            if self.debug:            
                print 'K for elem', el, '\n', local_k
                
            #Element K-Matrix converted into Global K-Matrix format 
            K[py.ix_([node1x,node1y,node2x,node2y,node3x,node3y],[node1x,node1y,node2x,node2y,node3x,node3y])] = K[py.ix_([node1x,node1y,node2x,node2y,node3x,node3y],[node1x,node1y,node2x,node2y,node3x,node3y])]+local_k

            #Adding contibution into Global Stiffness           
            self.k_global = self.k_global + K
            
        if self.debug:            
                print 'Global Stiffness','\n', self.k_global        

    def sample(self, model, evidence):
        g = evidence['g']
        h = evidence['h']
        C = evidence['C']
        z = evidence['z']
        shot_id = evidence['shot_id']
        noise_proportion = evidence['noise_proportion']
        observation_var_g = evidence['observation_var_g']
        observation_var_h = evidence['observation_var_h']

        canopy_cover = model.known_params['canopy_cover']
        z_min = model.known_params['z_min']
        z_max = model.known_params['z_max']

        prior_p = model.hyper_params['T']['p']

        N = len(z)
        T = zeros(N)
        noise_rv = stats.uniform(z_min, z_max - z_min)
        min_index = min(z.index)
        for i in shot_id.index:
            l = zeros(3)
            index = i-min_index
            shot_index = shot_id[i]-min(shot_id)
            l[0] = noise_proportion*noise_rv.pdf(z[i])
            g_norm = stats.norm(g[shot_index], sqrt(observation_var_g))
            C_i = canopy_cover[C[shot_index]]
            l[1] = (1-noise_proportion)*(1-C_i)*g_norm.pdf(z[i])
            h_norm = stats.norm(h[shot_index] + g[shot_index], sqrt(observation_var_h))
            if z[i] > g[shot_index]+3:
                l[2] = (1-noise_proportion)*(C_i)*h_norm.pdf(z[i])
            p = l/sum(l)
            T[index] = Categorical(p).rvs()

        return T

    def sample(self, T, g, g0=None):
        if g0==None:
            g0 = g

        v, h = self.v, self.h
        VH, HH, b_init = self

        V = zeros((T, v))
        H = zeros((T, h))
        B = zeros((T, h))
        
        VH_t = 1*VH

        VH_t[2] = VH[2] + b_init
        V[[0]], H_t_stoch = rbm.sample(VH_t, g0, 1, self.vis_gauss)
        H[[0]] = sigmoid(VH_t * V[[0]])
        if self.vis_gauss:
            V[[0]] = VH_t.T() * H_t_stoch
        else:
            V[[0]] = sigmoid(VH_t.T() * H_t_stoch)
        for t in range(1, T):
            B[[t]] = HH*H[[t-1]]

            VH_t[2] = VH[2] + B[t]
            V[[t]], H_t_stoch = rbm.sample(VH_t, g, 1, self.vis_gauss)

            H[[t]] = sigmoid(VH_t * V[[t]])
            if self.vis_gauss:
                V[[t]] = VH_t.T() * H_t_stoch
            else:
                V[[t]] = sigmoid(VH_t.T() * H_t_stoch)
        return V

def openRomb(integrand, a, b,eps=1e-6,jmax=14,k=5):
    """
    Returns the integral on the _open_interval_ (a,b).
    Integration is performed by Romberg's method of order 2k,
    where, e.g., k=2 is Simpson's rule.
    """
    jmaxp=jmax+1

    s = 0.*M.zeros(jmaxp)
    h = 0.*M.zeros(jmaxp+1)
    ss = 0.
    dss = 0.

    h[0]=1.0
    for j in range(0,jmax):
        s[j]=tripleInt(integrand,a,b,s[j],j)
        if j >= k:
            ss,dss = interpPoly(h[j-k:j],s[j-k:j],k,0.0)
            if M.fabs(dss) <= eps*M.fabs(ss):
                return ss
        s[j+1]=s[j]
        h[j+1]=h[j]/9.

    print 'Non-convergence in openRomb'
    return ss

def mk_image(galaxy):
    base = './../../images_v5/GS_2.5as_matched/gs_all_'

    i_img = pyf.getdata(base+str(galaxy)+'_I.fits')
    j_img = pyf.getdata(base+str(galaxy)+'_J.fits')
    h_img = pyf.getdata(base+str(galaxy)+'_H.fits')

    x = pyl.hstack(i_img)
    i_lim = scoreatpercentile(x,99)
    x = pyl.hstack(j_img)
    j_lim = scoreatpercentile(x,99)
    x = pyl.hstack(h_img)
    h_lim = scoreatpercentile(x,99)

    img = pyl.zeros((h_img.shape[0], h_img.shape[1], 3), dtype=float)
    img[:,:,0] = img_scale.asinh(h_img, scale_max=h_lim, non_linear=0.5)
    img[:,:,1] = img_scale.asinh(j_img, scale_max=j_lim, non_linear=0.5)
    img[:,:,2] = img_scale.asinh(i_img, scale_max=i_lim, non_linear=0.5)

    img = pyl.zeros((h_img.shape[0], h_img.shape[1], 3), dtype=float)
    img[:,:,0] = img_scale.asinh(h_img, scale_min=-0.1*h_lim, scale_max=h_lim,
                                non_linear=0.5)
    img[:,:,1] = img_scale.asinh(j_img, scale_min=-0.1*j_lim, scale_max=j_lim,
                                non_linear=0.5)
    img[:,:,2] = img_scale.asinh(i_img, scale_min=-0.1*i_lim, scale_max=i_lim,
                                non_linear=0.5)

    return img

    def _get_angles(steps,track_length):
        angles = pl.zeros(track_length-2)
        polar = pl.zeros(pl.shape(steps))
        for i in range(track_length-1):
            polar[i,0] = pl.norm(steps[i,:])
            polar[i,1] = pl.arctan(steps[i,0]/steps[i,1])

            if pl.isnan( polar[i,1]):
                polar[i,1] = 0

            if (steps[i,0] >= 0):
                if (steps[i,1] >= 0):
                    pass
                elif (steps[i,1] < 0):
                    polar[i,1] += 2.*pl.pi
            elif (steps[i,0] < 0):
                if (steps[i,1] >= 0):
                    polar[i,1] += pl.pi
                elif (steps[i,1] < 0):
                    polar[i,1] += pl.pi

        for i in range(track_length-2):
            angles[i] = polar[i+1,1] - polar[i,1]

        return angles

def spikefano(timestamps, start_time=0, zero_times=0, end_time=None, window_len=.1, subwindow_len=None):
  """Given the time stamps compute the fano factor with a jumping window.
  Inputs:
    timestamps - the spike timestamps
    window_len - length of window to look at ff (same units as timestamps). One window gets us one ff estimate
                 The fano factor is the LS fit of fano_windows (variance,mean) points
    subwindow_len - length of one spike count computation window

  Outputs:
    t   - time array
    ff  - fano factors
  """
  window_edges, windows, subwindows = window_spike_train(timestamps, start_time, zero_times, end_time, window_len=window_len, subwindow_len=subwindow_len)
  t = pylab.zeros(windows.shape[1])
  ff = pylab.zeros(windows.shape[1])
  for n in xrange(windows.shape[1]):
    spk_count = pylab.zeros(subwindows.shape[1]*subwindows.shape[2])
    for m in xrange(subwindows.shape[0]):#Epochs
      for l in xrange(subwindows.shape[2]):#Subwindows
        #FF computation
        sbw0 = subwindows[m,n,l,0]
        sbw1 = subwindows[m,n,l,1]
        spk_count[m*subwindows.shape[2]+l] = sbw1 - sbw0

    mean = spk_count.mean()
    std = spk_count.std()
    ff[n] = std**2/mean
    t[n] = window_len * (n+.5)

  return t, ff

    def SCpm(SC_0=SC_0, i=i, r=r, f=f, m_all_cause=m_all_cause, age_mesh=dm.get_param_age_mesh()):
        SC = pl.zeros([2, len(age_mesh)])
        p = pl.zeros(len(age_mesh))
        m = pl.zeros(len(age_mesh))

        SC[:, 0] = SC_0
        p[0] = SC_0[1] / (SC_0[0] + SC_0[1])
        m[0] = dismod3.utils.trim(
            m_all_cause[age_mesh[0]] - f[age_mesh[0]] * p[0],
            0.1 * m_all_cause[age_mesh[0]],
            1 - dismod3.settings.NEARLY_ZERO,
        )  # trim m[0] to avoid numerical instability

        for ii, a in enumerate(age_mesh[:-1]):
            A = pl.array([[-i[a] - m[ii], r[a]], [i[a], -r[a] - m[ii] - f[a]]]) * (age_mesh[ii + 1] - age_mesh[ii])

            SC[:, ii + 1] = pl.dot(scipy.linalg.expm(A), SC[:, ii])

            p[ii + 1] = dismod3.utils.trim(
                SC[1, ii + 1] / (SC[0, ii + 1] + SC[1, ii + 1]),
                dismod3.settings.NEARLY_ZERO,
                1 - dismod3.settings.NEARLY_ZERO,
            )
            m[ii + 1] = dismod3.utils.trim(
                m_all_cause[age_mesh[ii + 1]] - f[age_mesh[ii + 1]] * p[ii + 1],
                0.1 * m_all_cause[age_mesh[ii + 1]],
                pl.inf,
            )

        SCpm = pl.zeros([4, len(age_mesh)])
        SCpm[0:2, :] = SC
        SCpm[2, :] = p
        SCpm[3, :] = m
        return SCpm

 def grad(self,data,weightcost):
     grad = zeros(0)
     if type(data)!=type([]):
         data = [data]
     numcases = len(data)
     numscoretypes = len(self.scorefuncs)
     if not type(weightcost) == type([]):
         weightcost = [weightcost] * numscoretypes
     posgrad = [None]*numscoretypes
     neggrad = [None]*numscoretypes
     for k in range(numscoretypes):
         if isscalar(weightcost[k]):
             weightcost[k] = \
                   array([weightcost[k]]*len(self.scorefuncs[k].params))
         posgrad[k] = zeros(self.scorefuncs[k].params.shape,dtype=float)
         neggrad[k] = zeros(self.scorefuncs[k].params.shape,dtype=float)
     for i in range(numcases):
         poscliques = self.posdata(data[i])
         negcliques = self.negdata(data[i])
         for k in range(numscoretypes):
             for posclique in poscliques[k]:
                 posgrad[k] += self.scorefuncs[k].grad(*posclique)
             if self.normalizeacrosscliques:
                 posgrad[k] = posgrad[k]/double(len(poscliques[k]))
             for weighting, negclique in negcliques[k]:
                 for w, neginst in zip(weighting,negclique):
                     neggrad[k] += w * self.scorefuncs[k].grad(*neginst)
             if self.normalizeacrosscliques:
                 neggrad[k] = neggrad[k]/double(len(poscliques[k]))
     for k in range(numscoretypes):
         grad = concatenate((grad,(posgrad[k]-neggrad[k])/double(numcases)\
                                  -weightcost[k]*self.scorefuncs[k].params))
     return -grad

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 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 checkmodelgrad(model,e,RETURNGRADS=False,*args):
    from pylab import norm
    """Check the correctness of passed-in model in terms of cost-/gradient-
       computation, using gradient approximations with perturbances of 
       size e. 
    """
    def updatemodelparams(model, newparams):
        model.params *= 0.0
        model.params += newparams.copy()
    def cost(params,*args):
        paramsold = model.params.copy()
        updatemodelparams(model,params.copy().flatten())
        result = model.cost(*args) 
        updatemodelparams(model,paramsold.copy())
        return result
    def grad(params,*args):
        paramsold = model.params.copy()
        updatemodelparams(model, params.copy().flatten())
        result = model.grad(*args)
        updatemodelparams(model, paramsold.copy())
        return result
    dy = model.grad(*args)
    l = len(model.params)
    dh = zeros(l,dtype=float)
    for j in range(l):
        dx = zeros(l,dtype=float)
        dx[j] = e
        y2 = cost(model.params+dx,*args)
        y1 = cost(model.params-dx,*args)
        dh[j] = (y2 - y1)/(2*e)
    print "analytic: \n", dy
    print "approximation: \n", dh
    if RETURNGRADS: return dy,dh
    else: return norm(dh-dy)/norm(dh+dy)