Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Created for the CDS 110 course at Caltech but which others may find useful.
#!/usr/bin/env python
#
# Copyright (c) 2012, 2015 by California Institute of Technology
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions
# are met:
#
# 1. Redistributions of source code must retain the above copyright
# notice, this list of conditions and the following disclaimer.
#
# 2. Redistributions in binary form must reproduce the above copyright
# notice, this list of conditions and the following disclaimer in the
# documentation and/or other materials provided with the distribution.
#
# 3. Neither the name of the California Institute of Technology nor
# the names of its contributors may be used to endorse or promote
# products derived from this software without specific prior
# written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL CALTECH
# OR THE CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
# SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
# LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
# USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
# ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
# OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
# SUCH DAMAGE.
"""
CDS 110a; Fall 2012.
Operation "Python rising: CDS110a".
SCL <slivingston@cds.caltech.edu>
Original release: 2012-10-05
Changes:
2012-10-08 13:40: - minor bugfixes;
- elaborated comments.
2015-01-09 10:12 - print is a function;
- update versions under "Getting help on the Web";
- add link to IPython;
- update author's email address.
2015-02-05 20:24 - bugfix ODE integration example;
thanks to Sean Sanchez for reporting it.
"""
## Getting help on the Web
## (may need to change version numbers...)
#
# https://docs.python.org/2.7/
# http://ipython.org/documentation.html
# http://docs.scipy.org/doc/numpy/reference/
# http://docs.scipy.org/doc/scipy/reference/
# http://matplotlib.org/contents.html
# http://python-control.sourceforge.net/manual/
## From the Python prompt,
#
# Note that while Python has a tolerable built-in prompt, a much
# better tool for interactive computing is IPython <http://ipython.org>.
#
# you can usually use the "help" command invoked with the name of the
# function or class (or module, or any object with a docstring) to see
# handy documentation. E.g., to learn about the str class, you could type
#
# >>> help(str)
#
# The documentation viewer recognizes most commands taken by the
# "less" command-line utility; they are similar to those for movement
# in Vi, e.g., type "F" to move "forward" one screenful, "B" to move
# "backward", "J" to move down one line, "K" to move up one line.
#
# For NumPy and SciPy (and packages using the same documentation
# style), you should use numpy.info instead of the "help" command.
# For instance, for the "square root" function,
#
# >>> import numpy as np
# >>> np.info(np.sqrt)
## From the terminal,
#
# you can see documentation as would be returned by the "help" command
# (see above) using the "pydoc" utility. E.g., you should see the
# same as in the above string class example by entering
#
# $ pydoc str
from __future__ import print_function
# Standard imports
import numpy as np
import scipy as sp
import scipy.linalg as la
import matplotlib.pyplot as plt
# Some specifics for our needs
from scipy.integrate import odeint
if __name__ == "__main__":
# From here below is mostly ported from "MATLAB_Review.m", which
# was written by Robert Karol in Fall 2011. This is done in the
# hope of helping those of you already familiar with MATLAB and
# want to kick the habit.
# Scientific format also works.
print(1.743e5)
print(4.23e-2)
## Declaring and using Variables
# Variables in MATLAB do not need to be declared as in many computer
# languages. It can distinguish integers, floats, and strings.
# MATLAB is case sensitive.
a = 9
A = 1/3.
b = 'Hello'
print("a = "+str(a)) # This is the more modern approach.
print("a =", a) # Comma-separated list works, too.
print("a = %d; A = %f" % (a, A)) # printf()-style has fallen out of favor.
# Exponentiation
print("="*60)
print("Exponentiation")
print("-"*34)
print(a**A)
print(np.exp(a)) # e^a
print(np.log(100)) # uses natural log
print(np.log10(100)) # log base 10
# Trig
print("="*60)
print("Trig")
print("-"*34)
print(np.pi)
print(np.sin(np.pi)) # Uses radians
print(np.cos(np.pi))
# Imaginary numbers
print("="*60)
print("Imaginary numbers")
print("-"*34)
print(np.exp(3.1415926j))
print(np.real(np.exp(3.1415926j)))
print(np.imag(np.exp(3.1415926j)))
## Arrays
# Matrices
print("="*60)
print("Matrices")
print("-"*34)
w = np.array([1, 2, 3.]) # row vector. Commas separate elements.
# Notice decimal after "3"; this ensures
# floating point representation.
y = np.array([[1],
[2],
[3]], # column vector.
dtype=np.float64) # Explicit type declaration.
print(y)
# While the line-breaks are not necessary, they make matrices
# *much* easier to read.
# 3 by 3 matrix:
a = np.array([[1, 2, 3],
[4, 5, 6],
[7, 8, 9.]])
print(a)
b = np.arange(1, 11) # Create vector of values 1..10 (the last
# value is *not* included).
c = np.arange(1, 11, 2) # Last element is step size (default 1).
e = np.linspace(1,6,25) # Generates a vector with 25 elements equally spaced between 1 and 6 inclusive.
## Matrix Operations
print("="*60)
print("Matrix Operations")
print("-"*34)
A = np.arange(1, 10)
# Reshape command used to change dimensions of the matrix.
A = A.reshape(3, 3)
print("A =", A)
B = np.arange(9, 0, -1)
B = B.reshape(3,3)
print("B =", B)
print(A + B) # Matrices must be the same size, to be added
print(A + 3) # Adding a scalar to a matrix adds it at ever element.
print(A - B) # Subtraction works as addition. Same size or scalar.
print(np.dot(A, B)) # Use numpy.dot() for matrix multiplication
print(A * B) # Elementwise Multiplication, every element of A is
# multiplied by the corresponding element in B.
print(A / B) # A(i,j)/B(i,j) Same size or scalar
print(3 ** B) # Matrix Exponentiation
print(B ** 3)
print(A ** B) # A(i,j) to the B(i,j)
## Array initialization
A = np.ones((5,4)) # Matrix of all ones, size 5 by 4.
A = np.zeros(3) # N.B., this will make vector with 3 elements.
A = np.diag([1, 2, 3, 4, 5, 6, 7])
# NumPy has quite powerful "index slicing" for accessing parts of ndarrays.
x = A[0,:] # Extract the first row, with all columns
y = A[:,2]
a = A[:3,:]
b = A[0:5:2, 1:7:3] # 1,3,5 rows and 2, 5 columns
x[0]
y[-1] # Index -1 is similar to "end" in MATLAB.
## Data structures and operations: arrays
# Create random arrays.
print("="*60)
print("Create random arrays")
print("-"*34)
A = np.random.randn(2,2)
print(A)
B = np.random.rand(2,3)
print(B)
# Concatenate matrices horizontally.
print("="*60)
print("Concatenate matrices horizontally.")
print("-"*34)
C = np.hstack((A, B))
print(C)
C[1,:] = C[0,:] # Place row 1 into row 2
C[0,:] = np.array([1, 2, 3, 4, 5]) # Replace row 1 with the numbers 1 - 5
print(C.shape) # Outputs the dimensions of C
print("="*60)
print("Permute columns as [0, 4, 2, 1, 3]")
print("-"*34)
C = C[:, [0, 4, 2, 1, 3]] # Permute columns
print(C)
A = np.zeros((4,5)) # Create a matrix of zeros
A = np.ones((4,5)) # Create a matrix of ones
np.sum(A, axis=0) # Sum along columns
np.sum(A, axis=1) # Sum along rows
np.sum(A) # Sum all elements
# Advanced operations.
print("="*60)
print("\"Advanced operations\"")
print("-"*34)
A = np.random.randn(3,3)
print("A = ")
print(A)
print("Inverse", la.inv(A)) # Takes the inverse of a matrix
print("A^{-1} * A = ", np.dot(la.inv(A), A))
print("Matrix exponential; exp(A) = ", la.expm(A)) # Performs the matrix exponential, exp(A) would be elementwise
print("determinant of A is ", la.det(A)) # Determinant of A
print("eigenvalues are ", la.eig(A)[0]) # Outputs the eigenvalues of A
# W has eigenvalues, V has eigenvectors
(W,V) = la.eig(A)
print(np.dot(V, W)-np.dot(A, V))
# search
A = np.random.rand(10)
print(A)
print("Indices with value greater than .7 are...")
print(np.arange(len(A))[A >= .7])
## Inf, NaN, and rounding
a = np.Inf # Infinity
print(a)
print(1/a) # 1/infinity = 0
a = np.NaN # Not a number
print(a)
print(1/a)
print(np.isnan(a)) # returns 1 if a is not a number
f = 1.234;
# Round
print(np.round(f))
print(np.floor(f))
print(np.ceil(f))
## Plotting
# Vector plot
x = np.linspace(-1, 1, 80)
plt.plot(x, x**2, 'b.-') # Overwrites previous plot
# Labeling commands
plt.title('This is a title')
plt.xlabel('x')
plt.ylabel('y')
## First Order ODE
# Define the anonymous function y' = f(t,y)
# In this case y' = -y/2.0
# This is an example of a stable one-dimensional linear system.
yp = lambda y, t: -y/2.0
t = np.linspace(1, 10, 20)
y = odeint(yp, 2, t) # Solve the ODE yp from time 1 to 10, using y'(t0)=2.
plt.figure()
plt.plot(t, y) # Plot ODE solution.
plt.xlabel("time $t_{\mathrm{sampled}}$") # In-line LaTeX math support.
plt.ylabel("output", size=18) # Notice how font size is changed.
plt.show() # show() causes the script to pause and display figures.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment