#########################################################################################
#
# Scientific Python Course
#
# L04 - SciPy
#
#########################################################################################
#
# ARCHER, 2015
#
# http://www.archer.ac.uk
# support@archer.ac.uk
#
#
# Neelofer Banglawala : nbanglaw@epcc.ed.ac.uk
# Kevin Stratford     : kevin@epcc.ed.ac.uk
#
# EPCC, University of Ediburgh, EPSRC, NERC, CRAY
#
#
# ********************************************************************
# * Reusing this material :                                          *
# *                                                                  *
# * This work is licensed under a Creative Commons Attribution-      *
# * Non-Commercial-ShareAlike 4.0 Internatiional License.            *
# * http://creativecommons.org/licenses/by-nc-sa/4.0/deed.en_US      *
# *                                                                  *
# * This means you are free to copy and redistributhe the material   *
# * and adapt and build on the material under the following terms:   *
# * You must give appropriate credit, provide a link to the license  *
# * and indicate if changes were made. If you adapt or build on the  *
# * material you must distribute your work under the same license as *
# * the original.                                                    *
# *                                                                  *
# ********************************************************************
#
#
#
# 
#########################################################################################
# 
#  S5-1 [SciPy] Integration : find pi I
# 
#########################################################################################
# 
# # numerically solve the integral using dblquad
# import numpy as np;
# from scipy.integrate import dblquad;
# 
# def integrand(y,x): # integrand
#     return 1;
# 
# # integral for the area of a circle with radius r
# def integrate_to_pi(r): 
#     (area,err) = dblquad(integrand, -1*r,r, 
#                          lambda x: -1*np.sqrt(r*r-x*x), 
#                          lambda x: np.sqrt(r*r-x*x)  );
#     return area/(r*r) ;
# 
# 
#########################################################################################
# 
#  S5-2 [SciPy] Integration : find pi II
# 
#########################################################################################
# 
# # result!
# print integrate_to_pi(50.0);
# 
# # compare with numpy pi
# np.pi - integrate_to_pi(50.0);
# 
# # see timing routine from numpy lecture
# # %timeit integrate_to_pi(50.0) 
# 
# # Ex: calculate the double integral of f(x, y) = y * sin (x) 
# # for x = [0, pi/2], y = [0, 1]
# # 
# # ANSWER: 1/2 
# 
# 
#########################################################################################
# 
#  S5-4 [SciPy] Integration : coupled masses II
# 
#########################################################################################
# 
# 
# %matplotlib 
# import numpy as np;
# 
# def x1_t(t,x1,x2,k,m): # exact solution for mass 1 at time t
#     w=np.sqrt(k/m);    # initial velocity assumed zero
#     a1=(x1+x2)/2.0;
#     a2=(x1-x2)/2.0;
#     return a1*np.cos(w*t) + a2*np.cos(np.sqrt(3)*w*t);
# 
# def x2_t(t,x1,x2,k,m): # exact solution for mass 2 at time t
#     w=np.sqrt(k/m);    # initial velocity assumed zero
#     a1=(x1+x2)/2.0;
#     a2=(x1-x2)/2.0;
#     return a1*np.cos(w*t) - a2*np.cos(np.sqrt(3)*w*t);
# 
# # "vectorize" solutions to act on an array of times 
# x1_sol = np.vectorize(x1_t);
# x2_sol = np.vectorize(x2_t);
# 
# 
#########################################################################################
# 
#  S5-5 [SciPy] Integration : coupled masses III
# 
#########################################################################################
# 
# def vectorfield(w, t, p):
#     """Defines the differential equations for the coupled 
#     spring-mass system. Arguments:
#         w :  vector of the state variables: w = [x1,y1,x2,y2]
#         t :  time
#         p :  vector of the parameters: p = [m,k] 
#     """
#     x1, y1, x2, y2 = w;
#     m, k = p;
# 
#     # Create f = (x1',y1',x2',y2'):
#     f = [y1, (-k * x1 + k * (x2 - x1)) / m, 
#          y2, (-k * x2 - k * (x2 - x1)) / m];
#     return f;
# 
#########################################################################################
# 
#  S5-6 [SciPy] Integration : coupled masses IV
# 
#########################################################################################
# 
# # Use ODEINT to solve differential equations defined 
# # by vectorfield
# from scipy.integrate import odeint;
# 
# # Parameters and initial values
# m = 1.0; k = 1.0;     # mass m, spring constant k
# x01 = 0.5; x02 = 0.0; # Initial displacements
# y01 = 0.0; y02 = 0.0; # Initial velocities : LEAVE AS ZERO 
# 
# # ODE solver parameters
# abserr = 1.0e-8; relerr = 1.0e-6; 
# stoptime = 10.0; numpoints = 250;
# 
# # Create time samples for the output of the ODE solver
# t = [stoptime * float(i) / (numpoints - 1) 
#      for i in range(numpoints)];  
# 
# # contd...
# 
#########################################################################################
# 
#  S5-7 [SciPy] Integration : coupled masses IV
# 
#########################################################################################
# 
# # ... contd
# 
# # Pack up the parameters and initial conditions as lists/arrays:
# p =  [m, k];  w0 = [x01, y01, x02, y02];
# 
# # Call the ODE solver
# # Note: args is a tuple
# wsol = odeint(vectorfield, w0, t, args=(p,), atol=abserr, 
#               rtol=relerr);
# 
# # Print and save the solution
# with open('coupled_masses.dat', 'w') as f: 
#     for t1, w1 in zip(t, wsol):
#         print >> f, t1, w1[0], w1[1], w1[2], w1[3]
# 
#########################################################################################
# 
#  S5-8 [SciPy] Integration : coupled masses VI
# 
#########################################################################################
# 
# # import modules for plotting
# 
# import matplotlib.pyplot as plt;
# from matplotlib.font_manager import FontProperties;
# 
# # get saved values from saved file
# t, x1, y1, x2, y2 = np.loadtxt('coupled_masses.dat', unpack=True);
# 
# # contd...
# 
#########################################################################################
# 
#  S5-9 [SciPy] Integration : coupled masses VII
# 
#########################################################################################
# 
# # figure proporties
# plt.figure(1, figsize=(10, 3.5)); plt.xlabel('t'); 
# plt.ylabel('x'); plt.grid(True); plt.hold(True); 
# 
# # plot exact solutions
# time=np.linspace(0,stoptime,50);
# plt.plot(time, x1_sol(time,x01,x02,k,m), 'r*', linewidth=1);
# plt.plot(time, x2_sol(time,x01,x02,k,m), 'mo', linewidth=1);
# # plot numerical solutions
# plt.plot(t, x1, 'b-', linewidth=1); plt.plot(t, x2, 'g-', linewidth=1);
# 
# plt.legend(('$x_{1,sol}$', '$x_{2,sol}$', '$x_{1,num}$', 
#         '$x_{2,num}$'),prop=FontProperties(size=12));
# plt.title('Mass Displacements for the\nCoupled Spring-Mass System');
# plt.savefig('coupled_masses.png', dpi=100); # save figure
# 
#########################################################################################
# 
#  S7-1 [SciPy] Optimisation : maxima I
# 
#########################################################################################
# 
# from scipy import optimize, special;
# x = np.arange(0,10,0.01);
# 
# for k in np.arange(0.5,5.5):
#      y = special.jv(k,x); # Bessel function 
#      plt.plot(x,y);
#      # flip Bessel function about x-axis, turn maxima into minima
#      f = lambda x: -1*special.jv(k,x);
#      # find minimum between 4 and 10
#      x_max = optimize.fminbound(f,0,6);
#     plt.plot([x_max], [special.jv(k,x_max)],'ro') ;
#      plt.title('Different Bessel functions & their local maxima');
# 
# 
#########################################################################################
# 
#  S7-2 [SciPy] Optimsation : maxima II
# 
#########################################################################################
# 
# # Ex: find ** all ** maxima of Airy function in range x = [-8, 2]
# x = np.arange(-8,2,0.01); # x range
# ai =  special.airy(x)[0]; # Airy function, see online docs 
# 
# # What does the Airy function look like?
# plt.figure(figsize=(8,4)); plt.plot(x,ai);
# 
# 
#########################################################################################
# 
#  S7-3 [SciPy] Optimsation : maxima III
# 
#########################################################################################
# 
# # find each maximum (given you can plot the function)
# f = lambda x: -special.airy(x)[0]; # use negated of airy function
#  
# plt.plot(x,ai);
# x_max = optimize.fminbound(f,-8,-6);
# plt.plot([x_max], [special.airy(x_max)[0]],'ro');
# x_max = optimize.fminbound(f,-6,-3);  
# plt.plot([x_max], [special.airy(x_max)[0]],'ro');
# x_max =  optimize.fminbound(f,-3,2);
# plt.plot([x_max], [special.airy(x_max)[0]],'ro');
# plt.title('Different Airy functions and their local maxima');
# 
# 
#########################################################################################
# 
#  S7-4 [SciPy] Optimsation : maxima IV
# 
#########################################################################################
# 
# # ... contd
# # But what if you don't know what function looks like?
# # Could use a fixed 'window' within which to search 
# # BUT this could pick up a false maximum
# minw = -8; maxw = -6; # min, max values for fixed "window"
# plt.plot(x,ai);
# for w in np.arange(4):
#     x_max = optimize.fminbound(f,minw,maxw); # find minima / maxima
#     plt.plot([x_max], [special.airy(x_max)[0]],'ro');
#     minw+=2; maxw +=2; 
#     
# plt.title('Different Airy functions and their local maxima');
# 
# # Ex ++: how could you check you have found a genuine 
# # maximum or minimum?
# 
#########################################################################################
# 
#  S7-6 [SciPy] Optimisation : least-squares fit II
# 
#########################################################################################
# 
# 
# # set up true function and "measured" data
# x = np.arange(0, 6e-2, 6e-2 / 30);
# A, k, theta = 10, 1.0 / 3e-2, np.pi / 6;
# y_true = A * np.sin(2 * np.pi * k * x + theta);
# y_meas = y_true + 2*np.random.randn(len(x));
# 
# # residual function, e_i
# def residuals(p, y, x):
#     A, k, theta = p;
#     err = y - A * np.sin(2 * np.pi * k * x + theta);
#     return err;
# 
#########################################################################################
# 
#  S7-7 [SciPy] Optimisation : least-squares fit II
# 
#########################################################################################
# 
# # for easy evaluation of model function parameters [A, K, theta]
# def peval(x, p):
#     return p[0] * np.sin(2 * np.pi * p[1] * x + p[2]);
# 
# # starting values of A, k and theta
# p0 = [8, 1 / 2.3e-2, np.pi / 3];
# print(np.array(p0));
# 
# 
#########################################################################################
# 
#  S7-8 [SciPy] Optimisation : least-squares fit II
# 
#########################################################################################
# 
# # do least squares fitting
# from scipy.optimize import leastsq
# 
# plsq = leastsq(residuals, p0, args=(y_meas, x));
# print(plsq[0]); print(np.array([A, k, theta]));
# 
# # and plot the true function, measured (noisy) data 
# # and the model function with fitted parameters 
# plt.plot(x, peval(x, plsq[0]),x,y_meas,'o',x,y_true);
# 
# plt.title('Least-squares fit to noisy data');
# plt.legend(['Fit', 'Noisy', 'True']);
# 
# 
#########################################################################################
# 
#  S7-9 [SciPy] Convolution functions I
# 
#########################################################################################
# 
# # for example 
# import numpy as np
# from scipy import signal
# a = np.array([[1, 2, 3, 4], [5, 6, 7, 8],[9, 10, 11, 12]])
# mask = np.array([[1,1,1],[1,0,1],[1,1,1]])
# cS = signal.convolve(a, mask)
# print a; print; print mask; print; print cS
# 
#########################################################################################
# 
#  S7-11 [SciPy] Convolution functions III
# 
#########################################################################################
# 
# # Ex: Scipy has another convolve function, ndimage.convolve  
# # What does this produce, when used in 'constant' mode (below)?
# from scipy import ndimage
# cN = ndimage.convolve(a, mask, mode='constant')
# 
# # The default mode for ndimage.convolve is 'reflect'. Investigate. 
# 
# # ANSWER
# # ndimage.convolve only produces the inner dotted matrix as shown
# # in the green matrix in the figure i.e. convolution that fits 
# # the size of array a (not a true (commutative) convolution, like 
# # signal.convolve. With 'reflect' it replaces 'missing' values 
# # with adjacent values, e.g. 
# # cN[0,0]= 23= 5+6+2 (as constant mode)  
# #       +5+1 (from the left)  +1+1+2 (from the top)
# # cN[0,3]= 41= 3+7+8 (same as constant mode) 
# #       +8+4 (from the right) +3+4+4 (from the top)
# cN = ndimage.convolve(a, mask, mode='reflect')
# print cN
# 
#########################################################################################
# 
#  S9-1 [SciPy] Linear algebra : inverse matrix II
# 
#########################################################################################
# 
# # find inverse of matrix
# from scipy import linalg
# A = np.array([[1,3,5],[2,5,1],[2,3,8]])
# 
# invA = linalg.inv(A) # the inveser
# print invA; print;
# 
# # check inverse gives identity 
# identity=A.dot(linalg.inv(A))
# 
# # This really is the identity if we clean up the errors
# np.abs(np.around(identity,0)) # (round up small errors)
# 
#########################################################################################
# 
#  END
# 
#########################################################################################