3rdparty/alglib/src/alglibinternal.cpp

 /*************************************************************************
-ALGLIB 3.14.0 (source code generated 2018-06-16)
+ALGLIB 3.15.0 (source code generated 2019-02-20)
 Copyright (c) Sergey Bochkanov (ALGLIB project).
 
 >>> SOURCE LICENSE >>>
@@ -438,6 +438,26 @@ ae_int_t getspline2dserializationcode(ae_state *_state)
 }
 
 
+ae_int_t getidwserializationcode(ae_state *_state)
+{
+    ae_int_t result;
+
+
+    result = 7;
+    return result;
+}
+
+
+ae_int_t getknnserializationcode(ae_state *_state)
+{
+    ae_int_t result;
+
+
+    result = 108;
+    return result;
+}
+
+
 #endif
 #if defined(AE_COMPILE_APSERV) || !defined(AE_PARTIAL_BUILD)
 
@@ -974,6 +994,65 @@ void rmatrixsetlengthatleast(/* Real    */ ae_matrix* x,
 }
 
 
+/*************************************************************************
+Grows X, i.e. changes its size in such a way that:
+a) contents is preserved
+b) new size is at least N
+c) new size can be larger than N, so subsequent grow() calls can return
+   without reallocation
+
+  -- ALGLIB --
+     Copyright 20.03.2009 by Bochkanov Sergey
+*************************************************************************/
+void bvectorgrowto(/* Boolean */ ae_vector* x,
+     ae_int_t n,
+     ae_state *_state)
+{
+    ae_frame _frame_block;
+    ae_vector oldx;
+    ae_int_t i;
+    ae_int_t n2;
+
+    ae_frame_make(_state, &_frame_block);
+    memset(&oldx, 0, sizeof(oldx));
+    ae_vector_init(&oldx, 0, DT_BOOL, _state, ae_true);
+
+    
+    /*
+     * Enough place
+     */
+    if( x->cnt>=n )
+    {
+        ae_frame_leave(_state);
+        return;
+    }
+    
+    /*
+     * Choose new size
+     */
+    n = ae_maxint(n, ae_round(1.8*x->cnt+1, _state), _state);
+    
+    /*
+     * Grow
+     */
+    n2 = x->cnt;
+    ae_swap_vectors(x, &oldx);
+    ae_vector_set_length(x, n, _state);
+    for(i=0; i<=n-1; i++)
+    {
+        if( i<n2 )
+        {
+            x->ptr.p_bool[i] = oldx.ptr.p_bool[i];
+        }
+        else
+        {
+            x->ptr.p_bool[i] = ae_false;
+        }
+    }
+    ae_frame_leave(_state);
+}
+
+
 /*************************************************************************
 Grows X, i.e. changes its size in such a way that:
 a) contents is preserved
@@ -1098,6 +1177,71 @@ void rmatrixgrowrowsto(/* Real    */ ae_matrix* a,
 }
 
 
+/*************************************************************************
+Grows X, i.e. appends cols in such a way that:
+a) contents is preserved
+b) new col count is at least N
+c) new col count can be larger than N, so subsequent grow() calls can return
+   without reallocation
+d) new matrix has at least MinRows row (if less than specified amount
+   of rows is present, new rows are added with undefined contents);
+   MinRows can be 0 or negative value = ignored
+
+  -- ALGLIB --
+     Copyright 20.03.2009 by Bochkanov Sergey
+*************************************************************************/
+void rmatrixgrowcolsto(/* Real    */ ae_matrix* a,
+     ae_int_t n,
+     ae_int_t minrows,
+     ae_state *_state)
+{
+    ae_frame _frame_block;
+    ae_matrix olda;
+    ae_int_t i;
+    ae_int_t j;
+    ae_int_t n2;
+    ae_int_t m;
+
+    ae_frame_make(_state, &_frame_block);
+    memset(&olda, 0, sizeof(olda));
+    ae_matrix_init(&olda, 0, 0, DT_REAL, _state, ae_true);
+
+    
+    /*
+     * Enough place?
+     */
+    if( a->cols>=n&&a->rows>=minrows )
+    {
+        ae_frame_leave(_state);
+        return;
+    }
+    
+    /*
+     * Sizes and metrics
+     */
+    if( a->cols<n )
+    {
+        n = ae_maxint(n, ae_round(1.8*a->cols+1, _state), _state);
+    }
+    n2 = ae_minint(a->cols, n, _state);
+    m = a->rows;
+    
+    /*
+     * Grow
+     */
+    ae_swap_matrices(a, &olda);
+    ae_matrix_set_length(a, ae_maxint(m, minrows, _state), n, _state);
+    for(i=0; i<=m-1; i++)
+    {
+        for(j=0; j<=n2-1; j++)
+        {
+            a->ptr.pp_double[i][j] = olda.ptr.pp_double[i][j];
+        }
+    }
+    ae_frame_leave(_state);
+}
+
+
 /*************************************************************************
 Grows X, i.e. changes its size in such a way that:
 a) contents is preserved
@@ -1370,6 +1514,7 @@ ae_bool isfinitevector(/* Real    */ ae_vector* x,
      ae_state *_state)
 {
     ae_int_t i;
+    double v;
     ae_bool result;
 
 
@@ -1384,15 +1529,12 @@ ae_bool isfinitevector(/* Real    */ ae_vector* x,
         result = ae_false;
         return result;
     }
+    v = (double)(0);
     for(i=0; i<=n-1; i++)
     {
-        if( !ae_isfinite(x->ptr.p_double[i], _state) )
-        {
-            result = ae_false;
-            return result;
+        v = 0.01*v+x->ptr.p_double[i];
     }
-    }
-    result = ae_true;
+    result = ae_isfinite(v, _state);
     return result;
 }
 
@@ -2174,6 +2316,27 @@ void countdown(ae_int_t* v, ae_state *_state)
 }
 
 
+/*************************************************************************
+This function returns +1 or -1 depending on sign of X.
+x=0 results in +1 being returned.
+*************************************************************************/
+double possign(double x, ae_state *_state)
+{
+    double result;
+
+
+    if( ae_fp_greater_eq(x,(double)(0)) )
+    {
+        result = (double)(1);
+    }
+    else
+    {
+        result = (double)(-1);
+    }
+    return result;
+}
+
+
 /*************************************************************************
 This function returns product of two real numbers. It is convenient when
 you have to perform typecast-and-product of two INTEGERS.
@@ -2419,6 +2582,57 @@ double rboundval(double x, double b1, double b2, ae_state *_state)
 }
 
 
+/*************************************************************************
+Returns number of non-zeros
+*************************************************************************/
+ae_int_t countnz1(/* Real    */ ae_vector* v,
+     ae_int_t n,
+     ae_state *_state)
+{
+    ae_int_t i;
+    ae_int_t result;
+
+
+    result = 0;
+    for(i=0; i<=n-1; i++)
+    {
+        if( !(v->ptr.p_double[i]==0) )
+        {
+            result = result+1;
+        }
+    }
+    return result;
+}
+
+
+/*************************************************************************
+Returns number of non-zeros
+*************************************************************************/
+ae_int_t countnz2(/* Real    */ ae_matrix* v,
+     ae_int_t m,
+     ae_int_t n,
+     ae_state *_state)
+{
+    ae_int_t i;
+    ae_int_t j;
+    ae_int_t result;
+
+
+    result = 0;
+    for(i=0; i<=m-1; i++)
+    {
+        for(j=0; j<=n-1; j++)
+        {
+            if( !(v->ptr.pp_double[i][j]==0) )
+            {
+                result = result+1;
+            }
+        }
+    }
+    return result;
+}
+
+
 /*************************************************************************
 Allocation of serializer: complex value
 *************************************************************************/
@@ -2702,6 +2916,28 @@ void unserializerealmatrix(ae_serializer* s,
 }
 
 
+/*************************************************************************
+Copy boolean array
+*************************************************************************/
+void copybooleanarray(/* Boolean */ ae_vector* src,
+     /* Boolean */ ae_vector* dst,
+     ae_state *_state)
+{
+    ae_int_t i;
+
+    ae_vector_clear(dst);
+
+    if( src->cnt>0 )
+    {
+        ae_vector_set_length(dst, src->cnt, _state);
+        for(i=0; i<=src->cnt-1; i++)
+        {
+            dst->ptr.p_bool[i] = src->ptr.p_bool[i];
+        }
+    }
+}
+
+
 /*************************************************************************
 Copy integer array
 *************************************************************************/
@@ -3000,6 +3236,25 @@ ae_int_t chunkscount(ae_int_t tasksize,
 }
 
 
+/*************************************************************************
+Returns maximum density for level 2 sparse/dense functions. Density values
+below one returned by this function are better to handle via sparse Level 2
+functionality.
+
+  -- ALGLIB routine --
+     10.01.2019
+     Bochkanov Sergey
+*************************************************************************/
+double sparselevel2density(ae_state *_state)
+{
+    double result;
+
+
+    result = 0.1;
+    return result;
+}
+
+
 /*************************************************************************
 Returns A-tile size for a matrix.
 
@@ -3877,6 +4132,12 @@ void tagsortmiddleir(/* Integer */ ae_vector* a,
     ae_int_t t;
     ae_int_t tmp;
     double tmpr;
+    ae_int_t p0;
+    ae_int_t p1;
+    ae_int_t at;
+    ae_int_t ak;
+    ae_int_t ak1;
+    double bt;
 
 
     
