added discussion of algorithm

docs/demos/algorithm.ipynb

@@ -0,0 +1,94 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# How Does Powerbox Work?"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "It may be useful to understand the workings of powerbox to some extent -- either to diagnose performance issues or to understand its behaviour in certain contexts.\n",
+    "\n",
+    "The basic algorithm (for a Gaussian field) is the following:\n",
+    "\n",
+    "1. Given a box length $L$ (parameter ``boxlength``) and number of cells along a side, $N$ (parameter ``N``), as well as Fourier convention parameters $(a,b)$, determine wavenumbers along a side of the box: $k = 2\\pi j/(bL)$, for $j\\in (-N/2,..., N/2)$.\n",
+    "2. From these wavenumbers along each side, determine the *magnitude* of the wavenumbers at every point of the $d$-dimensional box, $k_j= \\sqrt{\\sum_{i=1}^d k_{i,j}^2}$.\n",
+    "3. Create an array, $G_j$, which assigns a complex number to each grid point. The complex number will have magnitude drawn from a standard normal, and phase distributed evenly on $(0,2\\pi)$.\n",
+    "4. Determine $\\delta_{k,j} = G_j \\sqrt{P(k_j)}$.\n",
+    "5. Determine $\\delta_x = V \\mathcal{F}^{-1}(\\delta_k)$, with $V = \\prod_{i=1}^{d} L_i$."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "For a Log-Normal field, the steps are slightly more complex, and involve determining the power spectrum that *would be* required on a Gaussian field to yield the same power spectrum for a log-normal field. The details of this approach can be found in [Coles and Jones (1991)](http://adsabs.harvard.edu/abs/1991MNRAS.248....1C) or [Beutler et al. (2011)](https://academic.oup.com/mnras/article/416/4/3017/976636)."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "One characteristic of this algorithm is that it contains *no information* below the resolution scale $L/N$. Thus, a good rule-of-thumb is to choose $N$ large enough to ensure that the smallest scale of interest is covered by a factor of 1.5, i.e., if the smallest length-scale of interest is $s$, then use $N = 1.5 L/s$.\n",
+    "\n",
+    "The range of $k$ used with this choice of $N$ also depends on the Fourier Convention used. For the default convention of $b=1$, the smallest scales are equivalent to $k = \\pi N/L$."
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python [conda env:powerbox]",
+   "language": "python",
+   "name": "conda-env-powerbox-py"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.6.5"
+  },
+  "latex_envs": {
+   "LaTeX_envs_menu_present": true,
+   "autoclose": false,
+   "autocomplete": true,
+   "bibliofile": "biblio.bib",
+   "cite_by": "apalike",
+   "current_citInitial": 1,
+   "eqLabelWithNumbers": true,
+   "eqNumInitial": 1,
+   "hotkeys": {
+    "equation": "Ctrl-E",
+    "itemize": "Ctrl-I"
+   },
+   "labels_anchors": false,
+   "latex_user_defs": false,
+   "report_style_numbering": false,
+   "user_envs_cfg": false
+  },
+  "toc": {
+   "base_numbering": 1,
+   "nav_menu": {},
+   "number_sections": true,
+   "sideBar": true,
+   "skip_h1_title": false,
+   "title_cell": "Table of Contents",
+   "title_sidebar": "Contents",
+   "toc_cell": false,
+   "toc_position": {},
+   "toc_section_display": true,
+   "toc_window_display": false
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}

docs/examples.rst

@@ -8,5 +8,6 @@ Other examples can be found in the :doc:`API documentation <api>` for each objec
    :maxdepth: 2
 
    demos/getting_started
+   demos/algorithm
    demos/cosmological_fields
    demos/dft

powerbox/powerbox.py

@@ -12,14 +12,17 @@
 
 try:
     from multiprocessing import cpu_count
+
     THREADS = cpu_count()
 
     from pyfftw import empty_aligned as empty
+
     HAVE_FFTW = True
 except ImportError:
     empty = np.empty
     HAVE_FFTW = False
 
+
 # TODO: add hankel-transform version of LogNormal
 
 
