def correlation(x_items, y_items):
    """Returns Pearson correlation between x and y, and its significance."""
    sum_x = sum_y = sum_x_sq = sum_y_sq = sum_xy = n = 0
    for x, y in zip(x_items, y_items):
        n += 1
        sum_x += x
        sum_x_sq += x * x
        sum_y += y
        sum_y_sq += y * y
        sum_xy += x * y
    try:
        r = 1.0 * ((n * sum_xy) - (sum_x * sum_y)) / \
           (sqrt((n * sum_x_sq)-(sum_x*sum_x))*sqrt((n*sum_y_sq)-(sum_y*sum_y)))
    except (ZeroDivisionError, ValueError): #no variation
        r = 0.0
    #check we didn't get a naughty value for r due to rounding error
    if r > 1.0:
        r = 1.0
    elif r < -1.0:
        r = -1.0
    if n < 3:
        prob = 1
    else:
        try:
            t = r/sqrt((1 - (r*r))/(n-2))
            prob = tprob(t, n-2)
        except ZeroDivisionError: #r was presumably 1
            prob = 0
    return (r, prob)

def txt_sqrt(norm, numeric=False):
    if numeric:
        return repr(sqrt(norm))
    else:
        if sqrt(norm) % 1 == 0:
            return str(sqrt(norm))
        else:
            return "sqrt(" + str(norm.nom) + ("./" + str(norm.denom)) * (norm.denom != 1) + ")"

 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 txt_sqrt(norm, numeric=False):
    if numeric:
        return repr(sqrt(norm))
    else:
        if sqrt(norm) % 1 == 0:
            return str(sqrt(norm))
        else:
            return 'sqrt(' + str(norm.nom) + \
                   ('./' + str(norm.denom)) * (norm.denom != 1) + ')'

    def _getinfo_TEST(self): # please leave in for debugging POV-Ray lmotor macro. mark 060324
        a = self.axen()
        xrot = -atan2(a[1], sqrt(1-a[1]*a[1]))*180/pi
        yrot = atan2(a[0], sqrt(1-a[0]*a[0]))*180/pi

        return  "[Object: Linear Motor] [Name: " + str(self.name) + "] " + \
                "[Force = " + str(self.force) + " pN] " + \
                "[Stiffness = " + str(self.stiffness) + " N/m] " + \
                "[Axis = " + str(self.axis[0]) + ", " +  str(self.axis[1]) + ", " +  str(self.axis[2]) + "]" + \
                "[xRotation = " + str(xrot) + ", yRotation = " + str(yrot) + "]"

