Minor changes in GPCC

sofianehaddad · Apr 5, 2024 · 270b53a · 270b53a
1 parent d95a00b
commit 270b53a
Show file tree

Hide file tree

Showing 2 changed files with 24 additions and 43 deletions.
diff --git a/lib/src/Uncertainty/Algorithm/MetaModel/Kriging/GaussianProcessConditionalCovariance.cxx b/lib/src/Uncertainty/Algorithm/MetaModel/Kriging/GaussianProcessConditionalCovariance.cxx
@@ -1,6 +1,6 @@
 //                                               -*- C++ -*-
 /**
- *  @brief The result of a kriging estimation
+ *  @brief The postprocessing of a GPR result (conditional covariance)
  *
  *  Copyright 2005-2024 Airbus-EDF-IMACS-ONERA-Phimeca
  *
@@ -77,13 +77,8 @@ void GaussianProcessConditionalCovariance::computePhi()
   if (F.getNbColumns() == 0) return;
   // Nothing to do if the design matrix has already been computed
   Matrix Q;
-  const UnsignedInteger method = result_.getLinearAlgebraMethod();
   LOGINFO("Solve linear system  L * phi= F");
-  Matrix phi;
-  if(method == 0)
-    phi = result_.getCholeskyFactor().solveLinearSystem(F);
-  else
-    phi = result_.getHMatCholeskyFactor().solveLower(F);
+  Matrix phi(solveTriangularSystem(F));
   // Compute QR decomposition of Phi_
   LOGINFO("Compute the QR decomposition of phi");
   Matrix G;
@@ -92,6 +87,21 @@ void GaussianProcessConditionalCovariance::computePhi()
   phiT_ = phi.transpose();
 }
 
+Matrix GaussianProcessConditionalCovariance::solveTriangularSystem(const Matrix & rhs) const
+{
+  Matrix result;
+  const UnsignedInteger method = result_.getLinearAlgebraMethod();
+  if (method == 0)
+  {
+    result = result_.getCholeskyFactor().solveLinearSystem(rhs);
+  }
+  else
+  {
+    result = result_.getHMatCholeskyFactor().solveLower(rhs);
+  }
+  return result;
+}
+
 /* Compute mean of new points conditionally to observations */
 Sample GaussianProcessConditionalCovariance::getConditionalMean(const Sample & xi) const
 {
@@ -150,18 +160,8 @@ CovarianceMatrix GaussianProcessConditionalCovariance::getConditionalCovariance(
   // 3) Compute r^t R^{-1} r'(x)
   // As we get the Cholesky factor L, we can solve triangular linear system
   // We define B = L^{-1} * r(x)
-  Matrix B;
-  const UnsignedInteger method = result_.getLinearAlgebraMethod();
-  if (method == 0)
-  {
-    LOGINFO("Solve L.B = SigmaYX");
-    B = result_.getCholeskyFactor().solveLinearSystem(crossCovariance);
-  }
-  else
-  {
-    LOGINFO("Solve L.B = SigmaYX (h-mat version)");
-    B = result_.getHMatCholeskyFactor().solveLower(crossCovariance);
-  }
+  LOGINFO("Solve L.B = SigmaYX");
+  const Matrix B(solveTriangularSystem(crossCovariance));
   // Use of gram to compute B^{t} * B
   // Default gram computes B*B^t
   // With transpose argument=true, it performs B^t*B
@@ -247,18 +247,8 @@ CovarianceMatrix GaussianProcessConditionalCovariance::getConditionalCovariance(
   // 3) Compute r^t R^{-1} r'(x)
   // As we get the Cholesky factor L, we can solve triangular linear system
   // We define B = L^{-1} * r(x)
-  Matrix B;
-  const UnsignedInteger method = result_.getLinearAlgebraMethod();
-  if (method == 0)
-  {
-    LOGINFO("Solve L.B = SigmaYX");
-    B = result_.getCholeskyFactor().solveLinearSystem(crossCovariance);
-  }
-  else
-  {
-    LOGINFO("Solve L.B = SigmaYX (h-mat version)");
-    B = result_.getHMatCholeskyFactor().solveLower(crossCovariance);
-  }
+  LOGINFO("Solve L.B = SigmaYX");
+  const Matrix B(solveTriangularSystem(crossCovariance));
   // Use of gram to compute B^{t} * B
   // Default gram computes B*B^t
   // With transpose argument=true, it performs B^t*B
@@ -396,18 +386,8 @@ Sample GaussianProcessConditionalCovariance::getConditionalMarginalVariance(cons
     // 3) Compute r^t R^{-1} r'(x)
     // As we get the Cholesky factor L, we can solve triangular linear system
     // We define B = L^{-1} * r(x)
-    Matrix B;
-    const UnsignedInteger method = result_.getLinearAlgebraMethod();
-    if (method == 0)
-    {
-      LOGINFO("Solve L.B = SigmaYX");
-      B = result_.getCholeskyFactor().solveLinearSystem(crossCovariance);
-    }
-    else
-    {
-      LOGINFO("Solve L.B = SigmaYX (h-mat version)");
-      B = result_.getHMatCholeskyFactor().solveLower(crossCovariance);
-    }
+    LOGINFO("Solve L.B = SigmaYX");
+    const Matrix B(solveTriangularSystem(crossCovariance));
     // We compute diag(B^t B)
     // We can notice that it corresponds to the sum of elements
     // for each column

diff --git a/...ncertainty/Algorithm/MetaModel/Kriging/openturns/GaussianProcessConditionalCovariance.hxx b/...ncertainty/Algorithm/MetaModel/Kriging/openturns/GaussianProcessConditionalCovariance.hxx
@@ -107,6 +107,7 @@ public:
 protected:
 
   void computePhi();
+  Matrix solveTriangularSystem(const Matrix & rhs) const;
 
 private: