# setup 
import numpy as np
import sympy as sp
import scipy
from scipy import linalg
from pprint import pprint
sp.init_printing(use_latex='mathjax')

import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = (12, 8)  # (width, height)
plt.rcParams['font.size'] = 14
plt.rcParams['legend.fontsize'] = 16
from matplotlib import patches

#get_ipython().magic('matplotlib')  # seperate window
get_ipython().magic('matplotlib inline') # inline plotting


"""
Explanation: Mechpy Tutorials
a mechanical engineering toolbox
source code - https://github.com/nagordon/mechpy
documentation - https://nagordon.github.io/mechpy/web/  

Neal Gordon
2017-02-20  

Math
Linear Algebra
Signal Processing  
plotting  

This is more or less my notes on using python for a variety of engineering tasks, including plotting, solving equations, linear algebra etc. Some of this work is my own, but many of it is from the scipy documentation.
Python Initilaization with module imports
End of explanation
"""


# using the range function to create a numpy array
x = np.array(range(10))
print('x = ',x, 'as a',x.dtype)

# using numpy function arange
x = np.arange(10)
print('x = ',x, 'as a',x.dtype)

x = np.linspace(0,9,10)
print('x = ',x, 'as a',x.dtype)

type(x)

# multiple array operations
x*3

# exponent array operations
x**2

# trig array operations
np.sin(x)

x.sum()

print([k for k in dir(x) if k[:2] != '__' ])

"""
Explanation: Math with Python
index
WHen using ipython, variables and modules can be stored in the console. Use %reset to remove any modules and variables from the kernal 
The modules used for math are numpy , scipy and matplotlib
All the these modules can be accessed with the module pylab, but is not updated as regular as the other packages.
the linear algebra system in python via numpy and scipy is 
a bit awkward at first, but once the nuances are understood
it has most of the functionality of MATLAB with zero the cost
The matrix class of scipy is handy but limits the type of data stored.
 As with MATLAB all "arrays" or "matrices" are just matrices. This is not true with python.
Thus, to avoid confusion, it is discouraged to use the matrix class
TO get started lets looks at numpy arrays, and then later we will use python powerful symbolic math package sympy for a different way to perform calculations
Linear Algebra with Python
index  
Python's numpy package allows python, a generic computing language to perform powerful mathematical calculations. Although python's math syntax is not as obvious as MATLAB's, the functionality is comparable. This document is designed to be an intro to that syntax 
Some references
http://nbviewer.ipython.org/github/carljv/cython_testing/blob/master/cython_linalg.ipynb
We can either use scipy, which includes numpy, 
http://docs.scipy.org/doc/
or use numpy directly
http://docs.scipy.org/doc/numpy/
Since there are many ways to solve linear algebra problems, (eg Octave/Matlab, julia, scipy, numpy)
I tend to prefer the most matlabesc approaches due to the ubiquity of Matlab and the simplicity of the syntax, which frankly, python suffers with.
The major difference between arrays and matrices in python is that arrays are n-dimensions, where matrices are only up to 2-dimensions
m
Numpy arrays
End of explanation
"""


# logical indexing
3 < x

# negate logical indexing
~(3 < x)

# logical indexing and 
(3 < x) & (x < 5)

# logical index or
(3 < x) | (x < 5)

# logical indexing and 
(3 < x) & (x < 5)

z = (3 < x) & (x < 5)
z

# if 0,1 is preferred, changethe datatype
z.astype(np.int)

"""
Explanation: Logical Indexing
End of explanation
"""


# createa 1x3 array
c = np.array([1,2,3])
print('c =',c,' of shape',c.shape)

#does not  convert to 3x1 since it is only a 1d array
# this is not what I would expected, so I will call it a nuance of arrays in python
c.transpose()

# we must reshape it
c = c.reshape((1,3))
print('c =',c,' of shape',c.shape)

c.transpose()

