Changeset 983

Timestamp:

10/24/08 15:58:30 (2 months ago)

Author:

lebrun

Message:

Promoted some NumericalPoint into NumericalPointWithDescription that were missed during the separation between NumericalPoint and Description into the getParameters() method of several distributions. This solves tickets #155.
Changed the return type of the getImportanceFactors() method of the QuadraticCumul class. This solves ticket #156.
Added a simplified constructor from a String to the class Description. It closes ticket #108.
Removed the dependence to rotRPackage for the Kolmogorov() method of the FittingTest class. It greatly improves both the performance and the generality of this method.
Improved the const-correctness of many methods.

Files:

• branches/lebrun/devel/lib/src/Base/Type/Description.cxx (modified) (1 diff)
• branches/lebrun/devel/lib/src/Base/Type/Description.hxx (modified) (1 diff)
• branches/lebrun/devel/lib/src/Base/Type/DescriptionImplementation.cxx (modified) (1 diff)
• branches/lebrun/devel/lib/src/Base/Type/DescriptionImplementation.hxx (modified) (1 diff)
• branches/lebrun/devel/lib/src/Uncertainty/Algorithm/Analytical/AnalyticalResult.cxx (modified) (17 diffs)
• branches/lebrun/devel/lib/src/Uncertainty/Algorithm/IsoProbabilisticTransformation/MarginalTransformation/InverseMarginalTransformationEvaluation.cxx (modified) (1 diff)
• branches/lebrun/devel/lib/src/Uncertainty/Algorithm/IsoProbabilisticTransformation/MarginalTransformation/MarginalTransformationEvaluation.cxx (modified) (3 diffs)
• branches/lebrun/devel/lib/src/Uncertainty/Algorithm/QuadraticCumul/QuadraticCumul.cxx (modified) (10 diffs)
• branches/lebrun/devel/lib/src/Uncertainty/Algorithm/QuadraticCumul/QuadraticCumul.hxx (modified) (1 diff)
• branches/lebrun/devel/lib/src/Uncertainty/Distribution/ComposedDistribution.cxx (modified) (5 diffs)
• branches/lebrun/devel/lib/src/Uncertainty/Distribution/Epanechnikov.cxx (modified) (1 diff)
• branches/lebrun/devel/lib/src/Uncertainty/Distribution/ExtraFunc/DistFunc.cxx (modified) (2 diffs)
• branches/lebrun/devel/lib/src/Uncertainty/Distribution/ExtraFunc/DistFunc.hxx (modified) (1 diff)
• branches/lebrun/devel/lib/src/Uncertainty/Distribution/IndependentCopula.cxx (modified) (1 diff)
• branches/lebrun/devel/lib/src/Uncertainty/Distribution/KernelMixture.cxx (modified) (1 diff)
• branches/lebrun/devel/lib/src/Uncertainty/Distribution/Mixture.cxx (modified) (3 diffs)
• branches/lebrun/devel/lib/src/Uncertainty/Distribution/MultiNomial.cxx (modified) (14 diffs)
• branches/lebrun/devel/lib/src/Uncertainty/Distribution/NormalCopula.cxx (modified) (16 diffs)
• branches/lebrun/devel/lib/src/Uncertainty/Distribution/RandomMixture.cxx (modified) (1 diff)
• branches/lebrun/devel/lib/src/Uncertainty/Distribution/StudentND.cxx (modified) (1 diff)
• branches/lebrun/devel/lib/src/Uncertainty/Model/DistributionImplementation.cxx (modified) (2 diffs)
• branches/lebrun/devel/lib/src/Uncertainty/Model/EllipticalDistribution.cxx (modified) (10 diffs)
• branches/lebrun/devel/lib/src/Uncertainty/StatTests/FittingTest.cxx (modified) (3 diffs)
• branches/lebrun/devel/lib/test/t_ComposedDistribution_std.at (modified) (2 diffs)
• branches/lebrun/devel/lib/test/t_FORM_sensitivity.at (modified) (1 diff)
• branches/lebrun/devel/lib/test/t_FittingTest_std.at (modified) (1 diff)
• branches/lebrun/devel/lib/test/t_FittingTest_std.cxx (modified) (3 diffs)
• branches/lebrun/devel/lib/test/t_Mixture_std.at (modified) (1 diff)
• branches/lebrun/devel/lib/test/t_QuadraticCumul_importanceFactors.at (modified) (1 diff)
• branches/lebrun/devel/python/test/t_ComposedDistribution_std.atpy (modified) (2 diffs)
• branches/lebrun/devel/python/test/t_FittingTest_std.atpy (modified) (1 diff)
• branches/lebrun/devel/python/test/t_FittingTest_std.py (modified) (4 diffs)
• branches/lebrun/devel/python/test/t_Mixture_std.atpy (modified) (1 diff)
• branches/lebrun/devel/python/test/t_QuadraticCumul_importanceFactors.atpy (modified) (1 diff)

