def __init__(self, path3D, matrix, shape, cap1=0, cap2=0, arrow=0,
                 larrow=3, warrow=2):
        """
        Constructor: Takes as arguments a path3D, matrix, 2D shape, and
        optional cap1, cap2, and arrow.  Calls getextrudeVertices() and
        getFaces() and stores the return values in self.vertices, self.vnormals
        , and self.faces.
        """
        self.path3D = path3D
        self.matrix = matrix
        self.shape = shape
        self.arrow = arrow
        if not isinstance(self.shape, Rectangle2D):
            self.arrow = 0
        self.larrow = larrow
        self.warrow = warrow
        self.cap1 = cap1
        self.cap2 = cap2
        if self.arrow: self.cap2 = 0    # no end cap if arrow
        self.norms = matrix[:,0]
        self.vertices, self.vnormals = self.getextrudeVertices()
        self.faces = self.getFaces()

        # Here need to create a numeric array for each properties:
        # colors
        self.colors = Numeric.ones((len(self.faces), 3),'f')
        # opacities
        self.opacities = Numeric.ones((len(self.faces),),'f')

 def test_constructor_shape(self):
     """__init__         -- make refCoords and resultCoords homogeneous"""
     n = len(self.random_points)
     ncoords = Ncoords( self.random_points) ### tested call ###
     # confirm shape to be nx4
     self.assertEqual( (n, 4), Numeric.shape(ncoords.resultCoords))
     self.assertEqual( (n, 4), Numeric.shape(ncoords.refCoords))
     # cofirm that the last column is all ones
     self.assertEqual(Numeric.ones(n).tolist(),
                      ncoords.resultCoords[:,3].tolist())
     self.assertEqual(Numeric.ones(n).tolist(),
                      ncoords.refCoords[:,3].tolist())

    def doit(self, atom1, atom2, angle, mov_atoms, returnVal=0):
        mol = atom1.top
        if mov_atoms is None:
            mov_atoms = mol.subTree(atom1, atom2, mol.allAtoms)
        assert len(mov_atoms)
        mov_coords = Numeric.array(mov_atoms.coords)
        lenCoords = len(mov_coords)
        x = Numeric.array(atom1.coords)
        y = Numeric.array(atom2.coords)
        rot = (angle * 3.14159/180.)%(2 * Numeric.pi)
        matrix = rotax(x, y, rot)
        _ones = Numeric.ones(lenCoords, 'f')
        _ones.shape = (lenCoords,1)
        mov_coords = Numeric.concatenate((mov_coords, _ones),1)
        newcoords = Numeric.dot(mov_coords, matrix)
        nc = newcoords[:,:3].astype('f')
        for i in range(lenCoords):
            at = mov_atoms[i]
            at._coords[at.conformation] = nc[i].tolist()

        event = EditAtomsEvent('coords', mov_atoms)
        self.vf.dispatchEvent(event)

        #have to return nc for setTorsionGC
        if returnVal:
            return nc

  def __make_rgb_matrices(self, rgb_matrices,shape,normalize,range_=False):
        """ 
        Sub-function of plot() that return the h,s,v matrices
        corresponding to the current matrices in
        sliced_matrices_dict. The shape of the matrices in the dict is
        passed, as well as the normalize boolean parameter.  The
        result specified a bitmap in hsv coordinate.
    
        Applies normalizing and cropping if required.
        """
        zero=zeros(shape,Float)
        one=ones(shape,Float)   

        r,g,b = rgb_matrices
        # Determine appropriate defaults for each matrix
        if r is None: r=zero 
        if g is None: g=zero 
        if b is None: b=zero 

        # CEBALERT: have I checked this works?
        if normalize!='None':
             r = self._normalize(r,range_=range_)
             g = self._normalize(g,range_=range_)
             b = self._normalize(b,range_=range_)

        return (r,g,b)

    def residusMaximus( self, atomValues, mask=None ):
        """
        Take list of value per atom, return list where all atoms of any
        residue are set to the highest value of any atom in that residue.
        (after applying mask)

        @param atomValues: list 1 x N, values per atom
        @type  atomValues: [ float ]
        @param mask: list 1 x N, 0|1, 'master' atoms of each residue
        @type  mask: [1|0]

        @return: Numpy array 1 x N of float
        @rtype: array
        """
        if mask is None:
            mask = N.ones( len( self.frames[0] ), N.int32 )

        ## eliminate all values that do not belong to the selected atoms
        masked = atomValues * mask

        result = []

        ## set all atoms of each residue to uniform value
        for res in range( 0, self.resMap()[-1]+1 ):

            ## get atom entries for this residue
            resAtoms = N.compress( N.equal( self.resMap(), res ), masked )

            ## get maximum value
            masterValue = max( resAtoms )

            result += resAtoms * 0.0 + masterValue

        return N.array( result )

