def fitData2(fileName):
    '''
    Models predictions using Terman's law model (cubic fit) and 
    Hooks Law (linear fit).

    Hook's Law functions up to the point where the spring reaches 
    it's elastic limit - when it stops behaving as a spring but instead 
    as a rope, etc (doesn't get longer b/c hang more weight on it)
    '''
    xVals, yVals = getData(fileName)
    extX = pylab.array(xVals + [1.05, 1.1, 1.15, 1.2, 1.25])
    xVals = pylab.array(xVals)
    yVals = pylab.array(yVals)
    xVals = xVals*9.81  # convert mass to force (F = mg)
    extX = extX*9.81    # convert mass to force (F = mg)
    pylab.plot(xVals, yVals, 'bo', label = 'Measured displacements')
    pylab.title('Measured Displacement of Spring')
    pylab.xlabel('|Force| (Newtons)')
    pylab.ylabel('Distance (meters)')
    a,b = pylab.polyfit(xVals, yVals, 1)
    estYVals = a*extX + b
    pylab.plot(extX, estYVals, label = 'Linear fit')
    a,b,c,d = pylab.polyfit(xVals, yVals, 3)
    estYVals = a*(extX**3) + b*extX**2 + c*extX + d
    pylab.plot(extX, estYVals, label = 'Cubic fit')
    pylab.legend(loc = 'best')

    def getCloneReplicates(self, clone, source, condition, applyFilter=False):
        '''Retrieve all growth curves for a clone+source+condition'''
        # Check if any other replicates should be returned
        # retArray is a 2xN multidimensional numpy array
        retArray = py.array([])
        first = True
        for i in xrange(1, self.numReplicates[clone] + 1):
            # Get replicate
            filterMe = self.dataHash[clone][i][source][condition]['filter']
            currCurve = self.dataHash[clone][i][source][condition]['od']

            # Check if filter is enabled and curve should be filtered
            if applyFilter and filterMe:
                continue

            # Create multidimensional array if first
            elif first:
                retArray = py.array([currCurve])
                first = False

            # Append to multidimensional array if not first
            else:
                retArray = py.concatenate((retArray,
                                           py.array([currCurve])))

        return retArray

    def t(self, k, cosTheta, pk, c):
        """
        Raw trispectrum

        Not recently tested
        """
        pk = c.pkInterp(k)
        f2term = (
            self.tf21(0, 1, 2, k, cosTheta, pk, c)
            + self.tf21(1, 2, 0, k, cosTheta, pk, c)
            + self.tf21(2, 0, 1, k, cosTheta, pk, c)
            + self.tf21(1, 2, 3, k, cosTheta, pk, c)
            + self.tf21(2, 3, 1, k, cosTheta, pk, c)
            + self.tf21(3, 1, 2, k, cosTheta, pk, c)
            + self.tf21(2, 3, 0, k, cosTheta, pk, c)
            + self.tf21(3, 0, 2, k, cosTheta, pk, c)
            + self.tf21(0, 2, 3, k, cosTheta, pk, c)
            + self.tf21(3, 0, 1, k, cosTheta, pk, c)
            + self.tf21(0, 1, 3, k, cosTheta, pk, c)
            + self.tf21(1, 3, 0, k, cosTheta, pk, c)
        ) * 4.0

        f3term = (
            self.tf31(M.array([0, 1, 2]), k, cosTheta, pk)
            + self.tf31(M.array([1, 2, 3]), k, cosTheta, pk)
            + self.tf31(M.array([2, 3, 1]), k, cosTheta, pk)
            + self.tf31(M.array([3, 1, 2]), k, cosTheta, pk)
        ) * 6.0

        # print cosTheta,f2term, f3term, ft2term+f3term
        return f2term + f3term