@@ -3891,76 +4152,167 @@ void tagsortmiddleir(/* Integer */ ae_vector* a,
     /*
      * General case, N>1: sort, update B
      */
-    i = 2;
-    do
+    for(i=2; i<=n; i++)
     {
         t = i;
         while(t!=1)
         {
             k = t/2;
-            if( a->ptr.p_int[offset+k-1]>=a->ptr.p_int[offset+t-1] )
+            p0 = offset+k-1;
+            p1 = offset+t-1;
+            ak = a->ptr.p_int[p0];
+            at = a->ptr.p_int[p1];
+            if( ak>=at )
             {
-                t = 1;
+                break;
             }
-            else
-            {
-                tmp = a->ptr.p_int[offset+k-1];
-                a->ptr.p_int[offset+k-1] = a->ptr.p_int[offset+t-1];
-                a->ptr.p_int[offset+t-1] = tmp;
-                tmpr = b->ptr.p_double[offset+k-1];
-                b->ptr.p_double[offset+k-1] = b->ptr.p_double[offset+t-1];
-                b->ptr.p_double[offset+t-1] = tmpr;
+            a->ptr.p_int[p0] = at;
+            a->ptr.p_int[p1] = ak;
+            tmpr = b->ptr.p_double[p0];
+            b->ptr.p_double[p0] = b->ptr.p_double[p1];
+            b->ptr.p_double[p1] = tmpr;
             t = k;
         }
     }
-        i = i+1;
+    for(i=n-1; i>=1; i--)
+    {
+        p0 = offset+0;
+        p1 = offset+i;
+        tmp = a->ptr.p_int[p1];
+        a->ptr.p_int[p1] = a->ptr.p_int[p0];
+        a->ptr.p_int[p0] = tmp;
+        at = tmp;
+        tmpr = b->ptr.p_double[p1];
+        b->ptr.p_double[p1] = b->ptr.p_double[p0];
+        b->ptr.p_double[p0] = tmpr;
+        bt = tmpr;
+        t = 0;
+        for(;;)
+        {
+            k = 2*t+1;
+            if( k+1>i )
+            {
+                break;
             }
-    while(i<=n);
-    i = n-1;
-    do
+            p0 = offset+t;
+            p1 = offset+k;
+            ak = a->ptr.p_int[p1];
+            if( k+1<i )
             {
-        tmp = a->ptr.p_int[offset+i];
-        a->ptr.p_int[offset+i] = a->ptr.p_int[offset+0];
-        a->ptr.p_int[offset+0] = tmp;
-        tmpr = b->ptr.p_double[offset+i];
-        b->ptr.p_double[offset+i] = b->ptr.p_double[offset+0];
-        b->ptr.p_double[offset+0] = tmpr;
-        t = 1;
-        while(t!=0)
+                ak1 = a->ptr.p_int[p1+1];
+                if( ak1>ak )
                 {
-            k = 2*t;
-            if( k>i )
+                    ak = ak1;
+                    p1 = p1+1;
+                    k = k+1;
+                }
+            }
+            if( at>=ak )
             {
-                t = 0;
+                break;
             }
-            else
+            a->ptr.p_int[p1] = at;
+            a->ptr.p_int[p0] = ak;
+            b->ptr.p_double[p0] = b->ptr.p_double[p1];
+            b->ptr.p_double[p1] = bt;
+            t = k;
+        }
+    }
+}
+
+
+/*************************************************************************
+Sorting function optimized for integer values (only keys, no labels),  can
+be used to sort middle of the array
+
+  -- ALGLIB --
+     Copyright 11.12.2008 by Bochkanov Sergey
+*************************************************************************/
+void sortmiddlei(/* Integer */ ae_vector* a,
+     ae_int_t offset,
+     ae_int_t n,
+     ae_state *_state)
+{
+    ae_int_t i;
+    ae_int_t k;
+    ae_int_t t;
+    ae_int_t tmp;
+    ae_int_t p0;
+    ae_int_t p1;
+    ae_int_t at;
+    ae_int_t ak;
+    ae_int_t ak1;
+
+
+    
+    /*
+     * Special cases
+     */
+    if( n<=1 )
     {
-                if( k<i )
+        return;
+    }
+    
+    /*
+     * General case, N>1: sort, update B
+     */
+    for(i=2; i<=n; i++)
     {
-                    if( a->ptr.p_int[offset+k]>a->ptr.p_int[offset+k-1] )
+        t = i;
+        while(t!=1)
         {
-                        k = k+1;
+            k = t/2;
+            p0 = offset+k-1;
+            p1 = offset+t-1;
+            ak = a->ptr.p_int[p0];
+            at = a->ptr.p_int[p1];
+            if( ak>=at )
+            {
+                break;
             }
+            a->ptr.p_int[p0] = at;
+            a->ptr.p_int[p1] = ak;
+            t = k;
         }
-                if( a->ptr.p_int[offset+t-1]>=a->ptr.p_int[offset+k-1] )
+    }
+    for(i=n-1; i>=1; i--)
     {
+        p0 = offset+0;
+        p1 = offset+i;
+        tmp = a->ptr.p_int[p1];
+        a->ptr.p_int[p1] = a->ptr.p_int[p0];
+        a->ptr.p_int[p0] = tmp;
+        at = tmp;
         t = 0;
+        for(;;)
+        {
+            k = 2*t+1;
+            if( k+1>i )
+            {
+                break;
             }
-                else
+            p0 = offset+t;
+            p1 = offset+k;
+            ak = a->ptr.p_int[p1];
+            if( k+1<i )
             {
-                    tmp = a->ptr.p_int[offset+k-1];
-                    a->ptr.p_int[offset+k-1] = a->ptr.p_int[offset+t-1];
-                    a->ptr.p_int[offset+t-1] = tmp;
-                    tmpr = b->ptr.p_double[offset+k-1];
-                    b->ptr.p_double[offset+k-1] = b->ptr.p_double[offset+t-1];
-                    b->ptr.p_double[offset+t-1] = tmpr;
-                    t = k;
+                ak1 = a->ptr.p_int[p1+1];
+                if( ak1>ak )
+                {
+                    ak = ak1;
+                    p1 = p1+1;
+                    k = k+1;
                 }
             }
+            if( at>=ak )
+            {
+                break;
+            }
+            a->ptr.p_int[p1] = at;
+            a->ptr.p_int[p0] = ak;
+            t = k;
         }
-        i = i-1;
     }
-    while(i>=1);
 }
 
 
@@ -4021,7 +4373,7 @@ void tagheappushi(/* Real    */ ae_vector* a,
     {
         k = (j-1)/2;
         v = a->ptr.p_double[k];
-        if( ae_fp_less(v,va) )
+        if( v<va )
         {
             
             /*
@@ -4115,7 +4467,7 @@ void tagheapreplacetopi(/* Real    */ ae_vector* a,
              * have no siblings due to heap structure)
              */
             v = a->ptr.p_double[k1];
-            if( ae_fp_greater(v,va) )
+            if( v>va )
             {
                 a->ptr.p_double[j] = v;
                 b->ptr.p_int[j] = b->ptr.p_int[k1];
@@ -4131,9 +4483,9 @@ void tagheapreplacetopi(/* Real    */ ae_vector* a,
              */
             v1 = a->ptr.p_double[k1];
             v2 = a->ptr.p_double[k2];
-            if( ae_fp_greater(v1,v2) )
+            if( v1>v2 )
             {
-                if( ae_fp_less(va,v1) )
+                if( va<v1 )
                 {
                     a->ptr.p_double[j] = v1;
                     b->ptr.p_int[j] = b->ptr.p_int[k1];
@@ -4146,7 +4498,7 @@ void tagheapreplacetopi(/* Real    */ ae_vector* a,
             }
             else
             {
-                if( ae_fp_less(va,v2) )
+                if( va<v2 )
                 {
                     a->ptr.p_double[j] = v2;
                     b->ptr.p_int[j] = b->ptr.p_int[k2];

3rdparty/alglib/src/alglibinternal.h

 /*************************************************************************
-ALGLIB 3.14.0 (source code generated 2018-06-16)
+ALGLIB 3.15.0 (source code generated 2019-02-20)
 Copyright (c) Sergey Bochkanov (ALGLIB project).
 
 >>> SOURCE LICENSE >>>
@@ -214,6 +214,8 @@ ae_int_t getmlpserializationcode(ae_state *_state);
 ae_int_t getmlpeserializationcode(ae_state *_state);
 ae_int_t getrbfserializationcode(ae_state *_state);
 ae_int_t getspline2dserializationcode(ae_state *_state);
+ae_int_t getidwserializationcode(ae_state *_state);
+ae_int_t getknnserializationcode(ae_state *_state);
 #endif
 #if defined(AE_COMPILE_APSERV) || !defined(AE_PARTIAL_BUILD)
 void seterrorflagdiff(ae_bool* flag,
@@ -272,6 +274,9 @@ void rmatrixsetlengthatleast(/* Real    */ ae_matrix* x,
      ae_int_t m,
      ae_int_t n,
      ae_state *_state);
+void bvectorgrowto(/* Boolean */ ae_vector* x,
+     ae_int_t n,
+     ae_state *_state);
 void ivectorgrowto(/* Integer */ ae_vector* x,
      ae_int_t n,
      ae_state *_state);
@@ -279,6 +284,10 @@ void rmatrixgrowrowsto(/* Real    */ ae_matrix* a,
      ae_int_t n,
      ae_int_t mincols,
      ae_state *_state);
+void rmatrixgrowcolsto(/* Real    */ ae_matrix* a,
+     ae_int_t n,
+     ae_int_t minrows,
+     ae_state *_state);
 void rvectorgrowto(/* Real    */ ae_vector* x,
      ae_int_t n,
      ae_state *_state);
@@ -361,6 +370,7 @@ void inc(ae_int_t* v, ae_state *_state);
 void dec(ae_int_t* v, ae_state *_state);
 void threadunsafeinc(ae_int_t* v, ae_state *_state);
 void countdown(ae_int_t* v, ae_state *_state);
+double possign(double x, ae_state *_state);
 double rmul2(double v0, double v1, ae_state *_state);
 double rmul3(double v0, double v1, double v2, ae_state *_state);
 ae_int_t idivup(ae_int_t a, ae_int_t b, ae_state *_state);
@@ -373,6 +383,13 @@ double rmaxabs3(double r0, double r1, double r2, ae_state *_state);
 double boundval(double x, double b1, double b2, ae_state *_state);
 ae_int_t iboundval(ae_int_t x, ae_int_t b1, ae_int_t b2, ae_state *_state);
 double rboundval(double x, double b1, double b2, ae_state *_state);
+ae_int_t countnz1(/* Real    */ ae_vector* v,
+     ae_int_t n,
+     ae_state *_state);
+ae_int_t countnz2(/* Real    */ ae_matrix* v,
+     ae_int_t m,
+     ae_int_t n,
+     ae_state *_state);
 void alloccomplex(ae_serializer* s, ae_complex v, ae_state *_state);
 void serializecomplex(ae_serializer* s, ae_complex v, ae_state *_state);
 ae_complex unserializecomplex(ae_serializer* s, ae_state *_state);
@@ -411,6 +428,9 @@ void serializerealmatrix(ae_serializer* s,
 void unserializerealmatrix(ae_serializer* s,
      /* Real    */ ae_matrix* v,
      ae_state *_state);
+void copybooleanarray(/* Boolean */ ae_vector* src,
+     /* Boolean */ ae_vector* dst,
+     ae_state *_state);
 void copyintegerarray(/* Integer */ ae_vector* src,
      /* Integer */ ae_vector* dst,
      ae_state *_state);
@@ -442,6 +462,7 @@ void splitlengtheven(ae_int_t tasksize,
 ae_int_t chunkscount(ae_int_t tasksize,
      ae_int_t chunksize,
      ae_state *_state);
+double sparselevel2density(ae_state *_state);
 ae_int_t matrixtilesizea(ae_state *_state);
 ae_int_t matrixtilesizeb(ae_state *_state);
 double smpactivationlevel(ae_state *_state);
@@ -521,6 +542,10 @@ void tagsortmiddleir(/* Integer */ ae_vector* a,
      ae_int_t offset,
      ae_int_t n,
      ae_state *_state);
+void sortmiddlei(/* Integer */ ae_vector* a,
+     ae_int_t offset,
+     ae_int_t n,
+     ae_state *_state);
 void tagheappushi(/* Real    */ ae_vector* a,
      /* Integer */ ae_vector* b,
      ae_int_t* n,

3rdparty/alglib/src/alglibmisc.h

 /*************************************************************************
-ALGLIB 3.14.0 (source code generated 2018-06-16)
+ALGLIB 3.15.0 (source code generated 2019-02-20)
 Copyright (c) Sergey Bochkanov (ALGLIB project).
 
 >>> SOURCE LICENSE >>>
@@ -95,7 +95,7 @@ namespace alglib
 Buffer object which is used to perform nearest neighbor  requests  in  the
 multithreaded mode (multiple threads working with same KD-tree object).
 
-This object should be created with KDTreeCreateBuffer().
+This object should be created with KDTreeCreateRequestBuffer().
 *************************************************************************/
 class _kdtreerequestbuffer_owner
 {
@@ -361,7 +361,7 @@ OUTPUT PARAMETERS
 
 
 IMPORTANT: KD-tree buffer should be used only with  KD-tree  object  which
-           was used to initialize buffer. Any attempt to use biffer   with
+           was used to initialize buffer. Any attempt to use buffer   with
            different object is dangerous - you  may  get  integrity  check
            failure (exception) because sizes of internal arrays do not fit
            to dimensions of KD-tree structure.
@@ -456,13 +456,18 @@ ae_int_t kdtreetsqueryknn(const kdtree &kdt, const kdtreerequestbuffer &buf, con
 
 
 /*************************************************************************
-R-NN query: all points within R-sphere centered at X
+R-NN query: all points within R-sphere centered at X, ordered by  distance
+between point and X (by ascending).
+
+NOTE: it is also possible to perform undordered queries performed by means
+      of kdtreequeryrnnu() and kdtreetsqueryrnnu() functions. Such queries
+      are faster because we do not have to use heap structure for sorting.
 
 IMPORTANT: this function can not be used in multithreaded code because  it
            uses internal temporary buffer of kd-tree object, which can not
            be shared between multiple threads.  If  you  want  to  perform
            parallel requests, use function  which  uses  external  request
-           buffer: KDTreeTsQueryRNN() ("Ts" stands for "thread-safe").
+           buffer: kdtreetsqueryrnn() ("Ts" stands for "thread-safe").
 
 INPUT PARAMETERS
     KDT         -   KD-tree
@@ -493,14 +498,60 @@ ae_int_t kdtreequeryrnn(const kdtree &kdt, const real_1d_array &x, const double
 ae_int_t kdtreequeryrnn(const kdtree &kdt, const real_1d_array &x, const double r, const xparams _xparams = alglib::xdefault);
 
 
+/*************************************************************************
+R-NN query: all points within R-sphere  centered  at  X,  no  ordering  by
+distance as undicated by "U" suffix (faster that ordered query, for  large
+queries - significantly faster).
+
+IMPORTANT: this function can not be used in multithreaded code because  it
+           uses internal temporary buffer of kd-tree object, which can not
+           be shared between multiple threads.  If  you  want  to  perform
+           parallel requests, use function  which  uses  external  request
+           buffer: kdtreetsqueryrnn() ("Ts" stands for "thread-safe").
+
+INPUT PARAMETERS
+    KDT         -   KD-tree
+    X           -   point, array[0..NX-1].
+    R           -   radius of sphere (in corresponding norm), R>0
+    SelfMatch   -   whether self-matches are allowed:
+                    * if True, nearest neighbor may be the point itself
+                      (if it exists in original dataset)
+                    * if False, then only points with non-zero distance
+                      are returned
+                    * if not given, considered True
+
+RESULT
+    number of neighbors found, >=0
+
+This  subroutine  performs  query  and  stores  its result in the internal
+structures of the KD-tree. You can use  following  subroutines  to  obtain
+actual results:
+* KDTreeQueryResultsX() to get X-values
+* KDTreeQueryResultsXY() to get X- and Y-values
+* KDTreeQueryResultsTags() to get tag values
+* KDTreeQueryResultsDistances() to get distances
+
+As indicated by "U" suffix, this function returns unordered results.
+
+  -- ALGLIB --
+     Copyright 01.11.2018 by Bochkanov Sergey
+*************************************************************************/
+ae_int_t kdtreequeryrnnu(const kdtree &kdt, const real_1d_array &x, const double r, const bool selfmatch, const xparams _xparams = alglib::xdefault);
+ae_int_t kdtreequeryrnnu(const kdtree &kdt, const real_1d_array &x, const double r, const xparams _xparams = alglib::xdefault);
+
+
 /*************************************************************************
 R-NN query: all points within  R-sphere  centered  at  X,  using  external
-thread-local buffer.
+thread-local buffer, sorted by distance between point and X (by ascending)
 
 You can call this function from multiple threads for same kd-tree instance,
 assuming that different instances of buffer object are passed to different
 threads.
 
+NOTE: it is also possible to perform undordered queries performed by means
+      of kdtreequeryrnnu() and kdtreetsqueryrnnu() functions. Such queries
+      are faster because we do not have to use heap structure for sorting.
+
 INPUT PARAMETERS
     KDT         -   KD-tree
     Buf         -   request buffer  object  created  for  this  particular
@@ -539,6 +590,55 @@ ae_int_t kdtreetsqueryrnn(const kdtree &kdt, const kdtreerequestbuffer &buf, con
 ae_int_t kdtreetsqueryrnn(const kdtree &kdt, const kdtreerequestbuffer &buf, const real_1d_array &x, const double r, const xparams _xparams = alglib::xdefault);
 
 
+/*************************************************************************
+R-NN query: all points within  R-sphere  centered  at  X,  using  external
+thread-local buffer, no ordering by distance as undicated  by  "U"  suffix
+(faster that ordered query, for large queries - significantly faster).
+
+You can call this function from multiple threads for same kd-tree instance,
+assuming that different instances of buffer object are passed to different
+threads.
+
+INPUT PARAMETERS
+    KDT         -   KD-tree
+    Buf         -   request buffer  object  created  for  this  particular
+                    instance of kd-tree structure with kdtreecreaterequestbuffer()
+                    function.
+    X           -   point, array[0..NX-1].
+    R           -   radius of sphere (in corresponding norm), R>0
+    SelfMatch   -   whether self-matches are allowed:
+                    * if True, nearest neighbor may be the point itself
+                      (if it exists in original dataset)
+                    * if False, then only points with non-zero distance
+                      are returned
+                    * if not given, considered True
+
+RESULT
+    number of neighbors found, >=0
+
+This  subroutine  performs  query  and  stores  its result in the internal
+structures  of  the  buffer object. You can use following  subroutines  to
+obtain these results (pay attention to "buf" in their names):
+* KDTreeTsQueryResultsX() to get X-values
+* KDTreeTsQueryResultsXY() to get X- and Y-values
+* KDTreeTsQueryResultsTags() to get tag values
+* KDTreeTsQueryResultsDistances() to get distances
+
+As indicated by "U" suffix, this function returns unordered results.
+
+IMPORTANT: kd-tree buffer should be used only with  KD-tree  object  which
+           was used to initialize buffer. Any attempt to use biffer   with
+           different object is dangerous - you  may  get  integrity  check
+           failure (exception) because sizes of internal arrays do not fit
+           to dimensions of KD-tree structure.
+
+  -- ALGLIB --
+     Copyright 18.03.2016 by Bochkanov Sergey
+*************************************************************************/
+ae_int_t kdtreetsqueryrnnu(const kdtree &kdt, const kdtreerequestbuffer &buf, const real_1d_array &x, const double r, const bool selfmatch, const xparams _xparams = alglib::xdefault);
+ae_int_t kdtreetsqueryrnnu(const kdtree &kdt, const kdtreerequestbuffer &buf, const real_1d_array &x, const double r, const xparams _xparams = alglib::xdefault);
+
+
 /*************************************************************************
 K-NN query: approximate K nearest neighbors
 
@@ -1696,12 +1796,23 @@ ae_int_t kdtreequeryrnn(kdtree* kdt,
      double r,
      ae_bool selfmatch,
      ae_state *_state);
+ae_int_t kdtreequeryrnnu(kdtree* kdt,
+     /* Real    */ ae_vector* x,
+     double r,
+     ae_bool selfmatch,
+     ae_state *_state);
 ae_int_t kdtreetsqueryrnn(kdtree* kdt,
      kdtreerequestbuffer* buf,
      /* Real    */ ae_vector* x,
      double r,
      ae_bool selfmatch,
      ae_state *_state);
+ae_int_t kdtreetsqueryrnnu(kdtree* kdt,
+     kdtreerequestbuffer* buf,
+     /* Real    */ ae_vector* x,
+     double r,
+     ae_bool selfmatch,
+     ae_state *_state);
 ae_int_t kdtreequeryaknn(kdtree* kdt,
      /* Real    */ ae_vector* x,
      ae_int_t k,

3rdparty/alglib/src/ap.cpp

 /*************************************************************************
-ALGLIB 3.14.0 (source code generated 2018-06-16)
+ALGLIB 3.15.0 (source code generated 2019-02-20)
 Copyright (c) Sergey Bochkanov (ALGLIB project).
 
 >>> SOURCE LICENSE >>>
@@ -791,6 +791,17 @@ void* aligned_malloc(size_t size, size_t alignment)
 #endif
 }
 
+void* aligned_extract_ptr(void *block)
+{
+#if AE_MALLOC==AE_BASIC_STATIC_MALLOC
+    return NULL;
+#else
+    if( block==NULL )
+        return NULL;
+    return *((void**)((char*)block-sizeof(void*)));
+#endif
+}
+
 void aligned_free(void *block)
 {
 #if AE_MALLOC==AE_BASIC_STATIC_MALLOC
@@ -799,7 +810,7 @@ void aligned_free(void *block)
     void *p;
     if( block==NULL )
         return;
-    p = *((void**)((char*)block-sizeof(void*)));
+    p = aligned_extract_ptr(block);
     free(p);
     if( _use_alloc_counter )
         ae_optional_atomic_sub_i(&_alloc_counter, 1);
@@ -891,11 +902,40 @@ on entry are correctly initialized by zeros.
 ************************************************************************/
 ae_bool ae_check_zeros(const void *ptr, ae_int_t n)
 {
-    const unsigned char *p = (const unsigned char*)ptr;
-    unsigned char c = 0x0;
-    ae_int_t i;
-    for(i=0; i<n; i++)
-        c |= p[i];
+    ae_int_t nu, nr, i;
+    unsigned long long c = 0x0;
+    
+    /*
+     * determine leading and trailing lengths
+     */
+    nu = n/sizeof(unsigned long long);
+    nr = n%sizeof(unsigned long long);
+    
+    /*
+     * handle leading nu long long elements
+     */
+    if( nu>0 )
+    {
+        const unsigned long long *p_ull;
+        p_ull = (const unsigned long long *)ptr;
+        for(i=0; i<nu; i++)
+            c |= p_ull[i];
+    }
+    
+    /*
+     * handle trailing nr char elements
+     */
+    if( nr>0 )
+    {
+        const unsigned char *p_uc;
+        p_uc  = ((const unsigned char *)ptr)+nu*sizeof(unsigned long long);
+        for(i=0; i<nr; i++)
+            c |= p_uc[i];
+    }
+    
+    /*
+     * done
+     */
     return c==0x0;
 }
 
@@ -1104,13 +1144,18 @@ void ae_db_init(ae_dyn_block *block, ae_int_t size, ae_state *state, ae_bool mak
      */
     ae_assert(size>=0, "ae_db_init(): negative size", state);
     block->ptr = NULL;
+    block->valgrind_hint = NULL;
     ae_touch_ptr(block->ptr);
+    ae_touch_ptr(block->valgrind_hint);
     if( make_automatic )
         ae_db_attach(block, state);
     else
         block->p_next = NULL;
     if( size!=0 )
+    {
         block->ptr = ae_malloc((size_t)size, state);
+        block->valgrind_hint = aligned_extract_ptr(block->ptr);
+    }
     block->deallocator = ae_free;
 }
 
@@ -1148,8 +1193,10 @@ void ae_db_realloc(ae_dyn_block *block, ae_int_t size, ae_state *state)
     {
         ((ae_deallocator)block->deallocator)(block->ptr);
         block->ptr = NULL;
+        block->valgrind_hint = NULL;
     }
     block->ptr = ae_malloc((size_t)size, state);
+    block->valgrind_hint = aligned_extract_ptr(block->ptr);
     block->deallocator = ae_free;
 }
 
@@ -1169,6 +1216,7 @@ void ae_db_free(ae_dyn_block *block)
     if( block->ptr!=NULL )
         ((ae_deallocator)block->deallocator)(block->ptr);
     block->ptr = NULL;
+    block->valgrind_hint = NULL;
     block->deallocator = ae_free;
 }
 
@@ -1184,11 +1232,18 @@ void ae_db_swap(ae_dyn_block *block1, ae_dyn_block *block2)
 {
     void (*deallocator)(void*) = NULL;
     void * volatile ptr;
+    void * valgrind_hint;
+    
     ptr = block1->ptr;
+    valgrind_hint = block1->valgrind_hint;
     deallocator = block1->deallocator;
+    
     block1->ptr = block2->ptr;
+    block1->valgrind_hint = block2->valgrind_hint;
     block1->deallocator = block2->deallocator;
+    
     block2->ptr = ptr;
+    block2->valgrind_hint = valgrind_hint;
     block2->deallocator = deallocator;
 }
 
@@ -1283,7 +1338,7 @@ void ae_vector_init_from_x(ae_vector *dst, x_vector *src, ae_state *state, ae_bo
     
     ae_vector_init(dst, (ae_int_t)src->cnt, (ae_datatype)src->datatype, state, make_automatic);
     if( src->cnt>0 )
-        memmove(dst->ptr.p_ptr, src->ptr, (size_t)(((ae_int_t)src->cnt)*ae_sizeof((ae_datatype)src->datatype)));
+        memmove(dst->ptr.p_ptr, src->x_ptr.p_ptr, (size_t)(((ae_int_t)src->cnt)*ae_sizeof((ae_datatype)src->datatype)));
 }
 
 /************************************************************************
@@ -1333,7 +1388,7 @@ void ae_vector_init_attach_to_x(ae_vector *dst, x_vector *src, ae_state *state,
     
     /* init */
     dst->cnt = cnt;
-    dst->ptr.p_ptr = src->ptr;
+    dst->ptr.p_ptr = src->x_ptr.p_ptr;
     dst->is_attached = ae_true;
 }
 
@@ -1539,7 +1594,7 @@ void ae_matrix_init_from_x(ae_matrix *dst, x_matrix *src, ae_state *state, ae_bo
     ae_matrix_init(dst, (ae_int_t)src->rows, (ae_int_t)src->cols, (ae_datatype)src->datatype, state, make_automatic);
     if( src->rows!=0 && src->cols!=0 )
     {
-        p_src_row = (char*)src->ptr;
+        p_src_row = (char*)src->x_ptr.p_ptr;
         p_dst_row = (char*)(dst->ptr.pp_void[0]);
         row_size = ae_sizeof((ae_datatype)src->datatype)*(ae_int_t)src->cols;
         for(i=0; i<src->rows; i++, p_src_row+=src->stride*ae_sizeof((ae_datatype)src->datatype), p_dst_row+=dst->stride*ae_sizeof((ae_datatype)src->datatype))
@@ -1610,7 +1665,7 @@ void ae_matrix_init_attach_to_x(ae_matrix *dst, x_matrix *src, ae_state *state,
         char *p_row;
         void **pp_ptr;
         
-        p_row = (char*)src->ptr;
+        p_row = (char*)src->x_ptr.p_ptr;
         rowsize = dst->stride*ae_sizeof(dst->datatype);
         pp_ptr  = (void**)dst->data.ptr;
         dst->ptr.pp_void = pp_ptr;
@@ -1908,7 +1963,7 @@ NOTES:
 ************************************************************************/
 void ae_x_set_vector(x_vector *dst, ae_vector *src, ae_state *state)
 {
-    if( src->ptr.p_ptr == dst->ptr )
+    if( src->ptr.p_ptr == dst->x_ptr.p_ptr )
     {
         /* src->ptr points to the beginning of dst, attached matrices, no need to copy */
         return;
@@ -1916,9 +1971,9 @@ void ae_x_set_vector(x_vector *dst, ae_vector *src, ae_state *state)
     if( dst->cnt!=src->cnt || dst->datatype!=src->datatype )
     {
         if( dst->owner==OWN_AE )
-            ae_free(dst->ptr);
-        dst->ptr = ae_malloc((size_t)(src->cnt*ae_sizeof(src->datatype)), state);
-        if( src->cnt!=0 && dst->ptr==NULL )
+            ae_free(dst->x_ptr.p_ptr);
+        dst->x_ptr.p_ptr = ae_malloc((size_t)(src->cnt*ae_sizeof(src->datatype)), state);
+        if( src->cnt!=0 && dst->x_ptr.p_ptr==NULL )
             ae_break(state, ERR_OUT_OF_MEMORY, "ae_malloc(): out of memory");
         dst->last_action = ACT_NEW_LOCATION;
         dst->cnt = src->cnt;
@@ -1937,7 +1992,7 @@ void ae_x_set_vector(x_vector *dst, ae_vector *src, ae_state *state)
             ae_assert(ae_false, "ALGLIB: internal error in ae_x_set_vector()", state);
     }
     if( src->cnt )
-        memmove(dst->ptr, src->ptr.p_ptr, (size_t)(src->cnt*ae_sizeof(src->datatype)));
+        memmove(dst->x_ptr.p_ptr, src->ptr.p_ptr, (size_t)(src->cnt*ae_sizeof(src->datatype)));
 }
 
 /************************************************************************
@@ -1972,7 +2027,7 @@ void ae_x_set_matrix(x_matrix *dst, ae_matrix *src, ae_state *state)
     char *p_dst_row;
     ae_int_t i;
     ae_int_t row_size;
-    if( src->ptr.pp_void!=NULL && src->ptr.pp_void[0] == dst->ptr )
+    if( src->ptr.pp_void!=NULL && src->ptr.pp_void[0] == dst->x_ptr.p_ptr )
     {
         /* src->ptr points to the beginning of dst, attached matrices, no need to copy */
         return;
@@ -1980,13 +2035,13 @@ void ae_x_set_matrix(x_matrix *dst, ae_matrix *src, ae_state *state)
     if( dst->rows!=src->rows || dst->cols!=src->cols || dst->datatype!=src->datatype )
     {
         if( dst->owner==OWN_AE )
-            ae_free(dst->ptr);
+            ae_free(dst->x_ptr.p_ptr);
         dst->rows = src->rows;
         dst->cols = src->cols;
         dst->stride = src->cols;
         dst->datatype = src->datatype;
-        dst->ptr = ae_malloc((size_t)(dst->rows*((ae_int_t)dst->stride)*ae_sizeof(src->datatype)), state);
-        if( dst->rows!=0 && dst->stride!=0 && dst->ptr==NULL )
+        dst->x_ptr.p_ptr = ae_malloc((size_t)(dst->rows*((ae_int_t)dst->stride)*ae_sizeof(src->datatype)), state);
+        if( dst->rows!=0 && dst->stride!=0 && dst->x_ptr.p_ptr==NULL )
             ae_break(state, ERR_OUT_OF_MEMORY, "ae_malloc(): out of memory");
         dst->last_action = ACT_NEW_LOCATION;
         dst->owner = OWN_AE;
@@ -2005,7 +2060,7 @@ void ae_x_set_matrix(x_matrix *dst, ae_matrix *src, ae_state *state)
     if( src->rows!=0 && src->cols!=0 )
     {
         p_src_row = (char*)(src->ptr.pp_void[0]);
-        p_dst_row = (char*)dst->ptr;
+        p_dst_row = (char*)dst->x_ptr.p_ptr;
         row_size = ae_sizeof(src->datatype)*src->cols;
         for(i=0; i<src->rows; i++, p_src_row+=src->stride*ae_sizeof(src->datatype), p_dst_row+=dst->stride*ae_sizeof(src->datatype))
             memmove(p_dst_row, p_src_row, (size_t)(row_size));
@@ -2030,8 +2085,8 @@ NOTES:
 void ae_x_attach_to_vector(x_vector *dst, ae_vector *src)
 {
     if( dst->owner==OWN_AE )
-        ae_free(dst->ptr);
-    dst->ptr = src->ptr.p_ptr;
+        ae_free(dst->x_ptr.p_ptr);
+    dst->x_ptr.p_ptr = src->ptr.p_ptr;
     dst->last_action = ACT_NEW_LOCATION;
     dst->cnt = src->cnt;
     dst->datatype = src->datatype;
@@ -2056,12 +2111,12 @@ NOTES:
 void ae_x_attach_to_matrix(x_matrix *dst, ae_matrix *src)
 {
     if( dst->owner==OWN_AE )
-            ae_free(dst->ptr);
+            ae_free(dst->x_ptr.p_ptr);
     dst->rows = src->rows;
     dst->cols = src->cols;
     dst->stride = src->stride;
     dst->datatype = src->datatype;
-    dst->ptr = &(src->ptr.pp_double[0][0]);
+    dst->x_ptr.p_ptr = &(src->ptr.pp_double[0][0]);
     dst->last_action = ACT_NEW_LOCATION;
     dst->owner = OWN_CALLER;
 }
@@ -2075,8 +2130,8 @@ dst                 vector
 void x_vector_clear(x_vector *dst)
 {
     if( dst->owner==OWN_AE )
-        aligned_free(dst->ptr);
-    dst->ptr = NULL;
+        aligned_free(dst->x_ptr.p_ptr);
+    dst->x_ptr.p_ptr = NULL;
     dst->cnt = 0;
 }
 
@@ -2638,8 +2693,8 @@ static void is_symmetric_rec_off_stat(x_matrix *a, ae_int_t offset0, ae_int_t of
         double v;
         ae_int_t i, j;
 
-        p1 = (double*)(a->ptr)+offset0*a->stride+offset1;
-        p2 = (double*)(a->ptr)+offset1*a->stride+offset0;
+        p1 = (double*)(a->x_ptr.p_ptr)+offset0*a->stride+offset1;
+        p2 = (double*)(a->x_ptr.p_ptr)+offset1*a->stride+offset0;
         for(i=0; i<len0; i++)
         {
             pcol = p2+i;
@@ -2697,7 +2752,7 @@ static void is_symmetric_rec_diag_stat(x_matrix *a, ae_int_t offset, ae_int_t le
     }
     
     /* base case */
-    p = (double*)(a->ptr)+offset*a->stride+offset;
+    p = (double*)(a->x_ptr.p_ptr)+offset*a->stride+offset;
     for(i=0; i<len; i++)
     {
         pcol = p+i;
@@ -2764,8 +2819,8 @@ static void is_hermitian_rec_off_stat(x_matrix *a, ae_int_t offset0, ae_int_t of
         double v;
         ae_int_t i, j;
 
-        p1 = (ae_complex*)(a->ptr)+offset0*a->stride+offset1;
-        p2 = (ae_complex*)(a->ptr)+offset1*a->stride+offset0;
+        p1 = (ae_complex*)(a->x_ptr.p_ptr)+offset0*a->stride+offset1;
+        p2 = (ae_complex*)(a->x_ptr.p_ptr)+offset1*a->stride+offset0;
         for(i=0; i<len0; i++)
         {
             pcol = p2+i;
@@ -2823,7 +2878,7 @@ static void is_hermitian_rec_diag_stat(x_matrix *a, ae_int_t offset, ae_int_t le
     }
     
     /* base case */
-    p = (ae_complex*)(a->ptr)+offset*a->stride+offset;
+    p = (ae_complex*)(a->x_ptr.p_ptr)+offset*a->stride+offset;
     for(i=0; i<len; i++)
     {
         pcol = p+i;
@@ -2894,8 +2949,8 @@ static void force_symmetric_rec_off_stat(x_matrix *a, ae_int_t offset0, ae_int_t
         double *p1, *p2, *prow, *pcol;
         ae_int_t i, j;
 
-        p1 = (double*)(a->ptr)+offset0*a->stride+offset1;
-        p2 = (double*)(a->ptr)+offset1*a->stride+offset0;
+        p1 = (double*)(a->x_ptr.p_ptr)+offset0*a->stride+offset1;
+        p2 = (double*)(a->x_ptr.p_ptr)+offset1*a->stride+offset0;
         for(i=0; i<len0; i++)
         {
             pcol = p2+i;
@@ -2936,7 +2991,7 @@ static void force_symmetric_rec_diag_stat(x_matrix *a, ae_int_t offset, ae_int_t
     }
     
     /* base case */
-    p = (double*)(a->ptr)+offset*a->stride+offset;
+    p = (double*)(a->x_ptr.p_ptr)+offset*a->stride+offset;
     for(i=0; i<len; i++)
     {
         pcol = p+i;
@@ -2981,8 +3036,8 @@ static void force_hermitian_rec_off_stat(x_matrix *a, ae_int_t offset0, ae_int_t
         ae_complex *p1, *p2, *prow, *pcol;
         ae_int_t i, j;
 
-        p1 = (ae_complex*)(a->ptr)+offset0*a->stride+offset1;
-        p2 = (ae_complex*)(a->ptr)+offset1*a->stride+offset0;
+        p1 = (ae_complex*)(a->x_ptr.p_ptr)+offset0*a->stride+offset1;
+        p2 = (ae_complex*)(a->x_ptr.p_ptr)+offset1*a->stride+offset0;
         for(i=0; i<len0; i++)
         {
             pcol = p2+i;
@@ -3023,7 +3078,7 @@ static void force_hermitian_rec_diag_stat(x_matrix *a, ae_int_t offset, ae_int_t
     }
     
     /* base case */
-    p = (ae_complex*)(a->ptr)+offset*a->stride+offset;
+    p = (ae_complex*)(a->x_ptr.p_ptr)+offset*a->stride+offset;
     for(i=0; i<len; i++)
     {
         pcol = p+i;
@@ -7487,7 +7542,7 @@ void alglib::real_1d_array::attach_to_ptr(ae_int_t iLen, double *pContent ) // T
     x.datatype = alglib_impl::DT_REAL;
     x.owner = alglib_impl::OWN_CALLER;
     x.last_action = alglib_impl::ACT_UNCHANGED;
-    x.ptr = pContent;
+    x.x_ptr.p_ptr = pContent;
     attach_to(&x, &_state);
     ae_state_clear(&_state);
 }
@@ -8070,7 +8125,7 @@ void alglib::real_2d_array::attach_to_ptr(ae_int_t irows, ae_int_t icols, double
     x.datatype = alglib_impl::DT_REAL;
     x.owner = alglib_impl::OWN_CALLER;
     x.last_action = alglib_impl::ACT_UNCHANGED;
-    x.ptr = pContent;
+    x.x_ptr.p_ptr = pContent;
     attach_to(&x, &_state);
     ae_state_clear(&_state);
 }

3rdparty/alglib/src/diffequations.cpp

 /*************************************************************************
-ALGLIB 3.14.0 (source code generated 2018-06-16)
+ALGLIB 3.15.0 (source code generated 2019-02-20)
 Copyright (c) Sergey Bochkanov (ALGLIB project).
 
 >>> SOURCE LICENSE >>>

3rdparty/alglib/src/diffequations.h

 /*************************************************************************
-ALGLIB 3.14.0 (source code generated 2018-06-16)
+ALGLIB 3.15.0 (source code generated 2019-02-20)
 Copyright (c) Sergey Bochkanov (ALGLIB project).
 
 >>> SOURCE LICENSE >>>

3rdparty/alglib/src/fasttransforms.cpp

 /*************************************************************************
-ALGLIB 3.14.0 (source code generated 2018-06-16)
+ALGLIB 3.15.0 (source code generated 2019-02-20)
 Copyright (c) Sergey Bochkanov (ALGLIB project).
 
 >>> SOURCE LICENSE >>>

3rdparty/alglib/src/fasttransforms.h

 /*************************************************************************
-ALGLIB 3.14.0 (source code generated 2018-06-16)
+ALGLIB 3.15.0 (source code generated 2019-02-20)
 Copyright (c) Sergey Bochkanov (ALGLIB project).
 
 >>> SOURCE LICENSE >>>

3rdparty/alglib/src/integration.cpp

 /*************************************************************************
-ALGLIB 3.14.0 (source code generated 2018-06-16)
+ALGLIB 3.15.0 (source code generated 2019-02-20)
 Copyright (c) Sergey Bochkanov (ALGLIB project).
 
 >>> SOURCE LICENSE >>>

3rdparty/alglib/src/integration.h

 /*************************************************************************
-ALGLIB 3.14.0 (source code generated 2018-06-16)
+ALGLIB 3.15.0 (source code generated 2019-02-20)
 Copyright (c) Sergey Bochkanov (ALGLIB project).
 
 >>> SOURCE LICENSE >>>

3rdparty/alglib/src/linalg.h

 /*************************************************************************
-ALGLIB 3.14.0 (source code generated 2018-06-16)
+ALGLIB 3.15.0 (source code generated 2019-02-20)
 Copyright (c) Sergey Bochkanov (ALGLIB project).
 
 >>> SOURCE LICENSE >>>
@@ -54,6 +54,64 @@ typedef struct
 #endif
 #if defined(AE_COMPILE_ABLAS) || !defined(AE_PARTIAL_BUILD)
 #endif
+#if defined(AE_COMPILE_DLU) || !defined(AE_PARTIAL_BUILD)
+#endif
+#if defined(AE_COMPILE_SPTRF) || !defined(AE_PARTIAL_BUILD)
+typedef struct
+{
+    ae_int_t nfixed;
+    ae_int_t ndynamic;
+    ae_vector idxfirst;
+    ae_vector strgidx;
+    ae_vector strgval;
+    ae_int_t nallocated;
+    ae_int_t nused;
+} sluv2list1matrix;
+typedef struct
+{
+    ae_int_t n;
+    ae_int_t k;
+    ae_vector nzc;
+    ae_int_t maxwrkcnt;
+    ae_int_t maxwrknz;
+    ae_int_t wrkcnt;
+    ae_vector wrkset;
+    ae_vector colid;
+    ae_vector isdensified;
+    ae_vector slscolptr;
+    ae_vector slsrowptr;
+    ae_vector slsidx;
+    ae_vector slsval;
+    ae_int_t slsused;
+    ae_vector tmp0;
+} sluv2sparsetrail;
+typedef struct
+{
+    ae_int_t n;
+    ae_int_t ndense;
+    ae_matrix d;
+    ae_vector did;
+} sluv2densetrail;
+typedef struct
+{
+    ae_int_t n;
+    sparsematrix sparsel;
+    sparsematrix sparseut;
+    sluv2list1matrix bleft;
+    sluv2list1matrix bupper;
+    sluv2sparsetrail strail;
+    sluv2densetrail dtrail;
+    ae_vector rowpermrawidx;
+    ae_matrix dbuf;
+    ae_vector v0i;
+    ae_vector v1i;
+    ae_vector v0r;
+    ae_vector v1r;
+    ae_vector tmp0;
+    ae_vector tmpi;
+    ae_vector tmpp;
+} sluv2buffer;
+#endif
 #if defined(AE_COMPILE_MATGEN) || !defined(AE_PARTIAL_BUILD)
 #endif
 #if defined(AE_COMPILE_TRFAC) || !defined(AE_PARTIAL_BUILD)
@@ -273,6 +331,14 @@ public:
 
 #endif
 
+#if defined(AE_COMPILE_DLU) || !defined(AE_PARTIAL_BUILD)
+
+#endif
+
+#if defined(AE_COMPILE_SPTRF) || !defined(AE_PARTIAL_BUILD)
+
+#endif
+
 #if defined(AE_COMPILE_MATGEN) || !defined(AE_PARTIAL_BUILD)
 
 #endif
@@ -1423,6 +1489,39 @@ NOTE: internal  temporary  copy  is  allocated   for   the   purposes   of
 void sparsetransposecrs(const sparsematrix &s, const xparams _xparams = alglib::xdefault);
 
 
+/*************************************************************************
+This function performs copying with transposition of CRS matrix.
+
+INPUT PARAMETERS
+    S0      -   sparse matrix in CRS format.
+
+OUTPUT PARAMETERS
+    S1      -   sparse matrix, transposed
+
+  -- ALGLIB PROJECT --
+     Copyright 23.07.2018 by Bochkanov Sergey
+*************************************************************************/
+void sparsecopytransposecrs(const sparsematrix &s0, sparsematrix &s1, const xparams _xparams = alglib::xdefault);
+
+
+/*************************************************************************
+This function performs copying with transposition of CRS matrix  (buffered
+version which reuses memory already allocated by  the  target as  much  as
+possible).
+
+INPUT PARAMETERS
+    S0      -   sparse matrix in CRS format.
+
+OUTPUT PARAMETERS
+    S1      -   sparse matrix, transposed; previously allocated memory  is
+                reused if possible.
+
+  -- ALGLIB PROJECT --
+     Copyright 23.07.2018 by Bochkanov Sergey
+*************************************************************************/
+void sparsecopytransposecrsbuf(const sparsematrix &s0, const sparsematrix &s1, const xparams _xparams = alglib::xdefault);
+
+
 /*************************************************************************
 This  function  performs  in-place  conversion  to  desired sparse storage
 format.
@@ -2494,6 +2593,14 @@ compatibility.
 void cmatrixsyrk(const ae_int_t n, const ae_int_t k, const double alpha, const complex_2d_array &a, const ae_int_t ia, const ae_int_t ja, const ae_int_t optypea, const double beta, const complex_2d_array &c, const ae_int_t ic, const ae_int_t jc, const bool isupper, const xparams _xparams = alglib::xdefault);
 #endif
 
+#if defined(AE_COMPILE_DLU) || !defined(AE_PARTIAL_BUILD)
+
+#endif
+
+#if defined(AE_COMPILE_SPTRF) || !defined(AE_PARTIAL_BUILD)
+
+#endif
+
 #if defined(AE_COMPILE_MATGEN) || !defined(AE_PARTIAL_BUILD)
 /*************************************************************************
 Generation of a random uniformly distributed (Haar) orthogonal matrix
@@ -3102,6 +3209,53 @@ OUTPUT PARAMETERS:
 void spdmatrixcholeskyupdatefixbuf(const real_2d_array &a, const ae_int_t n, const bool isupper, const boolean_1d_array &fix, real_1d_array &bufr, const xparams _xparams = alglib::xdefault);
 
 
+/*************************************************************************
+Sparse LU decomposition with column pivoting for sparsity and row pivoting
+for stability. Input must be square sparse matrix stored in CRS format.
+
+The algorithm  computes  LU  decomposition  of  a  general  square  matrix
+(rectangular ones are not supported). The result  of  an  algorithm  is  a
+representation of A as A = P*L*U*Q, where:
+* L is lower unitriangular matrix
+* U is upper triangular matrix
+* P = P0*P1*...*PK, K=N-1, Pi - permutation matrix for I and P[I]
+* Q = QK*...*Q1*Q0, K=N-1, Qi - permutation matrix for I and Q[I]
+
+This function pivots columns for higher sparsity, and then pivots rows for
+stability (larger element at the diagonal).
+
+INPUT PARAMETERS:
+    A       -   sparse NxN matrix in CRS format. An exception is generated
+                if matrix is non-CRS or non-square.
+    PivotType-  pivoting strategy:
+                * 0 for best pivoting available (2 in current version)
+                * 1 for row-only pivoting (NOT RECOMMENDED)
+                * 2 for complete pivoting which produces most sparse outputs
+
+OUTPUT PARAMETERS:
+    A       -   the result of factorization, matrices L and U stored in
+                compact form using CRS sparse storage format:
+                * lower unitriangular L is stored strictly under main diagonal
+                * upper triangilar U is stored ON and ABOVE main diagonal
+    P       -   row permutation matrix in compact form, array[N]
+    Q       -   col permutation matrix in compact form, array[N]
+
+This function always succeeds, i.e. it ALWAYS returns valid factorization,
+but for your convenience it also returns  boolean  value  which  helps  to
+detect symbolically degenerate matrices:
+* function returns TRUE, if the matrix was factorized AND symbolically
+  non-degenerate
+* function returns FALSE, if the matrix was factorized but U has strictly
+  zero elements at the diagonal (the factorization is returned anyway).
+
+
+  -- ALGLIB routine --
+     03.09.2018
+     Bochkanov Sergey
+*************************************************************************/
+bool sparselu(const sparsematrix &a, const ae_int_t pivottype, integer_1d_array &p, integer_1d_array &q, const xparams _xparams = alglib::xdefault);
+
+
 /*************************************************************************
 Sparse Cholesky decomposition for skyline matrixm using in-place algorithm
 without allocating additional storage.
@@ -6693,6 +6847,7 @@ void sparsetrsv(sparsematrix* s,
      /* Real    */ ae_vector* x,
      ae_state *_state);
 void sparseresizematrix(sparsematrix* s, ae_state *_state);
+void sparseinitduidx(sparsematrix* s, ae_state *_state);
 double sparsegetaveragelengthofchain(sparsematrix* s, ae_state *_state);
 ae_bool sparseenumerate(sparsematrix* s,
      ae_int_t* t0,
@@ -6718,6 +6873,12 @@ void sparsegetcompressedrow(sparsematrix* s,
      ae_state *_state);
 void sparsetransposesks(sparsematrix* s, ae_state *_state);
 void sparsetransposecrs(sparsematrix* s, ae_state *_state);
+void sparsecopytransposecrs(sparsematrix* s0,
+     sparsematrix* s1,
+     ae_state *_state);
+void sparsecopytransposecrsbuf(sparsematrix* s0,
+     sparsematrix* s1,
+     ae_state *_state);
 void sparseconvertto(sparsematrix* s0, ae_int_t fmt, ae_state *_state);
 void sparsecopytobuf(sparsematrix* s0,
      ae_int_t fmt,
@@ -7156,6 +7317,60 @@ void cmatrixsyrk(ae_int_t n,
      ae_bool isupper,
      ae_state *_state);
 #endif
+#if defined(AE_COMPILE_DLU) || !defined(AE_PARTIAL_BUILD)
+void cmatrixluprec(/* Complex */ ae_matrix* a,
+     ae_int_t offs,
+     ae_int_t m,
+     ae_int_t n,
+     /* Integer */ ae_vector* pivots,
+     /* Complex */ ae_vector* tmp,
+     ae_state *_state);
+void rmatrixluprec(/* Real    */ ae_matrix* a,
+     ae_int_t offs,
+     ae_int_t m,
+     ae_int_t n,
+     /* Integer */ ae_vector* pivots,
+     /* Real    */ ae_vector* tmp,
+     ae_state *_state);
+void cmatrixplurec(/* Complex */ ae_matrix* a,
+     ae_int_t offs,
+     ae_int_t m,
+     ae_int_t n,
+     /* Integer */ ae_vector* pivots,
+     /* Complex */ ae_vector* tmp,
+     ae_state *_state);
+void rmatrixplurec(/* Real    */ ae_matrix* a,
+     ae_int_t offs,
+     ae_int_t m,
+     ae_int_t n,
+     /* Integer */ ae_vector* pivots,
+     /* Real    */ ae_vector* tmp,
+     ae_state *_state);
+#endif
+#if defined(AE_COMPILE_SPTRF) || !defined(AE_PARTIAL_BUILD)
+ae_bool sptrflu(sparsematrix* a,
+     ae_int_t pivottype,
+     /* Integer */ ae_vector* pr,
+     /* Integer */ ae_vector* pc,
+     sluv2buffer* buf,
+     ae_state *_state);
+void _sluv2list1matrix_init(void* _p, ae_state *_state, ae_bool make_automatic);
+void _sluv2list1matrix_init_copy(void* _dst, void* _src, ae_state *_state, ae_bool make_automatic);
+void _sluv2list1matrix_clear(void* _p);
+void _sluv2list1matrix_destroy(void* _p);
+void _sluv2sparsetrail_init(void* _p, ae_state *_state, ae_bool make_automatic);
+void _sluv2sparsetrail_init_copy(void* _dst, void* _src, ae_state *_state, ae_bool make_automatic);
+void _sluv2sparsetrail_clear(void* _p);
+void _sluv2sparsetrail_destroy(void* _p);
+void _sluv2densetrail_init(void* _p, ae_state *_state, ae_bool make_automatic);
+void _sluv2densetrail_init_copy(void* _dst, void* _src, ae_state *_state, ae_bool make_automatic);
+void _sluv2densetrail_clear(void* _p);
+void _sluv2densetrail_destroy(void* _p);
+void _sluv2buffer_init(void* _p, ae_state *_state, ae_bool make_automatic);
+void _sluv2buffer_init_copy(void* _dst, void* _src, ae_state *_state, ae_bool make_automatic);
+void _sluv2buffer_clear(void* _p);
+void _sluv2buffer_destroy(void* _p);
+#endif
 #if defined(AE_COMPILE_MATGEN) || !defined(AE_PARTIAL_BUILD)
 void rmatrixrndorthogonal(ae_int_t n,
      /* Real    */ ae_matrix* a,
@@ -7251,6 +7466,11 @@ void spdmatrixcholeskyupdatefixbuf(/* Real    */ ae_matrix* a,
      /* Boolean */ ae_vector* fix,
      /* Real    */ ae_vector* bufr,
      ae_state *_state);
+ae_bool sparselu(sparsematrix* a,
+     ae_int_t pivottype,
+     /* Integer */ ae_vector* p,
+     /* Integer */ ae_vector* q,
+     ae_state *_state);
 ae_bool sparsecholeskyskyline(sparsematrix* a,
      ae_int_t n,
      ae_bool isupper,