def array2DToImage( array2D, cmap, width=None, 
	height=None, numComponents=4,texture=1, maxi=None, mini=None):
    #build the image:
    if width is None:
        width=array2D.shape[0]
    if height is None:
        height=array2D.shape[1]
    #texture sets the image dimensions to the smallest power of two
    if texture:
        dim1 = dim2=1
        while dim1< width: dim1=dim1<<1
        while dim2< height: dim2=dim2<<1

    colors = Map(array2D.ravel(), cmap,mini=mini, maxi=maxi )

    ###7/19: because these are Numeric arrays???!!!!
    ###7/19colors.shape = (width, height,3)
    colors.shape = (height, width,3)
    colors = colors*255
    colors = colors.astype('B')
    tex2Dimage = Numeric.ones((dim2,dim1,numComponents), 'B')
    ###7/19: because these are Numeric arrays???!!!!
    ###7/19tex2Dimage = Numeric.ones((dim1,dim2,numComponents), 'B')
    ###7/19 tex2Dimage[:width,:height,:3] = colors
    tex2Dimage[:height,:width,:3] = colors
    return tex2Dimage

    def __make_hsv_matrices(self,hsc_matrices,shape,normalize,range_=False):
        """
        Sub-function of plot() that return the h,s,v matrices corresponding
        to the current matrices in sliced_matrices_dict. The shape of the matrices
        in the dict is passed, as well as the normalize boolean parameter.
        The result specified a bitmap in hsv coordinate.

        Applies normalizing and cropping if required.
        """
        zero=zeros(shape,Float)
        one=ones(shape,Float)

        s,h,c = hsc_matrices
        # Determine appropriate defaults for each matrix
        if s is None: s=one # Treat as full strength by default
        if c is None: c=one # Treat as full confidence by default
        if h is None:       # No color, gray-scale plot.
            h=zero
            c=zero

        # If normalizing, offset the matrix so that the minimum
        # value is 0.0 and then scale to make the maximum 1.0
        if normalize!='None':
             s=self._normalize(s,range_=range_)
             # CEBALERT: I meant False, right?
             c=self._normalize(c,range_=False)


        # This translation from SHC to HSV is valid only for black backgrounds;
        # it will need to be extended also to support white backgrounds.
        hue,sat,val=h,c,s
        return (hue,sat,val)

    def group( self, a_indices, maxPerCenter ):
        """
        Group a bunch of integers (atom indices in PDBModel) so that each
        group has at most maxPerCenter items.
        
        @param a_indices: atom indices
        @type  a_indices: [int]
        @param maxPerCenter: max entries per group
        @type  maxPerCenter: int
        
        @return: list of lists of int
        @rtype: [[int],[int]..]
        """
        ## how many groups are necessary?
        n_centers = len( a_indices ) / maxPerCenter
        if len( a_indices ) % maxPerCenter:
            n_centers += 1

        ## how many items/atoms go into each group?
        nAtoms = N.ones(n_centers, N.Int) * int(len( a_indices ) / n_centers)
        i=0
        while N.sum(nAtoms) != len( a_indices ):
            nAtoms[i] += 1
            i += 1

        ## distribute atom indices into groups
        result = []
        pos = 0
        for n in nAtoms:
            result += [ N.take( a_indices, N.arange(n) + pos) ]
            pos += n

        return result

    def doit(self,mobAtoms):
        """
mobAtoms: AtomSet that is being frozen.
Assuming the AtomSet are from same molecule
        """

        # fixme: need checking for mat (matrix) and mobAtoms
        geomContainer = mobAtoms[0].top.geomContainer
        mGeom = geomContainer.masterGeom
        mat = mGeom.rotation
        mat = Numeric.reshape(mat, (4,4))

        # update coords
        mobAtoms = mobAtoms.findType(Atom)
        coords = mobAtoms.coords
        hCoords = Numeric.concatenate((coords,Numeric.ones((len(coords),1),\
                                                            'd')), 1)
        tCoords = Numeric.dot(hCoords, mat)[:,:3]
        tCoords = tCoords.tolist()
        mobAtoms.updateCoords(tCoords, 0) # overwritten the original coords
        
        # reset the rotation matrix of masterGeom
        identity = Numeric.identity(4,'f').ravel()
        mGeom.SetRotation(Numeric.identity(4,'f').ravel() ) 

        event = EditAtomsEvent('coords', mobAtoms)
        self.vf.dispatchEvent(event)
        
        mGeom.viewer.Redraw()
        
        return

 def __init__(self, cntrl, uknots, vknots):
     self._bezier = None
     cntrl = numerix.asarray(cntrl, numerix.Float)
     (dim, nu, nv) = cntrl.shape
     if dim < 2 or dim > 4:
         raise NURBSError, 'Illegal control point format'
     elif dim < 4:
         self.cntrl = numerix.zeros((4, nu, nv), numerix.Float)
         self.cntrl[0:dim,:,:] = cntrl
         self.cntrl[-1,:,:] = numerix.ones((nu,nv), numerix.Float)
     else:
         self.cntrl = cntrl
         
     # Force the u knot sequence to be a vector in ascending order
     # and normalise between [0.0,1.0]
     uknots = numerix.sort(numerix.asarray(uknots, numerix.Float))
     nku = uknots.shape[0]
     uknots = (uknots - uknots[0])/(uknots[-1] - uknots[0])
     if uknots[0] == uknots[-1]:
         raise NURBSError, 'Illegal uknots sequence'
     self.uknots = uknots
     
     # Force the v knot sequence to be a vector in ascending order
     # and normalise between [0.0,1.0]  
     vknots = -numerix.sort(-numerix.asarray(vknots, numerix.Float))
     nkv = vknots.shape[0]
     vknots = (vknots - vknots[0])/(vknots[-1] - vknots[0])
     if vknots[0] == vknots[-1]:
         raise NURBSError, 'Illegal vknots sequence'
     self.vknots = vknots
     
     # Spline Degree
     self.degree = [nku-nu-1, nkv-nv-1]
     if self.degree[0] < 0 or self.degree[1] < 0:
         raise NURBSError, 'NURBS order must be a positive integer'

    def castHmmDic( self, hmmDic, repete, hmmGap, key ):
        """
        Blow up hmmDic to the number of repetes of the profile used.
        Correct scores for possible deletions in the search sequence.

        @param hmmDic: dictionary from L{getHmmProfile}
        @type  hmmDic: dict
        @param repete: repete information from L{align}
        @type  repete: int
        @param hmmGap: information about gaps from L{align}
        @type  hmmGap: [int]
        @param key: name of scoring method to adjust for gaps and repetes
        @type  key: str
        
        @return: dictionary with information about the profile
        @rtype: dict        
        """
        s = hmmDic[key]

        for i in range( repete ):
            mask = N.ones( len(s) )
            N.put( mask, hmmGap[i], 0 )
            if i == 0:
                score = N.compress( mask, s, 0 )
            if i > 0:
                score = N.concatenate( ( N.compress( mask, s, 0 ), score ) )

        hmmDic[key] = score

        return hmmDic

    def Reset(self):
        if __debug__:
         if hasattr(DejaVu, 'functionName'): DejaVu.functionName()
	"""Reset members to default values"""

	self.pattern = Numeric.ones( (128,)) * 0xFF
	self.pattern = self.pattern.astype('B')

    def addDensity( self, radius=6, minasa=None, profName='density' ):
        """
        Count the number of heavy atoms within the given radius.
        Values are only collected for atoms with |minasa| accessible surface
        area.

        @param minasa: relative exposed surface - 0 to 100%
        @type  minasa: float
        @param radius: in Angstrom
        @type  radius: float
        """
        mHeavy = self.m.maskHeavy()

        xyz = N.compress( mHeavy, self.m.getXyz(), 0 )

        if minasa and self.m.profile( 'relAS', 0 ) == 0:
            self.addASA()

        if minasa:
            mSurf = self.m.profile2mask( 'relAS', minasa )
        else:
            mSurf = N.ones( self.m.lenAtoms() )

        ## loop over all surface atoms
        surf_pos = N.nonzero( mSurf )
        contacts = []

        for i in surf_pos:
            dist = N.sum(( xyz - self.m.xyz[i])**2, 1)
            contacts += [ N.sum( N.less(dist, radius**2 )) -1]

        self.m.atoms.set( profName, contacts, mSurf, default=-1,
                          comment='atom density radius %3.1fA' % radius,
                          version= T.dateString() + ' ' + self.version() )

