dff

import numpy
import pylab
import matplotlib.pyplot
import scipy.optimize
from scipy.optimize import curve_fit

''' A Program That Determines The Reduced Chi Squared Value For Various Theoretical Models.'''
'''The Best Fit Parameters Are Derived Using Levenberg-Marquardt Algorithm Which Solves The Non-Linear Least Squares Problem.'''

#Recorded Observed Data
xdata = numpy.array([ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
ydata = numpy.array([-10, -7, 1, 5, 7, 7, 6, 7, 10, 15])
xerror = numpy.array([0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3])
yerror = numpy.array([0.1, 0.1, 0.4, 0.3, 0.2, 0.1, 0.1, 0.2, 0.1, 0.1])


    #Theoretical Models To Be Tested
    def test_linear_model():
        print 'Testing linear model', '\n'

        #Scipy Optimise cuve_fit Model Produces Expected Values Using A Linear Model
        def func(x, a, b):
            return a*x + b

        #Constants Of Theoretical Model
        popt, pcov = curve_fit(func, xdata, ydata)
        print 'Constant Ax is',("%.2f" %popt[0])
        print 'Constant B is',("%.2f" %popt[1]), '\n'
        return func, popt

    def test_quadratic_model():
        print 'Testing quadratic model', '\n'

        #Scipy Optimise cuve_fit Model Produces Expected Values Using A Quadratic Model Polynomial Model
        def func(x, a, b, c):
            return a*x**2 + b*x + c

            #Constants Of Theoretical Model
        popt, pcov = curve_fit(func, xdata, ydata)
        print 'Constant Ax^2 is', ("%.2f" %popt[0])
        print 'Constant Bx is', ("%.2f" %popt[1])
        print 'Constant C is', ("%.2f" %popt[2]), '\n'

    def test_cubic_model():
        print 'Testing cubic model', '\n'

        #Scipy Optimise cuve_fit Model Produces Expected Values Using A Cubic Model Polynomial Model
        def func(x, a, b, c, d):
            return a*x**3 + b*x**2 + c*x + d

        #Constants Of Theoretical Model
        popt, pcov = curve_fit(func, xdata, ydata)
        print 'Constant Ax^3 is', ("%.2f" %popt[0])
        print 'Constant Bx^2 is', ("%.2f" %popt[1])
        print 'Constant Cx is', ("%.2f" %popt[2])
        print 'Constant D is', ("%.2f" %popt[3]), '\n'

    def test_quartic_model():
        print 'Testing quartic model', '\n'

        #Scipy Optimise cuve_fit Model Produces Expected Values Using A Quartic Polynomial Model
        def func(x, a, b, c, d, e):
            return a*x**4 + b*x**3 + c*x**2 + d*x + e

        #Constants Of Theoretical Model
        popt, pcov = curve_fit(func, xdata, ydata)
        print 'Constant Ax^4 is', ("%.2f" %popt[0])
        print 'Constant Bx^3 is', ("%.2f" %popt[1])
        print 'Constant Cx^2 is', ("%.2f" %popt[2])
        print 'Constant Dx is', ("%.2f" %popt[3])
        print 'Constant E is', ("%.2f" %popt[4]), '\n'

    def test_quintic_model():
        print 'Testing quintic model', '\n'

        #Scipy Optimise cuve_fit Model Produces Expected Values Using A Quintic Polynomial Model
        def func(x, a, b, c, d, e, f):
            return a*x**5 + b*x**4 + c*x**3 + d*x**2 + e*x + f

        #Constants Of Theoretical Model
        popt, pcov = curve_fit(func, xdata, ydata)
        print 'Constant Ax^5 is', ("%.2f" %popt[0])
        print 'Constant Bx^4 is', ("%.2f" %popt[1])
        print 'Constant Cx^3 is', ("%.2f" %popt[2])
        print 'Constant Dx^2 is', ("%.2f" %popt[3])
        print 'Constant Ex is', ("%.2f" %popt[4])
        print 'Constant F is', ("%.2f" %popt[5]), '\n'

    def test_logarithmic_model():
        print 'Testing logarithmic model', '\n'

        #Scipy Optimise cuve_fit Model Produces Expected Values Using A Logarithmic Model
        def func(x, a, b):
            return a*numpy.log(x)+ b

        #Constants Of Least Squared Theoretical Model
        popt, pcov = curve_fit(func, xdata, ydata)
        print 'Constant A is', ("%.2f" %popt[0])
        print 'Constant Bx is', ("%.2f" %popt[1]), '\n'

    def test_powerlaw_model():
        print 'Testing power law model', '\n'

        #Scipy Optimise cuve_fit Model Produces Expected Values Using A Power Law Model
        def func(x, a, b):
            return a*(x**b)

        #Constants Of Least Squared Theoretical Model
        popt, pcov = curve_fit(func, xdata, ydata)
        print 'Constant A is', ("%.2f" %popt[0])
        print 'Constant B is', ("%.2f" %popt[1]), '\n'

    def test_exponential_model():
        print 'Testing exponential model', '\n'

        #Scipy Optimise cuve_fit Model Produces Expected Values Using A Exponential Model
        def func(x, a, b, c):
            return a*numpy.exp(b*x + c)

        #Constants Of Least Squared Theoretical Model
        popt, pcov = curve_fit(func, xdata, ydata)
        print 'Constant A is', ("%.2f" %popt[0])
        print 'Constant Bx is', ("%.2f" %popt[1])
        print 'Constant C is', ("%.2f" %popt[2]), '\n'

    def test_exponential_model():
        print 'Testing exponetial offset model', '\n'

        #Scipy Optimise cuve_fit Model Produces Expected Values Using A Exponential Offset Model
        def func(x, a, b, c, d):
            return a*numpy.exp(b*x + c) + d

        #Constants Of Least Squared Theoretical Model
        popt, pcov = curve_fit(func, xdata, ydata)
        print 'Constant A is', ("%.2f" %popt[0])
        print 'Constant Bx is', ("%.2f" %popt[1])
        print 'Constant C is', ("%.2f" %popt[2])
        print 'Constant D is', ("%.2f" %popt[3]), '\n'

def main():
    commands = {
        1: test_linear_model,
        2: test_quadratic_model,
        3: test_cubic_model,
        4: test_quartic_model,
        5: test_quintic_model,
        6: test_logarithmic_model,
        7: test_powerlaw_model,
        8: test_exponential_model,
        9: test_exponentialoffset_model,
    }

    # User Prompted Input Values
    test = int(raw_input("Enter A Number Assigned Theoretical Model To Test: "))
    func, popt = commands[test]()


    #Derived Chi Squared Value For This Model
    chi_squared = numpy.sum(((func(xdata, *popt) - ydata) / xerror) ** 2)
    reduced_chi_squared = chi_squared / (len(xdata) - len(popt))
    print 'The degrees of freedom for this test is', len(xdata) - len(popt)
    print 'The chi squared value is: ', ("%.2f" % chi_squared)
    print 'The reduced chi squared value is: ', ("%.2f" % reduced_chi_squared)


    #Observed Values Are Plotted With Expected Values
    matplotlib.pyplot.figure()
    matplotlib.pyplot.scatter(xdata, ydata, s=0)
    matplotlib.pyplot.plot(xdata,func(xdata, *popt), label='Theoretical Model')
    matplotlib.pyplot.errorbar(xdata, ydata, xerr=xerror, yerr=yerror, linestyle="", color='red')
    matplotlib.pyplot.xlabel(' Observed Values For x ')
    matplotlib.pyplot.ylabel(' f(x)')
    matplotlib.pyplot.legend(loc='lower right')
    matplotlib.pyplot.show()

if __name__ == '__main__':
    main()