LeastSquares.h

#ifndef LeastSquares_h
#define LeastSquares_h
/* SPDX-License-Identifier: CC0-1.0 */
#include <vector>
 
#include "tnt/tnt_array2d.h"
 
#include <armadillo>
 
#include "BasisFunction.h"
 
namespace Isis {
  class LeastSquares {
    public:
 
      LeastSquares(Isis::BasisFunction &basis, bool sparse = false,
                   int sparseRows = 0, int sparseCols = 0, bool jigsaw = false);
      ~LeastSquares();
      void AddKnown(const std::vector<double> &input, double expected,
                    double weight = 1.0);
 
      std::vector<double> GetInput(int row) const;
      double GetExpected(int row) const;
      int Rows() const;
 
      enum SolveMethod { SVD,     
                         QRD,     
                         SPARSE   
                       };
 
      int Solve(Isis::LeastSquares::SolveMethod method = SVD);
      double Evaluate(const std::vector<double> &input);
      std::vector<double> Residuals() const;
      double Residual(int i) const;
      void Weight(int index, double weight);
 
      int Knowns() const {
        return p_expected.size();
      };
 
      double GetSigma0() { return p_sigma0; }
      int GetDegreesOfFreedom() { return p_degreesOfFreedom; }
      void Reset ();
 
      void ResetSparse() { Reset(); }
      std::vector<double> GetEpsilons () const { return p_epsilonsSparse; }
      void SetParameterWeights(const std::vector<double> weights) { p_parameterWeights = weights; }
      void SetNumberOfConstrainedParameters(int n) { p_constrainedParameters = n; }
 
    private:
      void SolveSVD();
      void SolveQRD();
      void SolveCholesky () {}
 
      int SolveSparse();
      void FillSparseA(const std::vector<double> &data);
      bool ApplyParameterWeights();
 
      arma::mat p_xSparse;          
      std::vector<double> p_epsilonsSparse;   
      std::vector<double> p_parameterWeights; 
      arma::SpMat<double> p_sparseA; 
      arma::SpMat<double> p_normals; 
      arma::mat p_ATb;                   
      arma::mat p_SLU_Factor;          
      bool p_jigsaw;
      bool p_sparse;
      bool p_solved;  
      int p_currentFillRow;
      int p_sparseRows;
      int p_sparseCols;
      int p_constrainedParameters; 
      int p_degreesOfFreedom;      
      double p_sigma0;               
      std::vector<std::vector<double> > p_input; 
      std::vector<double> p_expected;            
      std::vector<double> p_sqrtWeight;          
      std::vector<double> p_residuals;           
      Isis::BasisFunction *p_basis; 
 
 
 
  };
};
 
 
#endif