Types.hxx

Go to the documentation of this file.

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *\
 * Copyright (c) 2025, Davide Stocco and Enrico Bertolazzi.                                      *
 *                                                                                               *
 * The Pipal project is distributed under the MIT License.                                       *
 *                                                                                               *
 * Davide Stocco                                                               Enrico Bertolazzi *
 * University of Trento                                                     University of Trento *
 * e-mail: davide.stocco@unitn.it                             e-mail: enrico.bertolazzi@unitn.it *
\* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 
#pragma once
 
#ifndef INCLUDE_PIPAL_DEFINES_HXX
#define INCLUDE_PIPAL_DEFINES_HXX
 
// Standard libraries
#include <string>
#include <iostream>
#include <iomanip>
#include <ostream>
 
// STL
#include <cmath>
#include <limits>
#include <algorithm>
#include <functional>
#include <type_traits>
#include <numeric>
#include <memory>
 
// Eigen library
#ifndef PIPAL_EIGEN_EXTERNAL
#include <Eigen/Dense>
#include <Eigen/Sparse>
#include <Eigen/SparseCholesky>
#endif
 
// Print Pipal errors
#ifndef PIPAL_ERROR
#define PIPAL_ERROR(MSG)                \
  {                                     \
    std::ostringstream os;              \
    os << MSG;                          \
    throw std::runtime_error(os.str()); \
  }
#endif
 
// Assert for Pipal
#ifndef PIPAL_ASSERT
#define PIPAL_ASSERT(COND, MSG) \
  if (!(COND))                  \
  {                             \
    PIPAL_ERROR(MSG);           \
  }
#endif
 
// Warning for Pipal
#ifndef PIPAL_WARNING
#define PIPAL_WARNING(MSG)         \
  {                                \
    std::cout << MSG << std::endl; \
  }
#endif
 
// Warning assert for Pipal
#ifndef PIPAL_ASSERT_WARNING
#define PIPAL_ASSERT_WARNING(COND, MSG) \
  if (!(COND))                          \
  {                                     \
    PIPAL_WARNING(MSG);                 \
  }
#endif
 
#ifndef PIPAL_DEFAULT_INTEGER_TYPE
#define PIPAL_DEFAULT_INTEGER_TYPE int
#endif

namespace Pipal
{
  using Integer = PIPAL_DEFAULT_INTEGER_TYPE;
 
  template<typename Real> using Vector       = Eigen::Vector<Real, Eigen::Dynamic>;
  template<typename Real> using Matrix       = Eigen::Matrix<Real, Eigen::Dynamic, Eigen::Dynamic>;
  template<typename Real> using SparseMatrix = Eigen::SparseMatrix<Real>;
  template<typename Real> using Array        = Eigen::Array<Real, Eigen::Dynamic, 1>;
 
  using Indices = Eigen::Array<Integer, Eigen::Dynamic, 1>;
  using Mask    = Eigen::Array<bool, Eigen::Dynamic, 1>;

  using Algorithm = enum class Algorithm : Integer {CONSERVATIVE = 0, ADAPTIVE = 1};

  static Indices find(Mask const & mask)
  {
    Indices out(mask.count());
    for (Integer i{0}, j{0}; i < mask.size(); ++i) {if (mask[i]) {out[j++] = i;}}
    return out;
  }

  template <typename Real>
  static void insert_block(SparseMatrix<Real> & mat, SparseMatrix<Real> const & blk,
    Integer const row_offset, Integer const col_offset)
  {
    #define CMD "Pipal::Solver::insert_block(...): "
 
    // Get sizes
    const Integer mat_cols{static_cast<Integer>(mat.cols())}, mat_rows{static_cast<Integer>(mat.rows())};
    const Integer blk_cols{static_cast<Integer>(blk.cols())}, blk_rows{static_cast<Integer>(blk.rows())};
 
    // Check dimensions
    PIPAL_ASSERT(row_offset >= 0,
      CMD "row offset must be non-negative.");
    PIPAL_ASSERT(col_offset >= 0,
      CMD "column offset must be non-negative.");
    PIPAL_ASSERT(row_offset + blk_rows <= mat_rows,
      CMD "inserting block exceeds matrix row dimensions.");
    PIPAL_ASSERT(col_offset + blk_cols <= mat_cols,
      CMD "inserting block exceeds matrix column dimensions.");
 
    // Insert block
    for (Integer r{0}; r < blk_rows; ++r) {
      for (Integer c{0}; c < blk_cols; ++c) {
        if (blk.coeff(r, c) != 0) {mat.coeffRef(r + row_offset, c + col_offset) = blk.coeff(r, c);}
      }
    }
 
    #undef CMD
  }

