include optional seeding of random number generator in spherical harmonic noise function for reproducability and parallel noise production

markmbaum · markmbaum · commit 63bdeea938a2 · 2021-11-25T08:28:17.000-05:00
diff --git a/docs/_static/documentation_options.js b/docs/_static/documentation_options.js
@@ -1,6 +1,6 @@
 var DOCUMENTATION_OPTIONS = {
     URL_ROOT: document.getElementById("documentation_options").getAttribute('data-url_root'),
-    VERSION: '0.9.3',
+    VERSION: '0.9.4',
     LANGUAGE: 'None',
     COLLAPSE_INDEX: false,
     BUILDER: 'html',
diff --git a/docs/chebyshev.html b/docs/chebyshev.html
@@ -6,7 +6,7 @@
     <meta charset="utf-8" />
     <meta name="viewport" content="width=device-width, initial-scale=1.0" /><meta name="generator" content="Docutils 0.17.1: http://docutils.sourceforge.net/" />
 
-    <title>chebyshev &#8212; orthopoly 0.9.3 documentation</title>
+    <title>chebyshev &#8212; orthopoly 0.9.4 documentation</title>
     <link rel="stylesheet" type="text/css" href="_static/pygments.css" />
     <link rel="stylesheet" type="text/css" href="_static/haiku.css" />
     <script data-url_root="./" id="documentation_options" src="_static/documentation_options.js"></script>
@@ -20,7 +20,7 @@
     <link rel="prev" title="Welcome to orthopoly’s documentation!" href="index.html" /> 
   </head><body>
       <div class="header" role="banner"><h1 class="heading"><a href="index.html">
-          <span>orthopoly 0.9.3 documentation</span></a></h1>
+          <span>orthopoly 0.9.4 documentation</span></a></h1>
         <h2 class="heading"><span>chebyshev</span></h2>
       </div>
       <div class="topnav" role="navigation" aria-label="top navigation">
diff --git a/docs/gegenbauer.html b/docs/gegenbauer.html
@@ -6,7 +6,7 @@
     <meta charset="utf-8" />
     <meta name="viewport" content="width=device-width, initial-scale=1.0" /><meta name="generator" content="Docutils 0.17.1: http://docutils.sourceforge.net/" />
 
-    <title>gegenbauer &#8212; orthopoly 0.9.3 documentation</title>
+    <title>gegenbauer &#8212; orthopoly 0.9.4 documentation</title>
     <link rel="stylesheet" type="text/css" href="_static/pygments.css" />
     <link rel="stylesheet" type="text/css" href="_static/haiku.css" />
     <script data-url_root="./" id="documentation_options" src="_static/documentation_options.js"></script>
@@ -20,7 +20,7 @@
     <link rel="prev" title="chebyshev" href="chebyshev.html" /> 
   </head><body>
       <div class="header" role="banner"><h1 class="heading"><a href="index.html">
-          <span>orthopoly 0.9.3 documentation</span></a></h1>
+          <span>orthopoly 0.9.4 documentation</span></a></h1>
         <h2 class="heading"><span>gegenbauer</span></h2>
       </div>
       <div class="topnav" role="navigation" aria-label="top navigation">
diff --git a/docs/genindex.html b/docs/genindex.html
@@ -5,7 +5,7 @@
   <head>
     <meta charset="utf-8" />
     <meta name="viewport" content="width=device-width, initial-scale=1.0" />
-    <title>Index &#8212; orthopoly 0.9.3 documentation</title>
+    <title>Index &#8212; orthopoly 0.9.4 documentation</title>
     <link rel="stylesheet" type="text/css" href="_static/pygments.css" />
     <link rel="stylesheet" type="text/css" href="_static/haiku.css" />
     <script data-url_root="./" id="documentation_options" src="_static/documentation_options.js"></script>
@@ -16,7 +16,7 @@
     <link rel="search" title="Search" href="search.html" /> 
   </head><body>
       <div class="header" role="banner"><h1 class="heading"><a href="index.html">
-          <span>orthopoly 0.9.3 documentation</span></a></h1>
+          <span>orthopoly 0.9.4 documentation</span></a></h1>
         <h2 class="heading"><span>Index</span></h2>
       </div>
       <div class="topnav" role="navigation" aria-label="top navigation">
diff --git a/docs/index.html b/docs/index.html
@@ -6,7 +6,7 @@
     <meta charset="utf-8" />
     <meta name="viewport" content="width=device-width, initial-scale=1.0" /><meta name="generator" content="Docutils 0.17.1: http://docutils.sourceforge.net/" />
 
-    <title>Welcome to orthopoly’s documentation! &#8212; orthopoly 0.9.3 documentation</title>
+    <title>Welcome to orthopoly’s documentation! &#8212; orthopoly 0.9.4 documentation</title>
     <link rel="stylesheet" type="text/css" href="_static/pygments.css" />
     <link rel="stylesheet" type="text/css" href="_static/haiku.css" />
     <script data-url_root="./" id="documentation_options" src="_static/documentation_options.js"></script>
@@ -19,7 +19,7 @@
     <link rel="next" title="chebyshev" href="chebyshev.html" /> 
   </head><body>
       <div class="header" role="banner"><h1 class="heading"><a href="#">
-          <span>orthopoly 0.9.3 documentation</span></a></h1>
+          <span>orthopoly 0.9.4 documentation</span></a></h1>
         <h2 class="heading"><span>Welcome to orthopoly’s documentation!</span></h2>
       </div>
       <div class="topnav" role="navigation" aria-label="top navigation">
diff --git a/docs/legendre.html b/docs/legendre.html
@@ -6,7 +6,7 @@
     <meta charset="utf-8" />
     <meta name="viewport" content="width=device-width, initial-scale=1.0" /><meta name="generator" content="Docutils 0.17.1: http://docutils.sourceforge.net/" />
 
-    <title>legendre &#8212; orthopoly 0.9.3 documentation</title>
+    <title>legendre &#8212; orthopoly 0.9.4 documentation</title>
     <link rel="stylesheet" type="text/css" href="_static/pygments.css" />
     <link rel="stylesheet" type="text/css" href="_static/haiku.css" />
     <script data-url_root="./" id="documentation_options" src="_static/documentation_options.js"></script>
@@ -20,7 +20,7 @@
     <link rel="prev" title="gegenbauer" href="gegenbauer.html" /> 
   </head><body>
       <div class="header" role="banner"><h1 class="heading"><a href="index.html">
-          <span>orthopoly 0.9.3 documentation</span></a></h1>
+          <span>orthopoly 0.9.4 documentation</span></a></h1>
         <h2 class="heading"><span>legendre</span></h2>
       </div>
       <div class="topnav" role="navigation" aria-label="top navigation">
diff --git a/docs/py-modindex.html b/docs/py-modindex.html
@@ -5,7 +5,7 @@
   <head>
     <meta charset="utf-8" />
     <meta name="viewport" content="width=device-width, initial-scale=1.0" />
-    <title>Python Module Index &#8212; orthopoly 0.9.3 documentation</title>
+    <title>Python Module Index &#8212; orthopoly 0.9.4 documentation</title>
     <link rel="stylesheet" type="text/css" href="_static/pygments.css" />
     <link rel="stylesheet" type="text/css" href="_static/haiku.css" />
     <script data-url_root="./" id="documentation_options" src="_static/documentation_options.js"></script>
@@ -19,7 +19,7 @@
 
   </head><body>
       <div class="header" role="banner"><h1 class="heading"><a href="index.html">
-          <span>orthopoly 0.9.3 documentation</span></a></h1>
+          <span>orthopoly 0.9.4 documentation</span></a></h1>
         <h2 class="heading"><span>Python Module Index</span></h2>
       </div>
       <div class="topnav" role="navigation" aria-label="top navigation">
diff --git a/docs/search.html b/docs/search.html
@@ -5,7 +5,7 @@
   <head>
     <meta charset="utf-8" />
     <meta name="viewport" content="width=device-width, initial-scale=1.0" />
-    <title>Search &#8212; orthopoly 0.9.3 documentation</title>
+    <title>Search &#8212; orthopoly 0.9.4 documentation</title>
     <link rel="stylesheet" type="text/css" href="_static/pygments.css" />
     <link rel="stylesheet" type="text/css" href="_static/haiku.css" />
     
@@ -22,7 +22,7 @@
 
   </head><body>
       <div class="header" role="banner"><h1 class="heading"><a href="index.html">
-          <span>orthopoly 0.9.3 documentation</span></a></h1>
+          <span>orthopoly 0.9.4 documentation</span></a></h1>
         <h2 class="heading"><span>Search</span></h2>
       </div>
       <div class="topnav" role="navigation" aria-label="top navigation">
diff --git a/docs/searchindex.js b/docs/searchindex.js
diff --git a/docs/spherical_harmonic.html b/docs/spherical_harmonic.html
@@ -6,7 +6,7 @@
     <meta charset="utf-8" />
     <meta name="viewport" content="width=device-width, initial-scale=1.0" /><meta name="generator" content="Docutils 0.17.1: http://docutils.sourceforge.net/" />
 
-    <title>spherical_harmonic &#8212; orthopoly 0.9.3 documentation</title>
+    <title>spherical_harmonic &#8212; orthopoly 0.9.4 documentation</title>
     <link rel="stylesheet" type="text/css" href="_static/pygments.css" />
     <link rel="stylesheet" type="text/css" href="_static/haiku.css" />
     <script data-url_root="./" id="documentation_options" src="_static/documentation_options.js"></script>
@@ -20,7 +20,7 @@
     <link rel="prev" title="legendre" href="legendre.html" /> 
   </head><body>
       <div class="header" role="banner"><h1 class="heading"><a href="index.html">
-          <span>orthopoly 0.9.3 documentation</span></a></h1>
+          <span>orthopoly 0.9.4 documentation</span></a></h1>
         <h2 class="heading"><span>spherical_harmonic</span></h2>
       </div>
       <div class="topnav" role="navigation" aria-label="top navigation">
@@ -483,7 +483,7 @@ <h2 class="heading"><span>spherical_harmonic</span></h2>
 
 <dl class="py function">
 <dt class="sig sig-object py" id="orthopoly.spherical_harmonic.noise">
-<span class="sig-prename descclassname"><span class="pre">orthopoly.spherical_harmonic.</span></span><span class="sig-name descname"><span class="pre">noise</span></span><span class="sig-paren">(</span><em class="sig-param"><span class="n"><span class="pre">N</span></span></em>, <em class="sig-param"><span class="n"><span class="pre">p</span></span></em>, <em class="sig-param"><span class="n"><span class="pre">tol</span></span><span class="o"><span class="pre">=</span></span><span class="default_value"><span class="pre">1e-12</span></span></em><span class="sig-paren">)</span><a class="headerlink" href="#orthopoly.spherical_harmonic.noise" title="Permalink to this definition">¶</a></dt>
+<span class="sig-prename descclassname"><span class="pre">orthopoly.spherical_harmonic.</span></span><span class="sig-name descname"><span class="pre">noise</span></span><span class="sig-paren">(</span><em class="sig-param"><span class="n"><span class="pre">N</span></span></em>, <em class="sig-param"><span class="n"><span class="pre">p</span></span></em>, <em class="sig-param"><span class="n"><span class="pre">tol</span></span><span class="o"><span class="pre">=</span></span><span class="default_value"><span class="pre">1e-12</span></span></em>, <em class="sig-param"><span class="n"><span class="pre">seed</span></span><span class="o"><span class="pre">=</span></span><span class="default_value"><span class="pre">1</span></span></em><span class="sig-paren">)</span><a class="headerlink" href="#orthopoly.spherical_harmonic.noise" title="Permalink to this definition">¶</a></dt>
 <dd><p>Generates the coefficients for a random triangular spherical harmonic expansion with a specific relationship between the degree and the power spectral density (noise)</p>
 <dl class="field-list simple">
 <dt class="field-odd">Parameters</dt>
@@ -499,6 +499,7 @@ <h2 class="heading"><span>spherical_harmonic</span></h2>
 </ul>
 </p></li>
 <li><p><strong>tol</strong> (<em>float</em>) – bisection method relative tolerance when normalizing across a single degree’s coefficients for the total power</p></li>
+<li><p><strong>seed</strong> (<em>int</em>) – optional seed for random number generator</p></li>
 </ul>
 </dd>
 <dt class="field-even">Returns</dt>
diff --git a/docs/util.html b/docs/util.html
@@ -6,7 +6,7 @@
     <meta charset="utf-8" />
     <meta name="viewport" content="width=device-width, initial-scale=1.0" /><meta name="generator" content="Docutils 0.17.1: http://docutils.sourceforge.net/" />
 
-    <title>util &#8212; orthopoly 0.9.3 documentation</title>
+    <title>util &#8212; orthopoly 0.9.4 documentation</title>
     <link rel="stylesheet" type="text/css" href="_static/pygments.css" />
     <link rel="stylesheet" type="text/css" href="_static/haiku.css" />
     <script data-url_root="./" id="documentation_options" src="_static/documentation_options.js"></script>
@@ -19,7 +19,7 @@
     <link rel="prev" title="spherical_harmonic" href="spherical_harmonic.html" /> 
   </head><body>
       <div class="header" role="banner"><h1 class="heading"><a href="index.html">
-          <span>orthopoly 0.9.3 documentation</span></a></h1>
+          <span>orthopoly 0.9.4 documentation</span></a></h1>
         <h2 class="heading"><span>util</span></h2>
       </div>
       <div class="topnav" role="navigation" aria-label="top navigation">
diff --git a/orthopoly/spherical_harmonic.py b/orthopoly/spherical_harmonic.py
@@ -25,6 +25,7 @@
 import warnings
 import numpy as np
 from itertools import product
+from random import Random
 from scipy.spatial import SphericalVoronoi
 from numpy import *
 
@@ -684,7 +685,7 @@ def spectrum(a, yn, ym):
 
     return(ns, ps)
 
-def noise(N, p, tol=1e-12):
+def noise(N, p, tol=1e-12, seed=1):
     """Generates the coefficients for a random triangular spherical harmonic expansion with a specific relationship between the degree and the power spectral density (noise)
 
     :param int N: maximum degree in expansion
@@ -696,6 +697,7 @@ def noise(N, p, tol=1e-12):
         * 'blue':   p = 1
         * 'violet': p = 2
     :param float tol: bisection method relative tolerance when normalizing across a single degree's coefficients for the total power
+    :param int seed: optional seed for random number generator
 
     :return: :class:`Expansion` object with expansion coefficients for noise"""
 
@@ -709,6 +711,8 @@ def noise(N, p, tol=1e-12):
         assert isnum(p), 'p must be a number or a keyword'
     #get triangular degree and order arrays
     yn, ym = Tnm(N)
+    #intialize a random number generator
+    rng = Random(seed)
     #find expansion coefficients
     a = zeros((len(yn),))
     for n in set(yn[yn != 0]): #must ignore the zeroth degree (DC component)
@@ -718,7 +722,7 @@ def noise(N, p, tol=1e-12):
         sl = (yn == n)
         L = sum(sl)
         #random uniform distribution across the terms of degree n
-        r = 2.0*np.random.rand(L) - 1.0
+        r = array([2*rng.random() - 1 for _ in range(L)]) #2.0*np.random.rand(L) - 1.0)
         #find a factor f where sum(f*r**2) = d, using the bisection method
         flo = 0.0
         dlo = 0.0
diff --git a/setup.py b/setup.py
@@ -5,7 +5,7 @@
 
 setuptools.setup(
     name="orthopoly",
-    version="0.9.3",
+    version="0.9.4",
     author="Mark Baum <markmbaum@protonmail.com>",
     author_email="markmbaum@protonmail.com",
     description="Python functions for orthogonal polynomials and (real, 2D, orthonormal) spherical harmonics",
diff --git a/sphinx/conf.py b/sphinx/conf.py
@@ -21,7 +21,7 @@
 author = 'Mark M. Baum'
 
 # The full version, including alpha/beta/rc tags
-release = '0.9.3'
+release = '0.9.4'
 
 
 # -- General configuration ---------------------------------------------------
diff --git a/tests/spherical-harmonic/test_noise.py b/tests/spherical-harmonic/test_noise.py
@@ -85,7 +85,7 @@ def plot_noise(N, p='red', ntheta=151, nphi=301, ni=4, nj=5, cmap='RdBu'):
     for i in range(ni):
         for j in range(nj):
             ax = fig.add_subplot(ni, nj, c, projection='mollweide')
-            ex = noise(N, p)
+            ex = noise(N, p, seed=i*j)
             Y = M.dot(ex.a).reshape(ntheta, nphi)
             v = abs(Y).max()
             ax.pcolormesh(#P - pi, T - pi/2, Y,
@@ -141,5 +141,5 @@ def noise_means(N, npts, neval, p='red', cmap='RdBu'):
 plot_noise(10)
 plt.show()
 
-print('inspecting mean value of many instances of noise')
-noise_means(10, 201, 10000)
+#print('inspecting mean value of many instances of noise')
+#noise_means(10, 201, 10000)