def Map(values, colorMap, mini=None, maxi=None):
    """Get colors corresponding to values in a colormap"""

    values = Numeric.array(values)
    if len(values.shape)==2 and values.shape[1]==1:
	values.shape = ( values.shape[0], )
    elif len(values.shape) > 1:
	print 'ERROR: values array has bad shape'
	return None

    cmap = Numeric.array(colorMap)
    if len(cmap.shape) !=2 or cmap.shape[1] not in (3,4):
	print 'ERROR: colorMap array has bad shape'
	return None

    if mini is None: mini = min(values)
    else: values = Numeric.maximum(values, mini)
    if maxi is None: maxi = max(values)
    else: values = Numeric.minimum(values, maxi)
    valrange = maxi-mini
    if valrange < 0.0001:
	ind = Numeric.ones( values.shape )
    else:
	colrange = cmap.shape[0]-1
	ind = ((values-mini) * colrange) / valrange
    col = Numeric.take(colorMap, ind.astype(viewerConst.IPRECISION))
    return col

    def shuffledLists( self, n, lst, mask=None ):
        """
        shuffle order of a list n times, leaving masked(0) elements untouched

        @param n: number of times to shuffle the list
        @type  n: int
        @param lst: list to shuffle
        @type  lst: [any]
        @param mask: mask to be applied to lst
        @type  mask: [1|0]

        @return: list of shuffeled lists
        @rtype: [[any]]        
        """
        if not mask:
            mask = N.ones( len(lst)  )

        if type( lst ) == list:
            lst = N.array( lst )
        
        pos = N.nonzero( mask )

        rand_pos = N.array( [ self.__shuffleList( pos ) for i in range(n) ] )

        result = []
        for p in rand_pos:

            r = copy.copy( lst )
            N.put( r, p, N.take( lst, pos ) )
            result += [r]

        return result

    def cluster( self, n_clusters, weight=1.13, converged=1e-11,
                 aMask=None, force=0 ):
        """
        Calculate new clusters.

        @param n_clusters: number of clusters
        @type  n_clusters: int
        @param weight: fuzziness weigth
        @type  weight: float (default: 1.13)
        @param converged: stop iteration if min dist changes less than
                          converged (default: 1e-11)
        @type  converged: float
        @param aMask: atom mask applied before clustering
        @type  aMask: [1|0]
        @param force: re-calculate even if parameters haven't changed
                      (default:0)
        @type  force: 1|0
        """
        if aMask == None:
            aMask = N.ones( self.traj.getRef().lenAtoms() )

        if self.fc == None or force or self.fcWeight != weight \
           or self.n_clusters != n_clusters or N.any( self.aMask != aMask) \
           or self.fcConverged != converged:

            self.n_clusters = n_clusters
            self.fcWeight = weight
            self.aMask = aMask

            self.fc = FuzzyCluster( self.__raveled(), self.n_clusters,
                                    self.fcWeight )

            self.fcCenters = self.fc.go( self.fcConverged,
                                         1000, nstep=10,
                                         verbose=self.verbose )

    def smooth(self):
        average = (Numeric.ones( (3,3) , 'f')/9.).astype('f')
#        average[1][1] = .0
        print average
        self.redraw(filter=average)
        image = (self.imageAsArray()).astype('B')
        return image

    def SetMaterial(self, values, prop=1, tagModified=True):
	"""Set the materials
WARNING: when back face colors are set, two sided lighting has to be enabled
        
we set RGB values for all properties except for shininess and opacity
If an alpha value are specified they will be ignored
Since IndexedGeomDSPL requires alpha values for all properties
we set them automatically to 1.0 for all properties except for
diffuse. The alpha channel of the diffuse component will be set later
in fixOpacity which should be called after all properties of the
material have been updated.
"""
        if tagModified:
            self._modified = True

        if prop is None: prop = self.diff
        elif type(prop) is types.StringType:
            prop = getattr(self, prop)

        assert prop in (0,1,2,3,4,5)

        values = array( values, 'f' )
	if prop < self.shini:
            assert len(values.shape)==2 and values.shape[1] in [3,4]
            alpha = ones( (values.shape[0], 1), 'f' )
            values = concatenate( (values[:,:3], alpha), 1 )
        else:
            if len(values.shape) != 1:
                values = array([values], 'f' )
            
        self.prop[prop] = values

def triangularPut(m1d, upper=1, lower=0):
    """Returns 2D masked array with elements of the given 1D array in the strictly upper (lower) triangle.
    Elements of the 1D array should be ordered according to the upper triangular part of the 2D matrix.
    The lower triangular part (if requested) equals to the transposed upper triangular part.
    If upper == lower == 1 a symetric matrix is returned.
    """
    assert upper in [0,1] and lower in [0,1], "[0|1] expected for upper / lower"
    m1d = MA.asarray(m1d)
    assert MA.rank(m1d) == 1, "1D masked array expected"
    m2dShape0 = math.ceil(math.sqrt(2*m1d.shape[0]))
    assert m1d.shape[0] == m2dShape0*(m2dShape0-1)/2, "the length of m1d does not correspond to n(n-1)/2"
    if upper:
        if lower:
            mask = Numeric.fromfunction(lambda i,j: i==j, (m2dShape0, m2dShape0))
        else:
            mask = Numeric.fromfunction(lambda i,j: i>=j, (m2dShape0, m2dShape0))
    else:
        if lower:
            mask = Numeric.fromfunction(lambda i,j: i<=j, (m2dShape0, m2dShape0))
        else:
            mask = Numeric.ones((m2dShape0, m2dShape0))

    m2d = MA.ravel(MA.zeros((m2dShape0, m2dShape0), m1d.dtype.char))
    condUpperTriang = Numeric.fromfunction(lambda i,j: i<j, (m2dShape0, m2dShape0))
    putIndices = Numeric.compress(Numeric.ravel(condUpperTriang), Numeric.arange(0, m2dShape0**2, typecode=Numeric.Int))
    MA.put(m2d, putIndices, m1d)
    m2d = MA.reshape(m2d, (m2dShape0, m2dShape0))
    m2d = MA.where(condUpperTriang, m2d, MA.transpose(m2d))
    return MA.array(m2d, mask=Numeric.logical_or(mask, MA.getmaskarray(m2d)))

def scale(sxyz):
    ret = numerix.identity(4).astype(numerix.Float)
    s = numerix.ones(3, numerix.Float)
    s[0:len(sxyz)] = sxyz
    ret[0,0] = s[0]
    ret[1,1] = s[1]
    ret[2,2] = s[2]
    return ret