Created for the CDS 110 course at Caltech but which others may find useful.

cds110-tutorial.py

#!/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.