def t_two_sample (a, b, tails=None, exp_diff=0):
    """Returns t, prob for two INDEPENDENT samples of scores a, and b.  
    
    From Sokal and Rohlf, p 223.  

    Usage:   t, prob = t_two_sample(a,b, tails, exp_diff)

    t is a float; prob is a probability.
    a and b should be lists of observations (numbers) supporting Mean, Count,
    and Variance. Need not be equal.
    tails should be None (default), 'high', or 'low'.
    exp_diff should be the expected difference in means (a-b); 0 by default.
 """
    try:
        #see if we need to back off to the single-observation for single-item
        #groups
        n1 = a.Count
        if n1 < 2:
            return t_one_observation(a.Sum, b, tails, exp_diff)
        n2 = b.Count
        if n2 < 2:
            return t_one_observation(b.Sum, a, reverse_tails(tails), exp_diff)
        #otherwise, calculate things properly
        x1 = a.Mean
        x2 = b.Mean
        df = n1+n2-2
        svar = ((n1-1)*a.Variance + (n2-1)*b.Variance)/df
        t = (x1-x2-exp_diff)/sqrt(svar*(1/n1 + 1/n2))
    except (ZeroDivisionError, ValueError, AttributeError, TypeError):
        #bail out if the sample sizes are wrong, the values aren't numeric or
        #aren't present, etc.
        return (None, None)
    prob = t_tailed_prob(t, df, tails)
    return t, prob

    def Length(self):
	"""Length of vector

	Returns the length of the vector generated by the basis.
	"""
	from Numeric import sqrt
	return sqrt(self.InnerProduct(self))

 def findHandles_exact(self, p1, p2, cutoff = 0.0, backs_ok = 1, offset = V(0,0,0)):
     """
     return a list of (dist, handle) pairs, in arbitrary order,
     which includes, for each handle (spherical surface) hit by the ray from p1 thru p2,
     its front-surface intersection with the ray,
     unless that has dist < cutoff and backs_ok,
     in which case include its back-surface intersection
     (unless *that* has dist < cutoff).
     """
     #e For now, just be simple, don't worry about speed.
     # Someday we can preprocess self.handlpos using Numeric functions,
     # like in nearSinglets and/or findSinglets
     # (I have untested prototype code for this in extrude-outs.py).
     hh = self.handles
     res = []
     v = norm(p2-p1)
     ## is this modifying the vector in-place, causing a bug?? offset += self.origin # treat our handles' pos as relative to this
     ## I don't know, but one of the three instances of += was doing this!!! probably i was resetting the atom or mol pos....
     offset = offset + self.origin # treat our handles' pos as relative to this
     radius_multiplier = self.radius_multiplier
     for (pos,radius,info) in hh:
         ## bug in this? pos += offset
         pos = pos + offset
         radius *= radius_multiplier
         dist, wid = orthodist(p1, v, pos)
         if radius >= wid: # the ray hits the sphere
             delta = sqrt(radius*radius - wid*wid)
             front = dist - delta # depth from p1 of front surface of sphere, where it's hit
             if front >= cutoff:
                 res.append((front,(pos,radius,info)))
             elif backs_ok:
                 back = dist + delta
                 if back >= cutoff:
                     res.append((back,(pos,radius,info)))
     return res

    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 test_randomSequence(self):
        """randomSequence: 99% of new frequencies should be within 3*SD"""
        r_num, c_num = 100,20
        num_elements = r_num*c_num
        alpha = "ABCDEFGHIJKLMNOPQRSTUVWXYZ"
        r = random([r_num,c_num])
        p = Profile(r,alpha[:c_num])
        p.normalizePositions()
        d = p.Data
        n = 1000
        
        #Test only works on normalized profile, b/c of 1-d below
        means = n*d
        three_stds = sqrt(d*(1-d)*n)*3

        a = Alignment([p.randomSequence() for x in range(n)])

        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:
                    col = char_order.index(i)
                    res[row, col] = freq[i]
            return res

        ap = absoluteProfile(a,p.CharOrder)
        failure = abs(ap-means) > three_stds
        assert sum(sum(failure))/num_elements <= 0.01

def init_diamond():
    # a chunk of diamond grid, to be tiled out in 3d
    drawing_globals.sp0 = sp0 = 0.0
    #bruce 051102 replaced 1.52 with this constant (1.544),
    #  re bug 900 (partial fix.)
    drawing_globals.sp1 = sp1 = DIAMOND_BOND_LENGTH / sqrt(3.0)
    sp2 = 2.0*sp1
    sp3 = 3.0*sp1
    drawing_globals.sp4 = sp4 = 4.0*sp1

    digrid=[[[sp0, sp0, sp0], [sp1, sp1, sp1]],
            [[sp1, sp1, sp1], [sp2, sp2, sp0]],
            [[sp2, sp2, sp0], [sp3, sp3, sp1]],
            [[sp3, sp3, sp1], [sp4, sp4, sp0]],
            [[sp2, sp0, sp2], [sp3, sp1, sp3]],
            [[sp3, sp1, sp3], [sp4, sp2, sp2]],
            [[sp2, sp0, sp2], [sp1, sp1, sp1]],
            [[sp1, sp1, sp1], [sp0, sp2, sp2]],
            [[sp0, sp2, sp2], [sp1, sp3, sp3]],
            [[sp1, sp3, sp3], [sp2, sp4, sp2]],
            [[sp2, sp4, sp2], [sp3, sp3, sp1]],
            [[sp3, sp3, sp1], [sp4, sp2, sp2]],
            [[sp4, sp0, sp4], [sp3, sp1, sp3]],
            [[sp3, sp1, sp3], [sp2, sp2, sp4]],
            [[sp2, sp2, sp4], [sp1, sp3, sp3]],
            [[sp1, sp3, sp3], [sp0, sp4, sp4]]]
    drawing_globals.digrid = A(digrid)
    drawing_globals.DiGridSp = sp4
    return

    def test_euclidean_distance_unexpected(self):
        """euclidean_distance: works always when frames are aligned. UNEXPECTED!
        """
        a = array([3,4])
        b = array([8,5])
        c = array([[2,3],[4,5]])
        d = array([[1,5],[8,2]])
        e = array([[4,5],[4,5],[4,5]])
        f = array([1,1,1,1,1])
        self.assertFloatEqual(euclidean_distance(a,c),sqrt(4))
        self.assertFloatEqual(euclidean_distance(c,a),sqrt(4))
        self.assertFloatEqual(euclidean_distance(a,e),sqrt(6))

        #IT DOES RAISE AN ERROR WHEN THE FRAMES ARE NOT ALIGNED
        self.assertRaises(ValueError,euclidean_distance,c,e)
        self.assertRaises(ValueError,euclidean_distance,c,f)

 def __str__(self):
     n = self.norm
     sn = sqrt(n)
     if int(sn) == sn:
         string = repr(sn) + "/sqrt(pi)"
     else:
         string = "sqrt(" + repr(n.nom) + ("./" + repr(n.denom)) * (n.denom != 1) + "/pi)"
     return string

def min_dist(coord, surface):
    """
    Return minimum distance between coord
    and surface.
    """
    d=surface-coord
    d2=sum(d*d, 1)
    return sqrt(min(d2))

 def __str__(self):
     n = self.norm
     sn = sqrt(n)
     if int(sn) == sn:
         string = repr(sn) + '/sqrt(pi)'
     else:
         string = 'sqrt(' + repr(n.nom) + \
                  ('./' + repr(n.denom)) * (n.denom != 1) + '/pi)'
     return string

    def computeEndPointsFromChunk(self, chunk, update=True):
        """
        Derives and returns the endpoints and radius of a Peptide chunk.
        @param chunk: a Peptide chunk
        @type  chunk: Chunk
        @return: endPoint1, endPoint2 and radius
        @rtype: Point, Point and float

        @note: computing the endpoints works fine when n=m or m=0. Otherwise,
               the endpoints can be slightly off the central axis, especially
               if the Peptide is short.
        @attention: endPoint1 and endPoint2 may not be the original endpoints,
                    and they may be flipped (opposites of) the original
                    endpoints.
        """
        # Since chunk.axis is not always one of the vectors chunk.evecs
        # (actually chunk.poly_evals_evecs_axis[2]), it's best to just use
        # the axis and center, then recompute a bounding cylinder.
        if not chunk.atoms:
            return None

        axis = chunk.axis
        axis = norm(axis)  # needed
        center = chunk._get_center()
        points = chunk.atpos - center  # not sure if basepos points are already centered
        # compare following Numeric Python code to findAtomUnderMouse and its caller
        matrix = matrix_putting_axis_at_z(axis)
        v = dot(points, matrix)
        # compute xy distances-squared between axis line and atom centers
        r_xy_2 = v[:, 0] ** 2 + v[:, 1] ** 2

        # to get radius, take maximum -- not sure if max(r_xy_2) would use Numeric code, but this will for sure:
        i = argmax(r_xy_2)
        max_xy_2 = r_xy_2[i]
        radius = sqrt(max_xy_2)
        # to get limits along axis (since we won't assume center is centered between them), use min/max z:
        z = v[:, 2]
        min_z = z[argmin(z)]
        max_z = z[argmax(z)]

        # Adjust the endpoints such that the ladder rungs (rings) will fall
        # on the ring segments.
        # TO DO: Fix drawPeptideLadder() to offset the first ring, then I can
        # remove this adjustment. --Mark 2008-04-12
        z_adjust = self.getEndPointZOffset()
        min_z += z_adjust
        max_z -= z_adjust

        endpoint1 = center + min_z * axis
        endpoint2 = center + max_z * axis

        if update:
            # print "Original endpoints:", self.getEndPoints()
            self.setEndPoints(endpoint1, endpoint2)
            # print "New endpoints:", self.getEndPoints()

        return (endpoint1, endpoint2, radius)