def main():
    mu = pl.array([[0], [12], [24], [36]])
    Sigma = pl.array([[3.01602775,  1.02746769, -3.60224613, -2.08792829],
                      [1.02746769,  5.65146472, -3.98616664,  0.48723704],
                      [-3.60224613, -3.98616664, 13.04508284, -1.59255406],
                      [-2.08792829,  0.48723704, -1.59255406,  8.28742469]])

    # The data matrix is created for above mu and Sigma.
    d, U = pl.eig(Sigma)
    L = pl.diagflat(d)
    A = pl.dot(U, pl.sqrt(L))
    X = pl.randn(4, 1000)

    # Y is the data matrix of random samples.
    Y = pl.dot(A, X) + pl.tile(mu, 1000)

    pl.figure(1)
    pl.clf()
    pl.plot(X[0], Y[1], '+', color='#0000FF', label='i=0,j=1')
    pl.plot(X[0], Y[2], '+', color='#FF0000', label='i=0,j=2')
    pl.plot(X[0], Y[3], '+', color='#00FF00', label='i=0,j=3')
    pl.plot(X[1], Y[0], 'x', color='#FFFF00', label='i=1,j=0')
    pl.plot(X[1], Y[2], 'x', color='#00FFFF', label='i=1,j=2')
    pl.plot(X[1], Y[3], 'x', color='#444444', label='i=1,j=3')
    pl.plot(X[2], Y[0], '.', color='#774411', label='i=2,j=0')
    pl.plot(X[2], Y[1], '.', color='#222222', label='i=2,j=1')
    pl.plot(X[2], Y[3], '.', color='#AAAAAA', label='i=2,j=3')
    pl.plot(X[3], Y[0], '+', color='#FFAA22', label='i=3,j=0')
    pl.plot(X[3], Y[1], '+', color='#22AAFF', label='i=3,j=1')
    pl.plot(X[3], Y[2], '+', color='#FFDD00', label='i=3,j=2')
    pl.legend()
    pl.savefig('fig21.png')

def fitData3(fileName):

    # xVals is type 'numpy.ndarray'
    # xVals[0] will return the 1st item in array
    xVals, yVals = getData(fileName)
    xVals = pylab.array(xVals[:-6])
    yVals = pylab.array(yVals[:-6])
    xVals = xVals*9.81  # convert mass to force (F = mg)

    observed_data_variance = calcVariance(xVals)

    # need to grab the Y values from somewhere ??? maybe estYVals
    # to compare with observed data Y values
    # to calculate the variance of errors
    # errors_variance = calcVariance(xVals)

    coefficient_determination = 0


    pylab.plot(xVals, yVals, 'bo', label = 'Measured points')
    pylab.title('Measured Displacement of Spring')
    pylab.xlabel('Force (Newtons)')
    pylab.ylabel('Distance (meters)')
    a,b = pylab.polyfit(xVals, yVals, 1)  # fix y = ax + b

    # use line equation to graph predicted values
    estYVals = a*xVals + b
    k = 1/a
    pylab.plot(xVals, estYVals, label = 'Linear fit, k = '
               + str(round(k, 5)))
    pylab.legend(loc = 'best')