@@ -42,10 +45,10 @@ def _make_hermitian(mag, pha):
     kspace : array
         A complex hermitian array with normally distributed amplitudes.
     """
-    revidx = [slice(None, None, -1)]*len(mag.shape)
-    mag = (mag + mag[revidx])/np.sqrt(2)
-    pha = (pha - pha[revidx])/2 + np.pi
-    return mag*(np.cos(pha) + 1j*np.sin(pha))
+    revidx = [slice(None, None, -1)] * len(mag.shape)
+    mag = (mag + mag[revidx]) / np.sqrt(2)
+    pha = (pha - pha[revidx]) / 2 + np.pi
+    return mag * (np.cos(pha) + 1j * np.sin(pha))
 
 
 class PowerBox(object):
@@ -139,7 +142,7 @@ def __init__(self, N, pk, dim=2, boxlength=1.0, ensure_physical=False, a=1., b=1
         self.V = self.boxlength ** self.dim
 
         if self.vol_normalised_power:
-            self.pk = lambda k: pk(k)/self.V
+            self.pk = lambda k: pk(k) / self.V
         else:
             self.pk = pk
 
@@ -148,15 +151,15 @@ def __init__(self, N, pk, dim=2, boxlength=1.0, ensure_physical=False, a=1., b=1
 
         self.seed = seed
 
-        if N%2 == 0:
+        if N % 2 == 0:
             self._even = True
         else:
             self._even = False
 
         self.n = N + 1 if self._even else N
 
         # Get the grid-size for the final real-space box.
-        self.dx = float(boxlength)/N
+        self.dx = float(boxlength) / N
 
     def k(self):
         "The entire grid of wavenumber magitudes"
@@ -175,20 +178,20 @@ def r(self):
     @property
     def x(self):
         "The co-ordinates of the grid along a side"
-        return np.arange(-self.boxlength/2, self.boxlength/2, self.dx)[:self.N]
+        return np.arange(-self.boxlength / 2, self.boxlength / 2, self.dx)[:self.N]
 
     def gauss_hermitian(self):
         "A random array which has Gaussian magnitudes and Hermitian symmetry"
         if self.seed:
             np.random.seed(self.seed)
 
-        mag = np.random.normal(0, 1, size=[self.n]*self.dim)
-        pha = 2*np.pi*np.random.uniform(size=[self.n]*self.dim)
+        mag = np.random.normal(0, 1, size=[self.n] * self.dim)
+        pha = 2 * np.pi * np.random.uniform(size=[self.n] * self.dim)
 
         dk = _make_hermitian(mag, pha)
 
         if self._even:
-            cutidx = [slice(None, -1)]*self.dim
+            cutidx = [slice(None, -1)] * self.dim
             dk = dk[cutidx]
 
         return dk
@@ -205,20 +208,20 @@ def delta_k(self):
         "A realisation of the delta_k, i.e. the gaussianised square root of the power spectrum (i.e. the Fourier co-efficients)"
         p = self.power_array()
 
-        if np.any(p<0):
+        if np.any(p < 0):
             raise ValueError("The power spectrum function has returned negative values.")
 
         gh = self.gauss_hermitian()
-        gh[...] = np.sqrt(p)*gh
+        gh[...] = np.sqrt(p) * gh
         return gh
 
     def delta_x(self):
         "The realised field in real-space from the input power spectrum"
         # Here we multiply by V because the (inverse) fourier-transform of the (dimensionless) power has
         # units of 1/V and we require a unitless quantity for delta_x.
-        dk = empty((self.N,)*self.dim,dtype='complex128')
+        dk = empty((self.N,) * self.dim, dtype='complex128')
         dk[...] = self.delta_k()
-        dk[...] = self.V*dft.ifft(dk, L=self.boxlength, a=self.fourier_a, b=self.fourier_b)[0]
+        dk[...] = self.V * dft.ifft(dk, L=self.boxlength, a=self.fourier_a, b=self.fourier_b)[0]
         dk = np.real(dk)
 
         if self.ensure_physical:
@@ -252,22 +255,22 @@ def create_discrete_sample(self, nbar, randomise_in_cell=True, min_at_zero=False
             a single tracer's co-ordinates.
         """
         dx = self.delta_x()
-        dx = (dx + 1)*self.dx ** self.dim*nbar
+        dx = (dx + 1) * self.dx ** self.dim * nbar
         n = dx
         self.n_per_cell = np.random.poisson(n)
 
         # Get all source positions
-        args = [self.x]*self.dim
+        args = [self.x] * self.dim
         X = np.meshgrid(*args)
 
         tracer_positions = np.array([x.flatten() for x in X]).T
         tracer_positions = tracer_positions.repeat(self.n_per_cell.flatten(), axis=0)
 
         if randomise_in_cell:
-            tracer_positions += np.random.uniform(size=(np.sum(self.n_per_cell), self.dim))*self.dx
+            tracer_positions += np.random.uniform(size=(np.sum(self.n_per_cell), self.dim)) * self.dx
 
         if min_at_zero:
-            tracer_positions += self.boxlength/2.0
+            tracer_positions += self.boxlength / 2.0
 
         if store_pos:
             self.tracer_positions = tracer_positions
@@ -321,17 +324,17 @@ def __init__(self, *args, **kwargs):
 
     def correlation_array(self):
         "The correlation function from the input power, on the grid"
-        pa = empty((self.N,)*self.dim)
+        pa = empty((self.N,) * self.dim)
         pa[...] = self.power_array()
-        return self.V*np.real(dft.ifft(pa, L=self.boxlength, a=self.fourier_a, b=self.fourier_b)[0])
+        return self.V * np.real(dft.ifft(pa, L=self.boxlength, a=self.fourier_a, b=self.fourier_b)[0])
 
     def gaussian_correlation_array(self):
         "The correlation function required for a Gaussian field to produce the input power on a lognormal field"
         return np.log(1 + self.correlation_array())
 
     def gaussian_power_array(self):
         "The power spectrum required for a Gaussian field to produce the input power on a lognormal field"
-        gca = empty((self.N,)*self.dim)
+        gca = empty((self.N,) * self.dim)
         gca[...] = self.gaussian_correlation_array()
         gpa = np.abs(dft.fft(gca, L=self.boxlength, a=self.fourier_a, b=self.fourier_b))[0]
         gpa[self.k() == 0] = 0
@@ -344,15 +347,15 @@ def delta_k(self):
         """
         p = self.gaussian_power_array()
         gh = self.gauss_hermitian()
-        gh[...] = np.sqrt(p)*gh
+        gh[...] = np.sqrt(p) * gh
         return gh
 
     def delta_x(self):
         "The real-space over-density field, from the input power spectrum"
-        dk = empty((self.N,)*self.dim,dtype='complex128')
+        dk = empty((self.N,) * self.dim, dtype='complex128')
         dk[...] = self.delta_k()
-        dk[...] = np.sqrt(self.V)*dft.ifft(dk, L=self.boxlength, a=self.fourier_a, b=self.fourier_b)[0]
+        dk[...] = np.sqrt(self.V) * dft.ifft(dk, L=self.boxlength, a=self.fourier_a, b=self.fourier_b)[0]
         dk = np.real(dk)
 
         sg = np.var(dk)
-        return np.exp(dk - sg/2) - 1
+        return np.exp(dk - sg / 2) - 1