Merge pull request #25 from zoj613/clean

Author: zoj613

DESCRIPTION

@@ -1,5 +1,5 @@
 Package: htnorm
-Version: 0.2.0
+Version: 1.0.0
 Title: Fast Simulation of Hyperplane-Truncated Multivariate Normal Distributions
 Author: Zolisa Bleki
 Maintainer: Zolisa Bleki <zolisa.bleki@gmail.com>

Makefile

@@ -16,9 +16,9 @@ LIBS_DIR ?= /usr/lib
 override LIBS_DIR := -L$(LIBS_DIR)
 LIBS := -lm -lblas -llapack
 
-SRCFILES = src/dist.c src/htnorm.c src/rng.c
+SRCFILES = src/htnorm_distributions.c src/htnorm.c src/htnorm_rng.c
 
-OBJ = src/dist.o src/htnorm.o src/rng.o
+OBJ = src/htnorm_distributions.o src/htnorm.o src/htnorm_rng.o
 
 
 %.o: %.c

R/htnorm.r

@@ -1,7 +1,7 @@
 # Copyright (c) 2020, Zolisa Bleki
 # SPDX-License-Identifier: BSD-3-Clause
 
-matrix_type <- c(0, 1, 2)
+matrix_type <- c("regular", "diagonal", "identity")
 
 
 validate_output <- function(res) {
@@ -60,10 +60,20 @@ strprec_mvn <- function(rng, mean, a, phi, omega, str_mean, a_id, o_id, out) {
     )
 
     if (!is.element(a_id, matrix_type) || !is.element(o_id, matrix_type))
-        stop("`a_type` and `o_type` need to be one of {0, 1, 2}")
+        stop("`a_type` and `o_type` need to be one of {'regular', 'diagonal', 'identity'}")
 
-    a_id <- as.integer(a_id)
-    o_id <- as.integer(o_id)
+    a_id <- switch(
+        a_id,
+        "regular" = as.integer(0),
+        "diagonal" = as.integer(1),
+        "identity" = as.integer(2)
+    )
+    o_id <- switch(
+        o_id,
+        "regular" = as.integer(0),
+        "diagonal" = as.integer(1),
+        "identity" = as.integer(2)
+    )
 
     if (is.null(out))
         out <- rep(0, length(mean))
@@ -118,7 +128,7 @@ strprec_mvn <- function(rng, mean, a, phi, omega, str_mean, a_id, o_id, out) {
 #' phi <- eig$vectors
 #' omega <- diag(eig$values)
 #' a <- diag(runif(length(mean)))
-#' rng$structured_precision_mvnorm(mean, a, phi, omega, a_type = 1, out = out)
+#' rng$structured_precision_mvnorm(mean, a, phi, omega, a_type = "diagonal", out = out)
 HTNGenerator <- function(seed = NULL, gen = "xrs128p") {
 
     if (is.numeric(seed)) {
@@ -151,7 +161,7 @@ HTNGenerator <- function(seed = NULL, gen = "xrs128p") {
     }
 
     res$structured_precision_mvnorm <- function(
-        mean, a, phi, omega, str_mean = FALSE, a_type = 0, o_type = 0, out= NULL
+        mean, a, phi, omega, str_mean = FALSE, a_type = "regular", o_type = "regular", out= NULL
     ) {
         if (is.null(out)) {
             strprec_mvn(

README.md

@@ -109,7 +109,7 @@ interface is such that the code can be easily integrated into other existing lib
 Since `v1.0.0`, it supports passing a `numpy.random.Generator` instance as a parameter to aid reproducibility.
 
 ```python
-from pyhtnorm import hyperplane_truncated_mvnorm
+from pyhtnorm import hyperplane_truncated_mvnorm, structured_precision_mvnorm
 import numpy as np
 
 rng = np.random.default_rng()
@@ -166,6 +166,13 @@ sum(samples)
 out <- rep(0, 1000)
 rng$hyperplane_truncated_mvnorm(mean, cov, G, r, out = out)
 sum(out)  #verify
+
+out <- rep(0, 1000)
+eig <- eigen(cov)
+phi <- eig$vectors
+omega <- diag(eig$values)
+a <- diag(runif(length(mean)))
+rng$structured_precision_mvnorm(mean, a, phi, omega, a_type = "diagonal", out = out)
 ```
 
 ## Licensing
@@ -187,7 +194,7 @@ see the [LICENSE][6] file.
 [2]: https://www.pcg-random.org/
 [3]: https://en.wikipedia.org/wiki/Xoroshiro128%2B
 [4]: ./include/htnorm.h 
-[5]: ./include/rng.h
+[5]: ./include/htnorm_rng.h
 [6]: ./LICENSE
 [7]: https://python-poetry.org/docs/pyproject/
 [8]: https://www.sciencedirect.com/science/article/abs/pii/S2211675317301574 

build-wheels.sh

@@ -7,35 +7,38 @@ set -e -u -x
 cd $(dirname $0)
 
 bin_arr=(
-    /opt/python/cp36-cp36m/bin
-    /opt/python/cp37-cp37m/bin
+    #/opt/python/cp37-cp37m/bin
     /opt/python/cp38-cp38/bin
-    /opt/python/cp39-cp39/bin
+    #/opt/python/cp39-cp39/bin
 )
 
 # add  python to image's path
-export PATH=/opt/python/cp38-cp38/bin/:$PATH
+export PATH=/opt/python/cp37-cp37m/bin/:$PATH
 # 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
+python get-poetry.py -y --file poetry-1.1.5-linux.tar.gz
 # install openblas
 yum install -y openblas-devel
 
 function build_poetry_wheels
 {
-    # build wheels for 3.6-3.9 with poetry 
+    # build wheels for 3.7-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
+        "${BIN}/python" ${HOME}/.poetry/bin/poetry run pip install numpy==1.18.1
         BUILD_WHEELS=1 "${BIN}/python" ${HOME}/.poetry/bin/poetry build -f wheel
-        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
+        mkdir -p ./tmp
+        for whl in dist/*.whl; do
+            auditwheel repair "$whl" --plat $1 -w "./tmp" 
+            whlname="$(basename "$(echo ./tmp/*.whl)")"
+            "${BIN}/python" -m pip install ./tmp/"$whlname"
+            # test if installed wheel imports correctly
+            "${BIN}/python" -c "from pyhtnorm import *"
+            mv ./tmp/*.whl wheelhouse/
+            rm "$whl"
+        done
     done
 }
 

build.py

@@ -6,8 +6,8 @@
 
 source_files = [
     "pyhtnorm/_htnorm.c",
-    "src/rng.c",
-    "src/dist.c",
+    "src/htnorm_rng.c",
+    "src/htnorm_distributions.c",
     "src/htnorm.c"
 ]
 
@@ -32,7 +32,7 @@
     Extension(
         "_htnorm",
         source_files,
-        include_dirs=[np.get_include(), './include'],
+        include_dirs=[np.get_include(), './include', 'src'],
         library_dirs=library_dirs,
         libraries=libraries,
         define_macros=[('NPY_NO_DEPRECATED_API', 0)],

include/htnorm.h

@@ -28,7 +28,7 @@
 #include <stdbool.h>
 #include <stddef.h>
 
-#include "rng.h"
+#include "htnorm_rng.h"
 
 
 typedef enum {NORMAL, DIAGONAL, IDENTITY} type_t;

include/rng.h include/htnorm_rng.h

@@ -34,21 +34,68 @@ References
 
 """
 from cpython.pycapsule cimport PyCapsule_GetPointer
+from libc.stdint cimport uint64_t
 
-import numpy as np
 cimport numpy as np
 from numpy.random cimport BitGenerator, bitgen_t
 
-from pyhtnorm.c_htnorm cimport *
+from numpy.random import default_rng
 
 
 np.import_array()
 
-cdef extern from "../src/dist.h":
+cdef extern from "htnorm_distributions.h":
     int HTNORM_ALLOC_ERROR
 
 
-cdef inline void validate_return_info(info):
+cdef extern from "htnorm_rng.h" nogil:
+    ctypedef struct rng_t:
+        void* base
+        uint64_t (*next_uint64)(void* state)
+        double (*next_double)(void* state)
+
+
+cdef extern from "htnorm.h" nogil:
+    ctypedef enum mat_type "type_t":
+        NORMAL
+        DIAGONAL
+        IDENTITY
+
+    ctypedef struct ht_config_t:
+        size_t gnrow
+        size_t gncol
+        const double* mean
+        const double* cov
+        const double* g
+        const double* r
+        bint diag
+
+    ctypedef struct sp_config_t:
+        mat_type a_id
+        mat_type o_id
+        size_t pnrow
+        size_t pncol
+        const double* mean
+        const double* a
+        const double* phi
+        const double* omega
+        bint struct_mean
+
+    void init_ht_config(ht_config_t* conf, size_t gnrow, size_t gncol,
+                        const double* mean, const double* cov, const double* g,
+                        const double* r, bint diag)
+
+    void init_sp_config(sp_config_t* conf, size_t pnrow, size_t pncol,
+                        const double* mean, const double* a, const double* phi,
+                        const double* omega, bint struct_mean, mat_type a_id,
+                        mat_type o_id)
+
+    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)
+
+
+cdef inline void validate_return_info(int info):
     if info == HTNORM_ALLOC_ERROR:
         raise MemoryError("Not enough memory to allocate resources.")
     elif info < 0:
@@ -62,13 +109,14 @@ cdef inline void validate_return_info(info):
         )
 
 
-cdef set _VALID_MATRIX_TYPES = {NORMAL, DIAGONAL, IDENTITY}
+cdef dict MAT_TYPE = {"regular": NORMAL, "diagonal": DIAGONAL, "identity": IDENTITY}
+cdef const char* BITGEN_NAME = "BitGenerator"
 
 
-cdef inline void initialize_rng(BitGenerator bitgenerator, rng_t* htnorm_rng):
+cdef inline void initialize_rng(object bitgenerator, rng_t* htnorm_rng):
     cdef bitgen_t* bitgen
 
-    bitgen = <bitgen_t*>PyCapsule_GetPointer(bitgenerator.capsule, "BitGenerator")
+    bitgen = <bitgen_t*>PyCapsule_GetPointer(bitgenerator.capsule, BITGEN_NAME)
     htnorm_rng.base = bitgen.state
     htnorm_rng.next_uint64 = bitgen.next_uint64
     htnorm_rng.next_double = bitgen.next_double
@@ -135,7 +183,6 @@ def hyperplane_truncated_mvnorm(
         the algorithm could not successfully generate the samples.
 
     """
-    cdef BitGenerator bitgenerator
     cdef rng_t rng
     cdef ht_config_t config
     cdef int info
@@ -152,7 +199,7 @@ def hyperplane_truncated_mvnorm(
     init_ht_config(&config, g.shape[0], g.shape[1], &mean[0],
                    &cov[0, 0], &g[0, 0], &r[0], diag)
 
-    bitgenerator = np.random.default_rng(random_state)._bit_generator
+    bitgenerator = default_rng(random_state)._bit_generator
     initialize_rng(bitgenerator, &rng)
 
     with bitgenerator.lock, nogil:
@@ -169,14 +216,15 @@ def structured_precision_mvnorm(
     double[:,::1] phi,
     double[:,::1] omega,
     bint mean_structured=False,
-    int a_type=0,
-    int o_type=0,
+    str a_type="regular",
+    str o_type="regular",
     double[:] out=None,
     random_state=None
 ):
     """
     structured_precision_mvnorm(mean, a, phi, omega, mean_structured=False,
-                                a_type=0, o_type=0, out=None, random_state=None)
+                                a_type="regular", o_type="regular",
+                                out=None, random_state=None)
 
     Sample from a MVN with a structured precision matrix :math:`\Lambda`
     .. math::
@@ -197,9 +245,9 @@ def structured_precision_mvnorm(
         such than ``mean = (precision)^-1 * phi^T * omega * t``. If this
         is set to True, then the `mean` parameter is assumed to contain the
         array ``t``.
-    a_type : {0, 1, 2}, optional, default=0
+    a_type : {"regular", "diagonal", "identity"}, optional, default="regular"
         Whether `a` ia a normal, diagonal or identity matrix.
-    o_type : {0, 1, 2}, optional, default=0
+    o_type : {"regular", "diagonal", "identity"}, optional, default="regular"
         Whether `omega` ia a normal, diagonal or identity matrix.
     out : 1d array, optional, default=None
         An array of the same shape as `mean` to store the samples. If not
@@ -228,7 +276,6 @@ def structured_precision_mvnorm(
         the algorithm could not successfully generate the samples.
 
     """
-    cdef BitGenerator bitgenerator
     cdef rng_t rng
     cdef sp_config_t config
     cdef int info
@@ -238,8 +285,8 @@ def structured_precision_mvnorm(
         raise ValueError('`omega` and `a` both need to be square matrices')
     elif (phi.shape[0] != omega.shape[0]) or (phi.shape[1] != a.shape[0]):
         raise ValueError('Shapes of `phi`, `omega` and `a` are not consistent')
-    elif not {a_type, o_type}.issubset(_VALID_MATRIX_TYPES):
-        raise ValueError(f"`a_type` & `o_type` must be one of {_VALID_MATRIX_TYPES}")
+    elif not {a_type, o_type}.issubset(MAT_TYPE):
+        raise ValueError(f"`a_type` & `o_type` must be one of {set(MAT_TYPE)}")
     elif has_out and out.shape[0] != mean.shape[0]:
         raise ValueError("`out` must have the same size as the mean array.")
     elif not has_out:
@@ -248,9 +295,9 @@ def structured_precision_mvnorm(
 
     init_sp_config(&config, phi.shape[0], phi.shape[1], &mean[0], &a[0, 0],
                    &phi[0, 0], &omega[0, 0], mean_structured,
-                   <mat_type>a_type, <mat_type>o_type)
+                   MAT_TYPE[a_type], MAT_TYPE[o_type])
 
-    bitgenerator = np.random.default_rng(random_state)._bit_generator
+    bitgenerator = default_rng(random_state)._bit_generator
     initialize_rng(bitgenerator, &rng)
 
     with bitgenerator.lock, nogil: