Index: trunk/npstat/stat/randomTukeyDepth.icc
===================================================================
--- trunk/npstat/stat/randomTukeyDepth.icc	(revision 904)
+++ trunk/npstat/stat/randomTukeyDepth.icc	(revision 905)
@@ -1,177 +1,238 @@
 #include <utility>
 #include <vector>
 #include <cmath>
 #include <algorithm>
 
 #include "npstat/stat/DistributionsND.hh"
 #include "npstat/stat/StatUtils.hh"
+#include "npstat/stat/StatAccumulatorArr.hh"
 
 namespace npstat {
     template<typename Real>
     void randomTukeyDepth1(const Matrix<Real>& sample,
                            AbsRandomGenerator& rng,
                            const Matrix<double>& covmat,
                            const unsigned nRandom,
-                           double* depth, const unsigned nDepth)
+                           double* depth, const unsigned nDepth,
+                           const bool usePointDirections)
     {
         typedef std::pair<double,unsigned> Pair;
 
         const unsigned dim = sample.nColumns();
         assert(dim);
         const unsigned n = sample.nRows();
         assert(n);
         assert(nDepth >= n);
         assert(depth);
 
         std::vector<double> vec(dim, 0.0);
         const GaussND gauss(&vec[0], dim, covmat);
 
         std::vector<unsigned> minOrder(n, n);
         std::vector<Pair> ordering(n);
 
+        std::vector<double> sampleMean(dim, 0.0);
+        const bool reallyUsePointDirs = usePointDirections && dim > 1U && n > 2U;
+        if (reallyUsePointDirs)
+        {
+            StatAccumulatorArr acc(dim);
+            for (unsigned i=0; i<n; ++i)
+                acc.accumulate(sample[i], dim);
+            for (unsigned i=0; i<dim; ++i)
+                sampleMean[i] = acc.mean(i);
+        }
+
         const unsigned nOrderings = nRandom + dim;
-        for (unsigned iord=0; iord<nOrderings; ++iord)
+        const unsigned nOrdTotal = reallyUsePointDirs ? nOrderings + n : nOrderings;
+        for (unsigned iord=0; iord<nOrdTotal; ++iord)
         {
             if (iord < dim)
             {
                 for (unsigned ipt=0; ipt<n; ++ipt)
                     ordering[ipt] = Pair(sample[ipt][iord], ipt);
             }
             else
             {
-                // Generate random direction
-                gauss.random(rng, &vec[0], dim);
+                if (iord < nOrderings)
+                {
+                    // Generate random direction
+                    gauss.random(rng, &vec[0], dim);
+                }
+                else
+                {
+                    // Take the direction from the center
+                    // to the point
+                    const double* pt = sample[iord - nOrderings];
+                    for (unsigned i=0; i<dim; ++i)
+                        vec[i] = pt[i] - sampleMean[i];
+                }
+
+                // Normalize the direction vector
                 long double vsum = 0.0L;
                 for (unsigned i=0; i<dim; ++i)
                     vsum += vec[i]*vec[i];
                 const double norm = std::sqrt(static_cast<double>(vsum));
-                assert(norm > 0.0);
-                for (unsigned i=0; i<dim; ++i)
-                    vec[i] /= norm;
+                if (norm > 0.0)
+                {
+                    for (unsigned i=0; i<dim; ++i)
+                        vec[i] /= norm;
+                }
+                else
+                    vec[0] = 1.0;
 
                 // Project the points
                 for (unsigned ipt=0; ipt<n; ++ipt)
                 {
                     vsum = 0.0L;
                     const Real* row = sample[ipt];
                     for (unsigned i=0; i<dim; ++i)
                         vsum += vec[i] * *row++;
                     ordering[ipt] = Pair(vsum, ipt);
                 }
             }
 
             // Order the points and see if this time
             // we get minimal ordering for some points
             std::sort(ordering.begin(), ordering.end());
             for (unsigned ipt=0, ipt2=n-1U; ipt<n; ++ipt, --ipt2)
             {
                 const unsigned idx = ordering[ipt].second;
                 const unsigned minIpt = ipt < ipt2 ? ipt : ipt2;
                 if (minIpt < minOrder[idx])
                     minOrder[idx] = minIpt;
             }
         }
 
         if (n > 2U)
         {
             unsigned maxPossibleOrder = n/2U;
             if (n % 2U == 0U)
                 --maxPossibleOrder;
             const double med = maxPossibleOrder;
             for (unsigned ipt=0; ipt<n; ++ipt)
                 depth[ipt] = minOrder[ipt]/med;
         }
         else
         {
             for (unsigned ipt=0; ipt<n; ++ipt)
                 depth[ipt] = 0.0;
         }
     }
 
     template<typename Real1, typename Real2>
     void randomTukeyDepth2(const Matrix<Real1>& inSample,
                            const Matrix<Real2>& refSample,
                            AbsRandomGenerator& rng,
                            const Matrix<double>& covmat,
                            const unsigned nRandom,
-                           double* depth, const unsigned nDepth)
+                           double* depth, const unsigned nDepth,
+                           const bool usePointDirections)
     {
         const unsigned dim = inSample.nColumns();
         assert(dim);
         assert(dim == refSample.nColumns());
         const unsigned nIn = inSample.nRows();
         assert(nIn);
         assert(nDepth >= nIn);
         assert(depth);
         const unsigned nRef = refSample.nRows();
         assert(nRef);
 
         std::vector<double> vec(dim, 0.0);
         const GaussND gauss(&vec[0], dim, covmat);
 
         std::vector<double> minCdf(nIn, 2.0);
         std::vector<double> ordering(nRef);
 
+        std::vector<double> sampleMean(dim, 0.0);
+        const bool reallyUsePointDirs = usePointDirections && dim > 1U;
+        if (reallyUsePointDirs)
+        {
+            StatAccumulatorArr acc(dim);
+            for (unsigned i=0; i<nRef; ++i)
+                acc.accumulate(refSample[i], dim);
+            for (unsigned i=0; i<dim; ++i)
+                sampleMean[i] = acc.mean(i);
+        }
+
         const unsigned nOrderings = nRandom + dim;
-        for (unsigned iord=0; iord<nOrderings; ++iord)
+        const unsigned nOrdTotal = reallyUsePointDirs ? nOrderings + nIn : nOrderings;
+        for (unsigned iord=0; iord<nOrdTotal; ++iord)
         {
             if (iord < dim)
             {
                 for (unsigned ipt=0; ipt<nRef; ++ipt)
                     ordering[ipt] = refSample[ipt][iord];
             }
             else
             {
-                // Generate random direction
-                gauss.random(rng, &vec[0], dim);
+                if (iord < nOrderings)
+                {
+                    // Generate random direction
+                    gauss.random(rng, &vec[0], dim);
+                }
+                else
+                {
+                    // Take the direction from the center
+                    // to the point
+                    const double* pt = inSample[iord - nOrderings];
+                    for (unsigned i=0; i<dim; ++i)
+                        vec[i] = pt[i] - sampleMean[i];
+                }
+
+                // Normalize the direction vector
                 long double vsum = 0.0L;
                 for (unsigned i=0; i<dim; ++i)
                     vsum += vec[i]*vec[i];
                 const double norm = std::sqrt(static_cast<double>(vsum));
-                assert(norm > 0.0);
-                for (unsigned i=0; i<dim; ++i)
-                    vec[i] /= norm;
+                if (norm > 0.0)
+                {
+                    for (unsigned i=0; i<dim; ++i)
+                        vec[i] /= norm;
+                }
+                else
+                    vec[0] = 1.0;
 
                 // Project the reference points
                 for (unsigned ipt=0; ipt<nRef; ++ipt)
                 {
                     vsum = 0.0L;
                     const Real2* row = refSample[ipt];
                     for (unsigned i=0; i<dim; ++i)
                         vsum += vec[i] * *row++;
                     ordering[ipt] = vsum;
                 }
             }
 
             // Order the reference sample
             std::sort(ordering.begin(), ordering.end());
 
             // Go over the input sample an check what
             // would be the order of points in that sample
             for (unsigned ipt=0; ipt<nIn; ++ipt)
             {
                 double pointValue;
                 if (iord < dim)
                     pointValue = inSample[ipt][iord];
                 else
                 {
                     long double vsum = 0.0L;
                     const Real1* row = inSample[ipt];
                     for (unsigned i=0; i<dim; ++i)
                         vsum += vec[i] * *row++;
                     pointValue = vsum;
                 }
 
                 double cdf = empiricalCdf(ordering, pointValue);
                 assert(0.0 <= cdf && cdf <= 1.0);
                 if (1.0 - cdf < cdf)
                     cdf = 1.0 - cdf;
                 if (cdf < minCdf[ipt])
                     minCdf[ipt] = cdf;
             }
         }
 
         for (unsigned ipt=0; ipt<nIn; ++ipt)
             depth[ipt] = 2.0*minCdf[ipt];
     }
 }
