Introduction to SciPy Tutorial
This tutorial is an introduction SciPy library and its various functions and utilities. The objective of this tutorial is to give a brief idea about the usage of SciPy library for scientific computing problems in Python. The SciPy library has several toolboxes to solve common scientific computing problems. It contains submodules for applications like Integration, Interpolation, Image Processing, Optimization, Special Functions and Statistics, etc. SciPy works efficiently on NumPy arrays and is standard scientific computing library in Python. SciPy library is composed of sub-modules designed for specific tasks. A concise list of SciPy sub-modules is shown below:
Fast Fourier Transform
Optimization and Fits
Since NumPy is required for most of these sub-modules, this is how standard SciPy modules are imported:
>>> import numpy as np >>> from scipy import signal #Replace ‘signal’ with the name of other sub-modules
The core SciPy namespace share common functions with NumPy module, but SciPy offers an added advantage in terms of scientific computing functions.
File Input/Output in SciPy
- To load text files you can use
Sometimes you may require to import an image, you can use misc.imread() function in the following syntax to load the image.
>>>from scipy import misc >>>misc.imread(‘Image_Name.png’)
Linear Algebra Operations in SciPy
The module for standard Linear Algebra operations is known as scipy.linalg and you are required to import it well before any operation using the following command operation:
>>>from scipy import linalg
- To calculate the determinant of a square matrix, we will use scipy.linalg.det() function in the following way:
>>>mat = np.array([[2,1],[4.3]]) #For a square matrix ‘mat’ >>>linalg.det(mat) 2.0
- Scipy.linalg.inv() is used to compute inverse of a square matrix. The syntax for the function is mentioned below:
>>>linalg.inv(sqr_mat) #For a square matrix ‘sqr_mat’
- Scipy.linalg also supports complex matrix utilities like Singular-Value Decomposition, Cholesky Decomposition, QR, LU, etc.
The syntax for SVD is as follows:
>>> linalg.svd(mat) #Returns 3 arguments after Singular Value Decomposition
Special Functions in SciPy
Scipy.special module contains list of transcendental function, which are most often used in mathematical operation across various disciplines. Here is the syntax for some of the most used function from the scipy.special modules are:
# To calculate the area under a Guassian curve, we use erf() function like this: scipy.special.erf() # Syntax for Gamma function: scipy.special.gamma() # In order to calculate the log of Gamma, we use the following syntax: scipy.special.gammaln() # Elliptic Function scipy.special.eppilj() # Nth order Bessel Function scipy.special.jn()
Optimization Functions in SciPy
Optimization is a mathematical problem of estimating a numerical solution of variables that follow a certain equation. SciPy offers module which provides algorithms for function minimization, root finding, curve fitting, etc. This modules is known as scipy.optimize and can be imported using the following command:
>>> from scipy import optimize
The subsequent sections, will cover various examples of functions available in this module for mathematical optimizations.
- Evaluating minimum of a Scalar Function.
Let’s assume a function f(x), which will be used to estimate minima using functions in the scipy.optimize module.
Note: f(x) is an arbitrary function.
We can plot the above mentioned function using the commands mentioned below:
>>>def f(x): #Defining function … return x**3 + x**2 + np.sin(x) +np.cos(x) >>> x = np.arrange(-50, 50, 0.01) #First two arguments are the limit and the last argument is the interval >>> plt.plot(x, f(x)) >>> plt.show() #To show the graph of the defined function
One of the efficient ways to calculate minimum of this function is to use the BFGS algorithm with a starting point. This algorithm calculates a gradient descent of the function from the starting point given in the argument and outputs the minima with zero gradient and positive second order derivative. Here’s the syntax to use BFGS function:
>>>optimize.fmin_bfgs(f, o) #First argument is function’s name and second argument is the starting point of gradient descent
This algorithm is ideal for a function with single minima or if we wish to calculate local minima of the function. However in a situation when the range of the global minima is unknown, therefore a starting point can’t be picked for BFGS function, which will definitely give us the global minima. For situations like these, scipy.optimize module has basinhopping() function which combines a local optimizer with stochastic sampling of starting points for the local optimizer, hence giving a costlier global minimum. The syntax to use this function is as follows:
>>> optimize.basinhopping(f, 0)
Brute force method can also be used for global optimization, but it is less efficient. The syntax for brute force method will be:
So far we have talked about calculating global optimization, however SciPy also has function which enables us to find the local minimum within an interval for variables, using fminbound() function. The syntax for using fminbound() function is mentioned below:
>>>optimize.fminbound(f, -50, 50) #Argument is interval for the variable.
- Roots Estimation for Scalar Function
In the previous section we calculated the global and local minima of the function f(x). In this section we will learn the method to estimate the roots of a scalar function by finding the solution of the equation; f(x) = 0. The syntax of using fsolve() function is as follows:
>>>roots = optimize.fsolve(f)
- Curve fitting for data points
Let’s say you have a data sample and you need to estimate the curve/function which was used to create those sampled data points. Scipy.optimize modules has curve_fit() function, which doesn the job by estimating variables of the function using least squares curve fitting.
For the sake of this example, let’s use the following function g(x):
Now we create the sampled data set with some random noise, using the following arguments:
>>>x_data = np.linspace(-100, 100, 0.1) >>>y_data = g(x_data)+ np.random.randn(x_data.size)
Note: Functions like linspace and random noise generator are called from NumPy Library.
In order to do the curve fitting, we need to know the initial functional form, therefore let’s assume a function which is similar to g(x) but with unknown variables.
>>>def h(x, a, b): … return a*x + b*np.cos(x) # We will now use the curve_fit() function to estimate the values of a & b >>> initial_guess_ab = [1, 1] >>> variables, variables_covariance = optimize.curve_fit(h, x_data, y_data, initial_guess_ab)
The argument ‘variables’ contains an approximation of the variables of function g(x).
Fourier Transformation in SciPy
Fourier Transformation is computed on a time domain signal to check its behaviour in frequency domain. Fourier transformation finds its application in disciplines like signal and noise processing, image processing, audio signal processing, etc. SciPy offers the fftpack module which lets the user compute fast Fourier transforms. Here’s an example of a sine function which will be used to calculate Fourier transform using the fftpack module.
>>>signal = np.sin(x +10) >>>x = np.linspace(-100, 100, 0.1) >>>time_step = 0.01
Since the signal is a real function, its fourier transform will be symmetric. The fftfreq() function is required to estimate the sampling frequencies and the fft() function will generate the Fast Fourier transform of the signal. The syntax to compute FFT is as follows:
>>>from scipy import fftpack >>>sampling_frequency = fftpack.fftfreq(signal.size, d=time_step) >>>signal_fft = fftpack.fft(signal)
Similarly an inverse Fourier transform can be computed, as shown below:
>>>original_signal = fftpack.ifft(signal_fft)
Numerical Integration in SciPy
Numerical Integration or Definite Integration is an important utility in scientific computation. Scipy.integrate module offers various integration schemes such as quad(), quadrature(), fixed_quad(), romberg(). Let’s say we had to integrate sin function from 0 degrees to pi/2, here’s the syntax to do that:
>>>from scipy.integrate import quad >>>result, error = quad(np.sin, 0, np.pi/2)
This module also offers functions for integrating ODEs (Ordinary Differential Equations). The function, odeint() is a general purpose integrator which can solve multi-order differential equations. In this section, we will take the example of a first order ODE;
The initial condition for x(t=0) = -1 (Let’s assume).
Before we integrate this let’s define the function for which we are computing the derivative:
>>>def derivative_equation(x, time) … return 10*time
Here’s the syntax to solve the above mentioned ODE:
>>>from scipy.integrate import odeint >>>time_step = np.linspace(0, 10, 100) #Time step is 0.1 from [0, 100) >>>x_value, info = odeint(derivative_equation, 1, time_step)
The output of this function will be following equation; x = 5*t^2, integrated over the interval of 0 to 10.
Image Processing in SciPy
The module for Image processing in SciPy is known as scipy.ndimage. The routines offered in this modules are capable of application like geometrical transformation of images, including changes in resolution, orientation, shape, etc; Image filtering using Gaussian, Weiner, Median, and other such filters; erosion, dilation, opening, closing of images can also be done using this module. Here’s an extensive list of functions from scipy.ndimage module:
>>>from scipy import ndimage >>>ndimage.shift(image, (x, y)) #Shifting Image with (x, y) coordinate >>>ndimage.rotate(image, angle) #Rotating image to that angle >>>ndimage.zoom(image, magnitude) #Zooming image with the magnitude >>>ndimage.median_filter(image, argument) #Filtering image using Median filter >>>ndimage.gaussian_filter(image, argument) #Filtering image using Gaussian filter >>>ndimage.binary_erosion(image) #Binary Erosion on image >>>ndimage.binary_dilation(image) #Binary Dilation on image >>>ndimage.binary_opening(image) #Binary Opening on image >>>ndimage.binary_closing(opened_image) #Binary Closing on opened image
Signal Processing in SciPy
SciPy offers various function for signal processing applications like, Fourier Transform, Filtering, Window Function, Sampling, Curve Fitting, Removing Linear trend, etc. Some of these functions are mentioned below:
- Resampling using Fourier transform
Let’s say we have a sine function which needs to be resampled to n points within the time domain interval. Scipy.signal module has a function called, resample(), which uses FFT to do the same. The syntax to resample the function is mentioned below:
>>>t = np.linspace(-10, 10, 200) #Defining Time Interval >>>y = np.sin(t) >>>signal.resample(y, 100) # Number of required samples is 100
- Removing Linear Trend
Let’s say we have a function which have a linear element in the equation. Scipy.signal has detrend() function which can be used to eliminate the linear element from the signal and give us the transient solution. The syntax to use the detrend() function is mentioned below:
>>>t = np.linspace(-10, 10, 200) #Defining Time Interval >>>y = np.sin(t) + t >>>signal.detrend(y) # To remove the linear‘t’ variable in the equation