diff --git a/docs/examples/newparallel/rmt/rmt.ipynb b/docs/examples/newparallel/rmt/rmt.ipynb
new file mode 100644
index 00000000000..e5f8e60432f
--- /dev/null
+++ b/docs/examples/newparallel/rmt/rmt.ipynb
@@ -0,0 +1,1270 @@
+
+
+ rmt
+ 2
+
+
+
+
+
+
+
+
+
+
+ from rmtkernel import ensemble_diffs, normalize_diffs, GOE
+import numpy as np
+from IPython.parallel import Client
+ python
+ 1
+ 0
+
+
+
+
+
+
+
+
+
+ def wigner_dist(s):
+ """Returns (s, rho(s)) for the Wigner GOE distribution."""
+ return (np.pi*s/2.0) * np.exp(-np.pi*s**2/4.)
+ python
+ 2
+ 1
+
+
+
+ def generate_wigner_data():
+ s = np.linspace(0.0,4.0,400)
+ rhos = wigner_dist(s)
+ return s, rhos
+ python
+ 3
+ 1
+
+
+
+ s, rhos = generate_wigner_data()
+ python
+ 4
+ 0
+
+
+
+ plot(s, rhos)
+xlabel('Normalized level spacing s')
+ylabel('Probability $\rho(s)$')
+ python
+ 17
+ 0
+
+
+
+
+
+
+
+
+
+
+
+
+ def serial_diffs(num, N):
+ """Compute the nearest neighbor distribution for num NxX matrices."""
+ diffs = ensemble_diffs(num, N)
+ normalized_diffs = normalize_diffs(diffs)
+ return normalized_diffs
+ python
+ 6
+ 1
+
+
+
+ serial_nmats = 1000
+serial_matsize = 50
+ python
+ 7
+ 1
+
+
+
+ %timeit -r1 -n1 serial_diffs(serial_nmats, serial_matsize)
+ python
+ 8
+ 0
+
+
+
+
+
+ serial_diffs = serial_diffs(serial_nmats, serial_matsize)
+ python
+ 9
+ 0
+
+
+
+
+
+
+ hist_data = hist(serial_diffs, bins=30, normed=True)
+plot(s, rhos)
+xlabel('Normalized level spacing s')
+ylabel('Probability $P(s)$')
+ python
+ 10
+ 0
+
+
+
+
+
+
+
+
+
+
+
+
+ def parallel_diffs(rc, num, N):
+ nengines = len(rc.targets)
+ num_per_engine = num/nengines
+ print "Running with", num_per_engine, "per engine."
+ ar = rc.apply_async(ensemble_diffs, num_per_engine, N)
+ diffs = np.array(ar.get()).flatten()
+ normalized_diffs = normalize_diffs(diffs)
+ return normalized_diffs
+ python
+ 11
+ 1
+
+
+
+ client = Client()
+view = client[:]
+view.run('rmtkernel.py')
+view.block = False
+ python
+ 12
+ 1
+
+
+
+ parallel_nmats = 40*serial_nmats
+parallel_matsize = 50
+ python
+ 13
+ 1
+
+
+
+ %timeit -r1 -n1 parallel_diffs(view, parallel_nmats, parallel_matsize)
+ python
+ 14
+ 0
+
+
+
+
+
+
+ pdiffs = parallel_diffs(view, parallel_nmats, parallel_matsize)
+ python
+ 15
+ 0
+
+
+
+
+
+
+
+
+ hist_data = hist(pdiffs, bins=30, normed=True)
+plot(s, rhos)
+xlabel('Normalized level spacing s')
+ylabel('Probability $P(s)$')
+ python
+ 16
+ 0
+
+
+
+
+
+
+
+
+
diff --git a/docs/examples/newparallel/rmt/rmt.py b/docs/examples/newparallel/rmt/rmt.py
index d03142e95f6..a3292aa4d80 100644
--- a/docs/examples/newparallel/rmt/rmt.py
+++ b/docs/examples/newparallel/rmt/rmt.py
@@ -1,56 +1,146 @@
-#-------------------------------------------------------------------------------
-# Driver code that the client runs.
-#-------------------------------------------------------------------------------
-# To run this code start a controller and engines using:
-# ipcluster -n 2
-# Then run the scripts by doing irunner rmt.ipy or by starting ipython and
-# doing run rmt.ipy.
-
-from rmtkernel import *
-import numpy
+# 2
+
+#
+
+# # Eigenvalue distribution of Gaussian orthogonal random matrices
+
+#
+
+# The eigenvalues of random matrices obey certain statistical laws. Here we construct random matrices
+# from the Gaussian Orthogonal Ensemble (GOE), find their eigenvalues and then investigate the nearest
+# neighbor eigenvalue distribution $\rho(s)$.
+
+#
+
+from rmtkernel import ensemble_diffs, normalize_diffs, GOE
+import numpy as np
from IPython.parallel import Client
+#
-def wignerDistribution(s):
+# ## Wigner's nearest neighbor eigenvalue distribution
+
+#
+
+# The Wigner distribution gives the theoretical result for the nearest neighbor eigenvalue distribution
+# for the GOE:
+#
+# $$\rho(s) = \frac{\pi s}{2} \exp(-\pi s^2/4)$$
+
+#
+
+def wigner_dist(s):
"""Returns (s, rho(s)) for the Wigner GOE distribution."""
- return (numpy.pi*s/2.0) * numpy.exp(-numpy.pi*s**2/4.)
+ return (np.pi*s/2.0) * np.exp(-np.pi*s**2/4.)
+#
-def generateWignerData():
- s = numpy.linspace(0.0,4.0,400)
- rhos = wignerDistribution(s)
+def generate_wigner_data():
+ s = np.linspace(0.0,4.0,400)
+ rhos = wigner_dist(s)
return s, rhos
-
-def serialDiffs(num, N):
- diffs = ensembleDiffs(num, N)
- normalizedDiffs = normalizeDiffs(diffs)
- return normalizedDiffs
+#
+
+s, rhos = generate_wigner_data()
+
+#
+
+plot(s, rhos)
+xlabel('Normalized level spacing s')
+ylabel('Probability $\rho(s)$')
+
+#
+
+# ## Serial calculation of nearest neighbor eigenvalue distribution
+
+#
+
+# In this section we numerically construct and diagonalize a large number of GOE random matrices
+# and compute the nerest neighbor eigenvalue distribution. This comptation is done on a single core.
+
+#
+
+def serial_diffs(num, N):
+ """Compute the nearest neighbor distribution for num NxX matrices."""
+ diffs = ensemble_diffs(num, N)
+ normalized_diffs = normalize_diffs(diffs)
+ return normalized_diffs
+
+#
+
+serial_nmats = 1000
+serial_matsize = 50
+
+#
+
+%timeit -r1 -n1 serial_diffs(serial_nmats, serial_matsize)
+#
-def parallelDiffs(rc, num, N):
+serial_diffs = serial_diffs(serial_nmats, serial_matsize)
+
+#
+
+# The numerical computation agrees with the predictions of Wigner, but it would be nice to get more
+# statistics. For that we will do a parallel computation.
+
+#
+
+hist_data = hist(serial_diffs, bins=30, normed=True)
+plot(s, rhos)
+xlabel('Normalized level spacing s')
+ylabel('Probability $P(s)$')
+
+#
+
+# ## Parallel calculation of nearest neighbor eigenvalue distribution
+
+#
+
+# Here we perform a parallel computation, where each process constructs and diagonalizes a subset of
+# the overall set of random matrices.
+
+#
+
+def parallel_diffs(rc, num, N):
nengines = len(rc.targets)
num_per_engine = num/nengines
print "Running with", num_per_engine, "per engine."
- ar = rc.apply_async(ensembleDiffs, num_per_engine, N)
- return numpy.array(ar.get()).flatten()
-
-
-# Main code
-if __name__ == '__main__':
- rc = Client()
- view = rc[:]
- print "Distributing code to engines..."
- view.run('rmtkernel.py')
- view.block = False
-
- # Simulation parameters
- nmats = 100
- matsize = 30
- # tic = time.time()
- %timeit -r1 -n1 serialDiffs(nmats,matsize)
- %timeit -r1 -n1 parallelDiffs(view, nmats, matsize)
-
- # Uncomment these to plot the histogram
- # import pylab
- # pylab.hist(parallelDiffs(rc,matsize,matsize))
+ ar = rc.apply_async(ensemble_diffs, num_per_engine, N)
+ diffs = np.array(ar.get()).flatten()
+ normalized_diffs = normalize_diffs(diffs)
+ return normalized_diffs
+
+#
+
+client = Client()
+view = client[:]
+view.run('rmtkernel.py')
+view.block = False
+
+#
+
+parallel_nmats = 40*serial_nmats
+parallel_matsize = 50
+
+#
+
+%timeit -r1 -n1 parallel_diffs(view, parallel_nmats, parallel_matsize)
+
+#
+
+pdiffs = parallel_diffs(view, parallel_nmats, parallel_matsize)
+
+#
+
+# Again, the agreement with the Wigner distribution is excellent, but now we have better
+# statistics.
+
+#
+
+hist_data = hist(pdiffs, bins=30, normed=True)
+plot(s, rhos)
+xlabel('Normalized level spacing s')
+ylabel('Probability $P(s)$')
+
diff --git a/docs/examples/newparallel/rmt/rmtkernel.py b/docs/examples/newparallel/rmt/rmtkernel.py
index 2c913e20d78..a44b2b57bbf 100644
--- a/docs/examples/newparallel/rmt/rmtkernel.py
+++ b/docs/examples/newparallel/rmt/rmtkernel.py
@@ -2,43 +2,41 @@
# Core routines for computing properties of symmetric random matrices.
#-------------------------------------------------------------------------------
-import numpy
-ra = numpy.random
-la = numpy.linalg
+import numpy as np
+ra = np.random
+la = np.linalg
def GOE(N):
"""Creates an NxN element of the Gaussian Orthogonal Ensemble"""
m = ra.standard_normal((N,N))
m += m.T
- return m
+ return m/2
-def centerEigenvalueDiff(mat):
+def center_eigenvalue_diff(mat):
"""Compute the eigvals of mat and then find the center eigval difference."""
N = len(mat)
- evals = numpy.sort(la.eigvals(mat))
- diff = evals[N/2] - evals[N/2-1]
- return diff.real
+ evals = np.sort(la.eigvals(mat))
+ diff = np.abs(evals[N/2] - evals[N/2-1])
+ return diff
-def ensembleDiffs(num, N):
- """Return an array of num eigenvalue differences for the NxN GOE
- ensemble."""
- diffs = numpy.empty(num)
+def ensemble_diffs(num, N):
+ """Return num eigenvalue diffs for the NxN GOE ensemble."""
+ diffs = np.empty(num)
for i in xrange(num):
mat = GOE(N)
- diffs[i] = centerEigenvalueDiff(mat)
+ diffs[i] = center_eigenvalue_diff(mat)
return diffs
-def normalizeDiffs(diffs):
+def normalize_diffs(diffs):
"""Normalize an array of eigenvalue diffs."""
return diffs/diffs.mean()
-def normalizedEnsembleDiffs(num, N):
- """Return an array of num *normalized eigenvalue differences for the NxN
- GOE ensemble."""
- diffs = ensembleDiffs(num, N)
- return normalizeDiffs(diffs)
+def normalized_ensemble_diffs(num, N):
+ """Return num *normalized* eigenvalue diffs for the NxN GOE ensemble."""
+ diffs = ensemble_diffs(num, N)
+ return normalize_diffs(diffs)