def compare_models(db, stoch="itn coverage", stat_func=None, plot_type="", **kwargs):
    if stat_func == None:
        stat_func = lambda x: x

    X = {}
    for k in sorted(db.keys()):
        c = k.split("_")[2]
        X[c] = []

    for k in sorted(db.keys()):
        c = k.split("_")[2]
        X[c].append([stat_func(x_ki) for x_ki in db[k].__getattribute__(stoch).gettrace()])

    x = pl.array([pl.mean(xc[0]) for xc in X.values()])
    xerr = pl.array([pl.std(xc[0]) for xc in X.values()])
    y = pl.array([pl.mean(xc[1]) for xc in X.values()])
    yerr = pl.array([pl.std(xc[1]) for xc in X.values()])

    if plot_type == "scatter":
        default_args = {"fmt": "o", "ms": 10}
        default_args.update(kwargs)
        for c in X.keys():
            pl.text(pl.mean(X[c][0]), pl.mean(X[c][1]), " %s" % c, fontsize=8, alpha=0.4, zorder=-1)
        pl.errorbar(x, y, xerr=xerr, yerr=yerr, **default_args)
        pl.xlabel("First Model")
        pl.ylabel("Second Model")
        pl.plot([0, 1], [0, 1], alpha=0.5, linestyle="--", color="k", linewidth=2)

    elif plot_type == "rel_diff":
        d1 = sorted(100 * (x - y) / x)
        d2 = sorted(100 * (xerr - yerr) / xerr)
        pl.subplot(2, 1, 1)
        pl.title("Percent Model 2 deviates from Model 1")

        pl.plot(d1, "o")
        pl.xlabel("Countries sorted by deviation in mean")
        pl.ylabel("deviation in mean (%)")

        pl.subplot(2, 1, 2)
        pl.plot(d2, "o")
        pl.xlabel("Countries sorted by deviation in std err")
        pl.ylabel("deviation in std err (%)")
    elif plot_type == "abs_diff":
        d1 = sorted(x - y)
        d2 = sorted(xerr - yerr)
        pl.subplot(2, 1, 1)
        pl.title("Percent Model 2 deviates from Model 1")

        pl.plot(d1, "o")
        pl.xlabel("Countries sorted by deviation in mean")
        pl.ylabel("deviation in mean")

        pl.subplot(2, 1, 2)
        pl.plot(d2, "o")
        pl.xlabel("Countries sorted by deviation in std err")
        pl.ylabel("deviation in std err")
    else:
        assert 0, "plot_type must be abs_diff, rel_diff, or scatter"

    return pl.array([x, y, xerr, yerr])

 def read_doscar(self, fname="DOSCAR"):
     """Read a VASP DOSCAR file."""
     f = open(fname)
     natoms = int(f.readline().split()[0])
     [f.readline() for n in range(4)]  # Skip next 4 lines.
     dos = []
     for na in xrange(natoms + 1):
         try:
             line = f.readline()
             if line == "":
                 raise Exception
         except Exception, e:
             errstr = (
                 "Failed reading "
                 + str(na)
                 + ":th DOS block, probably "
                 + "this DOSCAR is from some old version of VASP that "
                 + "doesn't "
                 + "first produce a block with integrated DOS. Inserting "
                 + "empty 0:th block."
             )
             sys.stderr.write(errstr)
             dos.insert(0, pl.zeros((ndos, dos[1].shape[1])))
             continue
         try:
             ndos = int(line.split()[2])
         except:
             print "Error, line is: " + line + "ENDLINE"
         line = f.readline().split()
         cdos = pl.zeros((ndos, len(line)))
         cdos[0] = pl.array(line)
         for nd in xrange(1, ndos):
             line = f.readline().split()
             cdos[nd] = pl.array(line)
         dos.append(cdos)

def anim_update(i):
    """
    i: frame num
    """

    ## equivalent time = (i/tot)*(runtime/fct)
    t = (float(i)/float(tot_frames))*(anim_run_time/float(fct))
    # lnx, lny = line.get_data()
    # pos_x, pos_y = P.get_position(t)
    # head.set_data([pos_x], [pos_y])
    # lnx = pylab.array([k for k in lnx] + [pos_x])
    # lny = pylab.array([k for k in lny] + [pos_y])
    # line.set_data(lnx, lny)

    for i in range(len(Projs)):
        pos_x, pos_y = Projs[i].get_position(t)
        
        p_heads[i].set_data([pos_x], [pos_y])

        lnx, lny = p_lines[i].get_data()
        lnx = pylab.array([k for k in lnx] + [pos_x])
        lny = pylab.array([k for k in lny] + [pos_y])
        p_lines[i].set_data(lnx, lny)

    return p_lines + p_heads