Index: trunk/npstat/stat/randomTukeyDepth.hh
===================================================================
--- trunk/npstat/stat/randomTukeyDepth.hh	(revision 904)
+++ trunk/npstat/stat/randomTukeyDepth.hh	(revision 905)
@@ -1,89 +1,118 @@
 #ifndef NPSTAT_RANDOMTUKEYDEPTH_HH_
 #define NPSTAT_RANDOMTUKEYDEPTH_HH_
 
 /*!
 // \file randomTukeyDepth.hh
 //
-// \brief Implementation of Tukey depth using random directions
+// \brief Implementation of Tukey depth using random directions and/or
+//        directions from the center of the cloud to the sample points
 //
 // Author: I. Volobouev
 //
 // March 2023
 */
 
 #include "npstat/nm/Matrix.hh"
 #include "npstat/rng/AbsRandomGenerator.hh"
 
+// For a detailed description of the random Tukey depth concept, see
+// Cuesta-Albertos, J. A. and Nieto-Reyes, A., "The random Tukey depth",
+// Computational Statistics & Data Analysis 52, 11 (July 2008), 4979-4988.
 namespace npstat {
     /**
     // Random Tukey depth calculated within the sample.
     // Function arguments are as follows:
     //
     // sample   should be dimensioned n x d, where d is the
     //          point dimensionality and n is the sample size.
     //
     // rng      random number generator in 1 or d dimensions.
     //
     // covmat   should be dimensioned d x d. The random directions
     //          will be generated by a multivariate normal
     //          distribution with this covariance matrix.
     //
     // nRandom  is the number of random directions to generate.
     //          In addition to random directions, the ordering
     //          will also be perfomed according to the marginals.
     //
     // depth    this array will be filled with depth values
     //          on exit (in the same order of points as in
     //          the "sample" argument).
     //
     // nDepth   is the length of the "depth" array. Should be
     //          at least as large as the sample size.
+    //
+    // usePointDirections  if this argument is set to "true",
+    //                     directions from the center of the
+    //                     cloud to each sample point will also
+    //                     be utilized in the depth calculation.
+    //
+    // The CPU time used by this function will scale as
+    // O(nRandom*n*log(n)) if "usePointDirections" is false.
+    // If "usePointDirections" is true, another step will
+    // be added that scales as O(n^2*log(n)).
     */
     template<typename Real>
     void randomTukeyDepth1(const Matrix<Real>& sample,
                            AbsRandomGenerator& rng,
                            const Matrix<double>& covmat,
                            unsigned nRandom,
-                           double* depth, unsigned nDepth);
+                           double* depth, unsigned nDepth,
+                           bool usePointDirections = false);
 
     /**
     // Random Tukey depth for each point in a sample calculated
     // using another (reference) sample. Function arguments are
     // as follows:
     //
     // sample   should be dimensioned n x d, where d is the
     //          point dimensionality and n is the sample size.
     //
-    // refSample  the depth of each point in the "sample"
-    //            argument will be defined with respect
-    //            to this sample.
+    // refSample  the depth of each point in the "sample" argument
+    //            will be defined with respect to this "reference"
+    //            sample. The size of this sample should normally
+    //            be substantialy larger than n.
     //
     // rng      random number generator in 1 or d dimensions.
     //
     // covmat   should be dimensioned d x d. The random directions
     //          will be generated by a multivariate normal
     //          distribution with this covariance matrix.
     //
     // nRandom  is the number of random directions to generate.
     //          In addition to random directions, the ordering
     //          will also be perfomed according to the marginals.
     //
     // depth    this array will be filled with depth values
     //          on exit (in the same order of points as in
     //          the "sample" argument).
     //
     // nDepth   is the length of the "depth" array. Should be
     //          at least as large as the "sample" size.
+    //
+    // usePointDirections  if this argument is set to "true",
+    //                     directions from the center of the
+    //                     "refSample" cloud to each point of
+    //                     the "sample" will also be utilized
+    //                     in the depth calculation.
+    //
+    // Assuming that the size of the reference sample is m,
+    // the CPU time used by this function will scale as
+    // O(nRandom*(m+n)*log(m)) if "usePointDirections" is false.
+    // If "usePointDirections" is true, another step will
+    // be added that scales as O(n*(m+n)*log(m)).
     */
     template<typename Real1, typename Real2>
     void randomTukeyDepth2(const Matrix<Real1>& sample,
                            const Matrix<Real2>& refSample,
                            AbsRandomGenerator& rng,
                            const Matrix<double>& covmat,
                            unsigned nRandom,
-                           double* depth, unsigned nDepth);
+                           double* depth, unsigned nDepth,
+                           bool usePointDirections = false);
 }
 
 #include "npstat/stat/randomTukeyDepth.icc"
 
 #endif // NPSTAT_RANDOMTUKEYDEPTH_HH_