def estimateGaussian(trainData, trainLabels):
    numClasses = max(trainLabels) + 1
    N = [0]* numClasses
    mu = [0.0] * numClasses
    Slist = [zeros((trainData.shape[1], trainData.shape[1]))] * numClasses
    pList = [0.0] * numClasses
    
    #calculate N, and sum x's for mu
    for x,t in izip(trainData, trainLabels):
        N[t] += 1
        mu[t] += x
    #normalize mu
    for i in range(numClasses):
        mu[i] = mu[i] / float(N[i])
    #calculate the class probabilities
    for i in range(numClasses):
        pList[i] = float(N[i]) / sum(N)
    #calculate S0 and S1
    for x,t in izip(trainData, trainLabels):
        Slist[t] += outer(x - mu[t], x - mu[t])
        try:
            inv(Slist[t])
        except LinAlgError:
            Slist[t] += 0.1 * identity(Slist[t].shape[0], Float64)
        
    return (numClasses, N, mu, Slist, pList)

def kron(a,b):
    """kronecker product of a and b

    Kronecker product of two arrays is block array
    [[ a[ 0 ,0]*b, a[ 0 ,1]*b, ... , a[ 0 ,n-1]*b  ],
     [ ...                                   ...   ],
     [ a[m-1,0]*b, a[m-1,1]*b, ... , a[m-1,n-1]*b  ]]
    """
    wrapper = get_array_wrap(a, b)
    b = asanyarray(b)
    a = array(a,copy=False,subok=True,ndmin=b.ndim)
    ndb, nda = b.ndim, a.ndim
    if (nda == 0 or ndb == 0):
        return _nx.multiply(a,b)
    as_ = a.shape
    bs = b.shape
    if not a.flags.contiguous:
        a = reshape(a, as_)
    if not b.flags.contiguous:
        b = reshape(b, bs)
    nd = ndb
    if (ndb != nda):
        if (ndb > nda):
            as_ = (1,)*(ndb-nda) + as_
        else:
            bs = (1,)*(nda-ndb) + bs
            nd = nda
    result = outer(a,b).reshape(as_+bs)
    axis = nd-1
    for _ in xrange(nd):
        result = concatenate(result, axis=axis)
    if wrapper is not None:
        result = wrapper(result)
    return result

def fishersLinearDiscriminent(trainData, trainLabels, testData, testLabels):
    numClasses = max(trainLabels) + 1
    N = [0] * numClasses
    m = [0] * numClasses
    for x,t in izip(trainData,trainLabels):
        m[t] += x
        N[t] += 1
    for i in range(numClasses):
        m[i] /= N[i]
    Sw = zeros((trainData.shape[1], trainData.shape[1]))
    for x,t in izip(trainData, trainLabels):
        Sw += outer(x-m[t], x-m[t])
    try:
        inv(Sw)
    except LinAlgError:
        Sw += 0.1 * identity(Sw.shape[0], Float64)    
    
    w = dot(inv(Sw),(m[0] - m[1]))
    meanVect = (N[0]*m[0] + N[1]*m[1]) / sum(N)
    
    numCorrect = 0
    for x,t in izip(testData, testLabels):
        if dot(w, (x-meanVect)) > 0:
            if t == 1:
                numCorrect += 1
        else:
            if t == 0:
                numCorrect += 1
    return float(numCorrect) / float(len(testLabels))

    def __init__ (self, worldfile):

        with open (worldfile) as world:
            reader = csv.reader (world)
            names = {}
            for k, line in enumerate (reader):
                if not "#" in line [0]:
                    names [line [0]] = k

        data = loadtxt (worldfile, delimiter = ',', usecols = range (1, 7), dtype = float64)
        self.names = names
        self.radii = data [:, 0]
        self.masses = data [:, 1]
        self.positions = data [:, 2:4]
        self.velocities = data [:, 4:]
        self.accelerations = zeros ( (len (self.names),) + self.positions.shape, dtype = float64)

        self.mm = outer (self.masses, self.masses)
        self.diagind = tuple (range (0, len (self.accelerations)))

        self.count = len (self.names)
        self.time = 0

