.. _program_listing_file_src_solvers_PathSolver.cpp:

Program Listing for File PathSolver.cpp
=======================================

|exhale_lsh| :ref:`Return to documentation for file <file_src_solvers_PathSolver.cpp>` (``src/solvers/PathSolver.cpp``)

.. |exhale_lsh| unicode:: U+021B0 .. UPWARDS ARROW WITH TIP LEFTWARDS

.. code-block:: cpp

   extern "C" {
   #include "License.h"
   #include "MCP_Interface.h"
   #include "Macros.h"
   #include "Options.h"
   #include "Output.h"
   #include "Output_Interface.h"
   #include "Path.h"
   #include "PathOptions.h"
   }
   #include "solvers/PathSolver.h"
   
   void Solvers::PATH::sort(int rows, int cols, int elements, int *row, int *col, double *data) {
     double *m_data;
     int    *m_start;
     int    *m_len;
     int    *m_row;
   
     int i = 0, cs = 0, ce = 0;
   
     m_start = new int[cols + 1];
     m_len   = new int[cols + 1];
     m_row   = new int[elements + 1];
     m_data  = new double[elements + 1];
   
     for (i = 0; i < cols; i++) {
        m_len[i] = 0;
     }
   
     for (i = 0; i < elements; i++) {
        if ((col[i] < 1) || (col[i] > cols)) {
           throw("column incorrect.\n");
        }
   
        if ((row[i] < 1) || (row[i] > rows)) {
           throw("column incorrect.\n");
        }
   
        m_len[col[i] - 1]++;
     }
   
     m_start[0] = 0;
     for (i = 1; i < cols; i++) {
        m_start[i]   = m_start[i - 1] + m_len[i - 1];
        m_len[i - 1] = 0;
     }
     m_len[i - 1] = 0;
   
     for (i = 0; i < elements; i++) {
        cs         = col[i] - 1;
        ce         = m_start[cs] + m_len[cs];
        m_row[ce]  = row[i];
        m_data[ce] = data[i];
        m_len[cs]++;
     }
   
     elements = 0;
     for (i = 0; i < cols; i++) {
        cs = m_start[i];
        ce = cs + m_len[i];
   
        while (cs < ce) {
           row[elements]  = m_row[cs];
           col[elements]  = i + 1;
           data[elements] = m_data[cs];
           elements++;
           cs++;
        }
     }
   }
   
   int Solvers::PATH::CreateLMCP(int    n,
                                           int    m_nnz,
                                           int    m_i[],
                                           int    m_j[],
                                           double m_ij[],
                                           double q[],
                                           double lb[],
                                           double ub[],
                                           double x[],
                                           double z[],
                                           int    verbose,
                                           double timeLimit) {
     MCP_Termination termination;
     Information     information;
   
     Output_SetLog(stdout);
   
   
   
     auto o = Options_Create();
     Path_AddOptions(o);
     Options_Default(o);
   
   
   
     double *xSol, *zSol;
     double  dnnz;
     int     i;
   
   
     if (verbose == 0) {
        Options_SetBoolean(o, "output", False);
        Options_SetBoolean(o, "output_errors", False);
        Options_SetBoolean(o, "output_crash_iterations", False);
        Options_SetBoolean(o, "output_major_iterations", False);
        Options_SetBoolean(o, "output_minor_iterations", False);
        Options_SetBoolean(o, "output_options", False);
        Options_SetBoolean(o, "output_warnings", False);
     }
     if (timeLimit > 0)
        Options_SetDouble(o, "time_limit", timeLimit);
   
   
     Options_SetInt(o, "cumulative_iteration_limit", 4e6);
     Options_SetInt(o, "major_iteration_limit", 1e6);
     Options_SetInt(o, "minor_iteration_limit", 1e5);
     Options_SetInt(o, "crash_iteration_limit", 5e2);
     Options_SetInt(o, "nms_memory_size", 1e2);
     Options_SetInt(o, "lemke_rank_deficiency_iterations", 1e2);
     // Very important!
     Options_Set(o, "lemke_start_type advanced");
     Options_SetBoolean(o, "crash_perturb", True);
     Options_SetDouble(o, "proximal_perturbation", 1);
     Options_SetInt(o, "crash_nbchange_limit", 50);
     Options_SetDouble(o, "convergence_tolerance", 1e-6);
   
   
   
     if (n == 0) {
        Options_Destroy(o);
        this->status = ZEROStatus::Numerical;
        return MCP_Solved;
     }
   
     dnnz = MIN(1.0 * m_nnz, 1.0 * n * n);
   
     if (verbose != 0) {
        fprintf(stdout,
                   "%d row/cols, %f non-zeros, %3.2f%% dense.\n\n",
                   n,
                   dnnz,
                   100.0 * dnnz / (1.0 * n * n));
     }
   
     this->Problem.n       = n;
     this->Problem.nnz     = m_nnz;
     this->Problem.x       = x;
     this->Problem.q       = q;
     this->Problem.lb      = lb;
     this->Problem.ub      = ub;
     this->Problem.m_start = new int[n];
     this->Problem.m_len   = new int[n];
     this->Problem.m_row   = new int[m_nnz + 1];
     this->Problem.m_data  = new double[m_nnz + 1];
   
     sort(n, n, m_nnz, m_i, m_j, m_ij);
   
     int m_index = 0, m_count = 0;
     for (i = 0; i < n; i++) {
        this->Problem.m_start[i] = m_count + 1;
        this->Problem.m_len[i]   = 0;
   
        while ((m_index < m_nnz) && (m_j[m_index] <= i + 1)) {
           if (m_ij[m_index] != 0) {
             this->Problem.m_len[i]++;
             this->Problem.m_row[m_count]  = m_i[m_index];
             this->Problem.m_data[m_count] = m_ij[m_index];
             m_count++;
           }
           m_index++;
        }
     }
     this->Problem.nnz = m_count;
   
   
   
     License_SetString(
           "2830898829&Courtesy&&&USR&45321&5_1_2021&1000&PATH&GEN&31_12_2025&0_0_0&6000&0_0");
     auto m = MCP_Create(this->Problem.n, this->Problem.nnz + 1);
     MCP_SetInterface(m, &this->PATH_Interface);
     MCP_SetPresolveInterface(m, &this->PATH_Presolve);
     if (verbose != 0)
        Output_SetInterface(&this->PATH_Output);
     MCP_Jacobian_Structure_Constant(m, static_cast<Boolean>(true));
   
   
   
     Information info;
     if (verbose != 0)
        info.generate_output = Output_Log | Output_Status | Output_Listing;
     info.use_start  = False;
     info.use_basics = False;
   
     try {
        termination = Path_Solve(m, &info);
     } catch (...) {
        throw ZEROException(ZEROErrorCode::SolverError, "PATHWrapper threw an exception");
     }
   
   
     switch (termination) {
     case MCP_Solved: {
        this->status = ZEROStatus::Solved;
        xSol         = MCP_GetX(m);
        zSol         = MCP_GetF(m);
        for (i = 0; i < n; i++) {
           x[i] = xSol[i];
   
           z[i] = zSol[i];
        }
     } break;
     case MCP_TimeLimit:
        this->status = ZEROStatus::TimeLimit;
        break;
     case MCP_Infeasible:
        this->status = ZEROStatus::NotSolved;
     default:
        this->status = ZEROStatus::Numerical;
     }
   
   
   
     MCP_Destroy(m);
     Options_Destroy(o);
     Path_Destroy();
   
     return termination;
   }
   
   Solvers::PATH::PATH(const arma::sp_mat   &M,
                             const arma::vec      &q,
                             const perps          &Compl,
                             const VariableBounds &Bounds,
                             arma::vec            &z,
                             arma::vec            &x,
                             double                timeLimit,
                             bool                  verbose) {
   
     // Note that for path z and x are inverted (z are variables, x equations definitions).
     int n   = 0;
     int nnz = 0;
   
     // M.save("LCP_M.dat", arma::csv_ascii);
     // q.save("LCP_q.dat",arma::csv_ascii);
   
     this->Filled = false;
     std::vector<double> _q, _Mij, _lb, _ub, _xsol, _zsol;
     std::vector<int>    _Mi, _Mj, _xmap, _zmap;
     unsigned int        row = 1; // Fortran style, we start from 1
   
   
   
     for (const auto p : Compl) {
        // For each complementarity
        // z[p.first] \perp x[p.second]
   
   
        {
           /*if (M.row(p.first).n_nonzero >0) {
             LOG_S(WARNING) << "Solvers::PathLCP: Empty row associated with z_"
                                                 << std::to_string(p.first);
           }
   
           else {
               */
           int lb = Bounds.at(p.second).first;
           int ub = Bounds.at(p.second).second;
   
           // if (lb != ub)
           {
             ++n;
             // Avoid inserting fixed variables
   
             _lb.push_back(lb > 0 ? lb : 0);
             _ub.push_back(ub >= 0 ? ub : 1e20); // PATH will treat 1e20 as infinite. Note that negative UBs means infinity
   
             _q.push_back(q.at(p.first));
             _xmap.push_back(p.second);
             _zmap.push_back(p.first);
   
             _xsol.push_back(0);
             _zsol.push_back(0);
   
             ++row;
           } // end bounds
        }   // end empty row
     }     // end while
   
   
     // Sparse iterator for M
     for (arma::sp_mat::const_iterator it = M.begin(); it != M.end(); ++it) {
        _Mi.push_back(it.row() + 1);
        _Mj.push_back(it.col() + 1);
        _Mij.push_back(*it);
        ++nnz;
     }
   
     int stat = 0;
     /*
     _Mi.push_back(0);
     _Mj.push_back(0);
     _Mij.push_back(0);
     _q.push_back(0);
     _lb.push_back(0);
     _ub.push_back(0);
       */
   
     try {
        stat = this->CreateLMCP(n,
                                        nnz,
                                        &_Mi[0],
                                        &_Mj[0],
                                        &_Mij[0],
                                        &_q[0],
                                        &_lb[0],
                                        &_ub[0],
                                        &_xsol[0],
                                        &_zsol[0],
                                        verbose,
                                        timeLimit);
     } catch (...) {
        throw ZEROException(ZEROErrorCode::SolverError, "PATH threw an exception");
     }
   
   
     if (stat == 1) {
        LOG_S(1) << "Solvers::PathLCP: Found a solution";
        z.zeros(M.n_cols);
        x.zeros(M.n_rows);
        for (unsigned int i = 0; i < _xsol.size(); ++i) {
           z.at(_xmap.at(i)) = _xsol.at(i);
           x.at(_zmap.at(i)) = _zsol.at(i);
        }
        // z.print("z");
        // x.print("x");
        this->status = ZEROStatus::NashEqFound;
     } else {
        LOG_S(1) << "Solvers::PathLCP: No solution found (STATUS=" << std::to_string(status) << ")";
     }
   }
   
   
   void Solvers::PATH::C_bounds(int n, double *x, double *lb, double *ub) {
     int i;
   
     for (i = 0; i < n; i++) {
        x[i]  = this->Problem.x[i];
        lb[i] = this->Problem.lb[i];
        ub[i] = this->Problem.ub[i];
        // std::cout <<"\nx_"<<std::to_string(i)<<" in ["<<std::to_string(lb[i])<<",
        // "<<std::to_string(ub[i])<<"]";
     }
     // std::cout << "\n done\n";
   }
   
   
   int Solvers::PATH::C_function_evaluation(int n, double *x, double *f) {
     int    col, colStart, colEnd, row;
     double value;
   
     for (col = 0; col < n; col++) {
        f[col] = Problem.q[col];
     }
   
     for (col = 0; col < n; col++) {
        value = x[col];
   
        if (value != 0) {
           colStart = Problem.m_start[col] - 1;
           colEnd   = colStart + Problem.m_len[col];
   
           while (colStart < colEnd) {
             row = Problem.m_row[colStart] - 1;
             f[row] += Problem.m_data[colStart] * value;
             colStart++;
           }
        }
     }
   
     return 0;
   }
   
   int Solvers::PATH::C_jacobian_evaluation(int     n,
                                                         double *x,
                                                         int     wantf,
                                                         double *f,
                                                         int    *nnz,
                                                         int    *col_start,
                                                         int    *col_len,
                                                         int    *row,
                                                         double *data) {
     int element;
   
     if (wantf) {
        this->C_function_evaluation(n, x, f);
     }
   
     if (!Filled) {
        for (element = 0; element < Problem.n; element++) {
           col_start[element] = Problem.m_start[element];
           col_len[element]   = Problem.m_len[element];
        }
   
        for (element = 0; element < Problem.nnz; element++) {
           row[element]  = Problem.m_row[element];
           data[element] = Problem.m_data[element];
        }
   
        Filled = true;
     }
   
     *nnz = Problem.nnz;
     return 0;
   }
   
   void *Solvers::PATH::mcp_typ(void *dat, int nnz, int *typ) {
     int i;
   
     for (i = 0; i < nnz; i++) {
        typ[i] = PRESOLVE_LINEAR;
     }
     return dat;
   }
   
   void Solvers::PATH::C_problem_size(int *n, int *nnz) {
     *n   = this->Problem.n;
     *nnz = this->Problem.nnz + 1;
   }
   void Solvers::PATH::bounds(void *dat, int n, double *x, double *lb, double *ub) {
     auto *self = static_cast<Solvers::PATH *>(dat);
     self->C_bounds(n, x, lb, ub);
     return;
   }
   void Solvers::PATH::problem_size(void *dat, int *n, int *nnz) {
     auto *self = static_cast<Solvers::PATH *>(dat);
     self->C_problem_size(n, nnz);
   }
   
   void Solvers::PATH::start(void *dat) {
     auto *self   = static_cast<Solvers::PATH *>(dat);
     self->Filled = 0;
   }
   
   int Solvers::PATH::function_evaluation(void *dat, int n, double *x, double *f) {
     auto *self = static_cast<Solvers::PATH *>(dat);
     return self->C_function_evaluation(n, x, f);
   }
   
   int Solvers::PATH::jacobian_evaluation(void   *dat,
                                                       int     n,
                                                       double *x,
                                                       int     wantf,
                                                       double *f,
                                                       int    *nnz,
                                                       int    *col_start,
                                                       int    *col_len,
                                                       int    *row,
                                                       double *data) {
     auto *self = static_cast<Solvers::PATH *>(dat);
     return self->C_jacobian_evaluation(n, x, wantf, f, nnz, col_start, col_len, row, data);
   }
   
   void *Solvers::PATH::messageCB(void *dat, int mode, char *buf) {
     // std::cout << buf;
     return dat;
   }