# or when we defined it, add brackets in it to create a true 1x3 array
c = np.array([[1,2,3]])
print('c =',c,' of shape',c.shape)

# now we can transpose it.
c.transpose() 

print('c =',c,' of shape',c.shape)

c.shape

# lets creat some matrices (but recall use the array)
a = np.array([[1, 2], [3, 4]])
b = np.array([[1, 4], [5, 7]])

a

b

a * b

# to perform a matrix multiplecation on an array, use 
a @ b

# or 
np.dot(a,b)

# to solve ax=b, use 
linalg.solve(a,b)

# or slower and less accurate
linalg.inv(a) @ b

# or 
linalg.inv(a).dot(b)

"""
Explanation: numpy arrays as a matrix
End of explanation
"""


# 3 * x0 + x1 = 9 and x0 + 2 * x1 = 8
a = np.array([[3,1], [1,2]])
b = np.array([[9],[8]])
x = np.linalg.solve(a, b)
x

# lets create a 2x5 zeros matrix
c = np.zeros((2,5))
c

# With numpy, we can create arrays of numbers as well as functions that accept arrays
f = lambda x,y: np.sqrt(x**2+y**2)

f(3,4)

a = 3.1
b = np.arange(15,35, 0.75)
f(a,b)

# We can also perform numerical differentiation and integration
x = np.arange(-np.pi,np.pi, 0.01)
y = np.sin(x)*np.exp(-x)
plt.plot(x,y);
plt.xlabel('x')
plt.ylabel('$f(x)$')
plt.title('numerical integration example')
plt.show()

# Using the trapezoidal rule, we can numerically integrate the function
np.trapz(y,x)

"""
Explanation: $ 3 x_0 + x_1 = 9 $ 
$ x_0 + 2x_1 = 8 $
$ Ax=b  $
$ \left[\begin{matrix}3 & 1\1 & 2\end{matrix}\right]  \left[\begin{matrix}x_0\x_1\end{matrix}\right] = \left[\begin{matrix}9\8\end{matrix}\right] $
$ \left[\begin{matrix}x_0\x_1\end{matrix}\right] = \left[\begin{matrix}2\3\end{matrix}\right] $
End of explanation
"""


%reset

"""
Explanation: Symbolic Math with Sympy
in ipython use help(integrate) or integrate? to check syntax or learn about the command. To clear the previous modules use %reset to clear variables and imports.  It is also very import to not import over other modules or difficult errors to troubleshoot may appear. This is especially important with numpy and sympy
End of explanation
"""


import sympy as sp
import numpy as np

sp.init_printing()

x,y = sp.symbols('x y')

"""
Explanation: Now we can initialize the console to have very nice looking mathematical expressions. If we dont do this, the text will be simple ascii characters
End of explanation
"""


y = (sp.pi + x)**2
y

y.subs(x,1.3)

# or to get a number
sp.N(y.subs(x,1.3))

# calculus is also easy
sp.diff(y)


f = sp.Function('f')(x)

f = sp.sin(x)*sp.exp(-x)
f

F = sp.integrate(f)
F

F1 = sp.integrate(f,(x,-sp.pi,sp.pi ))
# as shown numerically, we can symbolically evaluate the function 
#  and get the same result
F1

float(F1)

sp.plot(f,(x,-sp.pi,sp.pi));

sp.plot(F,(x,-sp.pi,sp.pi));

# Other various math tasks can be much easier when using sympy

n = sp.symbols('n')

sum_fcn = sp.Sum(1/n**2, (n,1,10))
sum_fcn

# to evaluate the summation simply add the evalf method
sum_fcn.evalf()

sp.latex(sum_fcn)

"""
Explanation: I really enjoy using sympy because I can type my equations in, evaluate them, and it shows them very nice, which is easier to read then lines of text. Simple expressions are easy to make
End of explanation
"""


from IPython.display import Latex
Latex('$sum\_fcn=' + sp.latex(sum_fcn) + '$' )

