<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 <center>
    <img src="http
Step1:  <div id='intro' />
Introduction
Back to TOC
In this jupyter notebook we will briefly introduce part of ecosystem in Python for Scientific Computing with explanation and examples.
We strongly suggest to follow the explanation step by step, read the suggested documents/links, and play with the jupyter notebook changing it to learn!
<div id='numpy' />
NumPy, the core of numerical computing in Python
Back to TOC
We strongly suggests to take a look to the associated paper explaning Arrays in NumPy link.
Now, we will start with some examples provided in Fig. 1 from the previous reference. Also see https
Step2:  b) Indexing (view)
Step3:  c) Indexing (copy)
Step4:  d) Vectorization (THE HEART OF THE HEART OF NUMERICAL COMPUTING!!)
Step5:  e) Broadcasting
Step6:  f) Reduction
Step7:  g) Examples
Step8:  <div id='matplotlib' />
MatPlotLib, comprehensive 2-D plotting
Back to TOC
In this case, we will go straight to the examples. See https
Step9:  <div id='scipy' />
SciPy, fundamental library for scientific computing
Back to TOC
In this case, we will go straight to the examples. See https
Step10:  <div id='sympy' />
SymPy, symbolic mathematics
Back to TOC
In this case, we will go straight to the examples. See https
<ASSISTANT_TASK:>
Python Code:
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
import sympy as sym
from scipy.integrate import odeint
import scipy.linalg as spla
from scipy import optimize, interpolate
from scipy.linalg import solve_triangular, toeplitz, lu
from scipy.optimize import root
# pip install bitstring
import bitstring as bs
%matplotlib inline
sym.init_printing()
from ipywidgets import interact, widgets, fixed, IntSlider, interact_manual, Checkbox, RadioButtons
from mpl_toolkits.mplot3d import Axes3D
import matplotlib as mpl
mpl.rcParams['font.size'] = 14
mpl.rcParams['axes.labelsize'] = 20
mpl.rcParams['xtick.labelsize'] = 14
mpl.rcParams['ytick.labelsize'] = 14
from matplotlib import cm
# pip install memory_profiler
%load_ext memory_profiler
from sklearn import datasets
from colorama import Fore, Back, Style
# https://pypi.org/project/colorama/
# Fore: BLACK, RED, GREEN, YELLOW, BLUE, MAGENTA, CYAN, WHITE, RESET.
# Back: BLACK, RED, GREEN, YELLOW, BLUE, MAGENTA, CYAN, WHITE, RESET.
# Style: DIM, NORMAL, BRIGHT, RESET_ALL
textBold = lambda x: Style.BRIGHT+x+Style.RESET_ALL
textBoldH = lambda x: Style.BRIGHT+Back.YELLOW+x+Style.RESET_ALL
textBoldB = lambda x: Style.BRIGHT+Back.BLUE+Fore.BLACK+x+Style.RESET_ALL
textBoldR = lambda x: Style.BRIGHT+Back.RED+Fore.BLACK+x+Style.RESET_ALL
Explanation: <center>
    <img src="http://sct.inf.utfsm.cl/wp-content/uploads/2020/04/logo_di.png" style="width:60%">
    <h1> INF-285 - Computación Científica </h1>
    <h2> The Beginning </h2>
    <h2> <a href="#acknowledgements"> [S]cientific [C]omputing [T]eam </a> </h2>
    <h2> Version: 1.02 </h2>