  template <typename Real>
  static void insert_block(SparseMatrix<Real> & mat, SparseMatrix<Real> const & blk,
    Indices const row_index, Indices const col_index, Integer const row_offset, Integer const col_offset)
  {
    #define CMD "Pipal::Solver::insert_block(...): "
 
    // Get sizes
    const Integer mat_cols{static_cast<Integer>(mat.cols())}, mat_rows{static_cast<Integer>(mat.rows())};
    const Integer idx_rows{static_cast<Integer>(row_index.size())};
    const Integer idx_cols{static_cast<Integer>(col_index.size())};
 
    // Check dimensions
    PIPAL_ASSERT(row_offset >= 0,
      CMD "row offset must be non-negative.");
    PIPAL_ASSERT(col_offset >= 0,
      CMD "column offset must be non-negative.");
    PIPAL_ASSERT(row_offset + idx_rows <= mat_rows,
      CMD "inserting block exceeds matrix row dimensions.");
    PIPAL_ASSERT(col_offset + idx_cols <= mat_cols,
      CMD "inserting block exceeds matrix column dimensions.");
 
    // Insert block
    for (Integer r{0}; r < idx_rows; ++r) {
      for (Integer c{0}; c < idx_cols; ++c) {
        if (blk.coeff(row_index(r), col_index(c)) != 0) {
          mat.coeffRef(r + row_offset, c + col_offset) = blk.coeff(row_index(r), col_index(c));
        }
      }
    }
 
    #undef CMD
  }

  template<typename Real>
  struct Parameter
  {
    static constexpr Real    rhs_bnd{1.0e+18};       
    static constexpr Real    grad_max{1.0e+02};      
    static constexpr Real    infeas_max{1.0e+02};    
    static constexpr Real    nnz_max{2.0e+04};       
    static constexpr Integer opt_err_mem{6};         
    static constexpr Real    ls_factor{5.0e-01};     
    static constexpr Real    ls_thresh{1.0e-08};     
    static constexpr Real    ls_frac{1.0e-02};       
    static constexpr Real    slack_min{1.0e-20};     
    static constexpr Real    shift_min{1.0e-12};     
    static constexpr Real    shift_factor1{5.0e-01}; 
    static constexpr Real    shift_factor2{6.0e-01}; 
    static constexpr Real    shift_max{1.0e+08};     
    static constexpr Real    rho_init{1.0e-01};      
    static constexpr Real    rho_min{1.0e-12};       
    static constexpr Real    rho_factor{5.0e-01};    
    static constexpr Integer rho_trials{8};          
    static constexpr Real    mu_init{1.0e-01};       
    static constexpr Real    mu_min{1.0e-12};        
    static constexpr Real    mu_factor{1.0e-01};     
    static constexpr Real    mu_factor_exp{1.5};     
    static constexpr Integer mu_trials{4};           
    static constexpr Real    mu_max{1.0e-01};        
    static constexpr Real    mu_max_exp0{0.0};       
    static constexpr Real    update_con_1{1.0e-02};  
    static constexpr Real    update_con_2{1.0e-02};  
    static constexpr Real    update_con_3{1.01};     
 
    Real      opt_err_tol{1.0e-6};            
    Integer   iter_max{1000};                 
    Algorithm algorithm{Algorithm::ADAPTIVE}; 
    Real      mu_max_exp{0.0};                
 
    Integer bfgs_update_freq{20}; 

    Parameter() = default;

    Parameter(Parameter const &) = delete;

    Parameter & operator=(Parameter const &) = delete;
 
  }; // struct Parameter