# linear algebra

s11, n22, e33 = sp.symbols("sigma11 nu22 epsilon33")
A = sp.Matrix([[s11],[n22],[e33]])
A

m11, m12, m21, m22 = sp.symbols("m11, m12, m21, m22")
b1, b2 = sp.symbols("b1, b2")

A = sp.Matrix([[m11, m12],[m21, m22]])
A

b = sp.Matrix([[b1], [b2]])
b

A**2

A.det()

A.inv()

"""
Explanation: Since this is a windows machine, the latex command printed two backslashes, if we change that to single backslashes we can get the latex version of the summation 
$$
\sum_{n=1}^{10}\frac{1}{n^{2}}
$$
End of explanation
"""


x = sp.Symbol('x')   # x = var('x')
M = sp.Matrix([[2,x],[x,3]])
M

M.eigenvals()

M.eigenvects()

M.eigenvects()[1][0]

Mval = M.eigenvects()[1][0]
Mval.evalf(subs={x:3.14})

print(sp.latex(M))

from IPython.display import Latex
Latex('$M=' + sp.latex(M) + '$' )

"""
Explanation: mixed symbols and numerals linear algebra
End of explanation
"""


np.zeros((5,3))

np.array([[1,2],[3,4]])

np.array([[1,2],[3,4]])

# Matrix multiplication can be achieved using the dot method
i = [[1,0,0],[0,1,0],[0,0,1]]  # identiy matrix
a = [[4,3,1],[5,7,2],[2,2,2]]
np.dot(i,a)

#Or, matrix multiplication can be done if a matrix is explicitly defined
np.array(i) @ np.array(a)

# Notice, when arrays are mutliplied, we get the dot product 
m = np.array(i) * np.array(a)
m

m.T  # transpose

m**2

np.array(a)**2

m

m[:,2]

m[2,:]

m[:2,:2]

m[1:,1:]

"""
Explanation: Pythons list is a generic data storage object. it can be easily extended to a numpy array, which is specialized for numerical and scientific computation
End of explanation
"""


# import sympy library and initialize latex printing
import sympy as sp
#sp.init_printing()
#sp.init_printing(use_latex='matplotlib')
sp.init_printing(use_latex='mathjax')

# add a symbolic character
x = sp.Symbol('x')

sp.sqrt(x**2)

r = sp.Rational(11, 13)
r

float(r)

f = sp.Function('f')
f

f(x)

h = sp.Lambda(x,x**2)
h

w = 2*(x**2-x)-x*(x+1)
w

w.args

sp.simplify(w)

sp.factor(x**2-1)

#partial fractions
y = 1/(x**2+3*x+2)
y

sp.apart(y,x)

f = sp.Function('f')(x)
sp.diff(f,x)

y = sp.Symbol('y')
g = sp.Function('g')(x,y)
g.diff(x,y)

a,b,c,d = sp.symbols("a b c d")
M = sp.Matrix([[a,b],[c,d]])
M

M*M

# if ipython is to be used as a calculator initialize with 
from sympy import init_session
init_session() 

from sympy import oo, Function, dsolve, Eq, Derivative, sin,cos,symbols
from sympy.abc import x
import sympy as sp
import numpy as np
import matplotlib.pyplot as mp
get_ipython().magic('matplotlib inline')
# this will print output as unicode

# assign a sympy variable
x = sp.var('x')
x

#assign a function 
f =  sp.sin(6*x)*sp.exp(-x)
f

f.subs(x,3)

float(f.subs(x,3))

sp.plot(f,(x,-5,5))

# a onetime pretty print
sp.pprint(f)

#or we can print the latex rendering
sp.latex(f)

# first derivative
df = f.diff()
df

# differentaite f'' wrt x
sp.diff(f,x,1)

# substitute x with pi
f.subs(x,np.pi)

#%% Numeric Computation from the documentation
from sympy.abc import x

