Eigen.h

Go to the documentation of this file.

//
//     Minotaur -- It's only 1/2 bull
//
//     (C)opyright 2008 - 2025 The Minotaur Team.
//
 
#ifndef MINOTAUREIGEN_H
#define MINOTAUREIGEN_H
 
#include "LinearFunction.h"
#include "QuadraticFunction.h"
 
namespace Minotaur {
 
  class QuadraticFunction;
  typedef QuadraticFunction *QuadraticFunctionPtr;
  typedef const QuadraticFunction *ConstQuadraticFunctionPtr;
  //class EigenVector;
  //typedef LinearFunction EigenVector;
  //typedef LinearFunctionPtr EigenVectorPtr;
  typedef std::pair<double, LinearFunctionPtr> EigenPair;
  typedef std::vector<EigenPair>::const_iterator EigenPairConstIterator;
 
  class Eigen;
  typedef Eigen *EigenPtr;
 
  // /**
  // In Minotaur, EigenVectors are used only in the context of quadratic
  // functions. Thus an Eigen Vector is just a linear function.
  // */
  // class EigenVector {
  //   public:
  //     // /**
  //     // Default constructor
  //     // */
  //     EigenVector();
 
  //     // /**
  //     // Destroy
  //     // */
  //     ~EigenVector() {}
 
  //     // /**
  //     // Display the nonzero coefficients of the eigen vector.
  //     // */
  //     void write(std::ostream &s) const;
 
  //   private:
  //     // /**
  //     // The vector
  //     // */
  //     LinearFunctionPtr lf_;
  // };
 
 
  // /**
  // The EigenCalculator class is used to calculate eigen vectors and values
  // for matrices and other objects that have an associated matrix (like a
  // quadratic).
  // */
  class EigenCalculator {
  public:
    EigenCalculator();
 
    ~EigenCalculator() {};
 
    EigenPtr findValues(ConstQuadraticFunctionPtr qf, bool scaletol = false);
 
    EigenPtr findValues(int n, double **H);
 
    EigenPtr findVectors(ConstQuadraticFunctionPtr qf);
 
    // /**
    // Let qf = x'Ax, lf = cx. First find eigen vectors of the hessian of
    // qf. Then, x'Ax = x'QRERQ'x, where Q is orthogonal (QQ' = I). R is a
    // diagonal matrix with i-th diagonal element being the square root of
    // i-th eigen value. E has entries 1,-1 along the diagonal. Let y =
    // RQ'x. Then c'x = c'QR^(-1)y = b'y.  b = R^(-1)Q'c. This function
    // calculates Q, RER, b.
    //
    // The "x" vectors in qf and lf are not the same. e.g.  x0*x0 + x1*x1 +
    // 2x0*x1 + 2x0 + x2 + x3.  here lf does not contain x1 and has extra
    // variables x2 and x3. The vector 'c' thus corresponds to only 2x0 and
    // we get: 1*(x0 + x1)^2 + 0*(x0 - x1)^2 + 1*(x0+x1) + 1*(x0-x1) + x2 +
    // x3.  p_terms will thus have (x0 + x1 + 0.5)^2, n_terms will have
    // terms that have negative eigen values.  lin_terms will have
    // x0-x1+x2+x3 and cb = -0.25.
    // */
    void getSumOfSquares(std::vector<LinearFunctionPtr> &p_terms,
                         std::vector<LinearFunctionPtr> &n_terms,
                         std::vector<double> &p_const,
                         std::vector<double> &n_const,
                         LinearFunctionPtr &lin_terms, double &c,
                         ConstQuadraticFunctionPtr qf,
                         ConstLinearFunctionPtr lf);
 
  private:
    ConstQuadraticFunctionPtr qf_;
 
    size_t n_;
 
    double *A_;
 
    int m_;
 
    double abstol_;
 
    char findVectors_;
    
    double *w_;
 
    // /**
    // the first M columns of z contain the orthonormal eigenvectors of the
    // matrix A corresponding to the selected eigenvalues, with the i-th
    // column of Z holding the eigenvector associated with W(i). If JOBZ =
    // 'N', then Z is not referenced.  Note: the user must ensure that at
    // least max(1,M) columns are supplied in the array Z; if RANGE = 'V',
    // the exact value of M is not known in advance and an upper bound must
    // be used.  Supplying N columns is always safe.
    // */
    double *z_;
 
    // /**
    // The i-th eigenvector is nonzero only in elements
    // ISUPPZ( 2*i-1 ) through ISUPPZ( 2*i ).
    // */
    int *isuppz_;
 
    std::map<ConstVariablePtr, UInt, CompareVariablePtr> indices_;
 
    std::vector<ConstVariablePtr> vars_;
 
    double getDotProduct_(ConstLinearFunctionPtr lf1,
                          ConstLinearFunctionPtr lf2);
 
    void calculate_();
 
    void fillA_();
 
    EigenPtr getEigen_(bool scaletol);
 
    LinearFunctionPtr getLinearFunction_(const int i);
  };
 
  class Eigen {
  public:
    Eigen();
 
    void add(double value, LinearFunctionPtr e_vector, double tol);
 
    UInt numNegative() const;
 
    UInt numZero() const;
 
    UInt numPositive() const;
 
    EigenPairConstIterator begin() const;
 
    EigenPairConstIterator end() const;
 
    void write(std::ostream &s) const;
 
  private:
    std::vector<EigenPair> evPairs_;
 
    UInt neg_;
 
    UInt zero_;
 
    UInt pos_;
  };
 
}  //namespace Minotaur
 
#endif  // MINOTAUREIGEN_H