  using Counter = struct Counter
  {
    Integer f{0}; 
    Integer g{0}; 
    Integer H{0}; 
    Integer k{0}; 
    Integer M{0}; 

    Counter() = default;

    Counter(Counter const &) = delete;

    Counter & operator=(Counter const &) = delete;
 
  }; // struct Counter

  template<typename Real>
  struct Input
  {
    std::string name; 
    Indices       I1; 
    Indices       I2; 
    Indices       I3; 
    Indices       I4; 
    Indices       I5; 
    Indices       I6; 
    Indices       I7; 
    Indices       I8; 
    Indices       I9; 
    Vector<Real>  x0; 
    Vector<Real>  b2; 
    Vector<Real>  l3; 
    Vector<Real>  u4; 
    Vector<Real>  l5; 
    Vector<Real>  u5; 
    Vector<Real>  b6; 
    Vector<Real>  l7; 
    Vector<Real>  u8; 
    Vector<Real>  l9; 
    Vector<Real>  u9; 
    Integer       n0; 
    Integer       n1; 
    Integer       n2; 
    Integer       n3; 
    Integer       n4; 
    Integer       n5; 
    Integer       n6; 
    Integer       n7; 
    Integer       n8; 
    Integer       n9; 
    Integer       nV; 
    Integer       nI; 
    Integer       nE; 
    Integer       nA; 
    Integer       vi; 

    Input() = default;

    Input(Input const &) = delete;

    Input & operator=(Input const &) = delete;
 
  }; // struct Input

  template<typename Real>
  struct Iterate
  {
    using LDLT = Eigen::SimplicialLDLT<SparseMatrix<Real>, Eigen::Lower>;
 
    Vector<Real>       x;     
    Real               rho;   
    Real               rho_;  
    Real               mu;    
    Real               f;     
    Real               fu;    
    Vector<Real>       g;     
    Array<Real>        r1;    
    Array<Real>        r2;    
    Array<Real>        cE;    
    SparseMatrix<Real> JE;    
    Integer            JEnnz; 
    Array<Real>        lE;    
    Array<Real>        s1;    
    Array<Real>        s2;    
    Array<Real>        cI;    
    SparseMatrix<Real> JI;    
    Integer            JInnz; 
    Array<Real>        lI;    
    SparseMatrix<Real> H;     
    Integer            Hnnz;  
    Real               v;     
    Real               vu;    
    Real               v0;    
    Real               phi;   
    LDLT               ldlt;  
    Integer            Annz;  
    Real               shift; 
    Vector<Real>       b;     
    Vector<Real>       kkt;   
    Vector<Real>       kkt_;  
    Integer            err;   
 
    Real               fs;      
    Array<Real>        cEs;     
    Array<Real>        cEu;     
    Array<Real>        cIs;     
    Array<Real>        cIu;     
    SparseMatrix<Real> A;       
    Real               shift22; 
    Real               v_;      
    bool               cut_;    

    Iterate() = default;

    Iterate(Iterate const &) = delete;

    Iterate & operator=(Iterate const &) = delete;
 
  }; // struct Iterate

  template<typename Real>
  struct Direction
  {
    Vector<Real> x;       
    Real         x_norm;  
    Real         x_norm_; 
    Array<Real>  r1;      
    Array<Real>  r2;      
    Array<Real>  lE;      
    Array<Real>  s1;      
    Array<Real>  s2;      
    Array<Real>  lI;      
    Real         l_norm;  
    Real         lred0;   
    Real         ltred0;  
    Real         ltred;   
    Real         qtred;   
    Real         m;       

    Direction() = default;

    Direction(Direction const &) = delete;

    Direction & operator=(Direction const &) = delete;
 
  }; // struct Direction

  template<typename Real>
  struct Acceptance
  {
    Real p0{0.0};  
    Real p{0.0};   
    Real d{0.0};   
    bool s{false}; 

    Acceptance() = default;

    Acceptance(Acceptance const &) = delete;

    Acceptance & operator=(Acceptance const &) = delete;
 
  }; // struct Acceptance
 
} // namespace Pipal
 
#endif // INCLUDE_PIPAL_DEFINES_HXX