</center>
<div id='toc' />
Table of Contents
Introduction
NumPy
MatPlotLib
SciPy
SymPy
Libraries
Acknowledgements
End of explanation
# NumPy array
a1 = np.arange(12)
# NumPy ndarray
a2 = np.array([[1, 2, 3], [4, 5, 6]], float)
a3 = np.reshape(a1,(4,3))
print(textBoldB(' Printing their values '))
print(textBoldH('a1:\n'),a1)
print(textBoldH('a2:\n'),a2)
print(textBoldH('a3:\n'),a3)
print(textBoldB(' Printing their datatypes '))
print(textBoldH('a1:\n'),a1.dtype)
print(textBoldH('a2:\n'),a2.dtype)
print(textBoldH('a3:\n'),a3.dtype)
print(textBoldB(' Printing their type name '))
print(textBoldH('a1:\n'),type(a1).__name__)
print(textBoldH('a2:\n'),type(a2).__name__)
print(textBoldH('a3:\n'),type(a3).__name__)
Explanation: <div id='intro' />
Introduction
Back to TOC
In this jupyter notebook we will briefly introduce part of ecosystem in Python for Scientific Computing with explanation and examples.
We strongly suggest to follow the explanation step by step, read the suggested documents/links, and play with the jupyter notebook changing it to learn!
<div id='numpy' />
NumPy, the core of numerical computing in Python
Back to TOC
We strongly suggests to take a look to the associated paper explaning Arrays in NumPy link.
Now, we will start with some examples provided in Fig. 1 from the previous reference. Also see https://numpy.org.
a) Data structure: An array and a bi-dimensional array.
End of explanation
print(textBoldB(' Showing original variable a3 '))
print(textBoldH('a3\n'),a3)
print(textBoldB(' Showing a subset of elements (SLICES) of variable a3 \n the second and third columns '))
print(textBoldH('a3[:,1:]\n'),a3[:,1:])
print(textBoldB(' Storing the outcome in a4 '))
a4 = a3[:,1:]
print(textBoldB(' What does it means it is a "view"? '))
print(textBoldB(' Let\'s change the original a3 and see what happens to a4. '))
a3[0,1]=100
print(textBoldH('a4\n'), a4)
print(textBoldB(' a4 has changed too! So a4 is a "view" of part of a3. '))
print('\n')
print(textBoldB(' So, how do we copy it?'))
print(textBoldB(' Let\'s use np.copy to copy the "view" obtained before,'+
                '\n then modify the original variable and see what had happened'))
a5 = np.copy(a3[:,1:])
print(textBoldH('a3\n'), a3)
print(textBoldH('a4=a3[:,1:])\n'), a4)
print(textBoldH('a5=np.copy(a3[:,1:])\n'), a5)
a3[0,1]=1
print(textBoldB(' Changing a3 as "a3[0,1]=1"'))
print(textBoldH('a3\n'), a3)
print(textBoldH('a4\n'), a4)
print(textBoldH('a5\n'), a5)
print(textBoldB(' Thus, chaning a3 have changed a4 since it is a "view" and not changed a5 since it is a copy. '))
print(textBoldB(' Another way to extract slices is with steps: '))
a6 = a3[:,::2]
print(textBoldH('a3\n'), a3)
print(textBoldH('a6 = a3[:,::2]\n'), a6)
print(textBoldB(' In this case we obtained the first and third columns, \n '+
               ' the pattern consider ini:end:step, where "ini" is the first index considered,\n'+
               ' "end" defines that the pattern will be considere until the index "end-1",\n'+
               ' and "step" indicates the step for selecting the next index, it is 1 by default.'))
Explanation: b) Indexing (view)
End of explanation
print(textBoldB(' Now we will list the alternatives to "copy" the data from ndarrays'))
print(textBoldH('a3\n'), a3)
print('\n')
print(textBoldR('with SCALARS:'))
print(textBoldB('o1 = a3[1,2]'))
o1 = a3[1,2]
print(textBoldH('o1\n'), o1)
print(textBoldB('Changing o1: o1 = 111'))
o1 = 111
print(textBoldH('a3\n'), a3)
print(textBoldH('o1\n'), o1)
print('\n')
print(textBoldR('with MASKS:'))
print(textBoldB('o2 = a3[a3>9]'))
o2 = a3[a3>9]
print(textBoldH('o2\n'), o2)
print('\n')
print(textBoldR('with ARRAYS:'))
print(textBoldB('o3 = a3[[0,1],[1,2]]'))
o3 = a3[[0,1],[1,2]]
print(textBoldH('o3\n'), o3)
print(textBoldH('o3:\n'),o3.dtype)
print(textBoldH('o3:\n'),type(o3).__name__)
print('\n')
print(textBoldR('with ARRAYS and BROADCASTING (this will be review next) :'))
print(textBoldB('o4 = a3[[[1],[2]],[1,0]]'))
o4 = a3[[[1],[2]],[1,0]]
print(textBoldH('o4\n'), o4)
print(textBoldH('o4:\n'),o4.dtype)
print(textBoldH('o4:\n'),type(o4).__name__)
Explanation: c) Indexing (copy)
End of explanation
print(textBoldH('a3\n'), a3)
print(textBoldB(' We will define a new ndarray n1 = np.ones((4,3)) and add it to the array a3.'))
n1 = np.ones((4,3))
print(textBoldB(' The first implementation one may think of is the following:'))
def sum_two_arrays(a3,n1):
    out = np.zeros_like(a3) # What does this function do?
    for i in range(4):
        for j in range(3):
            out[i,j] = a3[i,j]+n1[i,j]
    return out
