From a435caeecb1f8dc0dd239bc6917e7052e75f7a70 Mon Sep 17 00:00:00 2001
From: Antonin Portelli <antonin.portelli@me.com>
Date: Mon, 29 Feb 2016 19:50:28 +0000
Subject: [PATCH] fits: covariance matrix is inverted using pseudoinverse, the
 interface allows to se the SVD tolerance

---
 lib/Chi2Function.cpp |  2 +-
 lib/FitInterface.cpp | 23 ++++++++++++++-----
 lib/FitInterface.hpp | 53 +++++++++++++++++++++++---------------------
 3 files changed, 46 insertions(+), 32 deletions(-)

diff --git a/lib/Chi2Function.cpp b/lib/Chi2Function.cpp
index fd7f51f..d189663 100644
--- a/lib/Chi2Function.cpp
+++ b/lib/Chi2Function.cpp
@@ -233,7 +233,7 @@ void Chi2Function::initBuffer(void) const
     lowerYX = upperYX.transpose();
 
     // inversion
-    buffer_->invVar = buffer_->invVar.inverse().eval();
+    buffer_->invVar = buffer_->invVar.pInverse(data_->getSvdTolerance());
     buffer_->isInit = true;
 }
 
diff --git a/lib/FitInterface.cpp b/lib/FitInterface.cpp
index 43de511..f99e1d1 100644
--- a/lib/FitInterface.cpp
+++ b/lib/FitInterface.cpp
@@ -128,16 +128,22 @@ Index FitInterface::getStatXDim(void) const
     return isXExact_.size() - isXExact_.sum();
 }
 
+double FitInterface::getSvdTolerance(void) const
+{
+    return svdTolerance_;
+}
+
 void FitInterface::setFitInterface(const FitInterface &fitInterface)
 {
     if (&fitInterface != this)
     {
-        isXExact_   = fitInterface.isXExact_;
-        isFitPoint_ = fitInterface.isFitPoint_;
-        isXXCorr_   = fitInterface.isXXCorr_;
-        isYYCorr_   = fitInterface.isYYCorr_;
-        isYXCorr_   = fitInterface.isYXCorr_;
-        isDataCorr_ = fitInterface.isDataCorr_;
+        isXExact_     = fitInterface.isXExact_;
+        isFitPoint_   = fitInterface.isFitPoint_;
+        isXXCorr_     = fitInterface.isXXCorr_;
+        isYYCorr_     = fitInterface.isYYCorr_;
+        isYXCorr_     = fitInterface.isYXCorr_;
+        isDataCorr_   = fitInterface.isDataCorr_;
+        svdTolerance_ = fitInterface.svdTolerance_;
     }
 }
 
@@ -156,6 +162,11 @@ void FitInterface::setYDim(const Index yDim)
     resize(getNData(), getXDim(), yDim);
 }
 
+void FitInterface::setSvdTolerance(const double tolerance)
+{
+    svdTolerance_ = tolerance;
+}
+
 void FitInterface::resize(const Index nData, const Index xDim, const Index yDim)
 {
     isXExact_.setConstant(xDim, 0);
diff --git a/lib/FitInterface.hpp b/lib/FitInterface.hpp
index 3860992..5605173 100644
--- a/lib/FitInterface.hpp
+++ b/lib/FitInterface.hpp
@@ -39,30 +39,32 @@ public:
     // destructor
     virtual ~FitInterface(void) = default;
     // access
-    void  assumeXExact(const Index i, const bool isExact = true);
-    void  assumeXXCorrelated(const Index i1, const Index i2,
-                             const bool isCorrelated = true);
-    void  assumeYYCorrelated(const Index j1, const Index j2,
-                             const bool isCorrelated = true);
-    void  assumeYXCorrelated(const Index j, const Index i,
-                             const bool isCorrelated = true);
-    void  assumeDataCorrelated(const Index k1, const Index k2,
-                               const bool isCorrelated = true);
-    void  assumeDataCorrelated(const bool isCorrelated = true);
-    void  fitPoint(const Index k, const bool isFitPoint = true);
-    void  fitPointRange(const Index k1, const Index k2,
+    void   assumeXExact(const Index i, const bool isExact = true);
+    void   assumeXXCorrelated(const Index i1, const Index i2,
+                              const bool isCorrelated = true);
+    void   assumeYYCorrelated(const Index j1, const Index j2,
+                              const bool isCorrelated = true);
+    void   assumeYXCorrelated(const Index j, const Index i,
+                              const bool isCorrelated = true);
+    void   assumeDataCorrelated(const Index k1, const Index k2,
+                                const bool isCorrelated = true);
+    void   assumeDataCorrelated(const bool isCorrelated = true);
+    void   fitPoint(const Index k, const bool isFitPoint = true);
+    void   fitPointRange(const Index k1, const Index k2,
                         const bool isFitPoint = true);
-    void  fitAllPoints(const bool isFitPoint = true);
-    Index getNData(void)     const;
-    Index getNFitPoint(void) const;
-    Index getXDim(void)      const;
-    Index getYDim(void)      const;
-    Index getStatXDim(void)  const;
-    void  setFitInterface(const FitInterface &fitInterface);
-    void  setNData(const Index nData);
-    void  setXDim(const Index xDim);
-    void  setYDim(const Index yDim);
-    void  resize(const Index nData, const Index xDim, const Index yDim);
+    void   fitAllPoints(const bool isFitPoint = true);
+    Index  getNData(void)        const;
+    Index  getNFitPoint(void)    const;
+    Index  getXDim(void)         const;
+    Index  getYDim(void)         const;
+    Index  getStatXDim(void)     const;
+    double getSvdTolerance(void) const;
+    void   setFitInterface(const FitInterface &fitInterface);
+    void   setNData(const Index nData);
+    void   setXDim(const Index xDim);
+    void   setYDim(const Index yDim);
+    void   setSvdTolerance(const double tolerance);
+    void   resize(const Index nData, const Index xDim, const Index yDim);
     // test
     bool isFitPoint(const Index k) const;
     bool isXExact(const Index i) const;
@@ -71,8 +73,9 @@ public:
     bool isYXCorrelated(const Index j, const Index i) const;
     bool isDataCorrelated(const Index k1, const Index k2) const;
 private:
-    IVec isXExact_, isFitPoint_;
-    IMat isXXCorr_, isYYCorr_, isYXCorr_, isDataCorr_;
+    IVec   isXExact_, isFitPoint_;
+    IMat   isXXCorr_, isYYCorr_, isYXCorr_, isDataCorr_;
+    double svdTolerance_{1.0e-10};
 };
 
 END_LATAN_NAMESPACE