def __call__(self,action=None):
        if action==None:
            self.angle = self.initial_angle
            self.velocity = self.initial_velocity
            return array((self.angle,self.velocity))
        else:
            if type(action) == ArrayType:
                action = action[0]
            self.last_action = action
            # setup constants
            m = self.mass
            g = self.g
            l = self.length
            dt = self.delta_t
            
            # Gravitational torque
            T_g = m * g * l * cos(pi/2-self.angle)
            # action torque
            T_a = action

            # rotational acceleration is Torque divided by moment-of-inertia
            a = (T_g + T_a)/(m*l)

            # update velocity from acceleration (including friction)
            self.velocity += (a * dt)
            self.velocity -= self.velocity * self.friction

            # update angle from velocity
            self.angle += (self.velocity * dt)/2

            return (array((normalize_angle(self.angle),self.velocity)),self.reward())

 def setData(self, *args):
     if len(args) == 3:
         self.locns, self.vals, self.sigs = args
         if len(locns) != len(vals) or len(vals) != len(sigs):
             raise DataError, 'Mismatch in lengths of arguments!'
     elif len(args) == 1:
         all = args[0]
         try:
             rows, cols = all.shape
             self.locns = all[:,0]
             self.vals = all[:,1]
             self.sigs = all[:,2]
         except AttributeError:
             # If *all* is not an array, treat it as a list.
             # Only the list format can handle >1-d locns using
             # a single setData argument.
             self.locns, self.vals, self.sigs = [], [], []
             for row in all:
                 self.locns.append(row[0])
                 self.vals.append(row[1])
                 self.sigs.append(row[2])
             self.locns = array(self.locns)
             self.vals = array(self.vals)
             self.sigs = array(self.sigs)
         except:
             raise DataError, 'Bad data array format!'
     self.ndata = len(vals)
     if len(self.locns.shape) == 1:
         self.dimen = 1
     elif len(self.locns.shape) ==2:
         self.dimen = self.locns.shape[1]
     else:
         raise DataError, 'Locations should be scalars or 1-d arrays!'

 def test__add_(self):
     """__add__: should not normalize input or output, just add"""
     p1 = Profile(array([[.3,.4,.1,0],[.1,.1,.1,.7]]),Alphabet="ABCD")
     p2 = Profile(array([[1,0,0,0],[1,0,0,1]]),Alphabet="ABCD")
     self.assertEqual((p1+p2).Data, array([[1.3,.4,.1,0],[1.1,.1,.1,1.7]]))
     self.assertRaises(ProfileError,self.empty.__add__, p1)
     self.assertEqual((self.empty + self.empty).Data.tolist(),[[]])

 def test_toOddsMatrix(self):
     """toOddsMatrix: should work on valid data or raise an error
     """
     p = Profile(array([[.1,.3,.5,.1],[.25,.25,.25,.25],\
         [.05,.8,.05,.1],[.7,.1,.1,.1],[.6,.15,.05,.2]]),\
         Alphabet="ACTG")
     p_exp = Profile(array([[.4, 1.2, 2, .4],[1,1,1,1],[.2,3.2,.2,.4],\
         [2.8,.4,.4,.4],[2.4,.6,.2,.8]]),Alphabet="ACTG")
     self.assertEqual(p.toOddsMatrix().Data,p_exp.Data)
     assert p.Alphabet is p.toOddsMatrix().Alphabet
     self.assertEqual(p.toOddsMatrix([.25,.25,.25,.25]).Data,p_exp.Data)
     
     #fails if symbol_freqs has wrong size
     self.assertRaises(ProfileError, p.toOddsMatrix,\
         [.25,.25,.25,.25,.25,.25])
     self.assertRaises(ProfileError, self.zero_entry.toOddsMatrix,\
         [.1,.2,.3])
     #works on empty profile
     self.assertEqual(self.empty.toOddsMatrix().Data.tolist(),[[]])
     #works with different input
     self.assertEqual(self.zero_entry.toOddsMatrix().Data,\
         array([[1.2,.8,0,2],[0,0,3.2,.8]]))
     self.assertFloatEqual(self.zero_entry.toOddsMatrix([.1,.2,.3,.4]).Data,\
         array([[3,1,0,1.25],[0,0,2.667,.5]]),1e-3)
     #fails when one of the background frequencies is 0
     self.assertRaises(ProfileError, self.zero_entry.toOddsMatrix,\
         [.1,.2,.3,0])

