Here you can find all of the code samples from A Student’s Guide to Python for Physical Modeling. You can copy and paste these fragments into a text editor or an IPython command line to run the code.
The entire code repository is also available on GitHub:
https://github.com/dr-kinder/code-samples
Some of the samples are drawn directly from the text. Others have been created to illustrate a point or technique. readme.txt contains a description of all files here.
These code samples are distributed under the BSD license. Be sure to read LICENSE.txt before using any of the code samples.
If you already know how to use the files, you can Jump to The Index.
Using the Code Samples
On this page, you will find the text of about 50 Python programs, which can be copied and pasted into the text editor of your choice or the IPython command prompt.
Just click on the link of the file name you are looking for. Then, use your mouse to select all of the text between lines of asterisks. (Do not copy the lines of *’s, though.) For example,
***************************************************************************
string_format.py
***************************************************************************
COPY ALL TEXT IN THIS REGION
***************************************************************************
string_percent.py
***************************************************************************
...After you copy and paste this text into your editor, you should be able to run the program in a Python interpreter. Your browser’s “Back” button will return you to the top of this page after you follow a link.
Long lines of code may be wrapped in your browser, but they should paste correctly into your text editor.
If you are ready to begin, you can Jump to The Index.
Code Archive
The entire code repository is available on GitHub:
https://github.com/dr-kinder/code-samples
You can download a zipped folder containing all of the code samples from GitHub. If you know how to use Git — or if you want to learn (see Appendix B of the Second Edition!) — you can fork or clone the repository to obtain your own copy of all of the code samples.
Troubleshooting
readme.txt contains a brief description of all the available files.
All of the .py files and the plain text code has been tested with Python 3.9. Most programs also run correctly with Python 2.7. However, if you are using Python 2, you should execute the following two commands before running the programs to ensure they perform as intended:
from __future__ import division, print_function
input = raw_inputTwo of the scripts, perrin.py and import_text.py, load data files. These data files are identified in the comments and included in the code repository and data_sets.zip. They should be copied into the current working directory before executing the scripts.
Text File
You can download a single text file that contains all of the code samples by following this link:
The link above will take you to a Google drive Web page where you can download a file called code_samples.txt. To save the file to your computer, click on the Download icon at the top of the page. This is a horizontal bar with a downward-pointing arrow above it.
code-samples.txt contains the text of about 50 Python programs, which can be copied and pasted into the text editor of your choice. (This is not a Python script or module. Trying to execute code_samples.txt with Python will result in an error.)
To use this file, first locate the file name you are looking for. Then, use your mouse to select all of the text between lines of asterisks. After you copy and paste this text into your editor, you should be able to run the program with a Python interpreter.
Program Index
LICENSE.txt
readme.txt
string_format.py
string_percent.py
for_loop.py
while_loop.py
vectorize.py
projectile.py
branching.py
nesting.py
import_text.py
save_load.py
print_write.py
simple_plot.py
graph_modifications.py
line3d.py
subplot.py
subplots.py
measurements.py
rotate.py
average.py
histogram.py
contour.py
matrix_inversion.py
quadrature.py
simple_oscillator.py
solve_ode.py
parametric_oscillator.py
ivp_comparison.py
vortex.py
gradient.py
streamlines.py
data_images.py
walker.py
html_movie.py
waves.py
convolution.py
perrin.py
sympy_examples.py
first_passage.py
data_dictionary.py
nd_random_walks.py
surprise.py
scope.py
name_collision.py
fancy_plot.py
legend.py
random_walk.py
surface.py
regression.py
bar3d.py
perrin.py
shading.py
Source Code
***************************************************************************
LICENSE.txt
***************************************************************************
Copyright (c) 2014, 2018, 2021, by Jesse M. Kinder and Philip C. Nelson
All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are
met:
* Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
* 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.
* Neither the name of the authors nor the names of other
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 JESSE M. KINDER OR
PHILIP C. NELSON 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.
================================================================================
Additional license information for html_movie.py:
html_movie.py was adapted from the scitools module created by Hans
Petter Langtangen. The scitools license is reproduced below.
================================================================================
Copyright (c) 2007-2009, Hans Petter Langtangen <hpl@simula.no> and
Simula Resarch Laboratory.
All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are
met:
* Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
* 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.
* Neither the name of Simula Research Laboratory 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 THE COPYRIGHT OWNER
OR 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.
***************************************************************************
readme.txt
***************************************************************************
CODE SAMPLES
============
-------------------------------------------------------------------------
These are code fragments for _A Student's Guide to Python for Physical
Modeling_ by Jesse M. Kinder and Philip Nelson.
www.amazon.com/Students-Guide-Python-Physical-Modeling/dp/0691223653
Use of these code fragments is subject to the terms of LICENSE.txt.
This file provides a brief description of each fragment. More details
can be found in the comments of the individual scripts and modules.
-------------------------------------------------------------------------
CONTENTS
--------
The order of appearance in the text (Second Edition) is as follows.
* Chapter 2
- string_format.py
- string_percent.py
- for_loop.py
- while_loop.py
* Chapter 3
- vectorize.py
- projectile.py
- branching.py
- nesting.py
* Chapter 4
- import_text.py
- save_load.py
- print_write.py
- simple_plot.py
- graph_modifications.py
- line3d.py
- subplot.py
- subplots.py
* Chapter 6
- measurements.py
- rotate.py
- average.py
- histogram.py
- contour.py
- matrix_inversion.py
- quadrature.py
- simple_oscillator.py
- solve_ode.py
- parametric_oscillator.py
- ivp_comparison.py
- vortex.py
- gradient.py
- streamlines.py
* Chapter 8
- data_images.py
- walker.py
- waves.py
* Chapter 9
- convolution.py
* Chapter 10
- first_passage.py
- data_dictionary.py
- nd_random_walks.py
* Epilogue
- surprise.py
* Appendix F
- scope.py
- name_collision.py
* Your Turn
- fancy_plot.py
- legend.py
- measurements.py
- random_walk.py
- surface.py
- regression.py
* Additional
- bar3d.py
- html_movie.py
- perrin.py
- shading.py
- sympy_examples.py
DESCRIPTIONS
------------
### string_format.py
Illustrate string formatting using the .format() method.
### string_percent.py
Illustrate string formatting using the % method.
### for_loop.py
Use a for loop to generate solutions to the quadratic equation.
This script illustrates the funamental form of a for loop. For
alternate solutions to the same problem, see while_loop.py and
vectorize.py.
### while_loop.py
Use a while loop to generate solutions to the quadratic equation.
For alternate solutions to the same problem, see for_loop.py and
vectorize.py.
### vectorize.py
Use vectorized code to generate solutions to the quadratic equation.
For alternate solutions to the same problem, see for_loop.py and
while_loop.py.
### projectile.py
Calculate how long an object is in the air when thrown from a spcified
height with a range of initial speeds assuming constant acceleration due
to gravity:
0.5 * g * t**2 - v0 * t - y0 = 0
This script illustrates good coding practice in the solution of a simple
problem: parameters with descriptive names, comments, whitespace, and
blocking with '#%%' for debugging in Spyder.
### branching.py
This script illustrates branching with the use of multiple conditional
statements:
if <condition1>:
...
elif <condition2>:
...
else:
...
### nesting.py
Use nested for loops to fill a two-dimensional array of values.
This script illustrates "nesting" --- one for loop inside of another.
### import_text.py
Load data from a text file by reading the file line by line.
This script reads in data from a text file and stores it in a NumPy
array. It can be adapted to load data from files that are difficult or
impossible to load with NumPy's np.loadtxt function.
### save_load.py
Save array data using NumPy's available methods, then load saved data.
This script demonstrates the various methods of saving and loading data
using NumPy arrays.
### print_write.py
Write same data to a file and print to display.
This script illustrates the similarites between writing to a text file
and printing to the screen.
### simple_plot.py
Create and display a basic plot.
### graph_modifications.py
This script creates a simple plot with two lines, then modifies several
features of the plot, including axis labels, data labels, legend, line
style, tick labels, and title.
### line3d.py
Create a three-dimensional parametric plot.
This script demonstrates how to create three-dimensional plots using the
Axes3D method from the mpl_toolkits.mplot3d module. See also surface.py.
### subplot.py
Create four plots in the same figure using plt.subplot().
This script demonstrates PyPlot's subplot method, which can be used to
display several plots side-by-side in the same figure. For another
approach, see subplots.py. The names are similar, but the behavior of
the plt.subplot and plt.subplots is different.
### subplots.py
Create four plots in the same figure using plt.subplots().
This script demonstrates PyPlot's subplots method, which can be used to
display several plots side-by-side in the same figure. For another
approach, see subplots.py. The names are similar, but the behavior of
the plt.subplot and plt.subplots is different.
### measurements.py
Functions to calculate distance between points using different metrics.
This script illustrates the fundamental form of user-defined functions
as well as keyword arguments and default values.
### rotate.py
Define function to rotate a vector in two dimensions.
### average.py
Compute and return the cummulative average of an array.
This script illustrates some principles of functional programming.
### histogram.py
Create histograms of random numbers.
This script illustrates how to use NumPy and PyPlot to create histograms
and bar plots.
### contour.py
Create a labeled contour plot.
This script illustrates how to generate a grid of coordinates for
contour and surface plots. It also demonstrates some options of
plt.contour and shows how to label contour lines.
### matrix_inversion.py
Invert a simple matrix to solve a system of linear equations.
This script illustrates the use of a special method from the SciPy linear
algebra library, scipy.linalg.
### quadrature.py
Integrate two functions using quad.
This script demonstrates numerical integration using the quad method of
scipy.integrate. The first function is a built-in NumPy funciton whose
integral can be computed with pencil and paper for comparison. The
second is a user-defined function.
### simple_oscillator.py
Define function to use in solution of differential equation for a simple
harmonic oscillator.
This script illustrates how to write a function that generates the array
required to integrate a second-order ordinary differential equation. It
is imported and used in solve_ode.py.
### solve_ode.py
Solution of ODE for harmonic oscillator.
This script imports the function F(y,t) in simple_oscillator.py then
uses the odeint method of scipy.integrate to solve the ordinary
differential equation defined by F(y,t).
See ivp_comparison.py for an example of the solve_ivp method, an
alternative to odeint.
### parametric_oscillator.py
Define a parametric function that accepts 4 parameters then integrate it
using odeint.
This script illustrates two methods for using scipy.integrate's odeint
methods to integrate a function that accepts more than two parameters.
### ivp_comparison.py
Compare different ODE solvers using solve_ivp.
solve_ivp offers an alternative to odeint for solving ordinary
differential equations. See solve_ode.py for an example of odeint.
### vortex.py
Create a quiver plot.
This script illustrates the use of PyPlot's quiver method.
### gradient.py
Calculate and display the gradient of a two-dimensional Gaussian.
This script illustrates the use of NumPy's gradient function and
demonstrates how to display a vector field. It displays the gradient as
a quiver plot superimposed on a filled contour plot of the Gaussian.
### streamlines.py
Create streamlines from a vector field.
This script demonstrates the use of PyPlot's streamplot method for
visualizing solutions to a differential equation defined by a vector
field.
### data_images.py
Illustrate differences between image and Cartesian coordinates.
The coordinates used in plotting functions and displaying images follow
different conventions. This script creates a figure that illustrates
the two conventions.
### walker.py
Make a movie out of the steps of a two-dimensional random walk.
This script demonstrates the use of the FuncAnimation method of
Matplotlib's animation module to create a movie. If ffmpeg or mencoder
is installed on this computer, the script will save the movie to an mp4
file.
### html_movie.py
Module to generate an HTML document from a collection of images. When
viewed in a Web browser, the document will display a movie whose frames
are the individual images.
This module is adapted from the scitools library developed by Hans
Petter Langtangen.
### waves.py
Create an HTML animation of a moving Gaussian waves.
This script illustrates a method for combining a series of plots into an
animation using HTML and Javascript. It uses the html_movie.py module,
which is adapted from the scitools library developed by Hans Petter
Langtangen.
### convolution.py
This script creates an eLoG (elongated Laplacian of Gaussian) filter
that emphasizes long, vertical lines in a figure. The effect of the
filter is demonstrated on a plus sign.
### sympy_examples.py
Demonstrate some useful methods available in the SymPy module.
Python may complain about undefined variables if you attempt to run this
script. init_session() defines several variables, but Python may not be
aware of this. It is better to run the commands one at a time from the
command line.
### first_passage.py
Define a function to simulate first passage of a random walker.
Simulate a random walker that starts at the origin and takes steps to
the right with probability p and to the left with probability 1-p.
Return the number of steps for the first passage of location x==L, or
give up after N steps.
### data_dictionary.py
Store input and data two different simulations in a dictionary.
This script requires first_passage.py to be in the same directory.
### nd_random_walks.py
Python class to simulate various random walks in N dimensions.
RandomWalk -- general random walk class
LatticeWalk -- random walk on a D-dimensional lattice
default is a cubic lattice in D dimensions
TriangularWalk -- 2D triangular lattice
HoneycombWalk -- 2D honeycomb lattice
DirectionalWalk -- random walks of variable step size in D dimensions
defaults to constant step length in random directions
UniformWalk -- step size drawn from uniform distribution
GaussianWalk -- step size drawn from normal distribution
ExponentialWalk -- step size drawn from exponential distribution
ParetoWalk -- step size drawn from power law distribution
### surprise.py
This script will create a familar but interesting image.
It may take about a minute to run.
### scope.py
Demonstrate Python's rules of scope --- i.e., how Python looks up
variable names.
### name_collision.py
Illustrate how Python's rules of scope prevent name collisions.
### fancy_plot.py
Add a title and axis labels to a simple plot.
### legend.py
Create a plot with a legend to distinguish multiple curves.
### random_walk.py
Monte Carlo simulation of a two-dimensional random walk.
This script illustrates the use of a random number generator to create a
time series for a random walk.
### surface.py
Create a three-dimensional surface plot.
This script demonstrates how to create three-dimensional plots using the
Axes3D method from the mpl_toolkits.mplot3d module. See also line3d.py.
### regression.py
Example of linear regression on data from the first passage problem.
first_passage.py must be in the working directory. The script could
take a while to finish if (samples x nmax) ~ 10**8 or more.
### bar3d.py
Create a 3D histogram ("Lego plot").
This script provides a visual example of the Central Limit Theorem.
### perrin.py
Generate figure displaying Perrin's experimental data on Brownian
motion. This script requires the data set 04brownian/g26perrindata.npy.
The script illustrates loading and plotting a data set and includes
LaTeX formatting of axis labels and grid lines.
### shading.py
Demonstrate shading of surface plots using Matplotlib's LightSource.
-------------------------------------------------------------------------
***************************************************************************
string_format.py
***************************************************************************
# string_format.py
# -------------------------------------------------------------------------
# Illustrate formatting strings using the .format() method.
# -------------------------------------------------------------------------
import numpy as np
print("The value of pi is approximately " + str(np.pi))
print("The value of {} is approximately {:.5f}".format('pi', np.pi))
s = "{1:d} plus {0:d} is {2:d}"
print(s.format(2,4,2+4))
print("Every {2} has its {3}.".format('dog','day','rose','thorn'))
print("The third element of the list is {0[2]:g}.".format(np.arange(10)))
***************************************************************************
string_percent.py
***************************************************************************
# string_percent.py
# -------------------------------------------------------------------------
# Illustrate string formatting using the % method.
# -------------------------------------------------------------------------
import numpy as np
print("The value of pi is approximately " + str(np.pi))
print("The value of %s is approximately %.5f" % ('pi', np.pi))
s = "%d plus %d is %d"
print(s % (2, 4, 2 + 4))
***************************************************************************
for_loop.py
***************************************************************************
# for_loop.py
# -------------------------------------------------------------------------
# Use a for loop to generate multiple solutions to the quadratic equation.
# -------------------------------------------------------------------------
import numpy as np
b, c = 2, -1
for a in np.arange(-1, 2, 0.3):
x = (-b + np.sqrt(b**2 - 4*a*c)) / (2*a)
print("a= {:.4f}, x= {:.4f}".format(a,x))
***************************************************************************
while_loop.py
***************************************************************************
# while_loop.py
# -------------------------------------------------------------------------
# Use a while loop to generate multiple solutions to the quadratic equation.
# -------------------------------------------------------------------------
import numpy as np
a, b, c = 2, 2, -1
while (b**2 - 4*a*c >= 0):
x = (-b + np.sqrt(b**2 - 4*a*c)) / (2*a)
print("a = {:.4f}, x = {:.4f}".format(a,x))
a = a - 0.3
print("done!")
***************************************************************************
vectorize.py
***************************************************************************
# vectorize.py
# -------------------------------------------------------------------------
# Use vectorized operations to generate multiple solutions to the
# quadratic equation.
# -------------------------------------------------------------------------
import numpy as np
b, c = 2, -1
a = np.arange(-1, 2, 0.3)
x = (-b + np.sqrt(b**2 - 4*a*c)) / (2*a)
***************************************************************************
projectile.py
***************************************************************************
# projectile.py
# -----------------------------------------------------------------------------
# Calculate how long an object is in the air when thrown from a specified height
# with a range of initial speeds assuming constant acceleration due to gravity:
# 0.5 * g * t**2 - v0 * t - y0 = 0
# -----------------------------------------------------------------------------
import numpy as np
#%% Initialization of variables.
initial_speed = 0.0 # v0 = initial vertical speed of ball in [m/s]
impact_time = 0.0 # t = time of impact in [s] (computed in loop)
#%% Initialization of parameters.
g = 9.8066 # gravitational acceleration in [m/s^2]
initial_height = 2.0 # y0 = height ball is thrown from in [m]
speed_increment = 5.0 # how much to increase speed in [m/s] for each iteration
cutoff_time = 10.0 # stop computing after impact time exceeds cutoff
#%% Calculate and display impact time. Increment initial speed each step.
# Repeat until impact time exceeds cutoff.
while impact_time < cutoff_time:
# Use quadratic equation to solve kinematic equation for impact time:
impact_time = (np.sqrt(initial_speed**2 + 2 * g * initial_height) + initial_speed) / g
print("speed= {} m/s; time= {:.1f} s".format(initial_speed, impact_time))
initial_speed += speed_increment
print("Calculation complete.")
***************************************************************************
branching.py
***************************************************************************
# branching.py
# -----------------------------------------------------------------------------
# This script illustrates branching.
# -----------------------------------------------------------------------------
import numpy as np
for trial in range(5):
userInput = input('Pick a number: ')
userNumber = float(userInput)
if userNumber < 0:
print('Square root is not real.')
else:
print('Square root of {} is {:.4f}.'.format(userNumber, np.sqrt(userNumber)))
userAgain = input('Another [y/n]? ')
if userAgain != 'y':
break
if trial == 4:
print('Sorry, only 5 per customer.')
elif userAgain == 'n':
print('Bye!')
else:
print('Sorry, I did not understand that.')
***************************************************************************
nesting.py
***************************************************************************
# nesting.py
# -------------------------------------------------------------------------
# Use nested for loops to fill a two-dimensional array of values.
# -------------------------------------------------------------------------
import numpy as np
# Set dimensions of array.
rows = 3
columns = 4
# Create empty array then fill with values.
A = np.zeros((rows, columns))
for i in range(rows):
for j in range(columns):
A[i, j] = i**2 + j**3
***************************************************************************
import_text.py
***************************************************************************
# import_text.py
# -------------------------------------------------------------------------
# Load data from a text file by reading the file line by line.
# This script requires the data set 01HIVseries/HIVseries.csv.
# -------------------------------------------------------------------------
import numpy as np
f = open("HIVseries.csv")
temp_data = []
for line in f:
print(line)
x,y = line.split(',')
temp_data += [ (float(x), float(y)) ]
f.close()
data_set = np.array(temp_data)
***************************************************************************
save_load.py
***************************************************************************
# save_load.py
# -------------------------------------------------------------------------
# Save array data using NumPy's available methods, then load saved data.
# -------------------------------------------------------------------------
import numpy as np
# Generate array data.
x = np.linspace(0, 1, 1001)
y = 3*np.sin(x)**3 - np.sin(x)
# Save data to text, .npy, and .npz files.
np.save('x_values', x)
np.save('y_values', y)
np.savetxt('x_values.dat', x)
np.savetxt('y_values.dat', y)
np.savez('xy_values', x_vals=x, y_vals=y)
# Load saved data.
x2 = np.load('x_values.npy')
y2 = np.loadtxt('y_values.dat')
w = np.load('xy_values.npz')
# Check equality of original and saved data.
print((x2 == x).all())
print((y2 == y).all())
print((w['x_vals'] == x2).all())
print((w['y_vals'] == y2).all())
***************************************************************************
print_write.py
***************************************************************************
# print_write.py
# -------------------------------------------------------------------------
# Write same data to a file and print to display.
# -------------------------------------------------------------------------
f = open('power.txt','w')
print(" N \t\t2**N\t\t3**N") # print labels for columns
f.write(" N \t\t2**N\t\t3**N\n") # write labels to file
print("---\t\t----\t\t----") # print separator
f.write("---\t\t----\t\t----\n") # write separator to file
#%% loop over integers from 0 to 10 and print/write results
for N in range(11):
print("{:d}\t\t{:d}\t\t{:d}".format(N, pow(2,N), pow(3,N)))
f.write("{:d}\t\t{:d}\t\t{:d}\n".format(N, pow(2,N), pow(3,N)))
f.close()
***************************************************************************
simple_plot.py
***************************************************************************
# simple_plot.py
# -------------------------------------------------------------------------
# Create and display a basic plot.
# -------------------------------------------------------------------------
import numpy as np
import matplotlib.pyplot as plt
num_points = 26
x_min, x_max = 0, 4
x_values = np.linspace(x_min, x_max, num_points)
y_values = x_values**2
plt.figure()
plt.plot(x_values, y_values)
plt.show()
***************************************************************************
graph_modifications.py
***************************************************************************
# graph_modifications.py
# -------------------------------------------------------------------------
# This script creates a simple plot with two lines, then modifies many
# features of the plot, including axis labels, labels and legend, line
# style, tick labels, and title.
# -------------------------------------------------------------------------
import numpy as np
import matplotlib.pyplot as plt
# Generate data for plot.
num_points = 26
x_min, x_max = 0, 4
x_values = np.linspace(x_min, x_max, num_points)
y_values = x_values**2
# Create empty figure.
plt.figure()
# Plot data.
plt.plot(x_values, y_values, label="Population 1") # label when plot is created
plt.plot(x_values, x_values**3, label="Population 2") # label when plot is created
plt.legend()
# Gain control of current Axes object.
ax = plt.gca()
# Give plot a title.
ax.set_title("My First Plot", family='monospace', size=24, weight='bold')
# Label the axes.
ax.set_xlabel("Time [days]")
ax.set_ylabel("Population")
# Change tick labels and font
ax.set_xticklabels(ax.get_xticks(), family='monospace', fontsize=10)
ax.set_yticklabels(ax.get_yticks(), family='monospace', fontsize=10)
# Change the legend
lines = ax.get_lines() # returns a list of line objects
lines[0].set_label("Infected Population") # change labels using line objects
lines[1].set_label("Cured Population") # change labels using line objects
ax.legend() # display legend in plot
# Make the first line a thick, red, dashed line.
plt.setp(lines[0], linestyle='--', linewidth=3, color='r')
# Change the legend again.
ax.legend(("Healthy", "Recovered")) # change labels using Axes object
plt.show()
***************************************************************************
line3d.py
***************************************************************************
# line3d.py
# -------------------------------------------------------------------------
# Create a three-dimensional parametric plot.
# -------------------------------------------------------------------------
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D # import 3D plotting tool
fig = plt.figure() # create a new figure
ax = Axes3D(fig) # create 3D plotting object attached to figure
t = np.linspace(0, 5*np.pi, 501) # define parameter for parametric plot
ax.plot(t * np.cos(t), t * np.sin(t), t) # generate 3D plot
plt.show()
***************************************************************************
subplot.py
***************************************************************************
# subplot.py
# -------------------------------------------------------------------------
# Create four plots in the same figure using plt.subplot().
# -------------------------------------------------------------------------
import numpy as np
import matplotlib.pyplot as plt
from numpy.random import random
t = np.linspace(0, 1, 101)
# Create figure.
plt.figure()
# Use plt.subplot to add subplots.
plt.subplot(2, 2, 1); plt.hist(random(20)) # Upper left
plt.subplot(2, 2, 2); plt.plot(t, t**2, t, t**3 - t) # Upper right
plt.subplot(2, 2, 3); plt.plot(random(20), random(20), 'r*') # Lower left
plt.subplot(2, 2, 4); plt.plot(t*np.cos(10*t), t*np.sin(10*t)) # Lower right
plt.suptitle("Data and Functions") # Entire plot
plt.show()
***************************************************************************
subplots.py
***************************************************************************
# subplot.py
# -------------------------------------------------------------------------
# Create four plots in the same figure using plt.subplots().
# -------------------------------------------------------------------------
import numpy as np
import matplotlib.pyplot as plt
from numpy.random import random
t = np.linspace(0, 1, 101)
# Create Figure and Axes objects.
fig, ax = plt.subplots(2,2)
# Use methods of Axes objects to add subplots.
ax[0,0].hist(random(20)) # Upper left
ax[0,1].plot(t, t**2, t, t**3 - t) # Upper right
ax[1,0].plot(random(20), random(20), 'r*') # Lower left
ax[1,1].plot(t*np.cos(10*t), t*np.sin(10*t)) # Lower right
fig.suptitle("Data and Functions") # Entire plot
plt.show()
***************************************************************************
measurements.py
***************************************************************************
# measurements.py
"""
Functions to calculate distance between points using different metrics.
"""
import numpy as np
def crow(pointA, pointB):
"""
Distance between points A and B "as the crow flies."
pointA = (x1, y1)
pointB = (x2, y2)
returns sqrt( (x2-x1)**2 + (y2-y1)**2 )
"""
interval = np.sqrt( (pointA[0] - pointB[0])**2 + \
(pointA[1] - pointB[1])**2 )
return interval
def taxicab(pointA, pointB):
"""
Distance between points A and B "as the cab drives."
pointA = (x1, y1)
pointB = (x2, y2)
returns |x2-x1| + |y2-y1|
"""
interval = abs(pointB[0] - pointA[0]) + \
abs(pointB[1] - pointA[1])
return interval
def distance(pointA, pointB=(0,0), metric='taxi'):
"""
Return distance between points A and B. If metric is 'taxi', use taxicab
metric. Otherwise, use Euclidean distance.
pointA = (x1, y1)
pointB = (x2, y2)
"""
if metric == 'taxi':
return taxicab(pointA, pointB)
else:
return crow(pointA, pointB)
***************************************************************************
rotate.py
***************************************************************************
# rotate.py
# -------------------------------------------------------------------------
# Define function to rotate a vector in two dimensions.
# -------------------------------------------------------------------------
import numpy as np
def rotate_vector(vector, angle):
"""
Rotate two-dimensional vector through given angle.
vector = (x,y)
angle = rotation angle in radians (counterclockwise)
Returns the image of vector under rotation as a NumPy array.
"""
rotation_matrix = np.array([[ np.cos(angle), -np.sin(angle) ],
[ np.sin(angle), np.cos(angle) ]])
return np.dot(rotation_matrix, vector)
***************************************************************************
average.py
***************************************************************************
# average.py
# -------------------------------------------------------------------------
# Compute and return the cummulative average of an array.
# -------------------------------------------------------------------------
import numpy as np
def running_average(x):
"""
Return cummulative average of an array.
"""
y = np.zeros(len(x)) # new array to store result
current_sum = 0.0 # running sum of elements of x
for i in range(len(x)):
current_sum += x[i] # increment sum
y[i] = current_sum / (i + 1.0) # update running average
return y
***************************************************************************
histogram.py
***************************************************************************
# histogram.py
# -------------------------------------------------------------------------
# Create histograms of random numbers.
# -------------------------------------------------------------------------
import numpy as np
import matplotlib.pyplot as plt
from numpy.random import random as rng
#%% Generate data and view PyPlot's default histogram.
data = rng(1000)
plt.figure()
plt.hist(data)
#%% Get binned data from NumPy and make a colorful histogram where the
# width of each bin is proportional to the number of elements in it.
counts, bin_edges = np.histogram(data)
bin_size = bin_edges[1] - bin_edges[0]
new_widths = bin_size * counts / counts.max()
plt.figure()
plt.bar(bin_edges[:-1], counts, width=new_widths, color=['r','g','b'])
#%% Provide logarithmically spaced bin edges rather than using defaults.
log2bins = np.logspace(-8, 0, num=9, base=2)
log2bins[0] = 0.0 # set first bin edge to zero instead of 1/256
plt.figure()
plt.hist(data, bins=log2bins)
plt.show()
***************************************************************************
contour.py
***************************************************************************
# contour.py
# -------------------------------------------------------------------------
# Create a labeled contour plot.
# -------------------------------------------------------------------------
import numpy as np
import matplotlib.pyplot as plt
#%% Create a grid of x and y coordinates
x_vals = np.linspace(-3, 3, 201)
y_vals = np.linspace(0, 10, 101)
X, Y = np.meshgrid(x_vals, y_vals)
#%% Generate function values.
Z = np.cos(X) * np.sin(Y)
#%% Plot and label contours.
plt.figure()
cs = plt.contour(X, Y, Z, 10, linewidths=3, colors='k')
plt.clabel(cs,fontsize=10)
plt.show()
***************************************************************************
matrix_inversion.py
***************************************************************************
# matrix_inversion.py
# -------------------------------------------------------------------------
# Invert a simple matrix to solve a system of linear equations.
# -------------------------------------------------------------------------
import numpy as np
from scipy.linalg import inv
#%% Set up and solve C.x = a
a = np.array([-1, 5])
C = np.array([[1, 3], [3, 4]])
x = np.dot(inv(C), a)
#%% Check solution.
error = np.dot(C,x) - a
***************************************************************************
quadrature.py
***************************************************************************
# quadrature.py
# -------------------------------------------------------------------------
# Integrate two functions using quad.
# -------------------------------------------------------------------------
import numpy as np
from scipy.integrate import quad
#%% Integrate a built-in function.
upper_limit = np.linspace(0,3*np.pi,16)
cos_integral = np.zeros(upper_limit.size)
for i in range(upper_limit.size):
cos_integral[i], error = quad(np.cos, 0, upper_limit[i])
#%% Now integrate a user-defined function.
def integrand(x): return np.exp(-x**2/2)
upper_limit = np.linspace(0, 5, 51)
gauss_integral = np.zeros(upper_limit.size)
for i in range(upper_limit.size):
gauss_integral[i], error = quad(integrand, 0, upper_limit[i])
***************************************************************************
simple_oscillator.py
***************************************************************************
# simple_oscillator.py
# -------------------------------------------------------------------------
# Define function to use in solution of differential equation for a simple
# harmonic oscillator.
# -------------------------------------------------------------------------
def F(y, t):
"""
Return derivatives for 2nd order ODE y'' = -y.
"""
dy = [0, 0] # preallocate list to store derivatives
dy[0] = y[1] # first derivative of y(t)
dy[1] = -y[0] # second derivative of y(t)
return dy
***************************************************************************
solve_ode.py
***************************************************************************
# solveODE.py
# -----------------------------------------------------------------------------
# Solution of ODE for harmonic oscillator.
# -----------------------------------------------------------------------------
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
# Import function to integrate:
from simple_oscillator import F
# array of time values to study
t_min = 0; t_max = 10; dt = 0.1
t = np.arange(t_min, t_max+dt, dt)
# two sets of initial conditions
initial_conditions = [ (1.0, 0.0), (0.0, 1.0) ]
plt.figure() # Create figure; then, add plots.
for y0 in initial_conditions:
y = odeint(F, y0, t)
plt.plot(t, y[:, 0], linewidth=2)
skip = 5
t_test = t[::skip] # compare at a subset of points
plt.plot(t_test, np.cos(t_test), 'bo') # exact solution for y0 = (1,0)
plt.plot(t_test, np.sin(t_test), 'go') # exact solution for y0 = (0,1)
plt.show()
***************************************************************************
parametric_oscillator.py
***************************************************************************
# parametric_oscillator.py
# -------------------------------------------------------------------------
# Define a parametric function that accepts 4 parameters then integrate it
# using odeint.
# -------------------------------------------------------------------------
import numpy as np
from scipy.integrate import odeint
def F(y, t, spring_constant=1.0, mass=1.0):
"""
Return derivatives for harmonic oscillator:
y'' = -(k/m) * y
y = displacement in [m]
k = spring_constant in [N/m]
m = mass in [kg]
"""
dy = [0, 0] # array to store derivatives
dy[0] = y[1]
dy[1] = -(spring_constant/mass) * y[0]
return dy
# -------------------------------------------------------------------------
# Integrate parametric function using two different methods.
# -------------------------------------------------------------------------
y0 = (1.0, 0.0) # initial conditions
t = np.linspace(0, 10, 101) # times at which y(t) will be evaluated
# Method 1 -- dummy function
def G(y, t): return F(y, t, 2.0, 0.5)
yA = odeint(G, y0, t)
# Method 2 -- keywords
yB = odeint(F, y0, t, args=(2.0, 0.5))
***************************************************************************
ivp_comparison.py
***************************************************************************
# ivp_comparison.py
# -------------------------------------------------------------------------
# Compare different ODE solvers using solve_ivp.
# -------------------------------------------------------------------------
import numpy as np, matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
# Define ODE to integrate: simple harmonic oscillator.
def f(t,y): return [ y[1], -y[0] ]
# Define time interval.
t_min = 0
t_max = 10
# Define initial conditions.
y0 = [1.0, 0.0]
# Integrate the ODE using defaults.
result = solve_ivp(f, (t_min, t_max), y0)
plt.plot(result.t, result.y[0], '^k', label='RK45')
# Specify time series and use different solver.
dt = 0.1
t_vals = np.arange(t_min, t_max + dt, dt)
result = solve_ivp(f, (t_min, t_max), y0, t_eval=t_vals, method='BDF')
plt.plot(result.t, result.y[0], '.r', label='BDF')
plt.legend()
plt.show()
***************************************************************************
vortex.py
***************************************************************************
# vortex.py
# -------------------------------------------------------------------------
# Create a quiver plot.
# -------------------------------------------------------------------------
import numpy as np
import matplotlib.pyplot as plt
coords = np.linspace(-1, 1, 11)
X, Y = np.meshgrid(coords, coords)
Vx, Vy = Y, -X
plt.figure()
# pivot='mid' places centers of arrows on grid points
plt.quiver(X, Y, Vx, Vy, pivot='mid')
plt.axis('square')
plt.show()
***************************************************************************
gradient.py
***************************************************************************
# gradient.py
# -------------------------------------------------------------------------
# Calculate and display the gradient of a two-dimensional Gaussian.
# -------------------------------------------------------------------------
import numpy as np
import matplotlib.pyplot as plt
coords = np.linspace(-2, 2, 101)
X, Y = np.meshgrid(coords[::5], coords[::5]) # coarse grid for vector field
R = np.sqrt(X**2 + Y**2)
Z = np.exp(-R**2)
x, y = np.meshgrid(coords, coords) # fine grid for contour plot
r = np.sqrt(x**2 + y**2)
z = np.exp(-r**2)
ds = coords[5] - coords[0] # coarse grid spacing
dX, dY = np.gradient(Z, ds) # calculate gradient
plt.figure()
plt.contourf(x, y, z, 25)
plt.set_cmap('coolwarm')
plt.quiver(X, Y, dX.transpose(), dY.transpose(), scale=25, color='k')
plt.axis('tight')
plt.axis('square')
plt.show()
***************************************************************************
streamlines.py
***************************************************************************
# streamlines.py
# -------------------------------------------------------------------------
# Create streamlines from a vector field.
# -------------------------------------------------------------------------
import numpy as np
import matplotlib.pyplot as plt
# generate grid of points
lower, upper, step = -2, 2, 0.1
coords = np.arange(lower, upper+step, step)
X, Y = np.meshgrid(coords, coords)
# define vector field
Vx, Vy = Y, -X
# display streamlines defined by vector field
plt.figure()
plt.streamplot(coords, coords, Vx, Vy)
plt.axis('square')
plt.show()
***************************************************************************
data_images.py
***************************************************************************
# data_images.py
# -------------------------------------------------------------------------
# Illustrate differences between image and Cartesian coordinates.
# -------------------------------------------------------------------------
#%% Setup
import numpy as np
import matplotlib.pyplot as plt
# Define a coordinate grid.
# Same spacing, but twice as many points along x-axis
x_max, y_max = 2, 1
x_num, y_num = 200, 100
# Create coordinate arrays.
x = np.linspace(0,x_max,x_num)
y = np.linspace(0,y_max,y_num)
#%% Generate data to plot.
# Assign function values to placeholder array.
z = np.zeros((x_num,y_num))
for i in range(x_num):
for j in range(y_num):
z[i][j] = (x[i] - 2*y[j])**2
# Use meshgrid to generate the same function values.
X,Y = np.meshgrid(x,y)
Z = (X-2*Y)**2
#%% Visualize results.
fig, ax = plt.subplots(2,3, figsize=(12,6))
fig.suptitle(r"Plots of $f(x,y) = (x-2y)^2$")
ax[0,0].imshow(z)
ax[0,0].set_title("Loop: Image Coordinates")
ax[0,1].imshow(z, origin="lower")
ax[0,1].set_title("Loop: Spatial Coordinates")
ax[0,2].imshow(Z)
ax[0,2].set_title("meshgrid: Image Coordinates")
ax[1,0].imshow(z.transpose(), origin="lower")
ax[1,0].set_title("Loop: Transpose + Spatial Coordinates")
ax[1,1].imshow(Z, origin="lower")
ax[1,1].set_title("meshgrid: Spatial Coordinates")
ax[1,2].pcolormesh(X, Y, Z, shading='auto')
ax[1,2].axis('image')
ax[1,2].set_title("pcolormesh")
plt.show()
***************************************************************************
walker.py
***************************************************************************
# walker.py
# -----------------------------------------------------------------------------
# Make a movie out of the steps of a two-dimensional random walk.
# -----------------------------------------------------------------------------
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation
from numpy.random import random as rand
# Set number of steps for each random walk.
num_steps = 100
# Create an empty figure of the desired size.
bound = 20
fig = plt.figure() # must have figure object for movie
ax = plt.axes(xlim=(-bound, bound), ylim=(-bound, bound))
ax.set_aspect('equal') # equal scaling for x and y axes
# Create a line and a point with no data. They will be updated during each
# frame of the animation.
(my_line,) = ax.plot([], [], lw=2) # line to show path
(my_point,) = ax.plot([], [], 'ro', ms=9) # dot to show current position
# Generate the random walk data.
x_steps = 2*(rand(num_steps) < 0.5) - 1 # generate random steps +/- 1
y_steps = 2*(rand(num_steps) < 0.5) - 1
x_coordinate = x_steps.cumsum() # sum steps to get position
y_coordinate = y_steps.cumsum()
# This function will generate each frame of the animation.
# It adds all of the data through frame n to a line
# and moves a point to the nth position of the walk.
def get_step(n, x, y, this_line, this_point):
this_line.set_data(x[:n+1], y[:n+1])
this_point.set_data(x[n], y[n])
# Call the animator and create the movie.
my_movie = animation.FuncAnimation(fig, get_step, frames=num_steps, \
fargs=(x_coordinate,y_coordinate,my_line,my_point) )
# Save the movie in the current directory.
# *** THIS WILL CAUSE AN ERROR UNLESS FFMPEG OR MENCODER IS INSTALLED. ***
my_movie.save('random_walk.mp4', fps=30)
***************************************************************************
html_movie.py
***************************************************************************
# html_movie.py
# -----------------------------------------------------------------------------
# This module provides a single function called movie that will create a Web
# page to display a sequence of images as a movie.
# -----------------------------------------------------------------------------
# This module is adapted from the scitools library developed by Hans Petter
# Langtangen. The source from which this module was derived may be found in the
# file <scitools/lib/scitools/easyviz/movie.py> of the current distribution of
# the scitools module (as of May 30, 2015), available at
#
# https://code.google.com/p/scitools/
#
# The MovieEncoder class and the movie function of <movie.py> from the scitools
# library have been merged into a single function called movie. The html_movie
# function has been renamed get_html.
#
# The movie function in the current module takes a collection of input files and
# generates an HTML file that displays the images as a movie. The output is
# nearly identical to that of scitools.std.movie called with the default 'html'
# encoder. This module is self-contained and does not require any of the other
# functions or libraries of the scitools module. The print statements have also
# been modified to work with both Python 2 and Python 3.
#
# scitools:
# ---------
# Copyright (c) 2007-2009, Hans Petter Langtangen <hpl@simula.no> and
# Simula Resarch Laboratory.
#
# All rights reserved.
#
# Both the current module and scitools are distributed under the BSD license.
# See LICENSE.txt for details.
# -----------------------------------------------------------------------------
"""
This module provides a single function called movie that will create a Web page
to display a sequence of images as a movie.
The input files can be given in a list or tuple of file names or specified by a
regular expression that will generate the filenames. I.e.,
input_files="frame_*.jpg" will use all .jpg files in the current directory with
a name of the form "frame_001.jpg", "frame_002.jpg", etc.
"""
# Modules that provide tools to search for file names using regular expressions.
import os
import glob
import re
# Specify types of image files that can be used to create a movie.
_legal_file_types = 'png gif jpg jpeg'.split()
# -------------------------------------------------------------------------
def movie(input_files, output_file="movie.html", fps=25):
"""
Take a list or tuple of image file names or a regular expression that will
generate those file names and then create an HTML file that uses JavaScript
to display those images as a movie.
Images must be PNG, JPEG, or GIF files, with extension .png, .jpg, .jpeg, or
.gif, respectively.
input_files: a list or tuple of file names, or regular expression
output_file: name of the output HTML movie file
fps: frame rate, in frames per second
The sequence of images will be determined by their order in input_files if
input_files is a list or tuple. If input_files is a regular expression like
"frame_*.jpg", the sequence of images will be the alphanumeric order of all
filenames that match the search pattern. Because "frame_10.jpg" comes
before "frame_2.jpg" in alphanumeric order, it is important to include
leading zeros in file names: "frame_002.jpg" comes before "frame_010.jpg".
"""
# Print blank lines before function output statements.
print('\n\n')
# Determine the file type of the input files.
if isinstance(input_files, (tuple,list)):
file_ = input_files[0]
elif isinstance(input_files, str):
file_ = input_files
else:
raise ValueError("The input files must be given as either a "\
"list or tuple of strings or a string, not '%s'" % \
type(input_files))
# Check that the input files do exist.
if isinstance(input_files, str):
# Are input_files on ffmpeg/mpeg2enc format or Unix wildcard format?
ffmpeg_format = r'(.+)%\d+d(\..+)'
m = re.search(ffmpeg_format, input_files, re.DOTALL)
if m:
wildcard_format = m.group(1) + '*' + m.group(2)
else:
wildcard_format = input_files
all_input_files = glob.glob(wildcard_format)
if not all_input_files:
print('No files of the form %s exist.' % input_files)
else:
print('Found %d files of the format %s.' % \
(len(all_input_files), input_files))
else:
# User provided a list or tuple of specific filenames.
all_input_files = input_files
error_encountered = False
for f in all_input_files:
if not os.path.isfile(f):
print('Input file %s does not exist.' % f)
error_encountered = True
if error_encountered:
raise IOError('Some input files were not found.')
fname, ext = os.path.splitext(file_)
if not ext:
raise ValueError("Unable to determine file type from file name.")
file_type = ext[1:] # remove the . (dot)
if not file_type in _legal_file_types:
raise TypeError("File type must be %s, not '%s'" % \
(_legal_file_types, file_type))
# Create an html file that can play image files
files = input_files
if isinstance(files, str):
files = glob.glob(files)
files.sort()
print('\nMaking HTML code for displaying ' + ', '.join(files))
# Turn frame rate into a pause, in milliseconds, between images.
interval_ms = 1000.0/fps
outf = output_file
if outf is None:
outf = 'movie.html'
else:
# Ensure .html extension
outf = os.path.splitext(outf)[0] + '.html'
casename = os.path.splitext(outf)[0]
# Get the necessary HTML text to create the movie.
header, jscode, form, footer, files = \
get_html(files, interval_ms, casename=casename)
# Write the HTML text to the output file.
f = open(outf, 'w')
f.write(header + jscode + form + footer)
f.close()
print("\n\nmovie in output file " + outf)
return None
# -------------------------------------------------------------------------
def get_html(plotfiles, interval_ms=300, width=800, height=600,
casename=None):
"""
Takes a list plotfiles, such as::
'frame00.png', 'frame01.png', ...
and creates javascript code for animating the frames as a movie in HTML.
The `plotfiles` argument can be of three types:
* A Python list of the names of the image files, sorted in correct
order. The names can be filenames of files reachable by the
HTML code, or the names can be URLs.
* A filename generator using Unix wildcard notation, e.g.,
``frame*.png`` (the files most be accessible for the HTML code).
* A filename generator using printf notation for frame numbering
and limits for the numbers. An example is ``frame%0d.png:0->92``,
which means ``frame00.png``, ``frame01.png``, ..., ``frame92.png``.
This specification of `plotfiles` also allows URLs, e.g.,
``http://mysite.net/files/frames/frame_%04d.png:0->320``.
If `casename` is None, a casename based on the full relative path of the
first plotfile is used as tag in the variables in the javascript code
such that the code for several movies can appear in the same file
(i.e., the various code blocks employ different variables because
the variable names differ).
The returned result is text strings that incorporate javascript to
loop through the plots one after another. The html text also features
buttons for controlling the movie.
The parameter `iterval_ms` is the time interval between loading
successive images and is in milliseconds.
The `width` and `height` parameters do not seem to have any effect
for reasons not understood.
The following strings are returned: header, javascript code, form
with movie and buttons, footer, and plotfiles::
header, jscode, form, footer, plotfiles = html_movie('frames*.png')
# Insert javascript code in some HTML file
htmlfile.write(jscode + form)
# Or write a new standalone file that act as movie player
filename = plotfiles[0][:-4] + '.html'
htmlfile = open(filename, 'w')
htmlfile.write(header + jscode + form + footer)
htmlfile.close
This function is based on code written by R. J. LeVeque, based on
a template from Alan McIntyre.
"""
# Alternative method:
# http://stackoverflow.com/questions/9486961/animated-image-with-javascript
# Start with expanding plotfiles if it is a filename generator
if not isinstance(plotfiles, (tuple,list)):
if not isinstance(plotfiles, (str,unicode)):
raise TypeError('plotfiles must be list or filename generator,\
not %s' % type(plotfiles))
filename_generator = plotfiles
if '*' in filename_generator:
# frame_*.png
if filename_generator.startswith('http'):
raise ValueError( 'Filename generator %s cannot contain *;\
must be like\
http://some.net/files/frame_%%04d.png:0->120'\
% filename_generator)
plotfiles = glob.glob(filename_generator)
if not plotfiles:
raise ValueError('No plotfiles on the form' %
filename_generator)
plotfiles.sort()
elif '->' in filename_generator:
# frame_%04d.png:0->120
# http://some.net/files/frame_%04d.png:0->120
p = filename_generator.split(':')
filename = ':'.join(p[:-1])
if not re.search(r'%0?\d+', filename):
raise ValueError( 'Filename generator %s has wrong syntax;\
missing printf specification as in\
frame_%%04d.png:0->120' % filename_generator)
if not re.search(r'\d+->\d+', p[-1]):
raise ValueError( 'Filename generator %s has wrong syntax;\
must be like frame_%%04d.png:0->120'\
% filename_generator)
p = p[-1].split('->')
lo, hi = int(p[0]), int(p[1])
plotfiles = [filename % i for i in range(lo,hi+1,1)]
# Check that the plot files really exist, if they are local on the computer
if not plotfiles[0].startswith('http'):
missing_files = [fname for fname in plotfiles
if not os.path.isfile(fname)]
if missing_files:
raise ValueError('Missing plot files: %s' %
str(missing_files)[1:-1])
if casename is None:
# Use plotfiles[0] as the casename, but remove illegal
# characters in variable names since the casename will be
# used as part of javascript variable names.
casename = os.path.splitext(plotfiles[0])[0]
# Use _ for invalid characters
casename = re.sub('[^0-9a-zA-Z_]', '_', casename)
# Remove leading illegal characters until we find a letter or underscore
casename = re.sub('^[^a-zA-Z_]+', '', casename)
filestem, ext = os.path.splitext(plotfiles[0])
if ext == '.png' or ext == '.jpg' or ext == '.jpeg' or ext == '.gif':
pass
else:
raise ValueError( 'Plotfiles (%s, ...) must be PNG, JPEG, or GIF\
files with extension .png, .jpg/.jpeg, or .gif'\
% plotfiles[0])
header = """\
<html>
<head>
</head>
<body>
"""
no_images = len(plotfiles)
jscode = """
<script language="Javascript">
<!---
var num_images_%(casename)s = %(no_images)d;
var img_width_%(casename)s = %(width)d;
var img_height_%(casename)s = %(height)d;
var interval_%(casename)s = %(interval_ms)d;
var images_%(casename)s = new Array();
function preload_images_%(casename)s()
{
t = document.getElementById("progress");
""" % vars()
i = 0
for fname in plotfiles:
jscode += """
t.innerHTML = "Preloading image ";
images_%(casename)s[%(i)s] = new Image(img_width_%(casename)s, img_height_%(casename)s);
images_%(casename)s[%(i)s].src = "%(fname)s";
""" % vars()
i = i+1
jscode += """
t.innerHTML = "";
}
function tick_%(casename)s()
{
if (frame_%(casename)s > num_images_%(casename)s - 1)
frame_%(casename)s = 0;
document.name_%(casename)s.src = images_%(casename)s[frame_%(casename)s].src;
frame_%(casename)s += 1;
tt = setTimeout("tick_%(casename)s()", interval_%(casename)s);
}
function startup_%(casename)s()
{
preload_images_%(casename)s();
frame_%(casename)s = 0;
setTimeout("tick_%(casename)s()", interval_%(casename)s);
}
function stopit_%(casename)s()
{ clearTimeout(tt); }
function restart_%(casename)s()
{ tt = setTimeout("tick_%(casename)s()", interval_%(casename)s); }
function slower_%(casename)s()
{ interval_%(casename)s = interval_%(casename)s/0.707; }
function faster_%(casename)s()
{ interval_%(casename)s = interval_%(casename)s*0.707; }
// --->
</script>
""" % vars()
plotfile0 = plotfiles[0]
form = """
<form>
<input type="button" value="Start movie" onClick="startup_%(casename)s()">
<input type="button" value="Pause movie" onClick="stopit_%(casename)s()">
<input type="button" value="Resume movie" onClick="restart_%(casename)s()">
<input type="button" value="Slower" onClick="slower_%(casename)s()">
<input type="button" value="Faster" onClick="faster_%(casename)s()">
</form>
<p><div ID="progress"></div></p>
<img src="%(plotfile0)s" name="name_%(casename)s" border=2/>
""" % vars()
footer = '\n</body>\n</html>\n'
return header, jscode, form, footer, plotfiles
***************************************************************************
waves.py
***************************************************************************
# waves.py
# -----------------------------------------------------------------------------
# Generate frames for an animation of a moving Gaussian waves.
# -----------------------------------------------------------------------------
import numpy as np
import matplotlib.pyplot as plt
# get html_movie module from physicalmodelingwithpython.blogspot.com/
from html_movie import movie
# Generate waves for each frame.
# Return a gaussian with specified center and spread using array s.
def gaussian(s, center=0.0, spread=1.0):
return np.exp(-2 * (s - center)**2 / spread**2)
# All lengths are in [m], all times are in [s], and all speeds are in [m/s].
# Define the range of values to display.
xmin, xmax = -4.0, 4.0
ymin, ymax = -3.0, 3.0
# Define array of positions.
dx = 0.01
x = np.arange(xmin, xmax + dx, dx)
# Define the duration and number of frames for the simulation.
tmin, tmax = 0.0, 4.0
num_frames = 100
t = np.linspace(tmin, tmax, num_frames)
# Define the initial position and speed of gaussian waves.
r_speed = 2.0 # speed of right-moving wave
r_0 = -4.0 # initial position of right-moving wave
l_speed = -2.0 # speed of left-moving wave
l_0 = 4.0 # initial position of left-moving wave
# Generate a figure and get access to its axes object.
plt.close('all')
fig = plt.figure(figsize=(10,10))
ax = plt.axes(xlim=(xmin, xmax), ylim=(ymin, ymax))
# Create three empty line objects and grab control.
# The loop below will update the lines in each frame.
ax.plot([], [], 'b--', lw=1) # line for right-moving wave
ax.plot([], [], 'r--', lw=1) # line for left-moving wave
ax.plot([], [], 'g-', lw=3) # line for sum of waves
lines = ax.get_lines() # get list of 3 line objects in plot
# It is essential that the frames be named in alphabetical order.
# {:03d} will display integers with three digits and insert zeros if needed:
# '000_movie.jpg', '001_movie.jpg', etc.
file_name = "{:03d}_movie.jpg"
# Generate frames and save each figure as a separate .jpg file.
for i in range(num_frames):
r_now = r_0 + r_speed * t[i] # update center of right-moving wave
l_now = l_0 + l_speed * t[i] # update center of left-moving wave
yR = gaussian(x, r_now) # get current data for right-moving wave
yL = -gaussian(x, l_now) # get current data for left-moving wave
lines[0].set_data(x,yR) # update right-moving wave
lines[1].set_data(x,yL) # update left-moving wave
lines[2].set_data(x,yR + yL) # update sum of waves
plt.savefig(file_name.format(i)) # save current plot
# Use html movie encoder adapted from scitools to create an HTML document that
# will display the frames as a movie. Open movie.html in Web browser to view.
movie(input_files='*.jpg', output_file='movie.html')
***************************************************************************
convolution.py
***************************************************************************
# convolution.py
# -------------------------------------------------------------------------
# This script creates an eLoG (elongated Laplacian of Gaussian) filter that
# emphasizes long, vertical lines in a figure. The effect of the filter is
# demonstrated on a plus sign.
# -------------------------------------------------------------------------
import numpy as np
import matplotlib.pyplot as plt
import scipy.ndimage as sim
#%% Create a grid of points for the Gaussian filter.
v = np.arange(-25, 26)
X, Y = np.meshgrid(v, v)
#%% Create elongated Gaussian filter, apply Laplacian.
gauss_filter = np.exp(-0.5*(X**2/2 + Y**2/45))
laplace_filter = np.array( [[0, -1, 0], [-1, 4, -1], [0, -1, 0]])
combined_filter = sim.convolve(gauss_filter, laplace_filter)
#%% Create a plus sign '+' to demonstrate effect of filter.
plus = np.zeros((51, 51))
plus[23:28, 25] = 1.0
plus[25, 23:28] = 1.0
plt.figure()
plt.imshow(plus)
plt.gray()
#%% Apply filter to '+' and display resulting image.
# Use vmin/vmax to emphasize features within a restricted range of intensity.
cplus = sim.convolve(plus, combined_filter)
plt.figure()
plt.imshow(cplus, vmin=0, vmax=0.5*cplus.max())
plt.gray()
plt.show()
***************************************************************************
perrin.py
***************************************************************************
# perrin.py
# -------------------------------------------------------------------------
# Generate figure displaying Perrin's experimental data on Brownian motion.
# This script requires the data set 04brownian/g26perrindata.npy.
# -------------------------------------------------------------------------
import numpy as np
import matplotlib.pyplot as plt
data = np.load('g26perrindata.npy')
plt.figure(figsize=(10,10))
plt.plot(data[:,0], data[:,1], 'b.', ms=16)
plt.axis([-20,20,-20,20])
plt.grid(b=True,which='major',ls='-',lw=1.5,c='g')
plt.xlabel(r'$\Delta x$ [$\mu$m]', fontsize=24)
plt.ylabel(r'$\Delta y$ [$\mu$m]', fontsize=24)
ax=plt.gca()
ax.set_xticklabels(ax.get_xticks(), weight='bold', fontsize=16)
ax.set_yticklabels(ax.get_yticks(), weight='bold', fontsize=16)
plt.show()
***************************************************************************
sympy_examples.py
***************************************************************************
# sympy_examples.py
# -------------------------------------------------------------------------
# Demonstrate some useful methods available in the SymPy module.
# Run each command separately at an IPython command line.
# -------------------------------------------------------------------------
from sympy import *
init_session()
expand( (x + y)**5 )
factor( x**6 - 1 )
pi.n(100)
plot(besselj(0, x), besselj(1, x), (x, 0, 10))
diff( x*sin(y), x, y )
integrate(cos(x)**2, x)
integrate(exp(-x**2), (x, -oo, oo))
Sum(k**3, (k, 0, m)).doit().factor()
dsolve(f(x).diff(x) + f(x) - x, f(x)).simplify()
***************************************************************************
first_passage.py
***************************************************************************
# first_passage.py
# -------------------------------------------------------------------------
# Define a function to simulate first passage of a random walker.
# -------------------------------------------------------------------------
import numpy as np
def first_passage(N, L, p=0.5, message=False):
"""
Simulate a random walker that starts at the origin and takes steps to the
right with probability p and to the left with probability 1-p.
Return the number of steps for the first passage of location x==L,
or give up after N steps and return np.nan.
Use message=True to display the results on screen.
"""
rng = np.random.default_rng() # create a random number generator
dx = 2*(rng.random(N) < p) - 1 # generate individual steps
x = np.cumsum(dx) # compute location after each step
at_target = np.nonzero(x==L)[0] # find indices where x == L
if at_target.size > 0:
n = at_target[0] + 1
if message:
print("First passage of x={} occurred after {} steps.".format(L, n))
return n
else:
if message:
print("Did not reach x={} after {} steps.".format(L, N))
return np.nan
***************************************************************************
data_dictionary.py
***************************************************************************
# data_dictionary.py
# -------------------------------------------------------------------------
# Store input and data two different simulations in a dictionary.
# This script requires first_passage.py to be in the same directory.
# -------------------------------------------------------------------------
from first_passage import first_passage
data = {} # empty dictionary to store all data
data['A'] = {} # empty dictionary within data to store Simulation A
data['B'] = {} # empty dictionary within data to store Simulation B
#%% Define and run simulations.
samples = 500
data['A']['input'] = dict(N=1000, L=10, p=0.5)
data['A']['results'] = \
[ first_passage(**data['A']['input']) for n in range(samples) ]
data['B']['input'] = dict(N=1000, L=20, p=0.5)
data['B']['results'] = \
[ first_passage(**data['B']['input']) for n in range(samples) ]
#%% Run more simulations. Use "+=" to append new list to old.
data['A']['results'] += \
[ first_passage(**data['A']['input']) for n in range(samples) ]
data['B']['results'] += \
[ first_passage(**data['B']['input']) for n in range(samples) ]
***************************************************************************
nd_random_walks.py
***************************************************************************
# nd_random_walks.py
# -------------------------------------------------------------------------
# Class to simulate various random walks in N dimensions.
# -------------------------------------------------------------------------
"""
Description
-----------
This module implements several N-dimensional random walks.
Each RandomWalk class has these standard methods:
- get_walk(N) --> N steps in a random walk
- get_endpoints(M,N) --> endpoints from M walks of N steps
- get_distances(M,N) --> final distances from M walks of N steps
The class hierarchy is as follows.
RandomWalk -- general random walk class
LatticeWalk -- random walk on a D-dimensional lattice
default is a cubic lattice in D dimensions
TriangularWalk -- 2D triangular lattice
HoneycombWalk -- 2D honeycomb lattice
DirectionalWalk -- random walks of variable step size in D dimensions
defaults to constant step length in random directions
UniformWalk -- step size drawn from uniform distribution
GaussianWalk -- step size drawn from normal distribution
ExponentialWalk -- step size drawn from exponential distribution
ParetoWalk -- step size drawn from power law distribution
"""
import numpy as np
class RandomWalk():
def __init__(self, dimension=1):
# Create a random number generator.
self.rng = np.random.default_rng()
# Store dimension and starting point internally.
self.D = dimension
self.r0 = np.zeros((self.D,1))
def _get_steps(self,N):
# Internal method. Generate individual steps of the random walk.
return np.zeros((self.D, N)) # placeholder for parent class
def get_walk(self,N):
return np.append(self.r0, np.cumsum(self._get_steps(N), 1), axis=1)
def get_endpoints(self,M,N):
return np.array([self.get_walk(N)[:,-1] for m in range(M)]).transpose()
def get_distances(self,M,N):
return np.sqrt(np.sum(self.get_endpoints(M,N)**2, axis=0))
class LatticeWalk(RandomWalk):
def __init__(self, dimension=1):
super().__init__(dimension)
# Create list of legal steps on D-dimensional cubic lattice.
M = np.eye(self.D)
self.basis = [M[:,n] for n in range(self.D)]
self.basis += [-M[:,n] for n in range(self.D)]
def _get_steps(self, N):
# Generate individual steps by randomly choosing from self.basis.
return self.rng.choice(self.basis, size=N).transpose()
class TriangularWalk(LatticeWalk):
def __init__(self, dimension=2):
super().__init__(dimension=2)
# Create list of legal steps on a triangular lattice:
# unit vector along x-axis, plus images under successive rotations.
cosTheta, sinTheta = np.cos(np.pi/3), np.sin(np.pi/3)
R = np.array([[cosTheta, -sinTheta], [sinTheta, cosTheta]])
v = np.array([0,1])
self.basis = []
for n in range(6):
self.basis += [v]
v = np.dot(R,v)
class HoneycombWalk(LatticeWalk):
def __init__(self, dimension=2):
super().__init__(dimension=2)
# Create list of legal steps on a triangular lattice:
# unit vector along x-axis, plus images under successive rotations
# generates on sublattice. _get_steps() accounts for the other.
cosTheta, sinTheta = np.cos(2*np.pi/3), np.sin(2*np.pi/3)
R = np.array([[cosTheta, -sinTheta], [sinTheta, cosTheta]])
v = np.array([0,1])
self.basis = []
for n in range(3):
self.basis += [v]
v = np.dot(R,v)
def _get_steps(self, N):
# Honeycomb lattice has two sublattices. (-1)**n accounts for this:
# Rotate step by 180 degrees if on sublattice B.
return self.rng.choice(self.basis, size=N).T * (-1)**np.arange(N)
class DirectionalWalk(RandomWalk):
def _direction(self, N):
# Use the Box-Muller transform to generate N random unit vectors.
vectors = self.rng.normal(size=(self.D,N))
lengths = np.sqrt(np.sum(vectors**2,0))
return vectors/lengths
def _magnitude(self, N):
# Only directions are random for this class. Step size is 1.
return np.ones(N)
def _get_steps(self, N):
# Return steps for random walk of N steps in D dimensions.
return self._direction(N) * self._magnitude(N)
class UniformWalk(DirectionalWalk):
def _magnitude(self, N):
# Draw N step sizes from the uniform distribution. Mean size is 1.
return 2*self.rng.random(size=N)
class GaussianWalk(DirectionalWalk):
def _magnitude(self, N):
# Draw N step sizes from chi distribution. Mean size is 1.
return np.sqrt(self.rng.chisquare(df=self.D, size=N)/self.D)
class ExponentialWalk(DirectionalWalk):
def _magnitude(self, N):
# Draw N steps from exponential distribution. Mean size is 1.
return self.rng.exponential(size=N)
class ParetoWalk(DirectionalWalk):
def __init__(self, dimension=1, nu=2):
super().__init__(dimension)
# Exponent and normalization factor for power law distribution.
self.nu = nu
self.norm = max(0.01, nu-1)
def _magnitude(self, N):
# Draw N step sizes from Pareto distribution. Mean is 1 if nu > 1.01.
return self.rng.pareto(a=self.nu, size=N) * self.norm
***************************************************************************
surprise.py
***************************************************************************
# surprise.py
# -------------------------------------------------------------------------
# This script will create a familar but interesting image.
# It may take about a minute to run.
# -------------------------------------------------------------------------
import numpy as np
import matplotlib.pyplot as plt
# This determines the number of colors used in the plot.
# The larger the value, the longer the script will take.
max_iterations = 32
# These parameters define the boundaries of the plot
x_min, x_max = -2.5, 1.5
y_min, y_max = -1.5, 1.5
# This parameter controls the grid spacing. A smaller value gives better
# resolution, but the script will take longer to run.
ds = 0.002
X = np.arange(x_min, x_max + ds, ds)
Y = np.arange(y_min, y_max + ds, ds)
data = np.zeros( (X.size, Y.size), dtype='uint')
for i in range(X.size):
for j in range(Y.size):
x0, y0 = X[i], Y[j]
x, y = x0, y0
count = 0
while count < max_iterations:
# Update x and y simultaneously.
x, y = (x0 + x*x - y*y, y0 + 2*x*y)
# Exit loop if (x,y) is too far from the origin.
if (x*x + y*y) > 4.0: break
count += 1
data[i, j] = max_iterations - count
plt.imshow(data.transpose(), interpolation='nearest', cmap='jet')
plt.axis('off')
plt.show()
***************************************************************************
scope.py
***************************************************************************
# scope.py
# -------------------------------------------------------------------------
# Demonstrate Python's rules of scope.
# -------------------------------------------------------------------------
def illustrate_scope():
s_enclosing = 'E'
def display():
s_local = 'L'
print( "Local --> {}".format(s_local) )
print( "Enclosing --> {}".format(s_enclosing) )
print( "Global --> {}".format(s_global) )
print( "Built-in --> {}".format(abs) )
display()
s_global = 'G'
illustrate_scope()
***************************************************************************
name_collision.py
***************************************************************************
# name_collision.py
# -------------------------------------------------------------------------
# Demonstrate how Python's rules of scope prevent name collisions.
# -------------------------------------------------------------------------
def name_collisions():
x, y = 'E', 'E'
def display():
x = 'L'
print( "Inside display() ..." )
print( "x= {}\ny= {}\nz= {}".format(x,y,z) )
display()
print( "Inside name_collision() ..." )
print( "x= {}\ny= {}\nz= {}".format(x,y,z) )
x, y, z = 'G', 'G', 'G'
name_collisions()
print( "Outside function ..." )
print( "x= {}\ny= {}\nz= {}".format(x,y,z) )
***************************************************************************
fancy_plot.py
***************************************************************************
# fancy_plot.py
# -------------------------------------------------------------------------
# Add a title and axis labels to a simple plot.
# -------------------------------------------------------------------------
import numpy as np
import matplotlib.pyplot as plt
x_min, x_max = -4, 4
num_points = 51
x_list = np.linspace(x_min, x_max, num_points)
y_list = x_list**2
plt.figure()
plt.plot(x_list, y_list, 'r', linewidth=3)
ax = plt.gca()
ax.set_title("A Second Order Polynomial", fontsize=16)
ax.set_xlabel("$x$", fontsize=24)
ax.set_ylabel("$y = x^2$", fontsize=24)
plt.show()
***************************************************************************
legend.py
***************************************************************************
# legend.py
# -------------------------------------------------------------------------
# Create a plot with a legend to distinguish multiple curves.
# -------------------------------------------------------------------------
import numpy as np
import matplotlib.pyplot as plt
num_curves = 3
x = np.linspace(0, 1, 201)
y = np.zeros((x.size, num_curves))
for n in range(num_curves):
y[:, n] = np.sin((n+1) * x * 2 * np.pi)
plt.figure()
plt.plot(x, y, linewidth=2)
ax = plt.gca()
ax.legend( ("sin(2$\\pi$x)", "sin(4$\\pi$x)", "sin(6$\\pi$x)") )
plt.show()
***************************************************************************
random_walk.py
***************************************************************************
# random_walk.py
# -------------------------------------------------------------------------
# Monte Carlo simulation of a two-dimensional random walk.
# -------------------------------------------------------------------------
import numpy as np
import matplotlib.pyplot as plt
from numpy.random import random as rng
num_steps = 500
x_step = rng(num_steps) > 0.5
y_step = rng(num_steps) > 0.5
x_step = 2*x_step - 1
y_step = 2*y_step - 1
x_position = np.cumsum(x_step)
y_position = np.cumsum(y_step)
plt.figure()
plt.plot(x_position, y_position)
plt.axis('equal')
# A more succinct alternative:
x = ( 2*(rng(num_steps) > 0.5) - 1 ).cumsum()
y = ( 2*(rng(num_steps) > 0.5) - 1 ).cumsum()
plt.figure()
plt.plot(x,y)
plt.axis('equal')
plt.show()
***************************************************************************
surface.py
***************************************************************************
# surface.py
# -------------------------------------------------------------------------
# Create a three-dimensional surface plot.
# -------------------------------------------------------------------------
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
# define grid of points
points = np.linspace(-1, 1, 101)
X, Y = np.meshgrid(points, points)
Z = X**2 + Y**2
# create and display surface plot
ax = Axes3D(plt.figure())
ax.plot_surface(X, Y, Z)
# Adjust rstride and cstride to use more or fewer points.
# The following command will use all points in X, Y, and Z:
# ax.plot_surface(X, Y, Z, rstride=1, cstride=1)
plt.show()
***************************************************************************
regression.py
***************************************************************************
# regression.py
# -------------------------------------------------------------------------
# Example of linear regression on data from first passage problem.
# first_passage.py must be in the working directory.
# Increase samples and nmax for better statistics. This script will take a
# several minutes to finish if (samples x nmax) ~ 10**9 or more.
# -------------------------------------------------------------------------
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from sklearn.linear_model import LinearRegression
from first_passage import first_passage
#%% Define and run simulations.
data = {} # empty dictionary to store all data
nmax = 10**5
parameters = dict(N=nmax, p=0.5) # common inputs
l_values = [10, 20, 50, 100, 200, 500]
samples = 10**3
for l in l_values:
data["L={}".format(l)] = \
[first_passage(L=l, **parameters) for n in range(samples)]
step_data = pd.DataFrame(data)
#%% Prepare data for modeling with sklearn.
model = LinearRegression()
X = np.array(l_values).reshape(-1,1)
Y = step_data.mean()
logX = np.log(X)
logY = np.log(Y)
xFit = np.linspace(1, X.max(), 201).reshape(-1,1)
#%% Plot data and fit models.
plt.figure()
plt.plot(X, Y, 'ko', mew=2, mfc='none', label="Data")
# Check for linear relationship.
model.fit(X, Y)
print("Linear Model: A + b*x")
print("A = {:.3f}, b = {:.3f}".format(model.intercept_, model.coef_[0]))
print("R2: ", model.score(X,Y))
yFit = model.predict(xFit)
plt.plot(xFit, yFit, 'r-', label="Linear Model")
# Check for exponential relationship.
model.fit(X, logY)
print("Exponential Model: A * exp(b*x)")
print("A = {:.3f}, b = {:.3f}".format(model.intercept_, model.coef_[0]))
print("R2: ", model.score(X,Y))
yFit = model.predict(xFit)
plt.plot(xFit, np.exp(yFit), 'g-', label="Exponential Model")
# Check for power-law relationship.
model.fit(logX, logY)
print("Power Law Model: A * x**b")
print("A = {:.3f}, b = {:.3f}".format(model.intercept_, model.coef_[0]))
print("R2: ", model.score(X,Y))
yFit = model.predict(np.log(xFit))
plt.plot(xFit, np.exp(yFit), 'b-', label="Power Law Model")
plt.xlabel("Target Distance")
plt.ylabel("Average Steps")
plt.legend()
plt.show()
***************************************************************************
bar3d.py
***************************************************************************
# bar3d.py
# -------------------------------------------------------------------------
# Create a 3D "Lego plot" to illustrate the Central Limit Theorem.
# -------------------------------------------------------------------------
"""
Description
-----------
This script demonstrates the creation of a 3D bar chart and it illustrates how
the distribution of a sum of random numbers differs from the distribution those
numbers are drawn from.
The script will bin num_points random (x,y) pairs in a 2D histogram, then
display a 3D bar chart of this data. Each x and y coordinate is the sum
of num_to_sum numbers drawn from the uniform random distribution on [0,1).
Adjust num_to_sum and use different random number generators to explore the
Central Limit Theorem. (The Central Limit Theorem implies that the
distribution of the *sum* of many random numbers should approach a normal
distribution.)
"""
# %% Setup
import numpy as np
import matplotlib.pyplot as plt
# Create a random number generator.
rng = np.random.default_rng() # create a random number generator object
rand = rng.random # assign its uniform distribution method to rand
# Import the 3D plotter.
from mpl_toolkits.mplot3d import Axes3D
# %% Generate data and histogram.
# Generate random (x,y) pairs to bin. Each coordinate is the sum of N numbers
# drawn from the uniform distribution on [0,1).
num_points = 10000
num_to_sum = 2
num_bins = 20
x_vals = rand((num_points, num_to_sum)).sum(1)
y_vals = rand((num_points, num_to_sum)).sum(1)
# Use histogram2d to get counts and bin locations.
counts, x_bins, y_bins = np.histogram2d(x_vals, y_vals, bins=num_bins)
# %% Transform data for plotting 3D bar chart.
# For 3D bar chart, we need starting (x,y,z) coordinates for each bar and
# (dx,dy,dz) to specify its shape. The following code uses meshgrid, diff, and
# flatten to construct these inputs from the output of histogram2d.
x_start, y_start = np.meshgrid(x_bins[:-1], y_bins[:-1])
z_start = np.zeros_like(x_start)
x_length, y_length = np.meshgrid(np.diff(x_bins), np.diff(y_bins))
z_length = counts
x = x_start.flatten()
y = y_start.flatten()
z = z_start.flatten()
dx = x_length.flatten()
dy = y_length.flatten()
dz = z_length.flatten()
# %% Plot transformed data.
fig = plt.figure()
ax3d = Axes3D(fig)
ax3d.bar3d(x, y, z, dx, dy, dz)
plt.show()
***************************************************************************
perrin.py
***************************************************************************
# perrin.py
# -------------------------------------------------------------------------
# Generate figure displaying Perrin's experimental data on Brownian motion.
# This script requires the data set 04brownian/g26perrindata.npy.
# -------------------------------------------------------------------------
import numpy as np
import matplotlib.pyplot as plt
data = np.load('g26perrindata.npy')
plt.figure(figsize=(10,10))
plt.plot(data[:,0], data[:,1], 'b.', ms=16)
plt.axis([-20,20,-20,20])
plt.grid(b=True,which='major',ls='-',lw=1.5,c='g')
plt.xlabel(r'$\Delta x$ [$\mu$m]', fontsize=24)
plt.ylabel(r'$\Delta y$ [$\mu$m]', fontsize=24)
ax=plt.gca()
ax.set_xticklabels(ax.get_xticks(), weight='bold', fontsize=16)
ax.set_yticklabels(ax.get_yticks(), weight='bold', fontsize=16)
plt.show()
***************************************************************************
shading.py
***************************************************************************
# shading.py
# -------------------------------------------------------------------------
# Demonstrate shading of surface plots using Matplotlib's LightSource.
# -------------------------------------------------------------------------
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
# Import Bessel function.
from scipy.special import jn
# Import colormaps.
from matplotlib import cm
# Import lighting object for shading surface plots.
from matplotlib.colors import LightSource
# Define grid of points.
points = np.linspace(-10, 10, 101)
X, Y = np.meshgrid(points, points)
R = np.sqrt(X**2 + Y**2)
Z = jn(0,R)
# Create an rgb array for single-color surfaces.
white = np.ones((Z.shape[0], Z.shape[1], 3))
red = white * np.array([1,0,0])
green = white * np.array([0,1,0])
blue = white * np.array([0,0,1])
# Set view parameters for all subplots.
azimuth = 45
altitude = 60
# Create empty figure.
fig = plt.figure(figsize=(18,12))
# -------------------------------------------------------------------------
# Generate first subplot.
# -------------------------------------------------------------------------
# Create a light source object for light from
# 0 degrees azimuth, 0 degrees elevation.
light = LightSource(0, 0)
# Generate face colors for a shaded surface using either
# a color map or the uniform rgb color specified above.
illuminated_surface = light.shade_rgb(red, Z)
# Create a subplot with 3d plotting capabilities.
# This command will fail if Axes3D was not imported.
ax = fig.add_subplot(2,2,1, projection='3d')
ax.view_init(altitude, azimuth)
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, linewidth=0,
antialiased=False, facecolors=illuminated_surface)
# -------------------------------------------------------------------------
# Repeat the commands above for the other three subplots, but use different
# illumination angles and colors.
# -------------------------------------------------------------------------
light = LightSource(90, 0)
illuminated_surface = light.shade_rgb(green, Z)
ax = fig.add_subplot(2,2,2, projection='3d')
ax.view_init(altitude, azimuth)
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, linewidth=0,
antialiased=False, facecolors=illuminated_surface)
# -------------------------------------------------------------------------
light = LightSource(90, 45)
illuminated_surface = light.shade_rgb(blue, Z)
ax = fig.add_subplot(2,2,3, projection='3d')
ax.view_init(altitude, azimuth)
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, linewidth=0,
antialiased=False, facecolors=illuminated_surface)
# -------------------------------------------------------------------------
light = LightSource(180, 45)
illuminated_surface = light.shade(Z, cmap=cm.coolwarm)
ax = fig.add_subplot(2,2,4, projection='3d')
ax.view_init(altitude, azimuth)
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, linewidth=0,
antialiased=False, facecolors=illuminated_surface)
# -------------------------------------------------------------------------
plt.tight_layout()
plt.savefig('shading.png')
Comments