IntelPython · oleksandr-pavlyk · Sep 5, 2019 · Sep 5, 2019 · Sep 5, 2019 · Sep 5, 2019
@@ -1,4 +1,4 @@
-{% set version = "1.0.2" %}
+{% set version = "1.1.0" %}
 {% set buildnumber = 0 %}
 
 ### If you change the iccver here, you must also set the path correctly in build.sh / bld.bat!!!
@@ -31,12 +31,12 @@ requirements:
     host:
       - python
       - setuptools
-      - mkl-devel
+      - mkl-devel >=2019.4
       - cython
       - numpy
     run:
       - python
-      - mkl
+      - mkl >=2019.4
       - {{ pin_compatible('numpy') }}
 
 test:

@@ -27,6 +27,7 @@
 from __future__ import division, absolute_import, print_function
 
 from .mklrand import *
+from ._version import __version__
 
 from numpy.testing.nosetester import _numpy_tester
 test = _numpy_tester().test

@@ -0,0 +1,2 @@
+__version__ = '1.1.0'
+
@@ -71,6 +71,7 @@ cdef extern from "randomkit.h":
         MRG32K3A = 6
         MCG59 = 7
         PHILOX4X32X10 = 8
+        NONDETERM = 9
 
     void irk_fill(void *buffer, size_t size, irk_state *state) nogil
 
@@ -130,6 +131,7 @@ cdef extern from "mkl_distributions.h":
     void irk_geometric_vec(irk_state *state, npy_intp len, int *res, double p) nogil
     void irk_negbinomial_vec(irk_state *state, npy_intp len, int *res, double a, double p) nogil
     void irk_binomial_vec(irk_state *state, npy_intp len, int *res, int n, double p) nogil
+    void irk_multinomial_vec(irk_state *state, npy_intp len, int *res, int n, int d, double *pvec) nogil
     void irk_hypergeometric_vec(irk_state *state, npy_intp len, int *res, int ls, int ss, int ms) nogil
 
     void irk_poisson_vec_PTPE(irk_state *state, npy_intp len, int *res, double lam) nogil