# lambdify using the math module, 10^2 faster than subs
expr = sp.sin(x)/x
f = sp.lambdify(x,expr)
f(3.14)

# lambdify using numpy
f = sp.lambdify(x,expr, "numpy")
f(np.linspace(1,3.14,20))

"""
Explanation: Symbolic mathematics with sympy
End of explanation
"""



## Signal Processing
#Page 174 Introduction for python for Science - David Pine

import numpy as np
from scipy import fftpack
import matplotlib.pyplot as plt
get_ipython().magic('matplotlib inline') # inline plotting
width = 2.0
freq = 0.5
t = np.linspace(-10, 10, 101) # linearly space time array
g = np.exp(-np.abs(t)/width)*np.sin(2.0 * np.pi * freq * t)
dt = t[1]-t[0] # increment between times in time array
G = fftpack.fft(g) # FFT of g
f = fftpack.fftfreq(g.size, d=dt) # frequenies f[i] of g[i]
f = fftpack.fftshift(f) # shift frequencies from min to max
G = fftpack.fftshift(G) # shift G order to coorespond to f
fig = plt.figure(1, figsize=(8,6), frameon=False)
ax1 = fig.add_subplot(211)
ax1.plot(t, g)
ax1.set_xlabel('t')
ax1.set_ylabel('g(t)')
ax2 = fig.add_subplot(212)
ax2.plot(f, np.real(G), color='dodgerblue', label='real part')
ax2.plot(f, np.imag(G), color='coral', label='imaginary part')
ax2.legend()
ax2.set_xlabel('f')
ax2.set_ylabel('G(f)')
plt.show()

from pylab import *
from scipy import fft
get_ipython().magic('matplotlib inline') # inline plotting
N = 2**9
F = 25
t = arange(N)/float(N)
x = cos(2*pi*t*F) + rand(len(t))*3
subplot(2,1,1)
plot(t,x)
ylabel('x []')
xlabel('t [seconds]')
title('A cosine wave')
grid()

subplot(2,1,2)
f = t*N
xf = fft(x)
plot(f,abs(xf))
title('Fourier transform of a cosine wave')
xlabel('xf []')
ylabel('xf []')
xlim([0,N])
grid()
show()

# note the spike at 25 hz and 512-25

"""
Explanation: Signal Processing
End of explanation
"""


# see here for example scripts
# C:\Users\Neal\Anaconda3\Lib\site-packages\sympy\mpmath\tests
from sympy import Function, dsolve, Eq, Derivative, sin,cos,symbols
from sympy.abc import x
import numpy as np
import sympy as sp
import matplotlib.pyplot as mp
import matplotlib.pyplot as plt

f = Function('f')
deq = dsolve(Derivative(f(x), x,x) + 9*f(x), f(x))
deq

diffeq1_sym = deq.args[1]
diffeq1_sym

diffeq1 = diffeq1_sym.subs({'C1':1, 'C2':0.5})
diffeq1

diffeq1_f = sp.lambdify(x,diffeq1, "numpy")
diffeq1_f

diffeq1arr = diffeq1_f(np.linspace(1,3.14,20))
diffeq1arr

# numpy plot
plt.plot(diffeq1_f(np.linspace(-10,10,2000)));
plt.title('plot of the numpy array');

# sympy plot
sp.plot(diffeq1, title='plot of the sympy function');

"""
Explanation: Differential Equations
End of explanation
"""


# quiver plot
f=lambda x: [x[0]**2 - 2*x[0] - x[1] + 0.5, x[0]**2 + 4*x[1]**2 - 4]
x,y=np.mgrid[-0.5:2.5:24j,-0.5:2.5:24j]
U,V=f([x,y])
plt.quiver(x,y,U,V,color='r', \
         linewidths=(0.2,), edgecolors=('k'), \
         headaxislength=5)
plt.show()

"""
Explanation: Differential Equations Quiver plot
End of explanation
"""