def rollingMeanScale(series, period, plotAxis=False):
    svr_rbf = SVR(kernel='rbf', C=1e4, gamma=.01, epsilon=.01)
    '''Fit Model to Data Series'''
    tS= numpy.array([series.index]).T
    y_rbf = svr_rbf.fit(tS, list(series))
    '''Up-sample to get rid of bias'''
    fFit = arange(series.index[0],series.index[-1]+.1,.25)
    trend = y_rbf.predict(numpy.array([fFit]).T)
    
    '''Take rolling mean over 1-day window'''
    shift = int(round(period/.5))
    rMean = pandas.rolling_mean(trend, shift*2)
    rMean = numpy.roll(rMean, -shift)
    rMean[:shift]=rMean[shift]
    rMean[-(shift+1):]=rMean[-(shift+1)]
    rMean = pandas.Series(rMean, index=fFit)
    
    '''Adjust Data Series by subtracting out trend'''
    series = series - array(rMean[array(series.index, dtype=float)])
    series = scaleMe(series)-.5
    
    if plotAxis:
        plotAxis.plot(fFit, trend, label='Series Trend')
        plotAxis.plot(fFit, rMean, label='Rolling Mean')
        plotAxis.set_title('Detrend the Data')
        plotAxis.legend(loc='lower left')

    return series

def main():
    f = open("final_position.txt","r")
    data = pl.genfromtxt(f,comments = "L")
    
    # need to get every other
    x = pl.array([])
    y = pl.array([])
    for i,j in enumerate(data[:-7,2]):
        if i%4 == 0:
            x = pl.append(x,data[i,4])
            y = pl.append(y,j)
    
    print(x)
    print(y)
    fit = np.polyfit(x,y,2)

    print(fit)
    
    #fited = fit[0]+fit[1]*x + fit[2]*x**2
    fited = np.poly1d(fit)
    print(fited)

    #pl.plot(pl.append(x,[.262,.264,.266]),fited(pl.append(x,[.262,.264,.266])),color="black")
    pl.scatter(x,y,color = "black")
    pl.xlabel("$A$",fontsize="30")
    pl.ylabel("$x$",fontsize="30")
    pl.savefig("fin_pts.png",transparent=True,dpi=300)
    
    os.system("open fin_pts.png")

def rotate_molecule(coords, rotp = m.array((0.,0.,0.)), phi = 0., \
        theta = 0., psi = 0.):
    """Rotate a molecule via Euler angles.

    See http://mathworld.wolfram.com/EulerAngles.html for definition.
    Input arguments:
    coords: Atom coordinates, as Nx3 2d pylab array.
    rotp: The point to rotate about, as a 1d 3-element pylab array
    phi: The 1st rotation angle around z axis.
    theta: Rotation around x axis.
    psi: 2nd rotation around z axis.

    """
# First move the molecule to the origin
# In contrast to MATLAB, numpy broadcasts the smaller array to the larger
# row-wise, so there is no need to play with the Kronecker product.
    rcoords = coords - rotp
# First Euler rotation about z in matrix form
    D = m.array(((m.cos(phi), m.sin(phi), 0.), (-m.sin(phi), m.cos(phi), 0.), \
            (0., 0., 1.)))
# Second Euler rotation about x:
    C = m.array(((1., 0., 0.), (0., m.cos(theta), m.sin(theta)), \
            (0., -m.sin(theta), m.cos(theta))))
# Third Euler rotation, 2nd rotation about z:
    B = m.array(((m.cos(psi), m.sin(psi), 0.), (-m.sin(psi), m.cos(psi), 0.), \
            (0., 0., 1.)))
# Total Euler rotation
    A = m.dot(B, m.dot(C, D))
# Do the rotation
    rcoords = m.dot(A, m.transpose(rcoords))
# Move back to the rotation point
    return m.transpose(rcoords) + rotp

def fitData(fileName):
    '''
    Using Pylab's polyfit to find equations of the line to best fit the data from Hooks Law
    experiment.

    Hooks Law represented with equation - y = ax + b 

    y - Measured distance
    x - Force 
    '''
    xVals, yVals = getData(fileName)
    xVals = pylab.array(xVals)
    yVals = pylab.array(yVals)
    xVals = xVals*9.81  # convert mass to force (F = mg)
    pylab.plot(xVals, yVals, 'bo', label = 'Measured points')
    pylab.title('Measured Displacement of Spring')
    pylab.xlabel('Force (Newtons)')
    pylab.ylabel('Distance (meters)')
    a,b = pylab.polyfit(xVals, yVals, 1)  # fit y = ax + b
    # use line equation to graph predicted values
    estYVals = a*xVals + b
    k = 1/a
    pylab.plot(xVals, estYVals, label = 'Linear fit, k = '
               + str(round(k, 5)))
    pylab.legend(loc = 'best')

