diff --git a/DESCRIPTION b/DESCRIPTION index bc760fc..e830028 100644 --- a/DESCRIPTION +++ b/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 diff --git a/Makefile b/Makefile index a2d8517..ff544d9 100644 --- a/Makefile +++ b/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 diff --git a/R/htnorm.r b/R/htnorm.r index 2f7c53f..558fa2d 100644 --- a/R/htnorm.r +++ b/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( diff --git a/README.md b/README.md index 556974a..8588b20 100644 --- a/README.md +++ b/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 diff --git a/build-wheels.sh b/build-wheels.sh index 74110e4..463e722 100755 --- a/build-wheels.sh +++ b/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 } diff --git a/build.py b/build.py index 7263e61..e6bfd64 100644 --- a/build.py +++ b/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)], diff --git a/include/htnorm.h b/include/htnorm.h index 11f0cd6..f9d7f37 100644 --- a/include/htnorm.h +++ b/include/htnorm.h @@ -28,7 +28,7 @@ #include #include -#include "rng.h" +#include "htnorm_rng.h" typedef enum {NORMAL, DIAGONAL, IDENTITY} type_t; diff --git a/include/rng.h b/include/htnorm_rng.h similarity index 100% rename from include/rng.h rename to include/htnorm_rng.h diff --git a/man/HTNGenerator.Rd b/man/HTNGenerator.Rd index 0e11fdb..7c36337 100644 --- a/man/HTNGenerator.Rd +++ b/man/HTNGenerator.Rd @@ -41,7 +41,7 @@ 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 = 1, out = out) +rng$structured_precision_mvnorm(mean, a, phi, omega, a_type = "diagonal", out = out) } \references{ Cong, Yulai; Chen, Bo; Zhou, Mingyuan. Fast Simulation of diff --git a/pyhtnorm/__init__.pxd b/pyhtnorm/__init__.pxd deleted file mode 100644 index e31493a..0000000 --- a/pyhtnorm/__init__.pxd +++ /dev/null @@ -1 +0,0 @@ -from pyhtnorm.c_htnorm cimport * diff --git a/pyhtnorm/_htnorm.pyx b/pyhtnorm/_htnorm.pyx index cc0143e..822043e 100644 --- a/pyhtnorm/_htnorm.pyx +++ b/pyhtnorm/_htnorm.pyx @@ -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 = PyCapsule_GetPointer(bitgenerator.capsule, "BitGenerator") + bitgen = 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, - a_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: diff --git a/pyhtnorm/c_htnorm.pxd b/pyhtnorm/c_htnorm.pxd deleted file mode 100644 index 815cce4..0000000 --- a/pyhtnorm/c_htnorm.pxd +++ /dev/null @@ -1,55 +0,0 @@ -# cython: language_level=3 -""" -Copyright (c) 2020-2021, Zolisa Bleki - -SPDX-License-Identifier: BSD-3-Clause */ -""" -from libc.stdint cimport uint64_t - - -cdef extern from "../include/rng.h": - ctypedef struct rng_t: - void* base - uint64_t (*next_uint64)(void* state) nogil - double (*next_double)(void* state) nogil - - -cdef extern from "../include/htnorm.h": - - 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 is_diag "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) nogil - - 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) nogil - - int htn_hyperplane_truncated_mvn(rng_t* rng, const ht_config_t* conf, double* out) nogil - - int htn_structured_precision_mvn(rng_t* rng, const sp_config_t* conf, double* out) nogil diff --git a/pyproject.toml b/pyproject.toml index 1992f50..617070d 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [tool.poetry] name = "pyhtnorm" -version = "1.0.0" +version = "2.0.0" description = "Fast Simulation of Hyperplane-Truncated Multivatiate Normal Distributions" authors = ["Zolisa Bleki"] license = "BSD-3-Clause" @@ -20,17 +20,15 @@ include = [ {path = "include", format = "sdist"}, {path = "src/*.c", format = "sdist"}, {path = "src/*.h", format = "sdist"}, - {path = "src/*.txt", format = "sdist"}, - - {path = "pyhtnorm/*.pxd"} + {path = "third_party_licenses"}, ] [tool.poetry.build] script = "build.py" [tool.poetry.dependencies] -python = "^3.6.1" -numpy = "^1.17.0" +python = "^3.7.0" +numpy = "^1.18.1" [tool.poetry.dev-dependencies] Cython = "^0.29.20" @@ -38,5 +36,5 @@ pytest = "*" toml = "^0.10.2" [build-system] -requires = ["poetry-core>=1.0.0a9", "wheel", "setuptools", "numpy>=1.17"] +requires = ["poetry-core>=1.0.0a9", "wheel", "setuptools", "numpy==1.18.1"] build-backend = "poetry.core.masonry.api" diff --git a/src/Makevars b/src/Makevars index 50f8b58..0304eb9 100644 --- a/src/Makevars +++ b/src/Makevars @@ -1,16 +1,16 @@ PKG_CPPFLAGS = -I../include -DHTNORM_COLMAJOR PKG_LIBS = $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) -OBJS = htnorm.o rng.o dist.o +OBJS = htnorm_distributions.o htnorm.o htnorm_rng.o all: @$(MAKE) $(SHLIB) @rm -f *.mod src/*.o *.o -$(SHLIB): r_wrapper.o +$(SHLIB): htnorm_r_wrapper.o -r_wrapper.o: $(OBJS) -htnorm.o: dist.o +htnorm_r_wrapper.o: $(OBJS) +htnorm.o: htnorm_distributions.o clean: @rm -rf *.o *.mod *.d *.rc *.so *.dylib *.dll *.a *.lib $(SHLIB) $(OBJECTS) diff --git a/src/htnorm.c b/src/htnorm.c index 78417db..19b0f0f 100644 --- a/src/htnorm.c +++ b/src/htnorm.c @@ -8,14 +8,13 @@ * macro for several lines per function because matrix manipulation with BLAS * / LAPACK in C heavily depends on the memory layout of the input. * */ -#include #include -#include "blas.h" -#include "dist.h" +#include "htnorm_blas.h" +#include "htnorm_distributions.h" -inline void +ALWAYS_INLINE(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, bool diag) { @@ -29,7 +28,7 @@ init_ht_config(ht_config_t* conf, size_t gnrow, size_t gncol, const double* mean } -inline void +ALWAYS_INLINE(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, bool struct_mean, type_t a_id, type_t o_id) diff --git a/src/always_inline.h b/src/htnorm_always_inline.h similarity index 76% rename from src/always_inline.h rename to src/htnorm_always_inline.h index f156b96..5d05797 100644 --- a/src/always_inline.h +++ b/src/htnorm_always_inline.h @@ -1,3 +1,10 @@ +/* Copyright (c) 2020, Zolisa Bleki + * + * SPDX-License-Identifier: BSD-3-Clause + * + * This file is derived from an example in the documentation of clang. + */ +#pragma once #ifndef HTNORM_ALWAYS_INLINE_H #define HTNORM_ALWAYS_INLINE_H diff --git a/src/blas.h b/src/htnorm_blas.h similarity index 100% rename from src/blas.h rename to src/htnorm_blas.h diff --git a/src/dist.c b/src/htnorm_distributions.c similarity index 96% rename from src/dist.c rename to src/htnorm_distributions.c index 5c05934..44782ae 100644 --- a/src/dist.c +++ b/src/htnorm_distributions.c @@ -5,9 +5,9 @@ #include #include -#include "blas.h" -#include "dist.h" -#include "zig_constants.h" +#include "htnorm_blas.h" +#include "htnorm_distributions.h" +#include "htnorm_ziggurat_constants.h" #define std_normal_rand_fill(rng, arr_size, arr) \ for (int ii = (arr_size); ii--;) (arr)[ii] = std_normal_rand((rng)) @@ -69,16 +69,17 @@ mvn_rand_cov(rng_t* rng, const double* mean, const double* cov, int nrow, if (factor == NULL) return HTNORM_ALLOC_ERROR; - // do cholesky factorization. int info; static const int incx = 1; static const double one = 1; memcpy(factor, cov, nrow * nrow * sizeof(*factor)); + // do cholesky factorization. POTRF(nrow, factor, nrow, info); if (!info) { std_normal_rand_fill(rng, nrow, out); // triangular matrix-vector product. L * z. TRMV(nrow, factor, nrow, out, incx); + // mean + L * z AXPY(nrow, one, mean, out); } diff --git a/src/dist.h b/src/htnorm_distributions.h similarity index 97% rename from src/dist.h rename to src/htnorm_distributions.h index f8b4e56..f731b3f 100644 --- a/src/dist.h +++ b/src/htnorm_distributions.h @@ -9,12 +9,11 @@ #ifndef HTNORM_DIST_H #define HTNORM_DIST_H -#include #include -#include "../include/rng.h" +#include "../include/htnorm_rng.h" #include "../include/htnorm.h" -#include "always_inline.h" +#include "htnorm_always_inline.h" // error number for failed memory allocation throughout the library #define HTNORM_ALLOC_ERROR -100 diff --git a/src/pcg32_minimal.h b/src/htnorm_pcg32_minimal.h similarity index 96% rename from src/pcg32_minimal.h rename to src/htnorm_pcg32_minimal.h index aba4634..e3a2f46 100644 --- a/src/pcg32_minimal.h +++ b/src/htnorm_pcg32_minimal.h @@ -8,7 +8,7 @@ #include -#include "always_inline.h" +#include "htnorm_always_inline.h" typedef struct pcg32_rand { uint64_t state; diff --git a/src/pcg64.h b/src/htnorm_pcg64.h similarity index 94% rename from src/pcg64.h rename to src/htnorm_pcg64.h index e53900e..1c4323b 100644 --- a/src/pcg64.h +++ b/src/htnorm_pcg64.h @@ -4,8 +4,8 @@ #ifndef HTNORM_PCG64_H #define HTNORM_PCG64_H -#include "pcg32_minimal.h" -#include "splitmax64.h" +#include "htnorm_pcg32_minimal.h" +#include "htnorm_splitmix64.h" typedef struct {pcg32_random_t gen[2];} pcg64_random_t; diff --git a/src/r_wrapper.c b/src/htnorm_r_wrapper.c similarity index 99% rename from src/r_wrapper.c rename to src/htnorm_r_wrapper.c index 129c107..2f72a0c 100644 --- a/src/r_wrapper.c +++ b/src/htnorm_r_wrapper.c @@ -4,7 +4,7 @@ #include #include -#include "../include/htnorm.h" +#include "htnorm.h" enum RNG_TYPE {Xoroshiro128p, PCG64}; diff --git a/src/rng.c b/src/htnorm_rng.c similarity index 93% rename from src/rng.c rename to src/htnorm_rng.c index 9794371..f55cb46 100644 --- a/src/rng.c +++ b/src/htnorm_rng.c @@ -4,10 +4,9 @@ #include #include -#include "pcg64.h" -#include "xoroshiro128p.h" - -#include "../include/rng.h" +#include "htnorm_pcg64.h" +#include "htnorm_xoroshiro128p.h" +#include "../include/htnorm_rng.h" ALWAYS_INLINE(void) diff --git a/src/splitmax64.h b/src/htnorm_splitmix64.h similarity index 95% rename from src/splitmax64.h rename to src/htnorm_splitmix64.h index 430a8ab..f6875d1 100644 --- a/src/splitmax64.h +++ b/src/htnorm_splitmix64.h @@ -11,7 +11,7 @@ See . */ #include -#include "always_inline.h" +#include "htnorm_always_inline.h" // generate a postive integer using the splitmix64 bit generator static ALWAYS_INLINE(uint64_t) diff --git a/src/xoroshiro128p.h b/src/htnorm_xoroshiro128p.h similarity index 96% rename from src/xoroshiro128p.h rename to src/htnorm_xoroshiro128p.h index 89929f2..3adcf3d 100644 --- a/src/xoroshiro128p.h +++ b/src/htnorm_xoroshiro128p.h @@ -9,7 +9,7 @@ #include -#include "always_inline.h" +#include "htnorm_always_inline.h" typedef struct {uint64_t s[2];} xrs128p_random_t; diff --git a/src/zig_constants.h b/src/htnorm_ziggurat_constants.h similarity index 100% rename from src/zig_constants.h rename to src/htnorm_ziggurat_constants.h diff --git a/tests/test_pyhtnorm.py b/tests/test_pyhtnorm.py index 602096d..a5029d1 100644 --- a/tests/test_pyhtnorm.py +++ b/tests/test_pyhtnorm.py @@ -82,7 +82,7 @@ def test_structured_mvn(structured_mvn_data): structured_precision_mvnorm(mean, a2, phi, omega2) # raise error if invalid matrix structure is specified with pytest.raises(ValueError): - structured_precision_mvnorm(mean, a, phi, omega, a_type=-1000) + structured_precision_mvnorm(mean, a, phi, omega, a_type="wrong type") # rest for non-numerical values a2 = a.copy() a2[0] = np.nan @@ -94,7 +94,7 @@ def test_structured_mvn(structured_mvn_data): arr1 = structured_precision_mvnorm(mean, a, phi, omega, random_state=rng) rng = np.random.default_rng(10) arr2 = structured_precision_mvnorm( - mean, a, phi, omega, o_type=1, a_type=1, random_state=rng + mean, a, phi, omega, o_type="diagonal", a_type="diagonal", random_state=rng ) assert np.allclose(arr1, arr2) # test results of passing output array through the `out` parameter diff --git a/tests/testthat/test_rhtnorm.r b/tests/testthat/test_rhtnorm.r index 32d12ec..aefe6a5 100644 --- a/tests/testthat/test_rhtnorm.r +++ b/tests/testthat/test_rhtnorm.r @@ -86,7 +86,7 @@ test_that("structured precision normal", { # test consistency of output when `a_type` or `o_type` is given gen1 <- HTNGenerator(10)$structured_precision_mvnorm(mean, a, phi, omega) gen2 <- HTNGenerator(10)$structured_precision_mvnorm( - mean, a, phi, omega, a_type = 1, o_type = 1 + mean, a, phi, omega, a_type = "diagonal", o_type = "diagonal" ) expect_equal(gen1, gen2) # test results of passing output array through the `out` parameter diff --git a/third_party_licenses/LICENSE-OPENBLAS.txt b/third_party_licenses/LICENSE-OPENBLAS.txt new file mode 100644 index 0000000..d8a6db2 --- /dev/null +++ b/third_party_licenses/LICENSE-OPENBLAS.txt @@ -0,0 +1,32 @@ +The wheel distribution of pyhtnorm contains OpenBLAS binaries. + + +Copyright (c) 2011-2014, The OpenBLAS Project +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are +met: + + 1. Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + 2. Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in + the documentation and/or other materials provided with the + distribution. + 3. Neither the name of the OpenBLAS project nor the names of + its contributors may be used to endorse or promote products + derived from this software without specific prior written + permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE +LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, +OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE +USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/src/LICENSE-PCG32.txt b/third_party_licenses/LICENSE-PCG32.txt similarity index 98% rename from src/LICENSE-PCG32.txt rename to third_party_licenses/LICENSE-PCG32.txt index 8f71f43..04619f9 100644 --- a/src/LICENSE-PCG32.txt +++ b/third_party_licenses/LICENSE-PCG32.txt @@ -1,3 +1,7 @@ +The code for the PCG32 bitgenerator at src/htnorm_pcg32_minimal.h was obtained +from M. E. O'Neill's pcg-random package. Website of original code: pcg-random.org + + Apache License Version 2.0, January 2004 http://www.apache.org/licenses/ diff --git a/src/LICENSE-ZIGGURAT.txt b/third_party_licenses/LICENSE-ZIGGURAT.txt similarity index 93% rename from src/LICENSE-ZIGGURAT.txt rename to third_party_licenses/LICENSE-ZIGGURAT.txt index 43a7928..5dc0667 100644 --- a/src/LICENSE-ZIGGURAT.txt +++ b/third_party_licenses/LICENSE-ZIGGURAT.txt @@ -1,3 +1,8 @@ +The ziggurat constants at src/htnorm_ziggurat_constants.h and well as the +standard normal distribution random number generator implementation where +obtained from the numpy package. + + Copyright (c) 2005-2017, NumPy Developers. All rights reserved.