Merge pull request #161 from stan-dev/feature/reduce-sum

Added initial attempt at docs for reduce-sum (design-doc pull request #17)

Author: mitzimorris

src/functions-reference/higher-order_functions.Rmd

@@ -10,7 +10,8 @@ if (knitr::is_html_output()) {
 cat(' * <a href="functions-algebraic-solver.html">Algebraic Equation Solver</a>\n')
 cat(' * <a href="functions-ode-solver.html">Ordinary Differential Equation (ODE) Solvers</a>\n')
 cat(' * <a href="functions-1d-integrator.html">1D Integrator</a>\n')
-cat(' * <a href="functions-map.html">Higher-Order Map</a>\n')
+cat(' * <a href="functions-reduce.html">Reduce-Sum</a>\n')
+cat(' * <a href="functions-map.html">Map-Rect</a>\n')
 }
 ```
 
@@ -134,7 +135,7 @@ package MINPACK-1 [@minpack:1980].
 
 The Jacobian of the solution with respect to auxiliary parameters is
 computed using the implicit function theorem. Intermediate Jacobians
-(of the the algebraic function's output with respect to the unknowns y
+(of the algebraic function's output with respect to the unknowns y
 and with respect to the auxiliary parameters theta) are computed using
 Stan's automatic differentiation.
 
@@ -382,8 +383,92 @@ Internally the 1D integrator uses the double-exponential methods in the Boost 1D
 
 The gradients of the integral are computed in accordance with the Leibniz integral rule. Gradients of the integrand are computed internally with Stan's automatic differentiation.
 
+## Reduce-Sum Function {#functions-reduce}
 
-## Higher-Order Map {#functions-map}
+Stan provides a higher-order reduce function for summation. A function
+which returns a scalar `g: U -> real` is mapped to every element of a
+list of type `U[]`, `{ x1, x2, ... }` and all the results are
+accumulated,
+
+`g(x1) + g(x2) + ...`
+
+For efficiency reasons the reduce function doesn't work with the
+element-wise evaluated function `g` itself, but instead works through
+evaluating partial sums, `f: U[] -> real`, where:
+
+```
+f({ x1 }) = g(x1)
+f({ x1, x2 }) = g(x1) + g(x2)
+f({ x1, x2, ... }) = g(x1) + g(x2) + ...
+```
+
+Mathematically the summation reduction is associative and forming
+arbitrary partial sums in an aribitrary order will not change the
+result. However, floating point numerics on computers only have
+a limited precision such that associativity does not hold
+exactly. This implies that the order of summation determines the exact
+numerical result. For this reason, the higher-order reduce function is
+available in two variants:
+
+* `reduce_sum`: Automatically choose partial sums partitioning based on a dynamic
+ scheduling algorithm.
+* `reduce_sum_static`: Compute the same sum as `reduce_sum`, but partition
+ the input in the same way for given data set (in `reduce_sum` this partitioning
+ might change depending on computer load). This should result in stable
+ numerical evaluations.
+
+### Specifying the Reduce-sum Function
+
+The higher-order reduce function takes a partial sum function `f`, an array argument `x`
+(with one array element for each term in the sum), a recommended
+`grainsize`, and a set of shared arguments. This representation allows
+parallelization of the resultant sum.
+
+<!-- real; reduce_sum; (F f, T[] x, int grainsize, T1 s1, T2 s2, ...); -->
+\index{{\tt \bfseries reduce\_sum }!{\tt (F f, T[] x, int grainsize, T1 s1, T2 s2, ...): real}|hyperpage}
+
+`real` **`reduce_sum`**`(F f, T[] x, int grainsize, T1 s1, T2 s2, ...)`<br>\newline
+`real` **`reduce_sum_static`**`(F f, T[] x, int grainsize, T1 s1, T2 s2, ...)`<br>\newline
+
+Returns the equivalent of `f(1, size(x), x, s1, s2, ...)`, but computes
+the result in parallel by breaking the array `x` into independent
+partial sums. `s1, s2, ...` are shared between all terms in the sum.
+
+* *`f`*: function literal referring to a function specifying the
+partial sum operation. Refer to the [partial sum function](#functions-partial-sum).
+* *`x`*: array of `T`, one for each term of the reduction, `T` can be any type,
+* *`grainsize`*: For `reduce_sum`, `grainsize` is the recommended size of the partial sum (`grainsize = 1` means pick totally automatically). For `reduce_sum_static`, `grainsize` determines the maximum size of the partial sums, type `int`,
+* *`s1`*: first (optional) shared argument, type `T1`, where `T1` can be any type
+* *`s2`*: second (optional) shared argument, type `T2`, where `T2` can be any type,
+* *`...`*: remainder of shared arguments, each of which can be any type.
+
+### The Partial sum Function {#functions-partial-sum}
+
+The partial sum function must have the following signature where the type `T`, and the
+types of all the shared arguments (`T1`, `T2`, ...) match those of the original
+`reduce_sum` (`reduce_sum_static`) call.
+
+```
+(int start, int end, T[] x_subset, T1 s1, T2 s2, ...):real
+```
+
+The partial sum function returns the sum of the `start` to `end` terms (inclusive) of the overall
+calculations. The arguments to the partial sum function are:
+
+*   *`start`*, the index of the first term of the partial sum, type `int`
+
+*   *`end`*, the index of the last term of the partial sum (inclusive), type `int`
+
+*   *`x_subset`*, the subset of `x` a given partial sum is responsible for computing, type `T[]`, where `T` matches the type of `x` in `reduce_sum` (`reduce_sum_static`)
+
+*   *`s1`*, first shared argument, type `T1`, matching type of `s1` in `reduce_sum` (`reduce_sum_static`)
+
+*   *`s2`*, second shared argument, type `T2`, matching type of `s2` in `reduce_sum` (`reduce_sum_static`)
+
+*   *`...`*, remainder of shared arguments, with types matching those in `reduce_sum` (`reduce_sum_static`)
+
+
+## Map-Rect Function {#functions-map}
 
 Stan provides a higher-order map function.  This allows map-reduce
 functionality to be coded in Stan as described in the user's guide.

src/reference-manual/expressions.Rmd

@@ -1118,27 +1118,33 @@ literals.*
 `integrate_1d`, | `real, real, real[]` | `real[], int[]` | `real`
 `integrate_ode_X`, | `real, real[], real[]` | `real[], int[]` | `real[]`
 `map_rect` | `vector, vector` | `real[], int[]` | `vector`
+`reduce_sum` | ```T[], T1, T2, ...``` | | `real`
 
-For example, the rectangular mapping function might be used in the
-following way to compute the log likelihood of a hierarchical model.
+`T`, `T1`, `T2`, and the types of `...` can be any Stan type.
+
+For example, the `integrate_ode_rk45` function can be used to integrate
+differential equations in Stan:
 
 ```stan
 functions {
-  vector foo_ll(vector phi, vector theta, real[] x_r, int[] x_i) {
+  real[] foo(real t, real[] y, real[] theta, real[] x_r, int[] x_i) {
     ...
 ...
-vector[11] phi;
-vector[2] thetas[N];
-real x_rs[N, 5];
-real x_is[N, 0];
+int<lower=1> T;
+real y0[2];
+real t0;
+real ts[T];
+real theta[1];
+real x_r[0];
+int x_i[0];
 ...
-target += sum(map_rect(foo_ll, phi, thetas, x_rs, x_is));
+real y_hat[T, 2] = integrate_ode_rk45(foo, y0, t0, ts, theta, x_r, x_i);
 ```
 
 The function argument is `foo`, the name of the user-defined
 function;  as shown in the [higher-order functions table](#higher-order-functions), `foo`
-takes two vectors, a real array, and an integer array as arguments and
-returns a vector.
+takes a real, three more real arrays, and an integer array as arguments and
+returns a real array.
 
 
 ### Functions Passed by Reference {-}

src/stan-users-guide/_bookdown.yml

@@ -29,7 +29,7 @@ rmd_files: [
   "problematic-posteriors.Rmd",
   "reparameterization.Rmd",
   "efficiency-tuning.Rmd",
-  "map-reduce.Rmd",
+  "parallelization.Rmd",
 
   "part-appendices.Rmd",
   "style-guide.Rmd",