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