Merge pull request #18 from zoj613/ziggurat

pyhtnorm/_htnorm.pyx

@@ -57,9 +57,7 @@ cdef class HTNGenerator:
         The name of the generator to use for random number generation. The
         value needs to be one of {'pcg', 'xrs'}, where 'pcg' is PCG64 and 'xrs'
         is the Xoroshiro128plus bit generator.
-
-    Methods
-    -------
+Methods -------
     hyperplane_truncated_mvnorm(mean, cov, g, r, diag=False, out=False)
     structured_precision_mvnorm(mean, a, phi, omega, mean_structured=False,
                                 a_type=0, o_type=0, out=None)
@@ -257,7 +255,7 @@ cdef class HTNGenerator:
         if (omega.shape[0] != omega.shape[1]) or (a.shape[0] != a.shape[1]):
             raise ValueError('`omega` and `a` both need to be square matrices')
 
-        if (phi.shape[1] != omega.shape[0]) or (phi.shape[0] != a.shape[1]):
+        if (phi.shape[0] != omega.shape[0]) or (phi.shape[1] != a.shape[0]):
             raise ValueError(
                 "Shapes of `phi`, `omega` and `a` are not consistent"
             )

src/LICENSE-ZIGGURAT.txt

@@ -0,0 +1,58 @@
+Copyright (c) 2005-2017, NumPy Developers.
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+
+* Redistributions of source code must retain the above copyright
+   notice, this list of conditions and the following disclaimer.
+
+* 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.
+
+* Neither the name of the NumPy Developers nor the names of any
+   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.
+
+
+
+The ziggurat methods were derived from Julia.
+
+Copyright (c) 2009-2019: Jeff Bezanson, Stefan Karpinski, Viral B. Shah,
+and other contributors:
+
+https://github.com/JuliaLang/julia/contributors
+
+Permission is hereby granted, free of charge, to any person obtaining
+a copy of this software and associated documentation files (the
+"Software"), to deal in the Software without restriction, including
+without limitation the rights to use, copy, modify, merge, publish,
+distribute, sublicense, and/or sell copies of the Software, and to
+permit persons to whom the Software is furnished to do so, subject to
+the following conditions:
+
+The above copyright notice and this permission notice shall be
+included in all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
+LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
+OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

src/dist.c

@@ -7,41 +7,53 @@
 
 #include "blas.h"
 #include "dist.h"
+#include "zig_constants.h"
 
 #define std_normal_rand_fill(rng_t, arr_size, arr) \
     for (size_t inc = (arr_size); inc--;) (arr)[inc] = std_normal_rand((rng_t))
 
 extern ALWAYS_INLINE(mvn_output_t*) mvn_output_new(size_t nrow, type_t factor_type);
 extern ALWAYS_INLINE(void) mvn_output_free(mvn_output_t* a);
 
-/* Generate a sample from the standard normal distribution using the 
- * Marsaglia-Polar method.
- *
- * TODO: Think about using a faster one, maybe the Ziggurat method?*/
+// Generate a sample from the standard normal distribution using the Ziggurat method.
+// This uses numpy's implementation of the algorithm.
 static ALWAYS_INLINE(double)
 std_normal_rand(rng_t* rng)
 {
-    double s, u, v, z;
-    static double x, y;
-    static bool cached = false;
+    uint64_t r, rabs;
+    int sign, idx;
+    double x, xx, yy;
 
-    if (cached) {
-        cached = false;
-        return y;
-    }
+    for (;;) {
+        r = rng->next_int(rng->base);
+        idx = r & 0xff;
+        r >>= 8;
+        sign = r & 0x1;
+        rabs = (r >> 1) & 0x000fffffffffffff;
+        x = rabs * wi_double[idx];
 
-    do {
-        u = 2.0 * rng->next_double(rng->base) - 1;
-        v = 2.0 * rng->next_double(rng->base) - 1;
-        s = u * u + v * v;
-    } while (!(s < 1.0));
+        if (sign & 0x1)
+            x = -x;
+        if (rabs < ki_double[idx])
+            return x; /* 99.3% of the time return here */
 
-    z = sqrt(-2 * log(s) / s);
-    y = v * z;
-    x = u * z;
-    cached = true;
-    return x;
-}
+        if (idx == 0) {
+            // use tail sampling method.
+            for (;;) {
+            xx = -ziggurat_nor_inv_r * log(1.0 - rng->next_double(rng->base));
+            yy = -log(1.0 - rng->next_double(rng->base));
+            if (yy + yy > xx * xx)
+                return ((rabs >> 8) & 0x1) ? -(ziggurat_nor_r + xx) : ziggurat_nor_r + xx;
+            }
+        }
+        else {
+            // try and get a sample at the wedge.
+            if (((fi_double[idx - 1] - fi_double[idx]) *
+                rng->next_double(rng->base) + fi_double[idx]) < exp(-0.5 * x * x))
+            return x;
+        }
+    }
+} 
 
 
 int