def tryFits1(fName):
    '''
    Calculate the coefficient of determination (R**2) to determine how
    well the model fits the data and ergo could make predictions.
    '''
    distances, heights = getTrajectoryData(fName)
    distances = pylab.array(distances)*36
    totHeights = pylab.array([0]*len(distances))
    for h in heights:
        totHeights = totHeights + pylab.array(h)
    pylab.title('Trajectory of Projectile (Mean of 4 Trials)')
    pylab.xlabel('Inches from Launch Point')
    pylab.ylabel('Inches Above Launch Point')
    meanHeights = totHeights/float(len(heights))
    pylab.plot(distances, meanHeights, 'bo')
    a,b = pylab.polyfit(distances, meanHeights, 1)
    altitudes = a*distances + b
    pylab.plot(distances, altitudes, 'r',
               label = 'Linear Fit' + ', R2 = '
               + str(round(rSquare(meanHeights, altitudes), 4)))
    a,b,c = pylab.polyfit(distances, meanHeights, 2)
    altitudes = a*(distances**2) + b*distances + c
    pylab.plot(distances, altitudes, 'g',
               label = 'Quadratic Fit' + ', R2 = '
               + str(round(rSquare(meanHeights, altitudes), 4)))
    pylab.legend()

def tryFits(fName):
    '''
    Linear fit does not fit the data. Not a logical assumption that the arrow
    flies in a straight line to the target.

    Quadratic fit mirrors a parabolic pathway.
    '''
    distances, heights = getTrajectoryData(fName)
    distances = pylab.array(distances)*36 # Convert yards to feet
    totHeights = pylab.array([0]*len(distances))
    for h in heights:
        totHeights = totHeights + pylab.array(h) # Get one avg measurement of height
    pylab.title('Trajectory of Projectile (Mean of 4 Trials)')
    pylab.xlabel('Inches from Launch Point')
    pylab.ylabel('Inches Above Launch Point')
    meanHeights = totHeights/float(len(heights))
    pylab.plot(distances, meanHeights, 'bo')
    a,b = pylab.polyfit(distances, meanHeights, 1)
    altitudes = a*distances + b
    pylab.plot(distances, altitudes, 'r',
               label = 'Linear Fit')
    a,b,c = pylab.polyfit(distances, meanHeights, 2)
    altitudes = a*(distances**2) + b*distances + c 
    pylab.plot(distances, altitudes, 'g',
               label = 'Quadratic Fit')
    pylab.legend()

def normal(name, pi, sigma, p, s):
    """ Generate PyMC objects for a normal model

    :Parameters:
      - `name` : str
      - `pi` : pymc.Node, expected values of rates
      - `sigma` : pymc.Node, dispersion parameters of rates
      - `p` : array, observed values of rates
      - `s` : array, standard error of rates

    :Results:
      - Returns dict of PyMC objects, including 'p_obs' and 'p_pred' the observed stochastic likelihood and data predicted stochastic

    """
    p = pl.array(p)
    s = pl.array(s)

    assert pl.all(s >= 0), "standard error must be non-negative"

    i_inf = pl.isinf(s)

    @mc.observed(name="p_obs_%s" % name)
    def p_obs(value=p, pi=pi, sigma=sigma, s=s):
        return mc.normal_like(value, pi, 1.0 / (sigma ** 2.0 + s ** 2.0))

    s_noninf = s.copy()
    s_noninf[i_inf] = 0.0

    @mc.deterministic(name="p_pred_%s" % name)
    def p_pred(pi=pi, sigma=sigma, s=s_noninf):
        return mc.rnormal(pi, 1.0 / (sigma ** 2.0 + s ** 2.0))

    return dict(p_obs=p_obs, p_pred=p_pred)

