Revision: 6270
http://matplotlib.svn.sourceforge.net/matplotlib/?rev=6270&view=rev
Author: fer_perez
Date: 2008-10-19 07:53:47 +0000 (Sun, 19 Oct 2008)
Log Message:
-----------
Cleanup and updates to many examples.
Added new tool to create skeletons and solutions from a single source
script.
Modified Paths:
--------------
trunk/py4science/examples/bessel.py
trunk/py4science/examples/butter_filter.py
trunk/py4science/examples/fft_imdenoise.py
trunk/py4science/examples/filtilt_demo.py
trunk/py4science/examples/fitting.py
trunk/py4science/examples/montecarlo_pi.py
trunk/py4science/examples/qsort.py
trunk/py4science/examples/quad_newton.py
trunk/py4science/examples/recarray_demo.py
trunk/py4science/examples/trapezoid.py
Added Paths:
-----------
trunk/py4science/examples/mkskel.py
trunk/py4science/examples/skel/fortran_wrap/test_example.py
trunk/py4science/examples/skel/fortran_wrap/test_fib.py
Modified: trunk/py4science/examples/bessel.py
===================================================================
--- trunk/py4science/examples/bessel.py 2008-10-19 07:48:51 UTC (rev 6269)
+++ trunk/py4science/examples/bessel.py 2008-10-19 07:53:47 UTC (rev 6270)
@@ -1,45 +1,46 @@
#!/usr/bin/env python
"""Plot some Bessel functions of integer order, using Scipy and pylab"""
-import scipy.special
+from scipy import special
-import numpy as N
-import scipy as S
-import pylab as P
+import numpy as np
+import matplotlib.pyplot as plt
-# shorthand
-special = S.special
-
def jn_asym(n,x):
"""Asymptotic form of jn(x) for x>>n"""
+
+ #@ The asymptotic formula is:
+ #@ j_n(x) ~ sqrt(2.0/pi/x)*cos(x-(n*pi/2.0+pi/4.0))
- return S.sqrt(2.0/S.pi/x)*S.cos(x-(n*S.pi/2.0+S.pi/4.0))
+ return np.sqrt(2.0/np.pi/x)*np.cos(x-(n*np.pi/2.0+np.pi/4.0)) #@
# build a range of values to plot in
-x = N.linspace(0,30,400)
+x = np.linspace(0,30,400)
# Start by plotting the well-known j0 and j1
-P.figure()
-P.plot(x,special.j0(x),label='$J_0$')
-P.plot(x,special.j1(x),label='$J_1$')
+plt.figure()
+plt.plot(x,special.j0(x),label='$J_0$') #@
+plt.plot(x,special.j1(x),label='$J_1$') #@
# Show a higher-order Bessel function
n = 5
-P.plot(x,special.jn(n,x),label='$J_%s$' % n)
+plt.plot(x,special.jn(n,x),label='$J_%s$' % n)
# and compute its asymptotic form (valid for x>>n, where n is the order). We
-# must first find the valid range of x where at least x>n:
-x_asym = S.compress(x>n,x)
-P.plot(x_asym,jn_asym(n,x_asym),label='$J_%s$ (asymptotic)' % n)
+# must first find the valid range of x where at least x>n.
+#@ Find where x>n and evaluate the asymptotic relation only there
+x_asym = x[x>n] #@
+plt.plot(x_asym,jn_asym(n,x_asym),label='$J_%s$ (asymptotic)' % n) #@
# Finish off the plot
-P.legend()
-P.title('Bessel Functions')
+plt.legend()
+plt.title('Bessel Functions')
# horizontal line at 0 to show x-axis, but after the legend
-P.axhline(0)
+plt.axhline(0)
-# EXERCISE: redo the above, for the asymptotic range 0<x<<n. The asymptotic
-# form in this regime is
+# Extra credit: redo the above, for the asymptotic range 0<x<<n. The
+# asymptotic form in this regime is:
+#
# J(n,x) = (1/gamma(n+1))(x/2)^n
# Now, let's verify numerically the recursion relation
@@ -47,17 +48,25 @@
jn = special.jn # just a shorthand
# Be careful to check only for x!=0, to avoid divisions by zero
-xp = S.compress(x>0.0,x) # positive x
+#@ xp contains the positive values of x
+xp = x[x>0.0] #@
# construct both sides of the recursion relation, these should be equal
+#@ Define j_np1 to hold j_(n+1) evaluated at the points xp
j_np1 = jn(n+1,xp)
+#@ Define j_np1_rec to express j_(n+1) via a recursion relation, at points xp
j_np1_rec = (2.0*n/xp)*jn(n,xp)-jn(n-1,xp)
# Now make a nice error plot of the difference, in a new figure
-P.figure()
-P.semilogy(xp,abs(j_np1-j_np1_rec),'r+-')
-P.title('Error in recursion for $J_%s$' % n)
-P.grid()
+plt.figure()
+#@ We now plot the difference between the two formulas above. Note that to
+#@ properly display the errors, we want to use a logarithmic y scale. Search
+#@ the matplotlib docs for the proper calls.
+plt.semilogy(xp,abs(j_np1-j_np1_rec),'r+-') #@
+
+plt.title('Error in recursion for $J_%s$' % n)
+plt.grid()
+
# Don't forget a show() call at the end of the script
-P.show()
+plt.show()
Modified: trunk/py4science/examples/butter_filter.py
===================================================================
--- trunk/py4science/examples/butter_filter.py 2008-10-19 07:48:51 UTC (rev
6269)
+++ trunk/py4science/examples/butter_filter.py 2008-10-19 07:53:47 UTC (rev
6270)
@@ -1,13 +1,13 @@
-import numpy as n
+import numpy as np
import scipy.signal as signal
-from pylab import figure, show
+from matplotlib.pyplot import figure, show
dt = 0.01
-t = n.arange(0, 2, dt)
-s = n.sin(2*n.pi*t)
+t = np.arange(0, 2, dt)
+s = np.sin(2*np.pi*t)
# sine corrupted wih gaussian white noise
-sn = s + 0.1*n.random.randn(len(s)) # noisy sine
+sn = s + 0.1*np.random.randn(len(s)) # noisy sine
# the nyquist frequency 1/(2dt) is the maximum frequency in a sampled
# signal
Modified: trunk/py4science/examples/fft_imdenoise.py
===================================================================
--- trunk/py4science/examples/fft_imdenoise.py 2008-10-19 07:48:51 UTC (rev
6269)
+++ trunk/py4science/examples/fft_imdenoise.py 2008-10-19 07:53:47 UTC (rev
6270)
@@ -1,62 +1,86 @@
#!/usr/bin/env python
-"""Image denoising example using 2-dimensional FFT."""
+"""Simple image denoising example using 2-dimensional FFT."""
import sys
import numpy as np
from matplotlib import pyplot as plt
-def mag_phase(F):
- """Return magnitude and phase components of spectrum F."""
-
- return (np.absolute(F), np.angle(F))
-
def plot_spectrum(F, amplify=1000):
"""Normalise, amplify and plot an amplitude spectrum."""
- M, Phase = mag_phase(F)
- M *= amplify/M.max()
- M[M > 1] = 1
- print M.shape, M.dtype
- plt.imshow(M, plt.cm.Blues)
+ # Note: the problem here is that we have a spectrum whose histogram is
+ # *very* sharply peaked at small values. To get a meaningful display, a
+ # simple strategy to improve the display quality consists of simply
+ # amplifying the values in the array and then clipping.
+ # Compute the magnitude of the input F (call it mag). Then, rescale mag by
+ # amplify/maximum_of_mag. Numpy arrays can be scaled in-place with ARR *=
+ # number. For the max of an array, look for its max method.
+ mag = abs(F) #@
+ mag *= amplify/mag.max() #@
+
+ # Next, clip all values larger than one to one. You can set all elements
+ # of an array which satisfy a given condition with array indexing syntax:
+ # ARR[ARR<VALUE] = NEWVALUE, for example.
+ mag[mag > 1] = 1 #@
+ # Display: this one already works, if you did everything right with mag
+ plt.imshow(mag, plt.cm.Blues)
+
if __name__ == '__main__':
try:
# Read in original image, convert to floating point for further
# manipulation; imread returns a MxNx4 RGBA image. Since the image is
- # grayscale, just extrac the 1st channel
- im = plt.imread('data/moonlanding.png').astype(float)[:,:,0]
+ # grayscale, just extract the 1st channel
+ #@ Hints:
+ #@ - use plt.imread() to load the file
+ #@ - convert to a float array with the .astype() method
+ #@ - extract all rows, all columns, 0-th plane to get the first
+ #@ channel
+ #@ - the resulting array should have 2 dimensions only
+ im = plt.imread('data/moonlanding.png').astype(float)[:,:,0] #@
+ print "Image shape:",im.shape
except:
print "Could not open image."
sys.exit(-1)
# Compute the 2d FFT of the input image
- F = np.fft.fft2(im)
+ #@ Hint: Look for a 2-d FFT in np.fft.
+ #@ Note: call this variable 'F', which is the name we'll be using below.
+ F = np.fft.fft2(im) #@
- # Now, make a copy of the original spectrum and truncate coefficients.
+ # In the lines following, we'll make a copy of the original spectrum and
+ # truncate coefficients. NO immediate code is to be written right here.
+
+ # Define the fraction of coefficients (in each direction) we keep
keep_fraction = 0.1
# Call ff a copy of the original transform. Numpy arrays have a copy
# method for this purpose.
- ff = F.copy()
+ ff = F.copy() #@
- # Set r and c to be the number of rows and columns of the array.
- r,c = ff.shape
+ # Set r and c to be the number of rows and columns of the array.
+ #@ Hint: use the array's shape attribute.
+ r,c = ff.shape #@
# Set to zero all rows with indices between r*keep_fraction and
# r*(1-keep_fraction):
- ff[r*keep_fraction:r*(1-keep_fraction)] = 0
+ ff[r*keep_fraction:r*(1-keep_fraction)] = 0 #@
# Similarly with the columns:
- ff[:, c*keep_fraction:c*(1-keep_fraction)] = 0
+ ff[:, c*keep_fraction:c*(1-keep_fraction)] = 0 #@
# Reconstruct the denoised image from the filtered spectrum, keep only the
- # real part for display
- im_new = np.fft.ifft2(ff).real
-
+ # real part for display.
+ #@ Hint: There's an inverse 2d fft in the np.fft module as well (don't
+ #@ forget that you only want the real part).
+ #@ Call the result im_new,
+ im_new = np.fft.ifft2(ff).real #@
+
# Show the results
+ #@ The code below already works, if you did everything above right.
plt.figure()
plt.subplot(221)
@@ -75,9 +99,6 @@
plt.title('Reconstructed Image')
plt.imshow(im_new, plt.cm.gray)
- plt.savefig('fft_imdenoise.png', dpi=150)
- plt.savefig('fft_imdenoise.eps')
-
# Adjust the spacing between subplots for readability
plt.subplots_adjust(hspace=0.32)
plt.show()
Modified: trunk/py4science/examples/filtilt_demo.py
===================================================================
--- trunk/py4science/examples/filtilt_demo.py 2008-10-19 07:48:51 UTC (rev
6269)
+++ trunk/py4science/examples/filtilt_demo.py 2008-10-19 07:53:47 UTC (rev
6270)
@@ -7,10 +7,11 @@
order. The function also computes the initial filter parameters in
order to provide a more stable response (via lfilter_zi). The
following code has been tested with Python 2.4.4 and Scipy 0.5.1.
+"""
-"""
-from numpy import vstack, hstack, eye, ones, zeros, linalg, \
-newaxis, r_, flipud, convolve, matrix, array
+from numpy import (vstack, hstack, eye, ones, zeros, linalg,
+ newaxis, r_, flipud, convolve, matrix, array)
+
from scipy.signal import lfilter
def lfilter_zi(b,a):
@@ -38,19 +39,18 @@
return array(zi_return)
-
-
def filtfilt(b,a,x):
#For now only accepting 1d arrays
ntaps=max(len(a),len(b))
edge=ntaps*3
if x.ndim != 1:
- raise ValueError, "Filiflit is only accepting 1 dimension arrays."
+ raise ValueError("Filiflit is only accepting 1 dimension arrays.")
- #x must be bigger than edge
+ # x must be bigger than edge
if x.size < edge:
- raise ValueError, "Input vector needs to be bigger than 3 *
max(len(a),len(b)."
+ e="Input vector needs to be bigger than 3 * max(len(a),len(b)."
+ raise ValueError(e)
if len(a) < ntaps:
a=r_[a,zeros(len(b)-len(a))]
Modified: trunk/py4science/examples/fitting.py
===================================================================
--- trunk/py4science/examples/fitting.py 2008-10-19 07:48:51 UTC (rev
6269)
+++ trunk/py4science/examples/fitting.py 2008-10-19 07:53:47 UTC (rev
6270)
@@ -7,10 +7,11 @@
from scipy.optimize import leastsq
from scipy.interpolate import splrep,splev
-import numpy as N
-import scipy as S
-import pylab as P
+import numpy as np
+import matplotlib as mpl
+import matplotlib.pyplot as plt
+
def func(pars):
a, alpha, k = pars
return a*exp(alpha*x_vals) + k
@@ -35,7 +36,7 @@
print 'Least-squares fit to the data'
print 'true', pars_true
print 'best', best
-print '|err|_l2 =',P.l2norm(pars_true-best)
+print '|err|_l2 =',np.linalg.norm(pars_true-best)
# scipy's splrep uses FITPACK's curfit (B-spline interpolation)
print
@@ -48,27 +49,27 @@
def plot_polyfit(x,y,n,fignum=None):
""" """
if fignum is None:
- fignum = P.figure().number
- P.plot(x,y,label='Data')
+ fignum = plt.figure().number
+ plt.plot(x,y,label='Data')
- fit_coefs = N.polyfit(x,y,n)
- fit_val = N.polyval(fit_coefs,x)
- P.plot(x,fit_val,label='Polynomial fit, $n=%d$' % n)
- P.legend()
+ fit_coefs = np.polyfit(x,y,n)
+ fit_val = np.polyval(fit_coefs,x)
+ plt.plot(x,fit_val,label='Polynomial fit, $n=%d$' % n)
+ plt.legend()
return fignum
# Now use pylab to plot
-P.figure()
-P.plot(x_vals,y_noisy,label='Noisy data')
-P.plot(x_vals,func(best),lw=2,label='Least-squares fit')
-P.legend()
-P.figure()
-P.plot(x_vals,y_noisy,label='Noisy data')
-P.plot(x_vals,smooth,lw=2,label='Spline-smoothing')
-P.legend()
+plt.figure()
+plt.plot(x_vals,y_noisy,label='Noisy data')
+plt.plot(x_vals,func(best),lw=2,label='Least-squares fit')
+plt.legend()
+plt.figure()
+plt.plot(x_vals,y_noisy,label='Noisy data')
+plt.plot(x_vals,smooth,lw=2,label='Spline-smoothing')
+plt.legend()
fignum = plot_polyfit(x_vals,y_noisy,1)
plot_polyfit(x_vals,y_noisy,2,fignum)
plot_polyfit(x_vals,y_noisy,3,fignum)
-P.show()
+plt.show()
Added: trunk/py4science/examples/mkskel.py
===================================================================
--- trunk/py4science/examples/mkskel.py (rev 0)
+++ trunk/py4science/examples/mkskel.py 2008-10-19 07:53:47 UTC (rev 6270)
@@ -0,0 +1,325 @@
+#!/usr/bin/env python
+"""Make skeletons out of Python scripts.
+
+Usage:
+
+ mkskel [--test] file1.py file2.py ....
+
+If --test is given, the test suite is run instead.
+
+For each input filename f.py, a pair of output files is generated, f_soln.py
+and f_skel.py.
+
+Source markup is very simple. The tests in the file show precisely how it
+works, but in summary:
+
+- Pure comment lines with the special marker (#@) are left in the skeleton
+ (only the marker is stripped, but they remain as valid comments). These are
+ typically used for hints.
+
+- Code lines terminated with the marker are:
+
+ - In the skeleton, replaced by a NotImplementedError call. Consecutive lines
+ are replaced by a single call.
+
+ - In the solution, kept as is but the marker is removed.
+"""
+
+from __future__ import with_statement
+
+#-----------------------------------------------------------------------------
+# Stdlib imports
+import os
+import re
+import shutil
+import sys
+
+# Third-party imports
+import nose
+
+# Constants
+MARKER = '#@'
+DEL_RE = re.compile(r'''^((\s*)(.*?))\s*%s\s*$''' % MARKER)
+HINT_RE = re.compile(r'''^(?P<space>\s*)%s\s+(?P<hint>.*)$''' % MARKER)
+
+
+#-----------------------------------------------------------------------------
+# Main code begins
+
+def src2soln(src):
+ """Remove markers from input source, leaving all else intact.
+
+ Inputs:
+ src : sequence of lines (file-like objects work out of the box)
+ """
+
+ out = []
+ addline = out.append
+ for line in src:
+ # Check for lines to delete and with hints
+ mdel = DEL_RE.match(line)
+ mhint = HINT_RE.match(line)
+
+ # All hints are unconditionally removed
+ if mhint is None:
+ if mdel:
+ # if marker is matched in code, strip it and leave the code
+ line = mdel.group(1)+'\n'
+ addline(line)
+
+ return ''.join(out)
+
+
+def src2skel(src):
+ """Remove markers from input source, replacing marked lines.
+
+ Marked lines are replaced with "raise NotImplementedError" calls that
+ summarize the total number of deleted lines.
+
+ Inputs:
+ src : sequence of lines (file-like objects work out of the box)
+ """
+
+ def flush_buffers(normal_lines,del_lines=0):
+ """Local function to reuse some common code"""
+
+ if state_cur == normal:
+ # add the normal lines
+ out.extend(normal_lines)
+ normal_lines[:] = []
+ else:
+ # Add the summary of 'raise' lines
+
+ # flush counter of code (disabled, we report static summary)
+ #msg = '1 line' if del_lines==1 else ('%s lines' % del_lines)
+ #exc = exc_tpl % msg
+ exc = exc_tpl
+
+ # Use the last value of 'space'
+ line = '%s%s' % (spaces[0],exc)
+ out.append(line)
+ del_lines = 0
+ spaces[:] = []
+
+ return del_lines
+
+ # used to report actual # of lines removed - disabled
+ #exc_tpl = "raise NotImplementedError('%s missing')\n"
+ exc_tpl = "raise NotImplementedError('insert missing code here')\n"
+
+ # states for state machine and other initialization
+ normal,delete = 0,1
+ state_cur = normal
+ del_lines = 0 # counter, in case we want to report # of deletions
+ spaces = []
+ normal_lines = []
+ out = []
+
+ # To remove multiple consecutive lines of input marked for deletion, we
+ # need a small state machine.
+ for line in src:
+ # Check for lines to delete and with hints
+ mdel = DEL_RE.match(line)
+ mhint = HINT_RE.match(line)
+
+ if mhint:
+ state_new = normal
+ hint = mhint.group('space')+'# ' + mhint.group('hint') +'\n'
+ normal_lines.append(hint)
+ else:
+ if mdel is None:
+ state_new = normal
+ normal_lines.append(line)
+ else:
+ state_new = delete
+ del_lines += 1
+ spaces.append(mdel.group(2))
+
+ # Flush output only when there's a change of state
+ if state_new != state_cur:
+ del_lines = flush_buffers(normal_lines,del_lines)
+
+ # Update state machine
+ state_cur = state_new
+
+ # Final buffer flush is unconditional
+ flush_buffers(normal_lines)
+
+ return ''.join(out)
+
+
+def transform_file(fpath,fname_skel,fname_soln):
+ """Run the cleanup routines for a given input, creating skel and soln.
+ """
+
+ # get the mode of the input so that we can create the output files with the
+ # same mode
+ fmode = os.stat(fpath).st_mode
+
+ with open(fpath) as infile:
+ # Generate the skeleton
+ skel = src2skel(infile)
+ with open(fname_skel,'w') as fskel:
+ fskel.write(skel)
+ os.chmod(fname_skel,fmode)
+
+ # Reset the input file pointer and generate the solution
+ infile.seek(0)
+ soln = src2soln(infile)
+ with open(fname_soln,'w') as fsoln:
+ fsoln.write(soln)
+ os.chmod(fname_soln,fmode)
+
+#-----------------------------------------------------------------------------
+# Main execution routines
+
+def copyforce(src,dest):
+ """Forceful file link/copy that overwrites destination files."""
+ try:
+ copy = os.link
+ except AttributeError:
+ copy = shutil.copy
+ if os.path.isfile(dest):
+ os.remove(dest)
+ copy(src,dest)
+
+
+def mvforce(src,dest):
+ """Forceful file copy that overwrites destination files."""
+ if os.path.isfile(dest):
+ os.remove(dest)
+ shutil.move(src,dest)
+
+
+def main(argv=None):
+ """Main entry point as a command line script for normal execution"""
+
+ if argv is None:
+ argv = sys.argv
+
+ # If there are subdirs called skel and soln, we populate them by moving the
+ # generated files there, otherwise they're left in the current dir.
+ skel_dir = 'skel'
+ soln_dir = 'soln'
+ has_skel_dir = os.path.isdir(skel_dir)
+ has_soln_dir = os.path.isdir(soln_dir)
+
+ # First, check that all files are present and abort immediately if any of
+ # them isn't there.
+ for fpath in argv[1:]:
+ if not os.path.isfile(fpath):
+ raise OSError("file %r not found" % fpath)
+
+ # If all files are there, then go ahead and process them unconditionally
+ for fpath in argv[1:]:
+ basename, ext = os.path.splitext(fpath)
+ fname_skel = basename + '_skel' + ext
+ fname_soln = basename + '_soln' + ext
+ transform_file(fpath,fname_skel,fname_soln)
+ # Move files over to final dirs if present
+ if has_skel_dir:
+ mvforce(fname_skel,os.path.join(skel_dir,fname_skel))
+ if has_soln_dir:
+ mvforce(fname_soln,os.path.join(soln_dir,fname_soln))
+
+#-----------------------------------------------------------------------------
+# Tests
+
+def str_match(s1,s2):
+ """Check that two strings are equal ignoring trailing whitespace."""
+ #print '***S1\n',s1,'\n***S2\n',s2 # dbg
+ nose.tools.assert_equal(s1.rstrip(),s2.rstrip())
+
+
+def test_simple():
+ src = """
+ first line
+ del line #@
+ second line
+ """
+ srclines = src.splitlines(True)
+
+ clean = """
+ first line
+ raise NotImplementedError('insert missing code here')
+ second line
+ """
+
+ cleaned = src2skel(srclines)
+ yield str_match,cleaned,clean
+
+ clean = """
+ first line
+ del line
+ second line
+ """
+ cleaned = src2soln(src.splitlines(True))
+ yield str_match,cleaned,clean
+
+
+def test_multi():
+ src = """
+ first line
+ #@ Hint: remember that
+ #@ idea we discussed before...
+ del line #@
+ del line2 #@
+ del line3 #@
+ second line:
+ del line4 #@
+ del line5 #@
+ third line
+
+ some indented code: #@
+ with more... #@
+ """
+ srclines = src.splitlines(True)
+
+ clean = """
+ first line
+ # Hint: remember that
+ # idea we discussed before...
+ raise NotImplementedError('insert missing code here')
+ second line:
+ raise NotImplementedError('insert missing code here')
+ third line
+
+ raise NotImplementedError('insert missing code here')
+ """
+ cleaned = src2skel(srclines)
+ yield str_match,cleaned,clean
+
+ clean = """
+ first line
+ del line
+ del line2
+ del line3
+ second line:
+ del line4
+ del line5
+ third line
+
+ some indented code:
+ with more...
+ """
+ cleaned = src2soln(srclines)
+ yield str_match,cleaned,clean
+
+
[EMAIL PROTECTED]
+def test():
+ """Simple self-contained test runner."""
+ nose.runmodule(__name__,exit=False,
+ argv=['--doctests',
+ #'-s',
+ #'--pdb-failures',
+ ])
+
+#-----------------------------------------------------------------------------
+# Execution from the command line.
+
+if __name__ == "__main__":
+ if '--test' in sys.argv:
+ test()
+ else:
+ main(sys.argv)
Property changes on: trunk/py4science/examples/mkskel.py
___________________________________________________________________
Added: svn:executable
+ *
Modified: trunk/py4science/examples/montecarlo_pi.py
===================================================================
--- trunk/py4science/examples/montecarlo_pi.py 2008-10-19 07:48:51 UTC (rev
6269)
+++ trunk/py4science/examples/montecarlo_pi.py 2008-10-19 07:53:47 UTC (rev
6270)
@@ -15,7 +15,7 @@
import math
import random
-import numpy as N
+import numpy as np
from scipy import weave
from scipy.weave import inline,converters
@@ -110,7 +110,7 @@
print 'pi - weave :',v2()
# make a simple 10x10 array
- a = N.arange(100)
+ a = np.arange(100)
a.shape = 10,10
# Print it using our printer
Modified: trunk/py4science/examples/qsort.py
===================================================================
--- trunk/py4science/examples/qsort.py 2008-10-19 07:48:51 UTC (rev 6269)
+++ trunk/py4science/examples/qsort.py 2008-10-19 07:53:47 UTC (rev 6270)
@@ -1,29 +1,60 @@
-"""
-Simple quicksort implementation.
+"""Simple quicksort implementation.
-See http://en.wikipedia.org/wiki/Quicksort for algorithm, pseudocode
-and C implementation for comparison
+From http://en.wikipedia.org/wiki/Quicksort we have this pseudocode (see also
+the C implementation for comparison).
+
+Note: what follows is NOT python code, it's meant to only illustrate the
+algorithm for you. Below you'll need to actually implement it in real Python.
+You may be surprised at how close a working Python implementation can be to
+this pseudocode.
+
+
+function quicksort(array)
+ var list less, greater
+ if length(array) <= 1
+ return array
+ select and remove a pivot value pivot from array
+ for each x in array
+ if x <= pivot then append x to less
+ else append x to greater
+ return concatenate(quicksort(less), pivot, quicksort(greater))
"""
def qsort(lst):
- """Return a sorted copy of the input list."""
+ """Return a sorted copy of the input list.
- if len(lst) <= 1:
- return lst
+ Input:
- # Select pivot and apply recursively
- pivot, rest = lst[0],lst[1:]
- less_than = [ lt for lt in rest if lt < pivot ]
- greater_equal = [ ge for ge in rest if ge >= pivot ]
-
- return qsort(less_than) + [pivot] + qsort(greater_equal)
+ lst : a list of elements which can be compared.
+ Examples:
+ >>> qsort([])
+ []
+
+ >>> qsort([3,2,5])
+ [2, 3, 5]
+ """
+
+ #@ Hint: remember that all recursive functions need an exit condition
+ if len(lst) <= 1: #@
+ return lst #@
+
+ #@ Select pivot and apply recursively
+ pivot, rest = lst[0],lst[1:] #@
+ less_than = [ lt for lt in rest if lt < pivot ] #@
+ greater_equal = [ ge for ge in rest if ge >= pivot ] #@
+
+ #@ Upon return, make sure to properly concatenate the output lists
+ return qsort(less_than) + [pivot] + qsort(greater_equal) #@
+
+
#-----------------------------------------------------------------------------
# Tests
#-----------------------------------------------------------------------------
import random
+import nose
import nose, nose.tools as nt
def test_sorted():
Modified: trunk/py4science/examples/quad_newton.py
===================================================================
--- trunk/py4science/examples/quad_newton.py 2008-10-19 07:48:51 UTC (rev
6269)
+++ trunk/py4science/examples/quad_newton.py 2008-10-19 07:53:47 UTC (rev
6270)
@@ -4,35 +4,35 @@
from math import sin
-import scipy, scipy.integrate, scipy.optimize
+from scipy.integrate import quad
+from scipy.optimize import newton
-quad = scipy.integrate.quad
-newton = scipy.optimize.newton
-
# test input function
def f(t):
- return t*(sin(t))**2
+ #@ f(t): t * sin^2(t)
+ return t*(sin(t))**2 #@
-# exact \int_0^t f(s) ds - u
def g(t):
+ "Exact form for g by integrating f(t)"
u = 0.25
return .25*(t**2-t*sin(2*t)+(sin(t))**2)-u
-# now let's construct g(t) via numerical integration
def gn(t):
+ "g(t) obtained by numerical integration"
u = 0.25
- return quad(f,0.0,t)[0] - u
+ #@ Hint: use quad, see its return value carefully.
+ return quad(f,0.0,t)[0] - u #@
# main
tguess = 10.0
print '"Exact" solution (knowing the analytical form of the integral)'
-t0 = newton(g,tguess,f)
+t0 = newton(g,tguess,f) #@
print "t0, g(t0) =",t0,g(t0)
print
print "Solution using the numerical integration technique"
-t1 = newton(gn,tguess,f)
+t1 = newton(gn,tguess,f) #@
print "t1, g(t1) =",t1,g(t1)
print
Modified: trunk/py4science/examples/recarray_demo.py
===================================================================
--- trunk/py4science/examples/recarray_demo.py 2008-10-19 07:48:51 UTC (rev
6269)
+++ trunk/py4science/examples/recarray_demo.py 2008-10-19 07:53:47 UTC (rev
6270)
@@ -1,67 +1,53 @@
"""
parse and load ge.csv into a record array
"""
-import time, datetime, csv
-import dateutil.parser
-import numpy
-if 0:
- # create a csv reader to parse the file
- fh = file('data/ge.csv')
- reader = csv.reader(fh)
- header = reader.next()
+import datetime
- # iterate over the remaining rows and convert the data to date
- # objects, ints, or floats as approriate
- rows = []
- for date, open, high, low, close, volume, adjclose in reader:
- date = dateutil.parser.parse(date).date()
- volume = int(volume)
- open, high, low, close, adjclose = map(float, (
- open, high, low, close, adjclose))
- rows.append((date, open, high, low, close, volume, adjclose))
+import numpy as np
- fh.close()
+import matplotlib
+from matplotlib import pyplot as plt
- # this is how you can use the function matplotlib.mlab.load to do the same
- #rows = load('data/ge.csv', skiprows=1, converters={0:datestr2num},
delimiter=',')
+# Disable latex support just in case, so we can rotate text
+matplotlib.rcParams['text.usetex'] = False
- r = numpy.rec.fromrecords(rows,
names='date,open,high,low,close,volume,adjclose')
+# Load a data file into a record array using the matplotlib csv2rec utility
+r = matplotlib.mlab.csv2rec('data/ge.csv')
+r.sort() #sort by date, the first column
# compute the average approximate dollars traded over the last 10 days
# hint: closing price * volume trades approx equals dollars trades
recent = r[-10:]
-dollars = numpy.mean(recent.close * recent.volume)
-print '$%1.2fM'%(dollars/1e6)
+dollars = np.mean(recent.close * recent.volume)
+print 'Total traded over last 10 days: $%1.2fM'%(dollars/1e6)
-# plot the adjusted closing price vs time since 2003 - hint, you must
-# use date2num to convert the date to a float for mpl. Make two axes,
-# one for price and one for volume. Use a bar chart for volume
-import matplotlib
-matplotlib.rcParams['usetex'] = False
-from matplotlib.dates import date2num
-import pylab
-mask = r.date>datetime.date(2003,1,1)
-date = date2num(r.date[mask])
-price = r.adjclose[mask]
+# plot the adjusted closing price vs time since 2003. Make two axes, one for
+# price and one for volume. Use a bar chart for volume
+
+# Record arrays are like 'mini-spreadsheets'
+mask = r.date > datetime.date(2003,1,1)
+date = r.date[mask]
+price = r.adj_close[mask]
volume = r.volume[mask]
-fig = pylab.figure()
+# Plotting
+fig = plt.figure()
fig.subplots_adjust(hspace=0)
ax1 = fig.add_subplot(211)
-ax1.plot(date, price, '-');
-ax1.xaxis_date()
+ax1.plot(date, price, '-')
+
ax1.grid()
for label in ax1.get_xticklabels():
label.set_visible(False)
ax2 = fig.add_subplot(212, sharex=ax1)
ax2.bar(date, volume);
-ax2.xaxis_date()
+
ax2.grid()
for label in ax2.get_xticklabels():
label.set_rotation(30)
label.set_horizontalalignment('right')
-
-pylab.show()
+fig.autofmt_xdate()
+plt.show()
Added: trunk/py4science/examples/skel/fortran_wrap/test_example.py
===================================================================
--- trunk/py4science/examples/skel/fortran_wrap/test_example.py
(rev 0)
+++ trunk/py4science/examples/skel/fortran_wrap/test_example.py 2008-10-19
07:53:47 UTC (rev 6270)
@@ -0,0 +1,9 @@
+import example
+import numpy
+
+x = numpy.arange(10.)
+
+yf = example.cumsum(x)
+yn = numpy.cumsum(x)
+
+numpy.testing.assert_array_almost_equal(yf, yn)
Added: trunk/py4science/examples/skel/fortran_wrap/test_fib.py
===================================================================
--- trunk/py4science/examples/skel/fortran_wrap/test_fib.py
(rev 0)
+++ trunk/py4science/examples/skel/fortran_wrap/test_fib.py 2008-10-19
07:53:47 UTC (rev 6270)
@@ -0,0 +1,13 @@
+import numpy as N
+
+import example
+
+n = 10
+fn = example.fib(n)
+
+print 'Fibonacci numbers:'
+print fn
+
+# Check validity
+assert N.alltrue(fn[2:]== fn[1:-1]+fn[:-2]), "Fibonacci mismatch"
+
Modified: trunk/py4science/examples/trapezoid.py
===================================================================
--- trunk/py4science/examples/trapezoid.py 2008-10-19 07:48:51 UTC (rev
6269)
+++ trunk/py4science/examples/trapezoid.py 2008-10-19 07:53:47 UTC (rev
6270)
@@ -7,7 +7,8 @@
"""Simple trapezoid integrator for sequence-based innput.
Inputs:
- - x,y: arrays of the same length.
+ - x,y: arrays of the same length (and more than one element). If the two
+ inputs have different lengths, a ValueError exception is raised.
Output:
- The result of applying the trapezoid rule to the input, assuming that
@@ -15,14 +16,20 @@
Minimally modified from matplotlib.mlab."""
- # Sanity checks
- if len(x)!=len(y):
- raise ValueError('x and y must have the same length')
- if len(x)<2:
- raise ValueError('x and y must have > 1 element')
+ # Sanity checks.
+ #@
+ #@ Hint: if the two inputs have mismatched lengths or less than 2
+ #@ elements, we raise ValueError with an explanatory message.
+ if len(x) != len(y): #@
+ raise ValueError('x and y must have the same length') #@
+ if len(x) < 2: #@
+ raise ValueError('x and y must have > 1 element') #@
# Efficient application of trapezoid rule via numpy
- return 0.5*((x[1:]-x[:-1])*(y[1:]+y[:-1])).sum()
+ #@
+ #@ Hint: think of using numpy slicing to compute the moving difference in
+ #@ the basic trapezoid formula.
+ return 0.5*((x[1:]-x[:-1])*(y[1:]+y[:-1])).sum() #@
def trapzf(f,a,b,npts=100):
"""Simple trapezoid-based integrator.
@@ -39,15 +46,25 @@
Output:
- The value of the trapezoid-rule approximation to the integral."""
- # Generate an equally spaced grid to sample the function at
- x = np.linspace(a,b,npts)
+ #@ Hint: you will need to apply the function f to easch element of the
+ #@ vector x. What are several ways to do this? Can you profile them to
see
+ #@ what differences in timings result for long vectors x?
+
+ # Generate an equally spaced grid to sample the function.
+ x = np.linspace(a,b,npts)#@
+
# For an equispaced grid, the x spacing can just be read off from the first
# two points and factored out of the summation.
- dx = x[1]-x[0]
+ dx = x[1]-x[0]#@
+
# Sample the input function at all values of x
- y = np.array(map(f,x))
+ #@
+ #@ Hint: you need to make an array out of the evaluations, and the python
+ #@ builtin 'map' function can come in handy.
+ y = np.array(map(f,x))#@
+
# Compute the trapezoid rule sum for the final result
- return 0.5*dx*(y[1:]+y[:-1]).sum()
+ return 0.5*dx*(y[1:]+y[:-1]).sum()#@
#-----------------------------------------------------------------------------
@@ -56,20 +73,25 @@
import nose, nose.tools as nt
import numpy.testing as nptest
+# A simple function for testing
def square(x): return x**2
def test_err():
+ """Test that mismatched inputs raise a ValueError exception."""
nt.assert_raises(ValueError,trapz,range(2),range(3))
def test_call():
+ "Test a direct call with equally spaced samples. "
x = np.linspace(0,1,100)
y = np.array(map(square,x))
nptest.assert_almost_equal(trapz(x,y),1./3,4)
def test_square():
+ "Test integrating the square() function."
nptest.assert_almost_equal(trapzf(square,0,1),1./3,4)
def test_square2():
+ "Another test integrating the square() function."
nptest.assert_almost_equal(trapzf(square,0,3,350),9.0,4)
This was sent by the SourceForge.net collaborative development platform, the
world's largest Open Source development site.
-------------------------------------------------------------------------
This SF.Net email is sponsored by the Moblin Your Move Developer's challenge
Build the coolest Linux based applications with Moblin SDK & win great prizes
Grand prize is a trip for two to an Open Source event anywhere in the world
http://moblin-contest.org/redirect.php?banner_id=100&url=/
_______________________________________________
Matplotlib-checkins mailing list
Matplotlib-checkins@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/matplotlib-checkins
