Attempt to include eigen-based LU decomposition.

src/core/CMakeLists.txt

@@ -114,6 +114,9 @@
     ${CMAKE_CURRENT_SOURCE_DIR}/matrices/codac2_Matrix.h
     ${CMAKE_CURRENT_SOURCE_DIR}/matrices/codac2_Row.h
     ${CMAKE_CURRENT_SOURCE_DIR}/matrices/codac2_Vector.h
+    ${CMAKE_CURRENT_SOURCE_DIR}/matrices/codac2_IntFullPivLU.h
+    ${CMAKE_CURRENT_SOURCE_DIR}/matrices/codac2_IntFullPivLU.cpp
+
 
     ${CMAKE_CURRENT_SOURCE_DIR}/paver/codac2_pave.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/paver/codac2_pave.h

src/core/matrices/codac2_IntFullPivLU.cpp

@@ -0,0 +1,380 @@
+/** 
+ *  \file codac2_IntFullPivLU.cpp
+ * ----------------------------------------------------------------------------
+ *  \date       2024
+ *  \author     Damien Massé
+ *  \copyright  Copyright 2024 Codac Team
+ *  \license    GNU Lesser General Public License (LGPL)
+ */
+
+#include <ostream>
+#include <cstdlib>
+#include "codac2_Matrix.h"
+#include "codac2_Interval.h"
+#include "codac2_IntervalMatrix.h"
+#include "codac2_IntervalVector.h"
+#include "codac2_BoolInterval.h"
+#include "codac2_IntFullPivLU.h"
+
+#include "gaol/gaol_fpu.h" /* for rounding settings */
+
+namespace codac2 {
+
+/* product of double values, keeping 0 */
+inline double prodkeepzero(double a, double b) {
+    if (a==0 || b==0) return 0;
+    else if (a==oo || b==oo) return oo;
+    else return a*b;
+}
+
+/* utility function : compute the bound using Gauss-Jordan successive
+   algorithm */
+IntervalMatrix IntFullPivLU::buildLUbounds(const IntervalMatrix &E) {
+ /* we need rounding up for floating-point operations */
+    GAOL_RND_PRESERVE(); /* save the rounding state */
+    round_upward(); /* round up */
+
+    const Index nCols = E.cols();
+    const Index nRows = E.rows();
+    const Index nC = std::min(nCols,nRows);
+
+    Matrix M = E.mag();
+    IntervalMatrix res = E;
+    for (Index k=0;k<nC;k++) {
+
+       double pivot = (M(k,k)>=1.0 ? oo : 1.0/(-(M(k,k)-1.0)));
+       Matrix mP = M.block(0,k,k,1) * M.block(k,0,1,k);
+       for (auto coef : mP.reshaped()) {
+           coef=prodkeepzero(coef,pivot);
+       }
+       M.block(0,0,k,k) += mP;
+       for (auto coef : M.block(0,k,k,1).reshaped()) {
+           coef=prodkeepzero(coef,pivot);
+       }
+       for (auto coef : M.block(k,0,1,k).reshaped()) {
+           coef=prodkeepzero(coef,pivot);
+       }
+       M(k,k) = prodkeepzero(M(k,k),pivot);
+
+       IntervalMatrix U = IntervalMatrix::Identity(k+1,k+1);
+       for (Index col=0;col<=k;col++) {
+         for (Index row=0;row<=k;row++) 
+              U(row,col)+=(M(row,col)!=oo ?
+			Interval(-1,1)*M(row,col) :
+                        Interval());
+       }
+       res.block(k+1,k,nRows-k-1,1) +=
+	  E.bottomLeftCorner(nRows-k-1,k+1)*U*E.block(0,k,k+1,1);
+       res.block(k+1,k+1,1,nCols-k-1) +=
+          E.block(k+1,0,1,k+1)*U*E.topRightCorner(k+1,nCols-k-1);
+    }
+
+ /* restore rounding */
+    GAOL_RND_RESTORE();
+    return res;
+}
+
+
+/************************************************************/
+
+
+/** constructor from Matrix of double 
+ */
+IntFullPivLU::IntFullPivLU(const Matrix &M) :
+   _LU(M), 
+   matrixLU_(M.rows(), M.cols())
+{
+   this->computeMatrixLU(IntervalMatrix(M), _LU.maxPivot()*_LU.threshold());
+}
+
+/** constructor from Matrix of Intervals
+ */
+IntFullPivLU::IntFullPivLU(const IntervalMatrix &M) :
+   _LU(M.mid()), 
+   matrixLU_(M.rows(), M.cols())
+{
+   this->computeMatrixLU(M,_LU.maxPivot()*_LU.threshold());
+}
+
+void IntFullPivLU::computeMatrixLU(const IntervalMatrix &M, double nonzero) {
+   /* specific case if _LU.maxPivot() is 0 (the matrix itself is 0) */
+    if (_LU.maxPivot()==0.0) {
+       nonzero=1.0; 
+	/* strong threshold. The most probable result is oo anyway */
+    }
+   /* 1) calculer la matrice d'erreurs
+          -> déterminer L'^{-1}, U'^{-1}
+      2) calcul les exposants de cette matrice
+      3) reconstuire L et U
+    */
+    int nRows=M.rows();
+    int nCols=M.cols();
+    int dim = std::min(nRows,nCols);
+    Matrix mLU = _LU.matrixLU();
+    /* modification of mLU if not full rank (check this) */
+    for (int i=0;i<dim;i++) {
+        if (std::fabs(mLU(i,i))<nonzero) {
+           mLU(i,i)=(mLU(i,i)<0.0 ? -nonzero : nonzero); 
+        }
+    }
+    /* "Inversion" of mLU (i.e. (pseudo)inversion of L and U) */
+    IntervalMatrix ImLU 
+		= IntervalMatrix::Zero(nRows,nCols);
+    /* first L */
+    for (int c=0;c<nCols;c++) {
+      if (c>=nRows-1) break;
+      for (int r=c+1;r<nRows;r++) {
+         Interval s(mLU(r,c));
+         for (int k=c+1;k<r;k++) {
+            if (k>=nCols) break;
+            s += mLU(r,k)*ImLU(k,c);
+         }
+         ImLU(r,c)=-s;
+      }
+    }
+    /* then U */
+    for (int c=0;c<nCols;c++) {
+      int r;
+      if (c<nRows) {
+         ImLU(c,c)=Interval(1.0)/mLU(c,c);
+         r=c-1;
+      } else {
+         r=nRows-1;
+      }
+      for (;r>=0;r--) {
+         Interval s(mLU(r,c));
+         if (c<nRows) s/=mLU(c,c);
+         for (int k=r+1;k<c;k++) {
+             if (k>=nRows) break;
+             s += mLU(r,k)*ImLU(k,c);
+         }
+         ImLU(r,c)=-s/mLU(r,r);
+      }
+    }
+
+    /* compute the error matrix */
+    IntervalMatrix InvL = 
+		IntervalMatrix::Identity(nRows,nRows);
+    if (nRows>nCols) 
+      InvL.leftCols(nCols).triangularView<Eigen::StrictlyLower>() = ImLU;
+    else 
+      InvL.triangularView<Eigen::StrictlyLower>() = ImLU.leftCols(nRows);
+    IntervalMatrix InvU = 
+		IntervalMatrix::Zero(nCols,nCols);
+    if (nCols>nRows) 
+      InvU.topRows(nRows).triangularView<Eigen::Upper>() = ImLU;
+    else 
+      InvU.triangularView<Eigen::Upper>() = ImLU.topRows(nCols);
+    IntervalMatrix error = 
+       InvL * (_LU.permutationP() * M * 
+       _LU.permutationQ()) * InvU - 
+	  IntervalMatrix::Identity(nRows,nCols);
+
+    IntervalMatrix eps = IntFullPivLU::buildLUbounds(error);
+
+    /* product */
+    for (Index c=0;c<nCols;c++) 
+    for (Index r=0;r<nRows;r++) {
+       matrixLU_(r,c)=mLU(r,c); /* for Id */
+       if (r>c) { /* low part */
+          matrixLU_(r,c)+=eps(r,c); /* Id for L part */
+          for (Index k=c+1;k<r;k++) {
+             if (k>=nCols) continue;
+             matrixLU_(r,c)+=mLU(r,k)*eps(k,c);
+          }
+       } else { /* high part */
+          for (Index k=r;k<=c;k++) {
+             if (k>=nRows) continue;
+             matrixLU_(r,c)+=eps(r,k)*mLU(k,c);
+          }
+       }
+    }
+}
+
+/* computing the possible rank for a square matrix is easy as
+   we have only to look at the diagonal of U
+   for non-square matrix, the remaining lines/columns _may_ have
+   the rank */
+Interval IntFullPivLU::rank() const {
+    Interval res(0.0);
+    /* first we compute the rank generated by the diagonal of U */
+    Index dim = std::min(matrixLU_.rows(),matrixLU_.cols());
+    for (Index i=0;i<dim;i++) {
+       if (matrixLU_(i,i).mig()>0.0) res+=1.0;
+       else 
+       if (matrixLU_(i,i).mag()>0.0) res+=Interval(0.0,1.0);
+    }
+    if (res.lb()==dim) return res; /* matrix if full rank */
+    /* complement with potential other lines / rows */
+    Index remain = std::max(matrixLU_.rows(),matrixLU_.cols())-dim;
+    if (remain==0.0) return res; /* square matrix */
+    res = min(dim,res+Interval(0.0,remain));
+    return res;
+}
+
+
+double IntFullPivLU::maxPivot() const {
+    double res=0.0;
+    /* first we compute the rank generated by the diagonal of U */
+    Index dim = std::min(matrixLU_.rows(),matrixLU_.cols());
+    for (Index i=0;i<dim;i++) {
+       res = std::max(matrixLU_(i,i).ub(),res);
+    }
+    return res;
+}
+
+/* possible dimension of kernel = number of columns - rank of the matrix */
+Interval IntFullPivLU::dimensionOfKernel() const {
+    Interval rk = this->rank();
+    return matrixLU_.cols()-rk;
+}
+
+BoolInterval IntFullPivLU::isSurjective() const {
+    if (matrixLU_.rows()>matrixLU_.cols()) return BoolInterval::FALSE;
+    Interval r=this->rank();
+    if (r.lb()==matrixLU_.rows()) return BoolInterval::TRUE;
+    if (r.ub()==matrixLU_.rows()) return BoolInterval::UNKNOWN;
+    return BoolInterval::FALSE;
+}
+
+BoolInterval IntFullPivLU::isInjective() const {
+    if (matrixLU_.rows()<matrixLU_.cols()) return BoolInterval::FALSE;
+    Interval r=this->rank();
+    if (r.lb()==matrixLU_.cols()) return BoolInterval::TRUE;
+    if (r.ub()==matrixLU_.cols()) return BoolInterval::UNKNOWN;
+    return BoolInterval::FALSE;
+}
+
+BoolInterval IntFullPivLU::isInvertible() const {
+    if (!matrixLU_.is_squared()) return BoolInterval::FALSE;
+    Interval r=this->rank();
+    if (r.lb()==matrixLU_.cols()) return BoolInterval::TRUE;
+    if (r.ub()==matrixLU_.cols()) return BoolInterval::UNKNOWN;
+    return BoolInterval::FALSE;
+}
+
+Interval IntFullPivLU::determinant() const {
+    assert_release(matrixLU_.is_squared());
+    Interval product(1.0);
+    for (Index i=0;i<matrixLU_.cols();i++) {
+       product *= matrixLU_(i,i);
+          /* if matrixLU_(i,i)=0.0, the determinant is 0 even if
+             a previous diagonal element was (-oo,oo) */ 
+    }
+    return product;
+}
+
+IntervalMatrix IntFullPivLU::kernel() const {
+    /* like the rank, the computation is complex , and 
+       may give strange results in some bases. We'll build the vectors
+       one by one and make the matrix afterwards */
+    std::vector<IntervalVector> kernel;
+    std::vector<bool> generating(matrixLU_.cols(),true);
+    for (Index c=0;c<matrixLU_.cols();c++) {
+      /* a potential kernel vector is built as follows :
+         we consider each column c and states that the vector must be 0
+         for all column > c , and 1 for c.
+         if we manage to build such a vector by filling dependant indices
+         < c, then this index is stated as "generating" and will be 0 for
+         the next vectors. Otherwise, it is states as "not generating" and 
+         and used in next vectors */
+      /* before anything, if c<matrixLU_.rows() and 
+         matrixLU_(c,c) does not contains 0, just states that 
+         c is not generating and go on */
+      if (c<matrixLU_.rows() && !matrixLU_(c,c).contains(0.0)) {
+          generating[c]=false;
+          continue;
+      }
+      IntervalVector vect = IntervalVector::Zero(matrixLU_.cols());
+      vect[c]=1.0;
+      for (Index c1=c-1;c1>=0;c1--) {
+         if (c1>=matrixLU_.rows())  {
+             /* we will consider two cases :
+                1) generating[c1] is true. then 
+                       we consider vect[c1]=0.0
+                2) generating[c1] is false. we did not manage to
+		   build a vector using an element outside the square
+		   part of U, which should almost never happens,
+		   but means that the decomposition badly failed.
+                   just consider vect[c1]=top in this case */
+             if (!generating[c1]) vect[c1]=Interval();
+             continue;
+         }
+         /* in the following, either 0.0 \notin matrixLU_(c1,c1)
+            and generating(c1) is false, which is fine, 
+            or 0.0 in matrixLU_(c1,c1), which is less fine, but we
+            go on nevertheless by checking the row */
+         Interval z(0.0);
+         for (Index c2=c1+1;c2<=c;c2++) {
+           z -= matrixLU_(c1,c2)*vect[c2];
+         }
+         /* vect[c1] * matrixLU_(c1,c1) = z */
+         vect[c1]=Interval();
+         Interval tmp = matrixLU_(c1,c1);
+         bwd_mul(z,vect[c1],tmp);
+         if (vect[c1].is_empty()) { 
+             /* means that matrixLU_(c1,c1)=0.0 
+		and z!=0.0... that's bad news */
+             generating[c]=false; break;
+         }
+         /* note: here, in theory if generating[c1] is 0 we should
+            be able to put vect[c1]=0.0 even if 0.0 is not in 
+            vect[c1], because the actual kernel is a combination
+            of found vectors. _however_, it means that the vector
+            we are constructing is _not_ in the kernel, just that it
+            can be used to build the kernel... but that may be the only
+            acceptable approach */
+         if (generating[c1]) vect[c1]=0.0;
+      }
+      if (!generating[c]) continue;
+      kernel.push_back(_LU.permutationQ()*vect);
+    }
+    /* build the matrix */
+    IntervalMatrix res(matrixLU_.cols(),kernel.size());
+    for (unsigned int c=0;c<kernel.size();c++) {
+       res.col(c)=kernel[c];
+    }
+    return res;
+}
+
+IntervalMatrix IntFullPivLU::solve(const IntervalMatrix &rhs) const {
+    IntervalMatrix prhs = _LU.permutationP()*rhs;
+    /* inverse L */
+    for (Index r = 0; r<matrixLU_.rows()-1; r++) {
+       for (Index r1=r+1;r1<matrixLU_.rows();r1++) {
+          prhs.row(r1) -= matrixLU_(r1,r)*prhs.row(r);
+       }
+    }
+    /* first check, if rows>cols, the 0 lines */
+    for (Index r=matrixLU_.cols(); r<matrixLU_.rows();r++) {
+        for (Index a=0;a<prhs.cols();a++) {
+            if (!prhs(r,a).contains(0.0)) {
+               Interval emp; emp.set_empty();
+               return 
+		IntervalMatrix::Constant(matrixLU_.cols(),rhs.cols(),emp);
+            }
+        }
+    }
+    /* then use the diagonal elements */
+    IntervalMatrix res = IntervalMatrix::Zero(matrixLU_.cols(),rhs.cols());
+    Index dim = std::min(matrixLU_.cols(),matrixLU_.rows());
+    for (Index r = dim-1;r>=0;r--) {
+        for (Index r1=r+1;r1<dim;r1++) {
+           prhs.row(r) -= matrixLU_(r,r1)*res.row(r1);
+        }
+        if (matrixLU_(r,r).contains(0.0)) {
+            for (Index a=0;a<prhs.cols();a++) {
+              if (!prhs(r,a).contains(0.0)) {
+                 for (auto coef : res.reshaped()) coef.set_empty();
+                 return res;
+              }
+            } 
+        } else {
+            res.row(r) = (1.0/matrixLU_(r,r))*prhs.row(r);
+        }
+    }
+    return _LU.permutationQ()*res;
+}
+
+
+}