n2 = sum_two_arrays(a3,n1)
# To quantify the computation time it takes let's measure it:
%timeit n2 = sum_two_arrays(a3,n1)
# The output it shows in the computer we run the experiments was:
# 33.9 µs ± 1.68 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
# This means it took about 34 microseconds.
# Now, let's try the NumPy way! i.e. the vectorized approach!
%timeit n2 = a3+n1
# The output it shows in the computer we run the experiments was:
# 1.05 µs ± 12.7 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)
# This means it was about 34 times faster!! 
# and also it only took "one" line of code where in the previous version took "7" lines of code!!
print(textBoldR(' Given the previous argument, it is crucial you get familir with VECTORIZATION computing with NumPy!! '))
Explanation: d) Vectorization (THE HEART OF THE HEART OF NUMERICAL COMPUTING!!)
End of explanation
# Creating a 1D array, this is the stringly suggested way to create "vectors" is NumPy.
# This can be considered as a row vector.
x1 = np.array([1,2,3])
print(textBoldH('x1:'),x1)
print(textBoldH('shape:'),x1.shape)
print(textBoldH('dtype:'),x1.dtype)
print(textBoldH('type name:'),type(x1).__name__)
# Creating a 2D array with one row.
x2 = x1[np.newaxis] # What does np.newaxis do?
print(textBoldH('x2:'),x2)
print(textBoldH('shape:'),x2.shape)
print(textBoldH('dtype:'),x2.dtype)
print(textBoldH('type name:'),type(x2).__name__)
# This emulates a column vector.
# Transpossing the 2D array. This emulates a column vector.
x2 = x1[np.newaxis].T
print(textBoldH('x2:'),x2)
print(textBoldH('shape:'),x2.shape)
print(textBoldH('dtype:'),x2.dtype)
print(textBoldH('type name:'),type(x2).__name__)
# Notice that transpossing a 1D array does not do anything (please try it by yourself!).
# Now that we have a 1D array (that emulates a row vector) and 
# a 1 column 2d array (that emulates a columns vectors), we can do a "BROADCAST" and see what happend.
x3=x1+x2
print(textBoldH('x3:'),x3)
print(textBoldH('shape:'),x3.shape)
print(textBoldH('dtype:'),x3.dtype)
print(textBoldH('type name:'),type(x3).__name__)
# In this cas we have operated elementwise and broadcasted the array.
# The previous operation is equivalent to the following operation.
x11 = np.tile(x1, (3,1)) # Repeat the 1d array 'vertically' 3 times
print(textBoldH('x11:'),x11)
x21 = np.tile(x2, 3) # Repeat the 2d array 'horizontally' 3 times 
print(textBoldH('x21:'),x21)
x31 = x11+x21
print(textBoldH('x31:'),x31)
print(textBoldH('shape:'),x31.shape)
print(textBoldH('dtype:'),x31.dtype)
print(textBoldH('type name:'),type(x31).__name__)
print(textBoldR(' As you may suspect, this way of doing the sum is not recommended.\n'+
               ' It was added just for comparison purposes.'))