def kron(a,b):
    """
    Kronecker product of two arrays.

    Computes the Kronecker product, a composite array made of blocks of the
    second array scaled by the first.

    Parameters
    ----------
    a, b : array_like

    Returns
    -------
    out : ndarray

    See Also
    --------

    outer : The outer product

    Notes
    -----

    The function assumes that the number of dimenensions of `a` and `b`
    are the same, if necessary prepending the smallest with ones.
    If `a.shape = (r0,r1,..,rN)` and `b.shape = (s0,s1,...,sN)`,
    the Kronecker product has shape `(r0*s0, r1*s1, ..., rN*SN)`.
    The elements are products of elements from `a` and `b`, organized
    explicitly by::

        kron(a,b)[k0,k1,...,kN] = a[i0,i1,...,iN] * b[j0,j1,...,jN]

    where::

        kt = it * st + jt,  t = 0,...,N

    In the common 2-D case (N=1), the block structure can be visualized::

        [[ a[0,0]*b,   a[0,1]*b,  ... , a[0,-1]*b  ],
         [  ...                              ...   ],
         [ a[-1,0]*b,  a[-1,1]*b, ... , a[-1,-1]*b ]]


    Examples
    --------
    >>> np.kron([1,10,100], [5,6,7])
    array([  5,   6,   7,  50,  60,  70, 500, 600, 700])
    >>> np.kron([5,6,7], [1,10,100])
    array([  5,  50, 500,   6,  60, 600,   7,  70, 700])

    >>> np.kron(np.eye(2), np.ones((2,2)))
    array([[ 1.,  1.,  0.,  0.],
           [ 1.,  1.,  0.,  0.],
           [ 0.,  0.,  1.,  1.],
           [ 0.,  0.,  1.,  1.]])

    >>> a = np.arange(100).reshape((2,5,2,5))
    >>> b = np.arange(24).reshape((2,3,4))
    >>> c = np.kron(a,b)
    >>> c.shape
    (2, 10, 6, 20)
    >>> I = (1,3,0,2)
    >>> J = (0,2,1)
    >>> J1 = (0,) + J             # extend to ndim=4
    >>> S1 = (1,) + b.shape
    >>> K = tuple(np.array(I) * np.array(S1) + np.array(J1))
    >>> c[K] == a[I]*b[J]
    True

    """
    b = asanyarray(b)
    a = array(a,copy=False,subok=True,ndmin=b.ndim)
    ndb, nda = b.ndim, a.ndim
    if (nda == 0 or ndb == 0):
        return _nx.multiply(a,b)
    as_ = a.shape
    bs = b.shape
    if not a.flags.contiguous:
        a = reshape(a, as_)
    if not b.flags.contiguous:
        b = reshape(b, bs)
    nd = ndb
    if (ndb != nda):
        if (ndb > nda):
            as_ = (1,)*(ndb-nda) + as_
        else:
            bs = (1,)*(nda-ndb) + bs
            nd = nda
    result = outer(a,b).reshape(as_+bs)
    axis = nd-1
    for _ in range(nd):
        result = concatenate(result, axis=axis)
    wrapper = get_array_prepare(a, b)
    if wrapper is not None:
        result = wrapper(result)
    wrapper = get_array_wrap(a, b)
    if wrapper is not None:
        result = wrapper(result)
    return result

