Merge pull request #10 from zoj613/documentation

DOC: Add documentation

Author: zoj613

README.md

@@ -1,3 +1,172 @@
-# Fast Simulation of Hyperplane-Truncated Multivatiate Normal Distributions
+# htnorm
 
-An implementation of a fast and exact simulation algorithm for a multivariate normal distribution truncated on the intersection of a set of hyperplanes.
+This repo provides a C implementation of a fast and exact sampler from a 
+multivariate normal distribution (MVN) truncated on a hyperplane as described [here][1]
+
+this repo implements the following from the paper:
+
+- efficient Sampling from a MVN truncated on a hyperplane: 
+
+    ![hptrunc](https://latex.codecogs.com/svg.latex?%5Cmathbf%7Bx%7D%20%5Csim%20%5Cmathcal%7BN%7D_%7B%5Cmathcal%7BS%7D%7D%28%5Cmathbf%7B%5Cmu%7D%2C%20%5Cmathbf%7B%5CSigma%7D%29%3B%20%5Chspace%7B2mm%7D%20%5Cmathcal%7BS%7D%20%3D%20%5C%7B%5Cmathbf%7Bx%7D%20%3A%20%5Cmathbf%7BG%7D%5Cmathbf%7Bx%7D%20%3D%20%5Cmathbf%7Br%7D%5C%7D%2C%20%5Cmathbf%7BG%7D%20%5Cin%20%5Cmathcal%7BR%7D%5E%7Bk_2%20%5Ctimes%20k%7D%2C%20rank%28%5Cmathbf%7BG%7D%29%20%3D%20k_2%20%3C%20k)
+
+- efficient sampling from a MVN with a stuctured precision matrix: 
+
+    ![struc](https://latex.codecogs.com/svg.latex?%5Cmathbf%7Bx%7D%20%5Csim%20%5Cmathcal%7BN%7D%5C%5B%5Cmathbf%7B%5Cmu%7D%2C%20%28%5Cmathbf%7BA%7D%20+%20%5Cmathbf%7B%5CPhi%7D%5ET%5Cmathbf%7B%5COmega%7D%5Cmathbf%7B%5CPhi%7D%29%5E%7B-1%7D%5C%5D%3B%20%5Chspace%7B2mm%7D%20%5Cmathbf%7B%5CPhi%7D%20%5Cin%20%5Cmathcal%7BR%7D%5E%7Bn%20%5Ctimes%20p%7D%2C%20%5Cmathbf%7B%5COmega%7D%20%5Cin%20%5Cmathcal%7BR%7D%5E%7Bn%20%5Ctimes%20n%7D%2C%20%5Cmathbf%7BA%7D%20%5Cin%20%5Cmathcal%7BR%7D%5E%7Bp%20%5Ctimes%20p%7D)
+
+- efficent sampling frfom a MVN with a structured precision and mean:
+
+    ![strucmean](https://latex.codecogs.com/svg.latex?%5Cmathbf%7Bx%7D%20%5Csim%20%5Cmathcal%7BN%7D%5CBig%5C%5B%28%5Cmathbf%7BA%7D%20+%20%5Cmathbf%7B%5CPhi%7D%5ET%5Cmathbf%7B%5COmega%7D%5Cmathbf%7B%5CPhi%7D%29%5E%7B-1%7D%5Cmathbf%7B%5CPhi%7D%5ET%5Cmathbf%7B%5COmega%7D%5Cmathbf%7Bt%7D%2C%20%28%5Cmathbf%7BA%7D%20+%20%5Cmathbf%7B%5CPhi%7D%5ET%5Cmathbf%7B%5COmega%7D%5Cmathbf%7B%5CPhi%7D%29%5E%7B-1%7D%5CBig%5C%5D%3B%20%5Chspace%7B2mm%7D%20%5Cmathbf%7B%5COmega%7D%20%5Cin%20%5Cmathcal%7BR%7D%5E%7Bn%20%5Ctimes%20n%7D%2C%20%5Cmathbf%7BA%7D%20%5Cin%20%5Cmathcal%7BR%7D%5E%7Bp%20%5Ctimes%20p%7D)
+
+The algorithms implemented have the following practical applications:
+- Topic models when unknown parameters can be interpreted as fractions.
+- Admixture models
+- discrete graphical models
+- Sampling from posterior distribution of an Intrinsic Conditional Autoregressive prior [icar][8]
+- Sampling from posterior conditional distributions of various bayesian regression problems.
+
+
+## Dependencies
+
+- a C compiler that supports the C99 standard or later
+- an installation of BLAS and LAPACK that exposes its C interface via the headers `<cblas.h>` and `<lapacke.h>`
+(e.g openBLAS).
+
+
+## Usage
+
+Building a shared library of `htnorm` can be done with the following:
+```bash
+$ cd src/
+# optionally set path to CBLAS and LAPACKE headers using INCLUDE_DIRS environmental variable
+$ export INCLUDE_DIRS="some/path/to/headers" 
+# optionally set path to BLAS installation shared library
+$ export LIBS_DIR="some/path/to/library/"
+# optionally set the linker flag for your BLAS installation (e.g -lopenblas)
+$ export LIBS=<flag here>
+$ make lib
+```
+Afterwards the shared library will be found in a `lib/` directory of the project root,
+and the library can be linked dynamically via `-lhtnorm`.
+
+The puplic API exposes the samplers through the function declarations
+- `int htn_hyperplane_truncated_mvn(rng_t* rng, const ht_config_t* conf, double* out)`
+- `int htn_structured_precision_mvn(rng_t* rng, const sp_config_t* conf, double* out)`
+
+The details of the parameters are documented in ther header files ["htnorm.h"][4].
+
+Random number generation is done using [PCG64][2] or [Xoroshiro128plus][3] bitgenerators. 
+The API allows using a custom generator, and the details are documented in the header file 
+["rng.h"][5].
+
+## Examples
+```C
+#include "htnorm.h"
+
+int main ()
+{
+    ...
+
+    // instantiate a random number generator
+    rng_t* rng = rng_new_pcg64();
+    ht_config_t config;
+    config.g = ...; // G matrix
+    config.gnrow = ...; // number of rows of G
+    config.gncol = ...; // number of columns of G
+    cofig.r = ...; // r array
+    config.mean = ...; // mean array
+    config.cov = ...; // the convariance matrix
+    confi.diag = ...; // whether covariance is diagonal
+
+    double* samples = ...; // array to store the samples
+    // now call the sampler
+    int res_info = htn_hyperplane_truncated_mvn(rng, &config, samples);
+
+    // res_info contains a number that indicates whether sampling failed or not.
+
+    ...
+
+    // finally free the RNG pointer at some point
+    rng_free(rng);
+
+    ...
+    return 0;
+}
+```
+
+## Python API
+
+A high level python interface to the library is also provided. Installing it from 
+source requires an installation of [poetry][7] and the following shell commands:
+
+```bash
+$ git clone https://github.com/zoj613/htnorm.git
+$ cd htnorm/
+$ poetry install
+# add htnorm to python's path
+$ export PYTHONPATH=$PWD:$PYTHONPATH
+```
+Below is an example of how to use htnorm in python to sample from a multivariate
+gaussian truncated on the hyperplane ![sumzero](https://latex.codecogs.com/svg.latex?%5Cmathbf%7B1%7D%5ET%5Cmathbf%7Bx%7D%20%3D%200) (i.e. making sure the sampled values sum to zero)
+
+```python
+from pyhtnorm import HTNGenerator
+import numpy as np
+
+rng = HTNGenerator()
+
+# generate example input
+k1 = 1000
+k2 = 1
+npy_rng = np.random.default_rng()
+temp = npy_rng.random((k1, k1))
+cov = a @ a.T + np.diag(npy_rng.random(k1))
+G = np.ones((k2, k1))
+r = np.zeros(k1)
+mean = npy_rng.random(k1)
+
+samples = rng.hyperplane_truncated_mvnorm(mean, cov, G, r)
+# verify if sampled values sum to zero
+print(sum(samples))
+
+# alternatively one can pass an array to store the results in
+out = np.empty(k1)
+rng.hyperplane_truncated_mvnorm(mean, cov, G, r, out=out)
+# verify
+print(out.sum())
+```
+
+For more details about the parameters of the `HTNGenerator` and its methods,
+see the docstrings via python's `help` function.
+
+The python API also exposes the `HTNGenerator` class as a Cython extension type
+that can be "cimported" in a cython script.
+
+
+## TODO
+
+- Add an `R` interface to the library.
+
+
+## Licensing
+
+`htnorm` is free software made available under the BSD-3 License. For details
+see the [LICENSE][6] file.
+
+
+## References
+- Cong, Yulai; Chen, Bo; Zhou, Mingyuan. Fast Simulation of Hyperplane-Truncated 
+   Multivariate Normal Distributions. Bayesian Anal. 12 (2017), no. 4, 1017--1037. 
+   doi:10.1214/17-BA1052.
+- Bhattacharya, A., Chakraborty, A., and Mallick, B. K. (2016). 
+  “Fast sampling with Gaussian scale mixture priors in high-dimensional regression.” 
+  Biometrika, 103(4):985. 
+
+
+[1]: https://projecteuclid.org/euclid.ba/1488337478
+[2]: https://www.pcg-random.org/
+[3]: https://en.wikipedia.org/wiki/Xoroshiro128%2B
+[4]: https://github.com/zoj613/htnorm/blob/main/include/htnorm.h 
+[5]: https://github.com/zoj613/htnorm/blob/main/include/rng.h
+[6]: https://github.com/zoj613/htnorm/blob/main/LICENSE
+[7]: https://python-poetry.org/docs/pyproject/
+[8]: https://www.sciencedirect.com/science/article/abs/pii/S1877584517301600

include/htnorm.h

@@ -1,6 +1,27 @@
 /* Copyright (c) 2020, Zolisa Bleki
  *
- * SPDX-License-Identifier: BSD-3-Clause */
+ * SPDX-License-Identifier: BSD-3-Clause 
+ *
+ * Fast Simulation of Hyperplane-Truncated Multivaariate Normal Distributions.
+ *
+ * This library implements the algorithms 2, 4 and the one in example 4 of
+ * Cong, Y., Chen, B., & Zhou, M. (2017).
+ *
+ * Algorithm 2 allows fast and exact sampling from a multivariate normal (MVN)
+ * distribution that is trancated on a hyperplane.
+ *
+ * Algorthm 4 allows efficient sampling from a MVN with a structured precision
+ * matrix.
+ *
+ * Algorithm described in example 4 allows efficient sampling from a MVN with
+ * a structured precision matrix and a structured mean dependent on the
+ * structure of the precision matrix.
+ *
+ * References:
+ * 1) Cong, Y., Chen, B., & Zhou, M. (2017). Fast simulation of
+ *    hyperplane-truncated multivariate normal distributions. Bayesian Analysis,
+ *    12(4), 1017-1037.
+ */
 #ifndef HTNORM_HTNORM_H
 #define HTNORM_HTNORM_H
 
@@ -11,6 +32,7 @@
 
 #define INIT_HT_CONFIG(x) (x) = {.diag = false}
 
+// initialize a pointer to a `ht_config_t` struct.
 #define INIT_SP_CONFIG(x) \
     (x) = {.struct_mean = false, .a_id = NORMAL, .o_id = NORMAL} 
 
@@ -54,8 +76,53 @@ typedef struct {
     bool struct_mean; 
 } sp_config_t;
 
-
+/* Sample from a multivariate normal distribution truncated on a hyperplane.
+ *
+ * Generate sample x from the distribution: N(mean, cov) truncated on the plane
+ * {x | Gx = r}, where the rank of G is less than that of the covariance.
+ *
+ * Paramaters
+ * ----------
+ *  rng:
+ *      A pointer to the `rng_t` struct (The random number generator).
+ *  conf:
+ *      A pointer to the input configuration struct `ht_config_t`.
+ *  out:
+ *      An array to store the generated samples.
+ *
+ *  Returns
+ *  -------
+ *  Zero if the sampling was successful. A positive integer is returned if the
+ *  sampled failed because the covariance is not positive definite or if a
+ *  factorization of the covariance was not successful. A negative integer is
+ *  returned if one of the inputs contains an illegal value (non-numerical/NaN).
+ */
 int htn_hyperplane_truncated_mvn(rng_t* rng, const ht_config_t* conf, double* out);
+
+
+/* Sample from a MVN with a structured precision matrix and/or structured mean.
+ *
+ * Sample from a MVN: N(mean, (A + phi^T * Omega * phi)^-1) or
+ * N((A + phi^T * Omega * phi)^-1 * phi^T * t, (A + phi^T * Omega * phi)^-1))
+ * This algorithm in very efficient and ideal when the dimension of matrix A
+ * is greater than that of matrix Omega.
+ *
+ * Parameters
+ * ----------
+ *  rng:
+ *      A pointer to the `rng_t` struct (The random number generator).
+ *  conf:
+ *      A pointer to the input configuration struct `sp_config_t`.
+ *  out:
+ *      An array to store the generated samples.
+ *
+ *  Returns
+ *  -------
+ *  Zero if the sampling was successful. A positive integer is returned if the
+ *  sampled failed because the covariance is not positive definite or if a
+ *  factorization of the covariance was not successful. A negative integer is
+ *  returned if one of the inputs contains an illegal value (non-numerical/NaN).
+ * */
 int htn_structured_precision_mvn(rng_t* rng, const sp_config_t* conf, double* out);
 
 #endif

include/rng.h

@@ -1,24 +1,50 @@
 /* Copyright (c) 2020, Zolisa Bleki
  *
- * SPDX-License-Identifier: BSD-3-Clause */
+ * SPDX-License-Identifier: BSD-3-Clause 
+ *
+ * The public API of the underlying random number generator. The `rng_t`
+ * struct is passed as a parameter to both sampling functions in the htnorm
+ * header.
+ */
 #ifndef HTNORM_RNG_H
 #define HTNORM_RNG_H
 
 #include <stdint.h>
 
+/* The bitgenerator interface. Users can use this struct in order to define
+ * their own bitgenerator, The minimum requirement is to provide a function
+ * that generate 64bit unsigned integers and doubles in the range (0, 1).
+ */
 typedef struct bitgen {
+    // the bsse bitgenerator
     void* base;
+    // a function pointer that takes `base` as an input an returns a positive intgger
     uint64_t (*next_int)(void* base);
+    // a function pointer that takes `base` as an input an returns a double in the range (0, 1)
     double (*next_double)(void* base);
 } rng_t;
 
+// Helper function to deallocate a `rng_t` pointer.
 void rng_free(rng_t* rng);
 
+/* The following functions provide an interface to easily generate random
+ * integers or standard uniform variables using either PCG64 or Xoroshiro128plus
+ * bit generators.
+ *
+ * For more details, see: 
+ * 1) https://www.pcg-random.org/
+ * 2) https://www.pcg-random.org/posts/visualizing-the-heart-of-some-prngs.html
+ * 3) https://www.pcg-random.org/posts/bounded-rands.html
+ */
 
+// Return a pointer to rng_t based on PCG64 bit generator that is randomly seeded.
 rng_t* rng_pcg64_new(void);
+// Return a pointer to rng_t based n PCG64 bit generator with a specified seed.
 rng_t* rng_pcg64_new_seeded(uint64_t seed);
 
+// Return a pointer to rng_t based on Xoroshiro128plus that is randomly seeded.
 rng_t* rng_xrs128p_new(void);
+// Return a pointer to rng_t based on Xoroshiro128plus with a specified seed.
 rng_t* rng_xrs128p_new_seeded(uint64_t seed);
 
 #endif