Explanation: e) Broadcasting
End of explanation
print(textBoldH('a3\n'), a3)
print('\n')
print(textBoldB(' First reduction, sum along axis 0:'))
print(textBoldB(' s1 = np.sum(a3,axis=0)'))
s1 = np.sum(a3,axis=0)
print(textBoldH('s1'), s1)
print(textBoldH('shape:'),s1.shape)
print(textBoldH('dtype:'),s1.dtype)
print(textBoldH('type name:'),type(s1).__name__)
print('\n')
print(textBoldB(' First reduction, sum along axis 1:'))
print(textBoldB(' s2 = np.sum(a3,axis=1)'))
s2 = np.sum(a3,axis=1)
print(textBoldH('s2'), s2)
print(textBoldH('shape:'),s2.shape)
print(textBoldH('dtype:'),s2.dtype)
print(textBoldH('type name:'),type(s2).__name__)
print('\n')
print(textBoldB(' First reduction, sum along axis 0 and 1:'))
print(textBoldB(' s3 = np.sum(a3)'))
s3 = np.sum(a3)
print(textBoldH('s3'), s3)
print(textBoldH('shape:'),s3.shape)
print(textBoldH('dtype:'),s3.dtype)
print(textBoldH('type name:'),type(s3).__name__)
print('\n')
print(textBoldR(' As you may suspect, the reduction reduces the dimensionality in respect of the original input.'))
Explanation: f) Reduction
End of explanation
# Create a random matrix and comput its determinant
np.random.seed(0) # Fixing seed for randomnumber generator.
n = 100
A = np.random.random((n,n))
print(textBoldH('|A|:'), np.linalg.det(A))
# Compute matrix vector product
x = np.random.random(n)
b = A @ x
print(textBoldH('A*x:'), b)
# Solving a linear system of equations A*x=b, i.e. assume we know A and b and we want x.
x_out = np.linalg.solve(A,b)
print(textBoldH('x_out=A^{-1} b:'), x_out)
# Computing the norm 2, between x and x_out
print(textBoldH('||x-x_out||_2:'), np.linalg.norm(x-x_out))
# It seems it is small, why is this?
# Computing dot product
v1 = np.random.random(n)
v2 = np.random.random(n)
print(textBoldH('dot product between v1 and v2:'), np.dot(v1,v2))
# Computing elementwise product
print(textBoldH('v1*v2:'), v1*v2)
# Squaring each element of v1
print(textBoldH('out:'), np.power(v1,2))
Explanation: g) Examples
End of explanation
def plot_trig_func(n=5):
    x = np.linspace(0,10,n) # Generates a linearly spaced set of n points from 0 to 10
    y1 = np.sin(x) # Vectorization in action!
    y2 = np.cos(x) # Vectorization in action!
    plt.figure(figsize=(5,5)) # Creating the figure
    plt.plot(x,y1,label=r'$\sin(x)$') # Plotting sin
    plt.plot(x,y2,'r.-',label=r'$\cos(x)$') # Plotting cos
    plt.grid(True) # Adding nice grid
    plt.legend(loc='best') # Adding legend
    plt.title('Plotting trigonometric functions') # Adding title
    plt.xlabel(r'$x$') # Adding x-label
    plt.ylabel(r'$y$') # Adding y-label
    plt.show() # Showing the plot
# The following is a iPython widget called interact, 
# it allows us to interact with the function that it is being called.
interact(plot_trig_func,n=(3,1000,1)) # Making this 'tunable'!
def plot_funct_2_var(theta=0):
    x = np.linspace(-2,2,20)
    y = x # Recall this is just a 'view'
    X, Y = np.meshgrid(x,y) # It generates two 2D arrays with with cartesian product of 'x' and 'y'.    
    # 'lambda' definition help us to organize better the code
    f = lambda x,y: np.sin(x*np.sin(theta)+y*np.cos(theta))*np.exp(-(x**2+y**2))
    Z = f(X,Y) # Evaluation a function of two variables
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10,5))
    fig.suptitle('Comparison of two ways of plotting $f(x,y)$')
    # This is useful for images.
    ax1.imshow(Z, cmap='hot', interpolation='nearest')
    ax1.set_ylabel('imshow')
    # This is useful for functions in two variables, i.e. z=f(x,y)
    ax2.pcolormesh(X,Y,Z,cmap='hot', shading='auto')
    ax2.axis('equal')
    ax2.set_ylabel('pcolormesh')
    plt.show()
