def __init__(self,**params):
        from plastk.rand import uniform
        from Numeric import zeros
        super(GNG,self).__init__(**params)

        N = self.initial_num_units
        
        self.weights = uniform(self.rmin,self.rmax,(N,self.dim))
        self.dists = zeros((N,1)) * 0.0
        self.error = zeros((N,1)) * 0.0

        self.connections = [{} for i in range(N)]
        
        self.last_input = zeros(self.dim)
        
        self.count = 0

        if self.initial_connections_per_unit > 0:
            for w in self.weights:
                self.present_input(w)
                ww = self.winners(self.initial_connections_per_unit+1)
                i = ww[0]
                for j in ww[1:]:
                    self.add_connection(i,j)

        self.nopickle += ['_activation_fn']
        self.unpickle()

    def test_vanishing_moments(self):
        """Test that coefficients in lp satisfy the
           vanishing moments condition
        """ 

        from daubfilt import daubfilt, number_of_filters

        for i in range(number_of_filters):
            D = 2*(i+1)

            P = D/2  # Number of vanishing moments
            N = P-1  # Dimension of nullspace of the matrix A
            R = P+1  # Rank of A, R = D-N = P+1 equations
        
            lp, hp = daubfilt(D)


            # Condition number of A grows with P, so we test only
            # the first 6 (and eps is slightly larger than machine precision)

            A    = zeros((R,D), Float)  # D unknowns, D-N equations
            b    = zeros((R,1), Float)  # Right hand side
            b[0] = sqrt(2)                
  
            A[0,:] = ones(D, Float)   # Coefficients must sum to sqrt(2)
            for p in range(min(P,6)): # the p'th vanishing moment (Cond Ap)
                for k in range(D):            
                    m=D-k;
                    A[p+1,k] = (-1)**m * k**p;

            assert allclose(b, mvmul(A,lp))         

    def analyze(self,observations):
        """use Viterbi algorithm to
        find the states corresponding to the observations"""
        B = self.B
        A = self.A
        T = len(observations)
        N = self.N
        Omega_X = self.omega_X
        obs = self._getObservationIndices(observations)
        # initialisation
        delta = []
        delta.append(B[obs[0]] * self.pi) # (32a)
        phi = [array([0]*N)]              # (32b)
        # recursion
        for O_t in obs[1:]:
            delta_t = zeros(N,Float)
            phi_t = zeros(N)
            for j in range(N):
                delta_t[j] = max(delta[-1]*A[:,j]*B[O_t][j]) # (33a)
                phi_t[j] = argmax(delta[-1]*A[:,j])          # (33b)
            delta.append(delta_t)
            phi.append(phi_t)

        # reconstruction
        i_star = [argmax(delta[-1])]                         # (34b)
        phi.reverse() # we start from the end
        for phi_t in phi[:-1]:
            i_star.append(phi_t[i_star[-1]])                 # (35)
        trajectory = [Omega_X[i] for i in i_star]
        trajectory.reverse() # put time back in the right direction
        return trajectory

    def test_pairs_to_array(self):
        """pairs_to_array should match hand-calculated results"""
        p2a = pairs_to_array
        p1 = [0, 1, 0.5]
        p2 = [2, 3, 0.9]
        p3 = [1, 2, 0.6]
        pairs = [p1, p2, p3]
        self.assertEqual(p2a(pairs), \
            array([[0,.5,0,0],[0,0,.6,0],[0,0,0,.9],[0,0,0,0]]))
        #try it without weights -- should assign 1
        new_pairs = [[0,1],[2,3],[1,2]]
        self.assertEqual(p2a(new_pairs), \
            array([[0,1,0,0],[0,0,1,0],[0,0,0,1],[0,0,0,0]]))
        #try it with explicit array size
        self.assertEqual(p2a(pairs, 5), \
            array([[0,.5,0,0,0],[0,0,.6,0,0],[0,0,0,.9,0],[0,0,0,0,0],\
            [0,0,0,0,0]]))
        #try it when we want to map the indices into gapped coords
        #we're effectively doing ABCD -> -A--BC-D-
        transform = array([1,4,5,7])
        result = p2a(pairs, transform=transform)
        self.assertEqual(result.shape, (8,8))
        exp = zeros((8,8), Float64)
        exp[1,4] = 0.5
        exp[4,5] = 0.6
        exp[5,7] = 0.9
        self.assertEqual(result, exp)

        result = p2a(pairs, num_items=9, transform=transform)
        self.assertEqual(result.shape, (9,9))
        exp = zeros((9,9), Float64)
        exp[1,4] = 0.5
        exp[4,5] = 0.6
        exp[5,7] = 0.9
        self.assertEqual(result, exp)