def norm(a):
    """Returns the norm of a matrix or vector

    Calculates the Euclidean norm of a vector.
    Applies the Frobenius norm function to a matrix 
    (a.k.a. Euclidian matrix norm)

    a = Numeric array
    """
    return sqrt(sum((a*a).flat))

def m2rotaxis(m):
    """
    Return angles, axis pair that corresponds to rotation matrix m.
    """
    # Angle always between 0 and pi
    # Sense of rotation is defined by axis orientation
    t=0.5*(trace(m)-1)
    t=max(-1, t)
    t=min(1, t)
    angle=acos(t)
    if angle<1e-15:
        # Angle is 0
        return 0.0, Vector(1,0,0)
    elif angle<pi:
        # Angle is smaller than pi
        x=m[2,1]-m[1,2]
        y=m[0,2]-m[2,0]
        z=m[1,0]-m[0,1]
        axis=Vector(x,y,z)
        axis.normalize()
        return angle, axis
    else:
        # Angle is pi - special case!
        m00=m[0,0]
        m11=m[1,1]
        m22=m[2,2]
        if m00>m11 and m00>m22:
            x=sqrt(m00-m11-m22+0.5)
            y=m[0,1]/(2*x)
            z=m[0,2]/(2*x)
        elif m11>m00 and m11>m22:
            y=sqrt(m11-m00-m22+0.5)
            x=m[0,1]/(2*y)
            z=m[1,2]/(2*y)
        else:
            z=sqrt(m22-m00-m11+0.5)
            x=m[0,2]/(2*z)
            y=m[1,2]/(2*z)
        axis=Vector(x,y,z)
        axis.normalize()
        return pi, axis

    def __sub__(self, other):
        """
        Calculate distance between two atoms.
        
        Example:
            >>> distance=atom1-atom2

        @param other: the other atom
        @type other: L{Atom}
        """
        diff=self.coord-other.coord
        return sqrt(sum(diff*diff))

    def writepov(self, file, dispdef):
        if self.hidden: return
        if self.is_disabled(): return #bruce 050421
        c = self.posn()
        a = self.axen()

        xrot = -atan2(a[1], sqrt(1-a[1]*a[1]))*180/pi
        yrot =  atan2(a[0], sqrt(1-a[0]*a[0]))*180/pi

        file.write("lmotor(" \
                   + povpoint([self.width *  0.5, self.width *  0.5, self.length *  0.5]) + "," \
                   + povpoint([self.width * -0.5, self.width * -0.5, self.length * -0.5]) + "," \
                   + "<0.0, " + str(yrot) + ", 0.0>," \
                   + "<" + str(xrot) + ", 0.0, 0.0>," \
                   + povpoint(c) + "," \
                   + "<" + str(self.color[0]) + "," + str(self.color[1]) + "," + str(self.color[2]) + ">)\n")

        for a in self.atoms:
            if vlen(c - a.posn()) > 0.001: #bruce 060808 add condition to see if this fixes bug 719 (two places in this file)
                file.write("spoke(" + povpoint(c) + "," + povpoint(a.posn())  + "," + str (self.sradius) +
                           ",<" + str(self.color[0]) + "," + str(self.color[1]) + "," + str(self.color[2]) + ">)\n")