interact(plot_funct_2_var,theta=(0,2*np.pi,2*np.pi/50))
def plot_funct_2_var_3D(theta=0):
    x = np.linspace(-2,2,40)
    y = x # Recall this is just a 'view'
    X, Y = np.meshgrid(x,y) # It generates two 2D arrays with with cartesian product of 'x' and 'y'.    
    # 'lambda' definition help us to organize better the code
    f = lambda x,y: np.sin(x*np.sin(theta)+y*np.cos(theta))*np.exp(-(x**2+y**2))
    Z = f(X,Y) # Evaluation a function of two variables
    fig, ax = plt.subplots(subplot_kw={"projection": "3d"},figsize=(10,10))
    # Plot the surface.
    surf = ax.plot_surface(X, Y, Z, cmap='hot', linewidth=2, antialiased=True)
    # Add a color bar which maps values to colors.
    fig.colorbar(surf, shrink=0.5, aspect=10)
    plt.show()
interact(plot_funct_2_var_3D,theta=(0,2*np.pi,2*np.pi/50))
Explanation: <div id='matplotlib' />
MatPlotLib, comprehensive 2-D plotting
Back to TOC
In this case, we will go straight to the examples. See https://matplotlib.org and https://github.com/matplotlib/cheatsheets#cheatsheets.
End of explanation
# Numerical Integration
import scipy.integrate as integrate
# Integrating $\int_0^1 \exp(-x^2)\,dx$
result = integrate.quad(lambda x: np.exp(-x**2), 0, 1)
print(result)
# Integrating $\int_0^1 1\,dx=1$
result = integrate.quad(lambda x: 1, 0, 1)
print(result)
# Minimization
from scipy.optimize import minimize
f = lambda x: np.sin(x[0])*(1+np.cos(x[1]))
x0 = np.array([0, 0])
res = minimize(f, x0, method='nelder-mead', options={'xatol': 1e-8, 'disp': True})
print(res.x)
# Plotting the previous function
x = np.linspace(-2*np.pi,2*np.pi)
y = x
X,Y = np.meshgrid(x,y)
f = lambda x,y: np.sin(x)*(1+np.cos(y))
Z = f(X,Y)
plt.figure(figsize=(8,8))
plt.pcolormesh(X,Y,Z,shading='auto')
plt.plot(res.x[0],res.x[1],'r.')
plt.colorbar()
plt.title('The red dot is the minimum found numerically')
plt.show()
# root finding
from scipy.optimize import root
f = lambda x: x - 0.5 * np.exp(x)
sol = root(f, 0.3)
print(sol)
Explanation: <div id='scipy' />
SciPy, fundamental library for scientific computing
Back to TOC
In this case, we will go straight to the examples. See https://www.scipy.org.
End of explanation
# Working with algebraic expression versus 'numerical' expressions
sym.sqrt(3)==np.sqrt(3)
sym.sqrt(3)>np.sqrt(3)
sym.sqrt(3)<np.sqrt(3)
# Defining symbolic variables
x = sym.symbols('x')
x
y=x**2+sym.cos(x)+(x**2+1)**2
y
# Expanding the terms
sym.expand(y)
# Computing derivatives
yp=sym.diff(y)
yp
# Factorization of expression
sym.factor(yp)
# Computing limits
sym.limit((1-sym.cos(x))/x, x, 0)
sym.limit((sym.sin(x))/x, x, 0)
# Solving algebraic equations
sym.solve(x**2-x+1,x)
Explanation: <div id='sympy' />
SymPy, symbolic mathematics
Back to TOC
In this case, we will go straight to the examples. See https://www.sympy.org/en/index.html.
End of explanation
<END_TASK>