DOCS: more of ariana's suggestions

doc/features.rst

@@ -3,7 +3,7 @@
 Fortran Simulation Features
 ###########################
 
-.. Why use wlcsim? Does it do what I need?
+Why use wlcsim? Does it do what I need?
 
 Semiflexable polymers
 =====================

wlcsim/bd/__init__.py

@@ -4,7 +4,7 @@
 A Brownian dynamics simulation requires specification of three things:
 
     1. **The system state**. For us, this is :math:`\vec{r}_j(t_i) \in
-       \mathbb{r}^{N\times 3}`, the positions of the :math:`N` discrete beads
+       \mathbb{R}^{N\times 3}`, the positions of the :math:`N` discrete beads
        representing the polymer at each time :math:`t_i`.
     2. **The forces**. Functions taking as input the system state and returning
        as output the forces on each bead. The main forces in any polymer

wlcsim/bd/homolog.py

@@ -3,8 +3,6 @@
 
 See `homolog_points_to_loops_list` for a description of how the homologous
 network structure is handled.
-
-A specialized
 """
 from ..plot import PolymerViewer
 from .rouse import recommended_dt
@@ -156,7 +154,7 @@ def split_homologs_X(X, N, loop_list):
     return X1, X2
 
 def split_homolog_x(x0, N, loop_list):
-    """Make an (N_tot,3) array from rouse_homologs into (2,N,3)."""
+    """Make an (N_tot,3) array from `rouse_homologs` into (2,N,3)."""
     N_tot, _ = x0.shape
     x1 = x0[:N,:].copy()
     x2 = np.zeros((N,3))

wlcsim/bd/rouse.py

@@ -92,13 +92,11 @@
     can be thought of as one large "effective" bead with a diffusivity
     :math:`N` times smaller than the individual beads, and will again exhibit
     regular diffusion (:math:`t^1` scaling).
-
     It is worthwhile to note that the case of the continuous Rouse polymer,
     there is no shortest time scale (because of its fractal nature, wherein at
     each time point, the Rouse chain's spatial configuration traces out a
     Brownian walk). In this case, the :math:`t^{1/2}` behavior simply stretches
     out to arbitrarily short time scales.
-
     For a more detailed discussion of the scaling behavior of Rouse MSD's, see
     the `measured_D_to_rouse` function.
 

wlcsim/bd/runge_kutta.py

@@ -20,11 +20,20 @@
 
 
 def rk4_thermal_lena(f, D, t, x0):
-    """x'(t) = f(x(t), t) + Xi(t), where Xi is thermal, diffusivity D.
+    r"""
+    Integrate using the method of Steph Weber (Phys. Rev. E., 82, 011913).
+
+    Solves the equation
+
+    .. math::
+
+        x'(t) = f(x(t), t) + \xi(t),
 
-    x0 is x(t[0]).
+    where :math:`\xi` represents thermal forces, and each dimension of ``x0``
+    has diffusivity :math:`D`.
 
-    :math:`f: R^n x R -> R^n`
+    ``x0`` is :math:`x(t=0)`, and :math:`f: \mathbb{R}^n \times \mathbb{R} ->
+    \mathbb{R}^n`
     """
     t = np.array(t)
     x0 = np.array(x0)
@@ -53,51 +62,22 @@ def rk4_thermal_lena(f, D, t, x0):
     return x
 
 
-def rk4_thermal_bruno(f, D, t, x0):
+def euler_maruyama(f, D, t, x0):
     r"""
-    Test new method: Attempt to keep :math:`\omega` constant.
+    The most well-known stochastic integrator.
 
-    WARNING: does not converge strongly (autocorrelation function seems
-    higher than should be for OU process...), as is...x'(t) = f(x(t), t) +
-    Xi(t), where Xi is thermal, diffusivity D
+    Solves the equation
 
-    x0 is x(t[0]).
+    .. math::
 
-    :math:`f: R^n x R -> R^n`
-    """
-    t = np.array(t)
-    x0 = np.array(x0)
-    xsize = x0.shape
-    x = np.zeros(t.shape + x0.shape)
-    dts = np.diff(t)
-    x[0] = x0
-    dxdt = np.zeros((4,) + x0.shape)  # one for each RK step
-    Fbrown = np.sqrt(2*D / ((t[1] - t[0])/2))
-    # at each step i, we use data (x,t)[i-1] to create (x,t)[i]
-    # in order to make it easy to pull into a new functin later, we'll call
-    # t[i-1] "t0", old x (x[i-1]) "x0", and t[i]-t[i-1] "h".
-    for i in range(1, len(t)):
-        h = dts[i-1]
-        t0 = t[i-1]
-        x0 = x[i-1]
-        x_est = x0
-        dxdt[0] = f(x0, t0) + Fbrown  # slope at beginning of time step
-        # random force estimate at midpoint
-        Fbrown = np.sqrt(2*D / ((t[i]-t[i-1])/2))*np.random.normal(size=xsize)
-        x_est = x0 + dxdt[0]*(h/2)  # first estimate at midpoint
-        dxdt[1] = f(x_est, t0 + (h/2)) + Fbrown  # estimated slope at midpoint
-        x_est = x0 + dxdt[1]*(h/2)  # second estimate at midpoint
-        dxdt[2] = f(x_est, t0 + (h/2)) + Fbrown  # second slope at midpoint
-        x_est = x0 + dxdt[2]*h  # first estimate at next time point
-        # random force estimate at endpoint (and next start point)
-        Fbrown = np.sqrt(2*D / ((t[i]-t[i-1])/2))*np.random.normal(size=xsize)
-        dxdt[3] = f(x_est, t0 + h) + Fbrown  # slope at end of time step
-        # final estimate is weighted average of slope estimates
-        x[i] = x0 + h*(dxdt[0] + 2*dxdt[1] + 2*dxdt[2] + dxdt[3])/6
-    return x
+        x'(t) = f(x(t), t) + \xi(t),
 
+    where :math:`\xi` represents thermal forces, and each dimension of ``x0``
+    has diffusivity :math:`D`.
 
-def euler_maruyama(f, D, t, x0):
+    ``x0`` is :math:`x(t=0)`, and :math:`f: \mathbb{R}^n \times \mathbb{R} ->
+    \mathbb{R}^n`
+    """
     t = np.array(t)
     x0 = np.array(x0)
     x = np.zeros(t.shape + x0.shape)
@@ -116,7 +96,8 @@ def euler_maruyama(f, D, t, x0):
 
 
 def srk1_roberts(f, D, t, x0):
-    r"""From wiki, from A. J. Roberts. Modify the improved Euler scheme to
+    r"""
+    From wiki, from A. J. Roberts. Modify the improved Euler scheme to
     integrate stochastic differential equations. [1], Oct 2012.
 
     If we have an Ito SDE given by
@@ -129,13 +110,22 @@ def srk1_roberts(f, D, t, x0):
 
     .. math::
 
-        \vec{K}_1 = h \vec{a}(t_k, \vec{X}_k) + (\Delta W_k - S_k\sqrt{h}) \vec{b}(t_k, \vec{X}_k)
-        \vec{K}_2 = h \vec{a}(t_{k+1}, \vec{X}_k + \vec{K}_1) + (\Delta W_k - S_k\sqrt{h}) \vec{b}(t_{k+1}, \vec{X}_k + \vec{K}_1)
+        \vec{K}_1 = h \vec{a}(t_k, \vec{X}_k) + (\Delta W_k - S_k\sqrt{h})
+            \vec{b}(t_k, \vec{X}_k)
+
+    .. math::
+
+        \vec{K}_2 = h \vec{a}(t_{k+1}, \vec{X}_k + \vec{K}_1)
+            + (\Delta W_k - S_k\sqrt{h}) \vec{b}(t_{k+1},\vec{X}_k + \vec{K}_1)
+
+    .. math::
+
         \vec{X}_{k+1} = \vec{X}_k + \frac{1}{2}(\vec{K}_1 + \vec{K}_2)
 
     where :math:`\Delta W_k = \sqrt{h} Z_k` for a normal random :math:`Z_k \sim
     N(0,1)`, and :math:`S_k=\pm1`, with the sign chosen uniformly at random
-    each time."""
+    each time.
+    """
     t = np.array(t)
     x0 = np.array(x0)
     x = np.zeros(t.shape + x0.shape)
@@ -164,8 +154,14 @@ def srk1_roberts(f, D, t, x0):
 
 # simple test case
 def ou(x0, t, k_over_xi, D, method=rk4_thermal_lena):
-    "simulate ornstein uhlenbeck process with theta=k_over_xi and sigma^2/2=D"
-    def f(x,t):
+    r"""
+    Test simulation: the Ornstein-Uhlenbeck process.
+
+    Simulates a standard Ornstein-Uhlenbeck process with ``theta=k_over_xi``
+    and ``sigma^2/2=D``. This is the appropriate parameterization to match the
+    limiting case of ``N=2`` in `wlcsim.bd.rouse.with_integrator`.
+    """
+    def f(x, t):
         return -k_over_xi*x
     return method(f, D=D, t=t, x0=x0)
 
@@ -179,7 +175,7 @@ def _get_scalar_corr(X):
     for i in range(num_samples):
         for j in range(num_t):
             for k in range(j, num_t):
-                corr[k-j] = corr[k-j] + X[i,k]*X[i,j]
+                corr[k-j] = corr[k-j] + X[i, k]*X[i, j]
                 count[k-j] = count[k-j] + 1
     return corr, count
 
@@ -193,7 +189,7 @@ def _get_vector_corr(X):
     for i in range(num_samples):
         for j in range(num_t):
             for k in range(j, num_t):
-                corr[k-j] = corr[k-j] + X[i,k]@X[i,j]
+                corr[k-j] = corr[k-j] + X[i, k]@X[i, j]
                 count[k-j] = count[k-j] + 1
     return corr, count
 
@@ -211,23 +207,23 @@ def _get_bead_msd(X, k=None):
     count = np.zeros((num_t,))
     for i in range(num_t):
         for j in range(i, num_t):
-            ta_msd[j-i] += (X[j,k] - X[i,k])@(X[j,k] - X[i,k])
+            ta_msd[j-i] += (X[j, k] - X[i, k])@(X[j, k] - X[i, k])
             count[j-i] += 1
     return ta_msd, count
 
 
 @jit(nopython=True)
 def _msd(x):
     result = np.zeros_like(x)
-    for delta in range(1,len(x)):
-        for i in range(delta,len(x)):
+    for delta in range(1, len(x)):
+        for i in range(delta, len(x)):
             result[delta] += (x[i] - x[i-delta])**2
         result[delta] = result[delta] / (len(x) - delta)
     return result
 
 
-# test different integrators below on simply OU process
 def test_ou_autocorr(method=srk1_roberts):
+    """Test an individual integrator visually using an OU process."""
     import scipy.stats
     import matplotlib.pyplot as plt
     k = 2