def meanVar(data,outMean=None):
    """
    each data[i] is a vector, return the average vector and the variance
    vector
 
    each vector is assumed to be the same length of course
 
    if a mean is available from outside, set outMean which is a vector
    of len(data[i]); then the mean returned is a residual
    """

    mean = []
    var = []
    n = len(data)
    print n,len(data[0])
    try:
        for i in range(n):
            data[i] -= outMean
    except TypeError:
        pass
    for i in range(len(data[0])):
        m = 0.0
        v = 0.0
        for j in range(n):
            m += data[j][i]
            v += data[j][i]*data[j][i]
        if outMean == None:
            v -= m*m/float(n)
        m /= float(n)
        v /= (n-1)
        mean.append(m)
        var.append(v)
    return (M.array(mean),M.array(var))

        def residuals(params,x,y,z):
            xo = params[0]
            xs = params[1]
            yo = params[2]
            ys = params[3]
            zo = params[4]
            zs = params[5]

            xc = empty(shape(x))
            for i in range(len(x)):
                xc[i] = (x[i] - xo) * xs

            yc = empty(shape(y))
            for i in range(len(y)):
                yc[i] = (y[i] - yo) * ys

            zc = empty(shape(z))
            for i in range(len(z)):
                zc[i] = (z[i] - zo) * zs

            res = []
            for i in range(len(xc)):
                norm = l2norm(array([xc[i],yc[i],zc[i]])) - 1.0
                res.append(norm)

            return array(res)

def createExpData(f,xVals):
    '''假设f是只有一个参数的指数函数， xVals是一个数组，其中的值可以作为参数传入f
       返回一个数组，包含xVals中的元素传入f生成的结果'''
    yVals=[]
    for i in range(len(xVals)):
        yVals.append(f(xVals[i]))
    return pylab.array(xVals),pylab.array(yVals)

def getApices(y):
    """ 
    returns the time (in frames) and position of initial and final apex
    height from a given trajectory y, which are obtained by fitting a cubic
    spline 
    
    ==========
    parameter:
    ==========
    y : *array* (1D)
        the trajectory. Should ideally start ~1 frame before an apex and end ~1
        frame behind an apex

    ========
    returns:
    ========

    [x0, xF], [y0, yF] : location (in frames) and value of the first and final
    apices

    """
    # the math behind here is: fitting a 2nd order polynomial and finding 
    # the root of its derivative. Here, only the results are applied, that's
    # why it appears like "magic" numbers
    c = dot(array([[.5, -1, .5], [-1.5, 2., -.5], [1., 0., 0.]]), y[:3])
    x0 = -1. * c[1] / (2. * c[0])
    y0 = polyval(c, x0)

    c = dot(array([[.5, -1, .5], [-1.5, 2., -.5], [1., 0., 0.]]), y[-3:])
    xF = -1. * c[1] / (2. * c[0])
    yF = polyval(c, xF)
    xF += len(y) - 3

    return [x0, xF], [y0, yF]

def findOrder(xVals, yVals, accuracy = 1.0e-1):
    X = pylab.array(xVals)
    Y = pylab.array(yVals)
    find = True
    n = 0
    while find:
        #n += 1
        Z = pylab.polyfit(xVals, yVals, n, None, True)
        '''
        est = 0.0
        estVal = []
        for x in xVals:
            for i in range(len(Z[0])):
                est += Z[0][i] * (x**(len(Z[0])-i-1))
            estVal.append(est)
        error = 0.0
        for j in range(len(yVals)):
            error += (estVal[j] - yVals[j]) ** 2
        if error <= accuracy:
            find = False
            print 'order ', n, 'is good enough'
        '''               
        if Z[1] <= accuracy:
            find = False
            print 'order ', n, 'is good enough'
        n += 1
    print (Z[0])
    return Z[0]