def pacbporflist2similarityarray(pacbps,queryorsbjct,length):
    """ """
    bea = zeros(length)
    for pacbporf in pacbps:
        spos  = pacbporf._get_original_alignment_pos_start()
        epos  = pacbporf._get_original_alignment_pos_end()
        q,m,s = pacbporf.get_unextended_aligned_protein_sequences()
        if queryorsbjct == 'query':
            start = spos.query_pos
            end   = epos.query_pos + 1
            seqa  = list(q)
            ma    = list(m)
        else:
            start = spos.sbjct_pos
            end   = epos.sbjct_pos + 1
            seqa  = list(s)
            ma    = list(m)

        for pos in range(len(seqa)-1,-1,-1):
            if seqa[pos] == '-':
                seqa.pop(pos)
                ma.pop(pos)
        # prepare replacement of match string into match score list
        matcharray = zeros(end-start)
        for pos in range(0,len(ma)):
            symbol = ma[pos]
            if symbol != ' ':
                matcharray[pos] = 1
        # update (binary) array
        bea[spos.query_pos:epos.query_pos+1] += matcharray

    # correct bea for values > 1
    bea = where(greater_equal(bea, 2), 1, bea)
    return bea

    def __init__(self, win):
        self.coords = zeros([30,3], Float)
        self.prev_coords = zeros([3,3], Float)

        self.peptide_mol = None
        self.length = 0
        self.prev_psi = 0

        PeptideGeneratorPropertyManager.__init__(self)
        GeneratorBaseClass.__init__(self, win)

 def accept(self):
   
   m = self.getPickledMesh()
   
   # Create some Numeric arrays for the mesh (much faster)
   triangleNumArray = m.getTriangulation()
   verticeNumArray = m.getMeshVertices()
   #QMessageBox.information(None, "DEBUG", 'Tri 0 has verts ' + str(self.triangleNumArray[0]) )
   #QMessageBox.information(None, "DEBUG", 'Vert 942 is at ' + str(self.verticeNumArray[942]) )
   triCnt = len(triangleNumArray)
   verCnt = len(verticeNumArray)
   triangleNumArray = array( triangleNumArray, Numeric.Int32 )
   #verticeNumArray = array( verticeNumArray, Numeric.Int32 )
   
   # Blank list of neighbors
   neighbors         = zeros( (triCnt, 3), Numeric.Int32 )
   linkedNeighborCnt = zeros( triCnt, Numeric.Int32 )
   
   self.progressBar.setMinimum(0)
   self.progressBar.setMaximum(triCnt)
   
   startTime = time.time()
   for T in xrange(triCnt):
     self.progressBar.setValue(T)
     if linkedNeighborCnt[T] < 3:
       me = triangleNumArray[T]
       for t in xrange(triCnt):
         if linkedNeighborCnt[t] < 3 and t != T:
           them = triangleNumArray[t]
           commonVerts = 0
           for i in xrange(3):
             if them[i] == me[0] or them[i] == me[1] or them[i] == me[2]:
               commonVerts += 1
           if commonVerts == 2:
             add = True
             for i in xrange(linkedNeighborCnt[T]):
               if neighbors[T][i] == t:
                 add = False
             if add:
               neighbors[T][linkedNeighborCnt[T]] = t
               linkedNeighborCnt[T] += 1
               """except:
                 QMessageBox.information(None, 'DEBUG', 'T is ' + str(T) )
                 QMessageBox.information(None, 'DEBUG', 'its neighbors are ' + str(neighbors[T]) )
                 QMessageBox.information(None, 'DEBUG', 'we tried to add ' + str(t) )"""
               neighbors[t][linkedNeighborCnt[t]] = T
               linkedNeighborCnt[t] += 1
   
   QMessageBox.information(None, 'DEBUG', 'Done in ' + str(time.time() - startTime) + ' seconds' )
   
   """badTriAreaCnt = 0

    def __init__(self, shp, ptlist, origin, selSense, **opts):
        """
        ptlist is a list of 3d points describing a selection.
        origin is the center of view, and normal gives the direction
        of the line of light. Form a structure for telling whether
        arbitrary points fall inside the curve from the point of view.
        """
        # bruce 041214 rewrote some of this method
        simple_shape_2d.__init__(self, shp, ptlist, origin, selSense, opts)

        # bounding rectangle, in integers (scaled 8 to the angstrom)
        ibbhi = array(map(int, ceil(8 * self.bboxhi) + 2))
        ibblo = array(map(int, floor(8 * self.bboxlo) - 2))
        bboxlo = self.bboxlo

        # draw the curve in these matrices and fill it
        # [bruce 041214 adds this comment: this might be correct but it's very
        # inefficient -- we should do it geometrically someday. #e]
        mat = zeros(ibbhi - ibblo)
        mat1 = zeros(ibbhi - ibblo)
        mat1[0, :] = 1
        mat1[-1, :] = 1
        mat1[:, 0] = 1
        mat1[:, -1] = 1
        pt2d = self.pt2d
        pt0 = pt2d[0]
        for pt in pt2d[1:]:
            l = ceil(vlen(pt - pt0) * 8)
            if l < 0.01:
                continue
            v = (pt - pt0) / l
            for i in range(1 + int(l)):
                ij = 2 + array(map(int, floor((pt0 + v * i - bboxlo) * 8)))
                mat[ij] = 1
            pt0 = pt
        mat1 += mat

        fill(mat1, array([1, 1]), 1)
        mat1 -= mat  # Which means boundary line is counted as inside the shape.
        # boolean raster of filled-in shape
        self.matrix = mat1  ## For any element inside the matrix, if it is 0, then it's inside.
        # where matrix[0, 0] is in x, y space
        self.matbase = ibblo

        # axes of the plane; only used for debugging
        self.x = self.right
        self.y = self.up
        self.z = self.normal

def getSegyTraceHeader(SH,THN='cdp',data='none'):
        """
        getSegyTraceHeader(SH,TraceHeaderName)
        """

        bps=getBytePerSample(SH)

        if (data=='none'):
                data = open(SH["filename"]).read()
                

        # MAKE SOME LOOKUP TABLE THAT HOLDS THE LOCATION OF HEADERS
#       THpos=TraceHeaderPos[THN]
        THpos=STH_def[THN]["pos"]
        THformat=STH_def[THN]["type"]
        ntraces=SH["ntraces"]
        thv = zeros(ntraces)
        for itrace in range(1,ntraces+1,1):
                i=itrace

                pos=THpos+3600+(SH["ns"]*bps+240)*(itrace-1);

                txt="getSegyTraceHeader : Reading trace header " + THN + " " + str(itrace)  + " of " + str(ntraces) + " " +str(pos)

                printverbose(txt,20);
                thv[itrace-1],index = getValue(data,pos,THformat,endian,1)
                txt="getSegyTraceHeader : " + THN + "=" + str(thv[itrace-1])
                printverbose(txt,30);
        

def writemb(index, data, dxsize, dysize, bands, mb_db):
    """
    Write raster 'data' (of the size 'dataxsize' x 'dataysize') read from
    'dataset' into the mbtiles document 'mb_db' with size 'tilesize' pixels.
    Later this should be replaced by new <TMS Tile Raster Driver> from GDAL.
    """
    if bands == 3 and tileformat == 'png':
        tmp = tempdriver.Create('', tilesize, tilesize, bands=4)
        alpha = tmp.GetRasterBand(4)
        #from Numeric import zeros
        alphaarray = (zeros((dysize, dxsize)) + 255).astype('b')
        alpha.WriteArray( alphaarray, 0, tilesize-dysize )
    else:
        tmp = tempdriver.Create('', tilesize, tilesize, bands=bands)

    tmp.WriteRaster( 0, tilesize-dysize, dxsize, dysize, data, band_list=range(1, bands+1))
    tiledriver.CreateCopy('tmp.png', tmp, strict=0)
    query = """insert into tiles 
        (zoom_level, tile_column, tile_row, tile_data) 
        values (%d, %d, %d, ?)""" % (index[0], index[1], index[2])
    cur = mb_db.cursor()
    d = open('tmp.png', 'rb').read()
    cur.execute(query, (sqlite3.Binary(d),))
    mb_db.commit()
    cur.close()
    return 0

 def insert_n_rows(self, i, n=1):
     " Insert `n` rows into each column at row `i`. "
     rows = list()
     for tc in self.typecodes:
         rows.append( zeros((n,),typecode=tc) )
     self.insert_rows(i, rows)
     self.update_rows()

    def epsilon_greedy(self,sensation,applicable_actions):
        """
        Given self.epsilon() and self.Q(), return a distribution over
        applicable_actions as an array where each element contains the
        a probability mass for the corresponding action.  I.e.  The
        action with the highest Q gets p = self.epsilon() and the
        others get the remainder of the mass, uniformly distributed.
        """
        Q = array([self.Q(sensation,action) for action in applicable_actions])

        # simple epsilon-greedy policy
        # get a vector with a 1 where each max element is, zero elsewhere
        mask = (Q == mmax(Q))

        num_maxes = len(nonzero(mask))
        num_others = len(mask) - num_maxes

        if num_others == 0: return mask
        
        e0 = self.epsilon()/num_maxes
        e1 = self.epsilon()/num_others

        result = zeros(len(mask))+0.0
        putmask(result,mask,1-e0)
        putmask(result,mask==0,e1)
        return result

 def absoluteProfile(alignment,char_order):
     f = a.columnFrequencies()
     res = zeros([len(f),len(char_order)])
     for row, freq in enumerate(f):
         for i in freq:
             res[row, i] = freq[i]
     return res

def pos_char_weights(alignment, order=DNA_ORDER):
    """Returns the contribution of each character at each position.

    alignment: Alignemnt object
    order: the order of characters in the profile (all observed chars
        in the alignment
    
    This function is used by the function position_based
    
    For example: 
    GYVGS
    GFDGF
    GYDGF
    GYQGG
    
        0       1       2       3       4       5   
    G   1/1*4                           1/1*4   1/3*1
    Y           1/2*3
    F           1/2*1                           1/3*2
    V                   1/3*1
    D                   1/3*2
    Q                   1/3*1
    S                                           1/3*1
    """
    counts = alignment.columnFrequencies()
    a = zeros([len(order), alignment.SeqLen],Float64)
    for col, c in enumerate(counts):
        for char in c:
            a[order.index(char),col] = 1/(len(c)*c[char])
    return Profile(a,Alphabet=order)

 def __call__(self,input,error=False):
     X = self.scale_input(input)
     if not error:
         Y = dot(self.w,X)
     else:
         Y = dot(self.w,X), zeros(self.num_outputs)
     return self.scale_output(Y)

def inertia_eigenvectors(basepos, already_centered = False):
    """
    Given basepos (an array of positions),
    compute and return (as a 2-tuple) the lists of eigenvalues and
    eigenvectors of the inertia tensor (computed as if all points had the same
    mass). These lists are always length 3, even for len(basepos) of 0,1, or 2,
    overlapping or colinear points, etc, though some evals will be 0 in these cases.
       Optional small speedup: if caller knows basepos is centered at the origin, it can say so.
    """
    #bruce 060119 split this out of shakedown_poly_evals_evecs_axis() in chunk.py
    basepos = A(basepos) # make sure it's a Numeric array
    if not already_centered and len(basepos):
        center = add.reduce(basepos)/len(basepos)
        basepos = basepos - center
    # compute inertia tensor
    tensor = zeros((3,3),Float)
    for p in basepos:
        rsq = dot(p, p)
        m= - multiply.outer(p, p)
        m[0,0] += rsq
        m[1,1] += rsq
        m[2,2] += rsq
        tensor += m
    evals, evecs = eigenvectors(tensor)
    assert len(evals) == len(evecs) == 3
    return evals, evecs

def fitsin(x,y,tau=1,delta=0,A=1,C=1,sig=1):
    """Fit data in x and y to a sin function"""
    data=zeros((len(x),3),typecode='d')
    data[:,2]=sig
    data[:,0]=x
    data[:,1]=y
    return LS(sinus,(tau,delta,A,C),data)

    def setIds(self, id_fun=lambda x: x.Data.split("_")[-1]):
        """
        Sets "LeafLabel", "LeafCts", and "ContainsAll" attributes

        id_fun: function that takes node and generate a unique id (label)
            for each node. By default will create a label consisting of 
            the string to the right of the last underscore in the data
            attribute. E.g. if the node has data label of 1234_HSA, the
            function will return a unique lable of "HSA". the idea being
            that if your tree has multiple human (HSA) sequences, the
            result of the function will be multiple nodes w/the same
            label. 

        The LeafLabel attribute is the the result of the id_fun function.

        The LeafCts attribute is an array with counts of the leaves with the 
            same label.

        The ContainsAll attribute is True when it contains every instance 
            of the LeafLabels of its terminal descendants. E.g. the set
            of LeafLabels of its terminal descendants occur nowhere else
            in the tree. 

        This is used by the uniqueIds function to remove duplicate species
        from the tree but can be used for any label you choose.
        """
        labels =  [id_fun(x)  for x in self.TerminalDescendants]
        u_labels = list(set(labels))
        len_u_labels = len(u_labels)
        labels_dict =  dict(zip(u_labels, range(len_u_labels)))
        all_cts = zeros(len(u_labels))

        for label in labels: 
            all_cts[labels_dict[label]] += 1
      
        for n in self.traverse(self_before=False, self_after=True):
            if not n.Children:
                setattr(n, "LeafLabel", id_fun(n))
                setattr(n, "LeafCts", zeros(len_u_labels))
                n.LeafCts[labels_dict[n.LeafLabel]] = 1
            else:
                n.LeafCts = zeros(len_u_labels)
                for c in n.Children:
                    n.LeafCts += c.LeafCts 
            nzero = nonzero(n.LeafCts)
            total = sum(take(all_cts, nzero)- take(n.LeafCts, nzero))
            setattr(n, "ContainsAll", (total == 0))

 def __call__(self,sensation,reward=None):
     if is_terminal(sensation):
         new_sensation = sensation
     else:
         new_sensation = zeros(self.num_features,'f')
         for f in sensation:
             new_sensation[f] = 1
     return super(LinearListAgent,self).__call__(new_sensation,reward)

def table_to_array(tbl, typecode='f'):
    shape = (tbl.ncols, tbl.nrows)
    
    a = zeros( (tbl.ncols, tbl.nrows), typecode)
    for j in range(tbl.ncols):
        a[j] = tbl[j].astype(typecode)

    return transpose(a)