#.........这里部分代码省略.........
    >>> z = np.polyfit(x, y, 3)
    >>> z
    array([ 0.08703704, -0.81349206,  1.69312169, -0.03968254])

    It is convenient to use `poly1d` objects for dealing with polynomials:

    >>> p = np.poly1d(z)
    >>> p(0.5)
    0.6143849206349179
    >>> p(3.5)
    -0.34732142857143039
    >>> p(10)
    22.579365079365115

    High-order polynomials may oscillate wildly:

    >>> p30 = np.poly1d(np.polyfit(x, y, 30))
    /... RankWarning: Polyfit may be poorly conditioned...
    >>> p30(4)
    -0.80000000000000204
    >>> p30(5)
    -0.99999999999999445
    >>> p30(4.5)
    -0.10547061179440398

    Illustration:

    >>> import matplotlib.pyplot as plt
    >>> xp = np.linspace(-2, 6, 100)
    >>> plt.plot(x, y, '.', xp, p(xp), '-', xp, p30(xp), '--')
    [<matplotlib.lines.Line2D object at 0x...>, <matplotlib.lines.Line2D object at 0x...>, <matplotlib.lines.Line2D object at 0x...>]
    >>> plt.ylim(-2,2)
    (-2, 2)
    >>> plt.show()

    """
    order = int(deg) + 1
    x = NX.asarray(x) + 0.0
    y = NX.asarray(y) + 0.0

    # check arguments.
    if deg < 0 :
        raise ValueError("expected deg >= 0")
    if x.ndim != 1:
        raise TypeError("expected 1D vector for x")
    if x.size == 0:
        raise TypeError("expected non-empty vector for x")
    if y.ndim < 1 or y.ndim > 2 :
        raise TypeError("expected 1D or 2D array for y")
    if x.shape[0] != y.shape[0] :
        raise TypeError("expected x and y to have same length")

    # set rcond
    if rcond is None :
        rcond = len(x)*finfo(x.dtype).eps

    # set up least squares equation for powers of x
    lhs = vander(x, order)
    rhs = y

    # apply weighting
    if w is not None:
        w = NX.asarray(w) + 0.0
        if w.ndim != 1:
            raise TypeError("expected a 1-d array for weights")
        if w.shape[0] != y.shape[0] :
            raise TypeError("expected w and y to have the same length")
        lhs *= w[:, NX.newaxis]
        if rhs.ndim == 2:
            rhs *= w[:, NX.newaxis]
        else:
            rhs *= w

    # scale lhs to improve condition number and solve
    scale = NX.sqrt((lhs*lhs).sum(axis=0))
    lhs /= scale
    c, resids, rank, s = lstsq(lhs, rhs, rcond)
    c = (c.T/scale).T  # broadcast scale coefficients

    # warn on rank reduction, which indicates an ill conditioned matrix
    if rank != order and not full:
        msg = "Polyfit may be poorly conditioned"
        warnings.warn(msg, RankWarning)

    if full :
        return c, resids, rank, s, rcond
    elif cov :
        Vbase = inv(dot(lhs.T,lhs))
        Vbase /= NX.outer(scale, scale)
        # Some literature ignores the extra -2.0 factor in the denominator, but
        #  it is included here because the covariance of Multivariate Student-T
        #  (which is implied by a Bayesian uncertainty analysis) includes it.
        #  Plus, it gives a slightly more conservative estimate of uncertainty.
        fac = resids / (len(x) - order - 2.0)
        if y.ndim == 1:
            return c, Vbase * fac
        else:
            return c, Vbase[:,:,NX.newaxis] * fac
    else :
        return c

 def project_ts(self, bv):  # TODO: test
     base, vel = bv
     P = np.eye(self.N) - outer(base, base)
     return base, np.dot(P, vel)

def polyfit(x, y, deg, rcond=None, full=False, w=None, cov=False):

    import numpy.core.numeric as NX
    from numpy.core import isscalar, abs, dot
    from numpy.lib.twodim_base import diag, vander
    from numpy.linalg import eigvals, lstsq, inv
    try:
        from numpy.core import finfo # 1.7
    except:
        from numpy.lib.getlimits import finfo # 1.3 support for cluster

    order = int(deg) + 1
    x = NX.asarray(x) + 0.0
    y = NX.asarray(y) + 0.0

    # check arguments.
    if deg < 0 :
        raise ValueError("expected deg >= 0")
    if x.ndim != 1:
        raise TypeError("expected 1D vector for x")
    if x.size == 0:
        raise TypeError("expected non-empty vector for x")
    if y.ndim < 1 or y.ndim > 2 :
        raise TypeError("expected 1D or 2D array for y")
    if x.shape[0] != y.shape[0] :
        raise TypeError("expected x and y to have same length")

    # set rcond
    if rcond is None :
        rcond = len(x)*finfo(x.dtype).eps

    # set up least squares equation for powers of x
    lhs = vander(x, order)
    rhs = y

    # apply weighting
    if w is not None:
        w = NX.asarray(w) + 0.0
        if w.ndim != 1:
            raise TypeError, "expected a 1-d array for weights"
        if w.shape[0] != y.shape[0] :
            raise TypeError, "expected w and y to have the same length"
        lhs *= w[:, NX.newaxis]
        if rhs.ndim == 2:
            rhs *= w[:, NX.newaxis]
        else:
            rhs *= w

    # scale lhs to improve condition number and solve
    scale = NX.sqrt((lhs*lhs).sum(axis=0))
    lhs /= scale
    c, resids, rank, s = lstsq(lhs, rhs, rcond)
    c = (c.T/scale).T  # broadcast scale coefficients

    # warn on rank reduction, which indicates an ill conditioned matrix
    if rank != order and not full:
        msg = "Polyfit may be poorly conditioned"
        warnings.warn(msg, RankWarning)

    if full :
        return c, resids, rank, s, rcond
    elif cov :
        Vbase = inv(dot(lhs.T,lhs))
        Vbase /= NX.outer(scale, scale)
        # Some literature ignores the extra -2.0 factor in the denominator, but
        #  it is included here because the covariance of Multivariate Student-T
        #  (which is implied by a Bayesian uncertainty analysis) includes it.
        #  Plus, it gives a slightly more conservative estimate of uncertainty.
        fac = resids / (len(x) - order - 2.0)
        if y.ndim == 1:
            return c, Vbase * fac
        else:
            return c, Vbase[:,:,NX.newaxis] * fac
    else :
        return c

'''
Created on Sep 23, 2012

@author: will
'''
from numpy.core.numeric import arange, outer
from numpy.ma.extras import dot

if __name__ == '__main__':
    x = arange(10)
    print x
    print dot(x,x)
    print outer(x,x)

def calcM(classKernelList, trainLabels):
    Mlist = []
    for (classKernel, label) in zip(classKernelList, unique(trainLabels)):
        Mlist.append(calcClassM(classKernel, trainLabels, label))
    Mdiff = Mlist[0] - Mlist[1]
    return outer(Mdiff, Mdiff)