@@ -867,7 +869,10 @@ _brng_dict = {
     'R250' : R250,
     'MRG32K3A' : MRG32K3A,
     'MCG59' : MCG59,
-    'PHILOX4X32X10' : PHILOX4X32X10
+    'PHILOX4X32X10' : PHILOX4X32X10,
+    'NONDETERM' : NONDETERM,
+    'NONDETERMINISTIC' : NONDETERM,
+    'NON_DETERMINISTIC' : NONDETERM,
 }
 
 _brng_dict_stream_max = {
@@ -879,7 +884,8 @@ _brng_dict_stream_max = {
     R250: 1,
     MRG32K3A: 1,
     MCG59: 1,
-    PHILOX4X32X10: 1
+    PHILOX4X32X10: 1,
+    NONDETERM: 1,
 }
 
 def _default_fallback_brng_token_(brng):
@@ -967,8 +973,8 @@ cdef class RandomState:
         If `seed` is ``None``, then `RandomState` will try to read data from
         ``/dev/urandom`` (or the Windows analogue) if available or seed from
         the clock otherwise.
-    brng : {'MT19937', 'SFMT19937', 'MT2203', 'R250',
-            'WH', 'MCG31', 'MCG59', 'MRG32K3A', 'PHILOX4X32X10'}, optional
+    brng : {'MT19937', 'SFMT19937', 'MT2203', 'R250', 'WH', 'MCG31',
+            'MCG59', 'MRG32K3A', 'PHILOX4X32X10', 'NONDETERM'}, optional
         Basic pseudo-random number generation algorithms, provided by
         Intel MKL. The default choice is 'MT19937' - the Mersenne Twister.
 
@@ -1016,8 +1022,8 @@ cdef class RandomState:
         seed : int or array_like, optional
             Seed for `RandomState`.
             Must be convertible to 32 bit unsigned integers.
-        brng : {'MT19937', 'SFMT19937', 'MT2203', 'R250',
-                'WH', 'MCG31', 'MCG59', 'MRG32K3A', 'PHILOX4X32X10'}, optional
+        brng : {'MT19937', 'SFMT19937', 'MT2203', 'R250', 'WH', 'MCG31',
+                'MCG59', 'MRG32K3A', 'PHILOX4X32X10', 'NONDETERM'}, optional
             Basic pseudo-random number generation algorithms, provided by
             Intel MKL. The default choice is 'MT19937' - the Mersenne Twister.
 
@@ -5539,27 +5545,13 @@ cdef class RandomState:
             raise ValueError("sum(pvals[:-1]) > 1.0")
 
         shape = _shape_from_size(size, d)
-
         multin = np.zeros(shape, np.int32)
+
         mnarr = <ndarray>multin
         mnix = <int*>PyArray_DATA(mnarr)
         sz = PyArray_SIZE(mnarr)
-        with self.lock, nogil, cython.cdivision(True):
-            i = 0
-            while i < sz:
-                Sum = 1.0
-                dn = n
-                for j from 0 <= j < d-1:
-#                    mnix[i+j] = irk_binomial(self.internal_state, dn, pix[j]/Sum)
-                    irk_binomial_vec(self.internal_state, 1, mnix + (i+j), dn, pix[j]/Sum)
-                    dn = dn - mnix[i+j]
-                    if dn <= 0:
-                        break
-                    Sum = Sum - pix[j]
-                if dn > 0:
-                    mnix[i+d-1] = dn
-
-                i = i + d
+
+        irk_multinomial_vec(self.internal_state, sz // d, mnix, n, d, pix)
 
         return multin
 

@@ -1136,6 +1136,31 @@ irk_binomial_vec(irk_state *state, npy_intp len, int *res, const int n, const do
     }
 }
 
+void
+irk_multinomial_vec(irk_state *state, npy_intp len, int *res, const int n, const int k, const double *pvec)
+{
+    int err;
+
+    if(len < 1)
+        return;
+
+    if (n==0) {
+        memset(res, 0, len*k*sizeof(int));
+    }
+    else {
+        while (len > MKL_INT_MAX) {
+            err = viRngMultinomial(VSL_RNG_METHOD_MULTINOMIAL_MULTPOISSON, state->stream, MKL_INT_MAX, res, n, k, pvec);
+            assert(err == VSL_STATUS_OK);
+            res += k*MKL_INT_MAX;
+            len -= k*MKL_INT_MAX;
+        }
+
+        err = viRngMultinomial(VSL_RNG_METHOD_MULTINOMIAL_MULTPOISSON, state->stream, len, res, n, k, pvec);
+        assert(err == VSL_STATUS_OK);
+    }
+}
+
+
 void
 irk_geometric_vec(irk_state *state, npy_intp len, int *res, const double p)
 {
@@ -1953,184 +1978,3 @@ static double irk_double(irk_state *state) {
 
     return res;
 }
-
-#ifdef COMPILE_UNUSED_CODE
-
-#ifndef min
-#define min(x,y) (((x)<(y))?(x):(y))
-#endif
-
-static long irk_binomial_btpe(irk_state *state, long n, double p)
-{
-    double r,q,fm,p1,xm,xl,xr,c,laml,lamr,p2,p3,p4;
-    double a,u,v,s,F,rho,t,A,nrq,x1,x2,f1,f2,z,z2,w,w2,x;
-    long m,y,k,i;
-
-    r = min(p, 1.0-p);
-    q = 1.0 - r;
-    fm = n*r+r;
-    m = (long)floor(fm);
-    p1 = floor(2.195*sqrt(n*r*q)-4.6*q) + 0.5;
-    xm = m + 0.5;
-    xl = xm - p1;
-    xr = xm + p1;
-    c = 0.134 + 20.5/(15.3 + m);
-    a = (fm - xl)/(fm-xl*r);
-    laml = a*(1.0 + a/2.0);
-    a = (xr - fm)/(xr*q);
-    lamr = a*(1.0 + a/2.0);
-    p2 = p1*(1.0 + 2.0*c);
-    p3 = p2 + c/laml;
-    p4 = p3 + c/lamr;
-
-  /* sigh ... */
-  Step10:
-    nrq = n*r*q;
-    u = irk_double(state)*p4;
-    v = irk_double(state);
-    if (u > p1) goto Step20;
-    y = (long)floor(xm - p1*v + u);
-    goto Step60;
-
-  Step20:
-    if (u > p2) goto Step30;
-    x = xl + (u - p1)/c;
-    v = v*c + 1.0 - fabs(m - x + 0.5)/p1;
-    if (v > 1.0) goto Step10;
-    y = (long)floor(x);
-    goto Step50;
-
-  Step30:
-    if (u > p3) goto Step40;
-    y = (long)floor(xl + log(v)/laml);
-    if (y < 0) goto Step10;
-    v = v*(u-p2)*laml;
-    goto Step50;
-
-  Step40:
-    y = (long)floor(xr - log(v)/lamr);
-    if (y > n) goto Step10;
-    v = v*(u-p3)*lamr;
-
-  Step50:
-    k = labs(y - m);
-    if ((k > 20) && (k < ((nrq)/2.0 - 1))) goto Step52;
-
-    s = r/q;
-    a = s*(n+1);
-    F = 1.0;
-    if (m < y)
-    {
-        for (i=m+1; i<=y; i++)
-        {
-            F *= (a/i - s);
-        }
-    }
-    else if (m > y)
-    {
-        for (i=y+1; i<=m; i++)
-        {
-            F /= (a/i - s);
-        }
-    }
-    if (v > F) goto Step10;
-    goto Step60;
-
-    Step52:
-    rho = (k/(nrq))*((k*(k/3.0 + 0.625) + 0.16666666666666666)/nrq + 0.5);
-    t = -k*k/(2*nrq);
-    A = log(v);
-    if (A < (t - rho)) goto Step60;
-    if (A > (t + rho)) goto Step10;
-
-    x1 = y+1;
-    f1 = m+1;
-    z = n+1-m;
-    w = n-y+1;
-    x2 = x1*x1;
-    f2 = f1*f1;
-    z2 = z*z;
-    w2 = w*w;
-    if (A > (xm*log(f1/x1)
-           + (n-m+0.5)*log(z/w)
-           + (y-m)*log(w*r/(x1*q))
-           + (13680.-(462.-(132.-(99.-140./f2)/f2)/f2)/f2)/f1/166320.
-           + (13680.-(462.-(132.-(99.-140./z2)/z2)/z2)/z2)/z/166320.
-           + (13680.-(462.-(132.-(99.-140./x2)/x2)/x2)/x2)/x1/166320.
-           + (13680.-(462.-(132.-(99.-140./w2)/w2)/w2)/w2)/w/166320.))
-    {
-        goto Step10;
-    }
-
-  Step60:
-    if (p > 0.5)
-    {
-        y = n - y;
-    }
-
-    return y;
-}
-
-static long irk_binomial_inversion(irk_state *state, long n, double p)
-{
-    double q, qn, np, px, U;
-    long X, bound;
-
-    q = 1.0 - p;
-    qn = exp(n * log(q));
-    np = n*p;
-    bound = min(n, np + 10.0*sqrt(np*q + 1));
-
-    X = 0;
-    px = qn;
-    U = irk_double(state);
-    while (U > px)
-    {
-        X++;
-        if (X > bound)
-        {
-            X = 0;
-            px = qn;
-            U = irk_double(state);
-        } else
-        {
-            U -= px;
-            px  = ((n-X+1) * p * px)/(X*q);
-        }
-    }
-    return X;
-}
-
-#undef min
-
-static long irk_binomial(irk_state *state, long n, double p)
-{
-    double q;
-
-    if (p <= 0.5)
-    {
-        if (p*n <= 30.0)
-        {
-            return irk_binomial_inversion(state, n, p);
-        }
-        else
-        {
-            return irk_binomial_btpe(state, n, p);
-        }
-    }
-    else
-    {
-        q = 1.0-p;
-        if (q*n <= 30.0)
-        {
-            return n - irk_binomial_inversion(state, n, q);
-        }
-        else
-        {
-            return n - irk_binomial_btpe(state, n, q);
-        }
-    }
-
-}
-
-#endif
@@ -88,6 +88,9 @@ extern void irk_noncentral_f_vec(irk_state *state, npy_intp len, double *res, do
 extern void irk_triangular_vec(irk_state *state, npy_intp len, double *res, double left, double mode, double right);
 
 extern void irk_binomial_vec(irk_state *state, npy_intp len, int *res, const int n, const double p);
+
+extern void irk_multinomial_vec(irk_state *state, npy_intp len, int *res, const int n, const int k, const double* pvec);
+
 extern void irk_geometric_vec(irk_state *state, npy_intp len, int *res, const double p);
 extern void irk_negbinomial_vec(irk_state *state, npy_intp len, int *res, const double a, const double p);
 extern void irk_hypergeometric_vec(irk_state *state, npy_intp len, int *res, const int ls, const int ss, const int ms);

@@ -118,7 +118,8 @@ const MKL_INT brng_list[BRNG_KINDS] = {
     VSL_BRNG_R250,
     VSL_BRNG_MRG32K3A,
     VSL_BRNG_MCG59,
-    VSL_BRNG_PHILOX4X32X10
+    VSL_BRNG_PHILOX4X32X10,
+    VSL_BRNG_NONDETERM
 };
 
 /* Mersenne-Twister 2203 algorithm and Wichmann-Hill algorithm 

@@ -45,7 +45,7 @@ typedef enum {
 } irk_error;
 
 /* if changing this, also adjust brng_list[BRNG_KINDS] in randomkit.c */
-#define BRNG_KINDS 9
+#define BRNG_KINDS 10
 
 typedef enum {
     MT19937       = 0,
@@ -56,7 +56,8 @@ typedef enum {
     R250          = 5,
     MRG32K3A      = 6,
     MCG59         = 7,
-    PHILOX4X32X10 = 8
+    PHILOX4X32X10 = 8,
+    NONDETERM     = 9
 } irk_brng_t;