Merge pull request #22 from zoj613/numpify

Author: zoj613

Makefile

@@ -1,6 +1,6 @@
 #Copyright (c) 2020, Zolisa Bleki
 #SPDX-License-Identifier: BSD-3-Clause */
-.PHONY: clean pkg test wheels cythonize lib
+.PHONY: clean pkg test wheels cythonize lib install
 
 NAME := htnorm
 CC := gcc
@@ -57,6 +57,9 @@ test:
 cythonize:
 	cythonize pyhtnorm/*.pyx
 
+install: clean cythonize
+	poetry install
+
 sdist: cythonize
 	poetry build -f sdist
 

README.md

@@ -77,16 +77,19 @@ int main (void)
 
 ## Python API
 
+### Dependencies
+- NumPy >= 1.17
+
 A high level python interface to the library is also provided. Linux users can 
 install it using wheels via pip (thus not needing to worry about availability of C libraries),
 ```bash
-pip install pyhtnorm
+pip install -U pyhtnorm
 ```
 Wheels are not provided for MacOS. To install via pip, one can run the following commands:
 ```bash
 #set the path to LAPACK shared library
 export LIBS_DIR=<some directory>
-pip install pyhtnorm
+pip install -U pyhtnorm
 ```
 Alternatively, one can install it from source. This requires an installation of [poetry][7] and the following shell commands:
 
@@ -99,37 +102,34 @@ $ 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)
+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). The python
+API is such that the code can be easily integrated into other existing libraries.
+Since `v1.0.0`, it supports passing a `numpy.random.Generator` instance.
+Thus, one can reuse the same generator without having to declare one specifically for `htnorm`.
 
 ```python
-from pyhtnorm import HTNGenerator
+from pyhtnorm import hyperplane_truncated_mvnorm
 import numpy as np
 
-rng = HTNGenerator()
+rng = np.random.default_rng()
 
 # generate example input
-k1 = 1000
-k2 = 1
-npy_rng = np.random.default_rng()
-temp = npy_rng.random((k1, k1))
-cov = temp @ temp.T + np.diag(npy_rng.random(k1))
+k1, k2 = 1000, 1
+temp = rng.random((k1, k1))
+cov = temp @ temp.T
 G = np.ones((k2, k1))
 r = np.zeros(k2)
-mean = npy_rng.random(k1)
+mean = rng.random(k1)
 
-samples = rng.hyperplane_truncated_mvnorm(mean, cov, G, r)
-print(sum(samples))  # verify if sampled values sum to zero
+# passing `random_state` is optional. If the argument is not used, a fresh
+# random generator state is instantiated internally using system entropy.
+o = hyperplane_truncated_mvnorm(mean, cov, G, r, random_state=rng)
+print(o.sum())  # verify if sampled values sum to zero
 # 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)
-print(out.sum())  # verify
+hyperplane_truncated_mvnorm(mean, cov, G, r, out=o)
 ```
 
-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.
+For more information about the function's arguments, refer to its docstring.
 
 A pure numpy implementation is demonstrated in this [example script][9].
 
@@ -183,9 +183,9 @@ see the [LICENSE][6] file.
 [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
+[4]: ./include/htnorm.h 
+[5]: ./include/rng.h
+[6]: ./LICENSE
 [7]: https://python-poetry.org/docs/pyproject/
 [8]: https://www.sciencedirect.com/science/article/abs/pii/S1877584517301600
-[9]: https://github.com/zoj613/htnorm/blob/main/examples/numpy_implementation.py
+[9]: ./examples/numpy_implementation.py

build-wheels.sh

@@ -10,27 +10,32 @@ bin_arr=(
     /opt/python/cp36-cp36m/bin
     /opt/python/cp37-cp37m/bin
     /opt/python/cp38-cp38/bin
+    /opt/python/cp39-cp39/bin
 )
+
 # add  python to image's path
 export PATH=/opt/python/cp38-cp38/bin/:$PATH
-# download && install poetry
-curl -sSL https://raw.githubusercontent.com/python-poetry/poetry/master/get-poetry.py | python
-
+# download install script
+curl -#sSL https://raw.githubusercontent.com/python-poetry/poetry/master/get-poetry.py > get-poetry.py
+# install using local archive
+python get-poetry.py -y --file poetry-1.1.4-linux.tar.gz
 # install openblas
 yum install -y openblas-devel
 
 function build_poetry_wheels
 {
-    # build wheels for 3.6-3.8 with poetry 
+    # build wheels for 3.6-3.9 with poetry 
     for BIN in "${bin_arr[@]}"; do
         rm -Rf build/*
+        # install build deps
+        "${BIN}/python" ${HOME}/.poetry/bin/poetry run pip install numpy
         BUILD_WHEELS=1 "${BIN}/python" ${HOME}/.poetry/bin/poetry build -f wheel
-    done
-
-    # add C libraries to wheels
-    for whl in dist/*.whl; do
-        auditwheel repair "$whl" --plat $1
-        rm "$whl"
+        auditwheel repair dist/*.whl --plat $1
+        whl="$(basename dist/*.whl)"
+        "${BIN}/python" -m pip install wheelhouse/"$whl"
+        # test if installed wheel imports correctly
+        "${BIN}/python" -c "from pyhtnorm import *"
+        rm dist/*.whl
     done
 }
 

build.py

@@ -1,6 +1,8 @@
 from distutils.core import Extension
 import os
 
+import numpy as np
+
 
 source_files = [
     "pyhtnorm/_htnorm.c",
@@ -28,11 +30,12 @@
 
 extensions = [
     Extension(
-        "pyhtnorm._htnorm",
+        "_htnorm",
         source_files,
-        include_dirs=['./include'],
+        include_dirs=[np.get_include(), './include'],
         library_dirs=library_dirs,
         libraries=libraries,
+        define_macros=[('NPY_NO_DEPRECATED_API', 0)],
         extra_compile_args=['-std=c99']
     )
 ]

examples/numpy_implementation.py

@@ -10,40 +10,42 @@
 import numpy as np
 
 
-class Generator(np.random.Generator):
-    def hyperplane_truncated_mvnorm(self, mean, cov, G, r, diag=False):
-        if diag:
-            cov_diag = np.diag(cov)
-            y = mean + cov_diag * self.standard_normal(mean.shape[0])
-            covg = cov_diag[:, None] * G.T
-        else:
-            y = self.multivariate_normal(mean, cov, method='cholesky')
-            covg = cov @ G.T
-        gcovg = G @ covg
-        alpha = np.linalg.solve(gcovg, r - G @ y)
-        return y + covg @ alpha
+def hyperplane_truncated_mvnorm(mean, cov, G, r, diag=False, random_state=None):
+    rng = np.random.default_rng(random_state)
+    if diag:
+        cov_diag = np.diag(cov)
+        y = mean + cov_diag * rng.standard_normal(mean.shape[0])
+        covg = cov_diag[:, None] * G.T
+    else:
+        y = rng.multivariate_normal(mean, cov, method='cholesky')
+        covg = cov @ G.T
+    gcovg = G @ covg
+    alpha = np.linalg.solve(gcovg, r - G @ y)
+    return y + covg @ alpha
 
-    def structured_precision_mvnorm(self, mean, a, phi, omega, a_type=0, o_type=0):
-        if a_type:
-            Ainv = 1 / np.diag(a)
-            y1 = self.standard_normal(a.shape[0]) / np.sqrt(Ainv)
-            ainv_phi = Ainv[:, None] * phi.T
-        else:
-            Ainv = np.linalg.inv(a)
-            y1 = self.multivariate_normal(
-                np.zeros(a.shape[0]), Ainv, method='cholesky'
-            )
-            ainv_phi = Ainv @ phi.T
-        if o_type:
-            omegainv = 1 / np.diag(omega)
-            y2 = self.standard_normal(omega.shape[0]) / np.sqrt(omegainv)
-            alpha = np.linalg.solve(
-                np.diag(omegainv) + phi @ ainv_phi, phi @ y1 + y2
-            )
-        else:
-            omegainv = np.linalg.inv(omega)
-            y2 = self.multivariate_normal(
-                np.zeros(omega.shape[0]), omegainv, method='cholesky'
-            )
-            alpha = np.linalg.solve(omegainv + phi @ ainv_phi, phi @ y1 + y2)
-        return mean + y1 - ainv_phi @ alpha
+
+def structured_precision_mvnorm(mean, a, phi, omega, a_type=0, o_type=0, random_state=None):
+    rng = np.random.default_rng(random_state)
+    if a_type:
+        Ainv = 1 / np.diag(a)
+        y1 = rng.standard_normal(a.shape[0]) / np.sqrt(Ainv)
+        ainv_phi = Ainv[:, None] * phi.T
+    else:
+        Ainv = np.linalg.inv(a)
+        y1 = rng.multivariate_normal(
+            np.zeros(a.shape[0]), Ainv, method='cholesky'
+        )
+        ainv_phi = Ainv @ phi.T
+    if o_type:
+        omegainv = 1 / np.diag(omega)
+        y2 = rng.standard_normal(omega.shape[0]) / np.sqrt(omegainv)
+        alpha = np.linalg.solve(
+            np.diag(omegainv) + phi @ ainv_phi, phi @ y1 + y2
+        )
+    else:
+        omegainv = np.linalg.inv(omega)
+        y2 = rng.multivariate_normal(
+            np.zeros(omega.shape[0]), omegainv, method='cholesky'
+        )
+        alpha = np.linalg.solve(omegainv + phi @ ainv_phi, phi @ y1 + y2)
+    return mean + y1 - ainv_phi @ alpha

include/rng.h

@@ -15,11 +15,11 @@
  * 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 {
+typedef struct {
     // 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);
+    uint64_t (*next_uint64)(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;