branches/lebrun/devel/lib/src/Base/Type/Description.cxx

@@ -53 +53 @@
-
-
+      /* Constructor with value */
+      Description::Description(const String & value)
+        : Common::TypedCollectionInterfaceObject<DescriptionImplementation>(new DescriptionImplementation(value))
+      {
+        // Nothing to do
+      }
+
+
-      /* Constructor with size and default value */
-      Description::Description(const UnsignedLong size,

branches/lebrun/devel/lib/src/Base/Type/Description.hxx

@@ -57 +57 @@
-        Description(const UnsignedLong size);
-
+        /** Constructor with value */
+        Description(const String & value);
+
-        /** Constructor with size and default value */
-        Description(const UnsignedLong size,

branches/lebrun/devel/lib/src/Base/Type/DescriptionImplementation.cxx

@@ -63 +63 @@
-
-
+      /* Constructor with size */
+      DescriptionImplementation::DescriptionImplementation(const String & value)
+        : PersistentCollection<String>(1, value)
+      {
+        // Nothing to do
+      }
+
+
-      /* Constructor with size and default value */
-      DescriptionImplementation::DescriptionImplementation(const UnsignedLong size,

branches/lebrun/devel/lib/src/Base/Type/DescriptionImplementation.hxx

@@ -59 +59 @@
-        /** Constructor with size */
-        DescriptionImplementation(const UnsignedLong size);     
+
+        /** Constructor with default value */
+        DescriptionImplementation(const String & value);        
-        
-        /** Constructor with size and default value */

branches/lebrun/devel/lib/src/Uncertainty/Algorithm/Analytical/AnalyticalResult.cxx

@@ -81 +81 @@
-                                         const String & name):
-        PersistentObject(name),
-NumericalPoint(standardSpaceDesignPoint.getDimension(), 0.0
+standardSpaceDesignPoint.getDimension(
-        limitStateVariable_(limitStateVariable),
-        isStandardPointOriginInFailureSpace_(isStandardPointOriginInFailureSpace),
-
-NumericalPoint(standardSpaceDesignPoint.getDimension(), 0.0
-Sensitivity(0)
+0
+standardSpaceDesignPoint.getDimension(
+0
-        isAlreadyComputedImportanceFactors_(false),
-        isAlreadyComputedHasoferReliabilityIndexSensitivity_(false)
@@ -97 +97 @@
-      AnalyticalResult::AnalyticalResult():
-        PersistentObject(),
-NumericalPoint(0)
-NumericalPoint(0)
+0
+0
-        // Fake event based on a fake 1D composite random vector, which requires a fake 1D NumericalMathFunction
-        limitStateVariable_(RandomVector(new CompositeRandomVector(NumericalMathFunction(Description(0),Description(1),Description(1)),
-                                                                                 new ConstantRandomVector(NumericalPoint(0)))),
-
+0
-        isStandardPointOriginInFailureSpace_(false),
-
-NumericalPoint(0)
-Sensitivity(0)
+0
+0
+0
-        isAlreadyComputedImportanceFactors_(false),
-        isAlreadyComputedHasoferReliabilityIndexSensitivity_(false)
@@ -179 +179 @@
-      void AnalyticalResult::computeImportanceFactors()
-      {
-
+const 
-        importanceFactors_ = NumericalPoint(dimension, -1.0);
-        /* First, check that the importance factors are well-defined */
@@ -185 +185 @@
-          {
-            /* get the input distribution */
-
+const 
-            /* get the input standard distribution */
-
+const 
-            /* get the marginal 1D from in the standard space where all marginals are identical */
-
+const 
-
-            /* evaluate the corresponding vector in the importance factors space */
@@ -195 +195 @@
-            /* if Rosenblatt  Transformation : Z = Phi^(-1)o F_{Xi}(Xi*) */
-            // for each marginals */
-marginalIndex++
-
-NumericalScalar y = inputDistribution.getMarginal(marginalIndex).computeCDF(NumericalPoint(1, physicalSpaceDesignPoint_[marginalIndex]
+++marginalIndex
+
+const NumericalScalar y(inputDistribution.getMarginal(marginalIndex).computeCDF(NumericalPoint(1, physicalSpaceDesignPoint_[marginalIndex])
-                importanceFactors_[marginalIndex] = standardMarginalDistribution.computeQuantile(y)[0];
-              }
-            /* compute  square norm of importanceFactors vector */
-
+const 
-            /* compute final importance factors */
-marginalIndex++
+++marginalIndex
-              {
-                importanceFactors_[marginalIndex] *= importanceFactors_[marginalIndex] * inorm2;
@@ -235 +235 @@
-            computeImportanceFactors();
-          }
-
+const 
-        NumericalSample data(dimension, 1);
-
-        /* build data for the pie */
-i++
+++i
-          {
-            data[i]=NumericalPoint(1, importanceFactors_[i]);
@@ -249 +249 @@
-        Description palette(dimension);
-        Description labels(dimension);
-
+const 
-        Description::const_iterator colorsIterator;
-        colorsIterator = colors.begin();
@@ -257 +257 @@
-          {
-            description = Description(dimension);
-i++
+++i
-              {
-                OSS oss;
@@ -264 +264 @@
-              }
-          }
-i++
+++i
-          {
-            OSS oss;
@@ -272 +272 @@
-            labels[i] = oss;
-            palette[i] = (*colorsIterator);
-colorsIterator++
+++colorsIterator
-            if (colorsIterator == colors.end())
-              {
@@ -325 +325 @@
-      {
-        /* get Set1 : parameters of the physical distribution */
-
-
+const 
+const 
-        /* get Set2 : parameters of the physical model */
-
-NumericalPoint
-
+const 
+const NumericalPointWithDescription
+const 
-
-        /* get SetIso : parameters included in the isoprobabilistic transformation which is a subset of Set1 */
-
-
+const 
+const 
-        /* scaling factor between sensitivity and gradients (ref doc : factor = -lambda/beta) */
-        /* gradient of the standard limite state function */
-
-
-
+const 
+const 
+const 
-        /* evaluate the gradients of the physical model with respect to Set2 (ref doc : K2) */
-        NumericalPoint physicalGradient;
@@ -350 +350 @@
-        if (setIso.getDimension() > 0)
-          {
-
+ 
-          }
-        /* associate to each element of Set1 the gradient value */
-
-        /* hasoferReliabilityIndexSensitivity is the collection Set1 + one other collection wich is Set2 */
-
-
+const 
+const 
-        hasoferReliabilityIndexSensitivity_ = Sensitivity(size);
-
-sensitivityIndex++
-
-
-
+++sensitivityIndex
+
+const 
+const 
-            NumericalPointWithDescription currentSensitivity(currentDimension);
-
+const 
-            // the currentSensitivity gets the description and the name from set1
-            currentSensitivity.setDescription(currentDescription);
-
+const 
-            currentSensitivity.setName(currentName);
-
-currentIndex++
-
-
+const 
+++currentIndex
+
+const 
-                /* if currentParameters[currentIndex] is into setIso, then we get the sensitivity value in the corresponding isoProbabilisticGradient */
-                if (position < setIso.getDimension())
@@ -396 +396 @@
-      UnsignedLong AnalyticalResult::computePosition(const String & marginalName, const String & marginalParameterName, const Description & parameterSetNames) const
-      {
-
-
-index++
+const 
+const 
+++index
-          {
-            if ( parameterSetNames[index] == fullName ) return index;
@@ -416 +416 @@
-            computeHasoferReliabilityIndexSensitivity();
-          }
-
-
+const 
+const 
-        // The first graph shows the sensitivities with respect to the marginal parameters if there exist
-        // in the cas where the distribution or some marginals are defined by user, it may not have parameters : we create a empty sensitivity graph
-        Sensitivity marginalSensitivity(dimension);
-i++
+++i
-          {
-            marginalSensitivity[i] = hasoferReliabilityIndexSensitivity_[i];
@@ -434 +434 @@
-          {
-            Sensitivity otherSensitivity(size - dimension);
-i++
+++i
-              {
-                otherSensitivity[i - dimension] = hasoferReliabilityIndexSensitivity_[i];
@@ -458 +458 @@
-        BarPlot sensitivityBarPlot(NumericalSample(0, 2), shift, "");
-
-
+const 
-        Description::const_iterator colorsIterator;
-        colorsIterator = colors.begin();
-
-        // Create the barplots
-
-collectionIndex++
+const 
+++collectionIndex
-          {
-            NumericalSample data(sensitivity[collectionIndex].getDimension(), 2);
-
-sensitivityIndex++
+const 
+++sensitivityIndex
-              {
-                data[sensitivityIndex][0] = width;
@@ -477 +477 @@
-            oss << sensitivity[collectionIndex].getName() << " " << sensitivity[collectionIndex].getDescription();
-            sensitivityGraph.addDrawable(BarPlot(data, shift, (*colorsIterator), "solid", "solid", oss));
-colorsIterator++
+++colorsIterator
-            if (colorsIterator == colors.end())
-              {

branches/lebrun/devel/lib/src/Uncertainty/Algorithm/IsoProbabilisticTransformation/MarginalTransformation/InverseMarginalTransformationEvaluation.cxx

@@ -90 +90 @@
-      InverseMarginalTransformationEvaluation::Matrix InverseMarginalTransformationEvaluation::parametersGradient(const NumericalPoint & in) const
-      {
-
-
-
+const 
+const 
+const 
-        Matrix result(parametersDimension, inputDimension);
-        UnsignedLong rowIndex(0);
-        for (UnsignedLong j = 0; j < inputDimension; ++j)
-          {
-
-
-
-
-
+const 
+const 
+const 
+const 
+const 
-            for (UnsignedLong i = 0; i < marginalParametersDimension; ++i)
-              {

branches/lebrun/devel/lib/src/Uncertainty/Algorithm/IsoProbabilisticTransformation/MarginalTransformation/MarginalTransformationEvaluation.cxx

@@ -42 +42 @@
-      {
-        Description description;
-
+const 
-        for (UnsignedLong i = 0; i < size; ++i)
-          {
@@ -68 +68 @@
-        throw (InvalidArgumentException, InternalException)
-      {
-
+const 
-        NumericalPoint result(dimension);
-        // Apply CDF over the components
@@ -81 +81 @@
-      MarginalTransformationEvaluation::Matrix MarginalTransformationEvaluation::parametersGradient(const NumericalPoint & in) const
-      {
-
-
-
+const 
+const 
+const 
-        Matrix result(parametersDimension, inputDimension);
-        UnsignedLong rowIndex(0);
-        for (UnsignedLong j = 0; j < inputDimension; ++j)
-          {
-
-
+const 
+const 
-            for (UnsignedLong i = 0; i < marginalParametersDimension; ++i)
-              {

branches/lebrun/devel/lib/src/Uncertainty/Algorithm/QuadraticCumul/QuadraticCumul.cxx

@@ -63 +63 @@
-        : Base::Common::PersistentObject(name),
-          limitStateVariable_(limitStateVariable),
-NumericalPoint(0)
-NumericalPoint(0)
-Matrix(0, 0)
-SymmetricTensor(0, 0)
+0
+0
+0, 0
+0, 0
-          isAlreadyComputedValue_(false),
-          isAlreadyComputedGradient_(false),
@@ -74 +74 @@
-          isAlreadyComputedCovariance_(false),
-          isAlreadyComputedImportanceFactors_(false),
-CovarianceMatrix(0)
-NumericalPoint(0)
-NumericalPoint(0)
-CovarianceMatrix(0)
-NumericalPoint(0)
+0
+0
+0
+0
+0
-      {
-        /* Check if the given random vector is a composite random vector, which is mandatory */
@@ -153 +153 @@
-
-      /* importance factors accessor */
-
+WithDescription
-      {
-        if(!isAlreadyComputedImportanceFactors_)
@@ -168 +168 @@
-      {
-        // To ensure that the importance factors are up to date
-importanceFactors_ = 
-
+
+const 
-        NumericalSample data(dimension, 1);
-
-        /* build data for the pie */
-i++
+++i
-          {
-            data[i]=NumericalPoint(1, importanceFactors_[i]);
@@ -186 +186 @@
-        Description::const_iterator colorsIterator;
-        colorsIterator = colors.begin();
-i++
+++i
-          {
-            OSS oss;
@@ -193 +193 @@
-            oss << 100.0 * data[i][0] << "%";
-            palette[i] = (*colorsIterator);
-colorsIterator++
+++colorsIterator
-            if (colorsIterator == colors.end())
-              {
@@ -254 +254 @@
-        covariance_ = CovarianceMatrix(gradientAtMean_.getNbColumns());
-        /* we unroll the matrix multiplications transpose(gradient).covariance.gradient */
-i++
-
-j++
+++i
+
+++j
-              {
-                covariance_(i, j) = 0;
-k++
+++k
-                  {
-l++
+++l
-                      {
-                        covariance_(i, j) += gradientAtMean_(l, i) * inputCovariance_(l, k) * gradientAtMean_(k, j);
@@ -283 +283 @@
-        /* we compute here the importance factors */
-        /* in this scalar case, gradientAtMean is a NumericalPoint */
-
+const 
-        
-        /* in this scalar case, importance factors is a NumericalPoint, defined as importanceFactors = gradient .* inputCovariance * gradient / outputCovariance, where .* means an element-wise multiplication between vectors */
-        importanceFactors_ = NumericalPoint(dimension, 0.0);
-        // This line looks strange, but it is here to ensure that the covariance has actually been computed
-covariance_ = 
-i++
-
-j++
+
+++i
+
+++j
-              {
-                importanceFactors_[i] += inputCovariance_(i, j) * gradientAtMean_(j, 0);
@@ -298 +298 @@
-          }
-        importanceFactors_.setDescription(limitStateVariable_.getImplementation()->getAntecedent()->getDescription());
-
-        isAlreadyComputedImportanceFactors_ = true;
-      } // QuadraticCumul::computeImportanceFactors()
@@ -323 +322 @@
-        /* developped formula */
-
-
+const 
-        /* i */
-
+const 
-        /* k */
-        meanSecondOrder_ = valueAtMean_;
-        /* loop on k */
-k++
+++k
-          {
-            NumericalScalar kSecondOrderContribution(0.0);
-            /* loop on i */
-i++
+++i
-              {
-                kSecondOrderContribution += 0.5 * inputCovariance_(i, i) * hessianAtMean_(i, i, k);
-                /* loop on j */
-j++
+++j
-                  {
-                    kSecondOrderContribution += inputCovariance_(i, j) * hessianAtMean_(i, j, k);

branches/lebrun/devel/lib/src/Uncertainty/Algorithm/QuadraticCumul/QuadraticCumul.hxx

@@ -104 +104 @@
-
-        /** importance factors accessor */
-
+WithDescription
-
-        /** ImportanceFactors graph */

branches/lebrun/devel/lib/src/Uncertainty/Distribution/ComposedDistribution.cxx

@@ -616 +616 @@
-                parameters[index] = marginalParameters[j];
-                description[index] = OSS() << marginalName << "_" << marginalDescription[j];
-index++
+++index
-              }
-          }
@@ -633 +633 @@
-        IsoProbabilisticTransformation right(p_rightFunction, p_rightGradient, p_rightHessian);
-        right.setParameters(parameters);
-
-        IsoProbabilisticTransformation transformation(isoprobabilistic, right);
-
@@ -664 +663 @@
-                parameters[index] = marginalParameters[j];
-                description[index] = OSS() << marginalName << "_" << marginalDescription[j];
-index++
+++index
-              }
-          }
@@ -700 +699 @@
-        const Description description(getDescription());
-        // First put the marginal parameters
-marginalIndex++
+++marginalIndex
-          {
-            // Each marginal distribution must output a collection of parameters of size 1, even if it contains an empty NumericalPoint
-            const NumericalPointWithDescriptionCollection marginalParameters(distributionCollection_[marginalIndex].getParametersCollection());
-
+WithDescription
-            point.setName(description[marginalIndex]);
-            parameters[marginalIndex] = point;
@@ -713 +712 @@
-            // Second put the dependence parameters
-            const NumericalPointWithDescriptionCollection copulaParameters(copula_.getParametersCollection());
-
+WithDescription
-            if (copulaParameters.getSize() != 0)
-              {

branches/lebrun/devel/lib/src/Uncertainty/Distribution/Epanechnikov.cxx

@@ -166 +166 @@
-      Epanechnikov::NumericalPointWithDescriptionCollection Epanechnikov::getParametersCollection() const
-      {
+        // No parameter
-        return NumericalPointWithDescriptionCollection(0);
-      }

branches/lebrun/devel/lib/src/Uncertainty/Distribution/ExtraFunc/DistFunc.cxx

@@ -25 +25 @@
- */
-#include <iostream>
+#include <cstdio>
+#include <cstdlib>
-#include <cmath>
-#include "DistFunc.hxx"
@@ -267 +269 @@
-      } // End of rGamma
-
+      /*********************************************/
+      /* Kolmogorov distribution                   */
+      /*********************************************/
+static double K(int n, double d);
+static void m_multiply(double *A, double *B, double *C, int m);
+static void m_power(double *A, int eA, double *V, int *eV, int m, int n);
+static double
+K(int n, double d)
+{
+    /* Compute Kolmogorov's distribution.
+       Code published in
+         George Marsaglia and Wai Wan Tsang and Jingbo Wang (2003),
+         "Evaluating Kolmogorov's distribution".
+         Journal of Statistical Software, Volume 8, 2003, Issue 18.
+         URL: http://www.jstatsoft.org/v08/i18/.
+    */
+
+   int k, m, i, j, g, eH, eQ;
+   double h, s, *H, *Q;
+
+   k = (int) (n * d) + 1;
+   m = 2 * k - 1;
+   h = k - n * d;
+   H = (double*) calloc(m * m, sizeof(double));
+   Q = (double*) calloc(m * m, sizeof(double));
+   for(i = 0; i < m; i++)
+       for(j = 0; j < m; j++)
+           if(i - j + 1 < 0)
+               H[i * m + j] = 0;
+           else
+               H[i * m + j] = 1;
+   for(i = 0; i < m; i++) {
+       H[i * m] -= pow(h, i + 1);
+       H[(m - 1) * m + i] -= pow(h, (m - i));
+   }
+   H[(m - 1) * m] += ((2 * h - 1 > 0) ? pow(2 * h - 1, m) : 0);
+   for(i = 0; i < m; i++)
+       for(j=0; j < m; j++)
+           if(i - j + 1 > 0)
+               for(g = 1; g <= i - j + 1; g++)
+                   H[i * m + j] /= g;
+   eH = 0;
+   m_power(H, eH, Q, &eQ, m, n);
+   s = Q[(k - 1) * m + k - 1];
+   for(i = 1; i <= n; i++) {
+       s = s * i / n;
+       if(s < 1e-140) {
+           s *= 1e140;
+           eQ -= 140;
+       }
+   }
+   s *= pow(10., eQ);
+   free(H);
+   free(Q);
+   return(s);
+}
+
+static void
+m_multiply(double *A, double *B, double *C, int m)
+{
+    /* Auxiliary routine used by K().
+       Matrix multiplication.
+    */
+    int i, j, k;
+    double s;
+    for(i = 0; i < m; i++)
+        for(j = 0; j < m; j++) {
+            s = 0.;
+            for(k = 0; k < m; k++)
+                s+= A[i * m + k] * B[k * m + j];
+            C[i * m + j] = s;
+        }
+}
+
+static void
+m_power(double *A, int eA, double *V, int *eV, int m, int n)
+{
+    /* Auxiliary routine used by K().
+       Matrix power.
+    */
+    double *B;
+    int eB , i;
+
+    if(n == 1) {
+        for(i = 0; i < m * m; i++)
+            V[i] = A[i];
+        *eV = eA;
+        return;
+    }
+    m_power(A, eA, V, eV, m, n / 2);
+    B = (double*) calloc(m * m, sizeof(double));
+    m_multiply(V, V, B, m);
+    eB = 2 * (*eV);
+    if((n % 2) == 0) {
+        for(i = 0; i < m * m; i++)
+            V[i] = B[i];
+        *eV = eB;
+    }
+    else {
+        m_multiply(A, B, V, m);
+        *eV = eA + eB;
+    }
+    if(V[(m / 2) * m + (m / 2)] > 1e140) {
+        for(i = 0; i < m * m; i++)
+            V[i] = V[i] * 1e-140;
+        *eV += 140;
+    }
+    free(B);
+}
+      NumericalScalar DistFunc::pKolmogorov(const UnsignedLong n, const NumericalScalar x, const Bool tail)
+      {
+        // If x <= 0, we are at the left of the support
+        if (x <= 0.0) return (tail ? 1.0 : 0.0);
+        // If x >= 1, we are at the right of the support
+        if (x >= 1.0) return (tail ? 0.0 : 1.0);
+//      // We use the summation formula of R. Von Mises
+//      const UnsignedLong kmax(floor(n * (1.0 - x)));
+//      NumericalScalar sum(pow(1.0 - x, n) / x);
+//      NumericalScalar logBinomial(log(n));
+//      for (UnsignedLong k = 1; k <= kmax; ++k)
+//        {
+//          const NumericalScalar term(x + k / static_cast<NumericalScalar>(n));
+//          sum += exp(logBinomial + (k - 1) * log(term) + (n - k) * log(1.0 - term));
+//          logBinomial += log(n - k) - log (k + 1);
+//        }
+//      NumericalScalar result((tail ? 2.0 * x * sum : 1.0 - 2.0 * x * sum));
+        NumericalScalar result2((tail ? 1.0 - K(n, x) : K(n, x)));
+        return result2;
+      }
+
+      NumericalScalar DistFunc::pKolmogorovAsymptotic(const UnsignedLong n, const NumericalScalar x, const Bool tail)
+      {
+        // If x <= 0, we are at the left of the support