def test4():
    test = HMM(['a','b'],['s1','s2','s3'],
               array([[.3,.5],[.7,.5]]),
               array([[.5,0],[.5,.5],[0,.5]]),
               array([.9,.1]))
    print test.simulate(10)
    print test.simulate(10,1)

    def test_score_profile(self):
        """score: should work correctly for Profile as input
        """
        p1 = Profile(array([[1,0,0,0],[0,1,0,0],[0,0,.5,.5],[0,0,0,1],\
            [.25,.25,.25,.25]]),"TCAG")
        p2 = Profile(array([[0,1,0,0],[.2,0,.8,0],[0,0,.5,.5],[1/3,1/3,0,1/3],\
            [.25,.25,.25,.25]]),"TCAG")
        p3 = Profile(array([[1,0,0,0],[0,1,0,0],[0,0,0,1]]),"TCAG")
        p4 = Profile(array([[1,0,0,0],[0,1,0,0]]),"TCAG")
        p5 = Profile(array([[1,0,0,0],[0,1,0,0],[0,0,0,1]]),"AGTC")

        #works on normal valid data
        self.assertFloatEqual(self.score2.score(p1,offset=0),\
            [.55,1.25,.45])
        self.assertFloatEqual(self.score2.score(p2,offset=0),
            [1.49,1.043,.483],1e-3)
        #works with different offset
        self.assertFloatEqual(self.score2.score(p1,offset=1),
            [1.25,0.45])
        self.assertFloatEqual(self.score2.score(p1,offset=2),
            [0.45])
        #raises error on invalid offset 
        self.assertRaises(ProfileError,self.score2.score,\
            p1,offset=3)
        #works on profile of minimal length
        self.assertFloatEqual(self.score2.score(p3,offset=0),
            [0.6])
        #raises error when profile is too short
        self.assertRaises(ProfileError, self.score2.score,p4,offset=0)
        #raises error on empty profile
        self.assertRaises(ProfileError,self.empty.score,p1)
        #raises error when character order doesn't match
        self.assertRaises(ProfileError,self.score2.score,p5) 

 def __init__(self, *args, **kwargs):
     Visualiser.__init__(self, *args, **kwargs)
     fin = NetCDFFile(self.vis_source, 'r')
     self.xPoints = array(fin.variables['x'], Float)
     self.yPoints = array(fin.variables['y'], Float)
     self.quantityCache = {}
     fin.close()

    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 distance_to_closest(alignment, distance_method=hamming_distance):
    """Returns vector of distances to closest neighbor for each s in alignment
    
    alignment: Alignment object
    distance_method: function used to calculate the distance between two seqs
    
    Function returns the closest distances according to the RowOrder in the
    alignment

    example:
    Alignment({1:'ABCD',2:'ABCC',3:'CBDD',4:'ACAA'},RowOrder=[3,2,1,4])
    [2,1,1,3]
    """
    #change sequences into arrays
    for item in alignment:
        alignment[item] = array(alignment[item])
    
    closest = []
    for key in alignment.RowOrder:
        seq = alignment[key]
        dist = None
        for other_key in alignment.RowOrder:
            if key == other_key:
                continue
            d = distance_method(seq,alignment[other_key])
            if dist:
                if d < dist:
                    dist = d
            else:
                dist = d
        closest.append(dist)
    return array(closest)

 def test_ungapped_to_gapped(self):
     """ungapped_to_gapped should match hand-calculated results"""
     u2g = ungapped_to_gapped
     self.assertEqual(u2g(self.s1, '-'), array([0,1,2,4]))
     self.assertEqual(u2g(self.s2, '-'), array([2,3]))
     self.assertEqual(u2g(self.s3, '-'), array([0,1]))
     self.assertEqual(u2g(self.s4, '-'), array([0,1]))
     self.assertEqual(u2g(self.s5, '-'), array([]))

 def test_unmasked_to_masked(self):
     """unmasked_to_masked should match hand-calculated results"""
     u2m = unmasked_to_masked
     self.assertEqual(u2m(self.m1), array([0,1,2,4]))
     self.assertEqual(u2m(self.m2), array([2,3]))
     self.assertEqual(u2m(self.m3), array([0,1]))
     self.assertEqual(u2m(self.m4), array([0,1]))
     self.assertEqual(u2m(self.m5), array([]))

    def distance(self,other):
        """Calculates the distance between two BaseUsages.

        Distance is measured in three directions, CG-content, CU-content, and
        GU-content.
        """
        return euclidean_distance(array(self.toCartesian()),\
            array(other.toCartesian()))

 def test_center_of_mass_two_array(self):
     """center_of_mass_two_array should behave correctly"""
     com2 = center_of_mass_two_array
     coor = take(self.square_odd,(0,1),1)
     weights = take(self.square_odd,(2,),1)
     self.assertEqual(com2(coor,weights), array([2,2]))
     weights = weights.flat
     self.assertEqual(com2(coor,weights), array([2,2]))

 def setUp(self):
     """setUp for all CenterOfMass tests"""
     self.simple = array([[1,1,1],[3,1,1],[2,3,2]])
     self.simple_list = [[1,1,1],[3,1,1],[2,3,2]]
     self.more_weight = array([[1,1,3],[3,1,3],[2,3,50]])
     self.square = array([[1,1,25],[3,1,25],[3,3,25],[1,3,25]])
     self.square_odd = array([[1,1,25],[3,1,4],[3,3,25],[1,3,4]])
     self.sec_weight = array([[1,25,1],[3,25,1],[3,25,3],[1,25,3]])

 def __init__(self, stacks, position):
     """
     Pass a reference to the stacks object that holds all the blocks
     and the index of the initial stack of the block
     """
     Block.__init__(self, stacks, position)
     self.box = box(pos=(position-(stacks.blockCount/2), 0, 0), size=(.9,.9,.9), color=color.blue)
     self.label = label(pos=array(self.box.pos) + array([0,0,1]), text=str(position), opacity=0, box=0, line=0)

 def testAngleBetween(self):
     vector1 = array((1, 0, 0))
     vector2 = array((0, 1, 0))
     angle = angleBetween(vector1, vector2)
     assert angle == 90, "Fails sanity check"
     assert alltrue(vector1 == array((1, 0, 0))) and \
            alltrue(vector2 == array((0, 1, 0))), \
            "Arguments were modified (recurrence of bug ####)"

 def optimize(self):
     """Optimize the objective function with respect to varying params."""
     p, d = [], []
     for param in self.varyingParams:
         p.append(param.getValue())
         d.append(param.delta)
     p = array(p)
     d = identity(len(d)) * array(d,Float)
     return powell(self._objective, p, drxns=d)

def hermite(x0, x1, v0, v1, t):
    a = array([1, t, t**2, t**3])
    M = array([[1, 0, 0, 0],
               [1, 1, 1, 1],
               [0, 1, 0, 0],
               [0, 1, 2, 3]])
    p = transpose([[x0, x1, v0, v1]])

    return dot(dot(a, inverse(M)), p)

 def test_euclidean_distance(self):
     """euclidean_distance: should return dist between 2 vectors or matrices
     """
     a = array([3,4])
     b = array([8,5])
     c = array([[2,3],[4,5]])
     d = array([[1,5],[8,2]])
     self.assertFloatEqual(euclidean_distance(a,b),sqrt(26))
     self.assertFloatEqual(euclidean_distance(c,d),sqrt(30))

 def __init__(self, x, y=None, z=None):
     if y is None and z is None:
         # Array, list, tuple...
         if len(x)!=3:
             raise "Vector: x is not a list/tuple/array of 3 numbers"
         self._ar=array(x, 'd')
     else:
         # Three numbers
         self._ar=array((x, y, z), 'd')