Line data    Source code

       1             : #pragma once
       2             : 
       3             : #include <panoc-alm/decl/alm.hpp>
       4             : #include <stdexcept>
       5             : 
       6             : namespace pa::detail {
       7             : 
       8          30 : inline void project_y(rvec y,     // inout
       9             :                       crvec z_lb, // in
      10             :                       crvec z_ub, // in
      11             :                       real_t M    // in
      12             : ) {
      13             :     // TODO: Handle NaN correctly
      14         155 :     auto max_lb = [M](real_t y, real_t z_lb) {
      15         125 :         real_t y_lb = z_lb == -inf ? 0 : -M;
      16         125 :         return std::max(y, y_lb);
      17         125 :     };
      18         155 :     auto min_ub = [M](real_t y, real_t z_ub) {
      19         125 :         real_t y_ub = z_ub == inf ? 0 : M;
      20         125 :         return std::min(y, y_ub);
      21         125 :     };
      22          30 :     y = y.binaryExpr(z_lb, max_lb).binaryExpr(z_ub, min_ub);
      23          30 : }
      24             : 
      25          26 : inline void update_penalty_weights(const ALMParams &params, real_t Δ,
      26             :                                    bool first_iter, rvec e, rvec old_e,
      27             :                                    real_t norm_e, real_t old_norm_e,
      28             :                                    crvec Σ_old, rvec Σ) {
      29          26 :     if (norm_e <= params.δ) {
      30           3 :         Σ = Σ_old;
      31           3 :         return;
      32             :     }
      33          23 :     if (params.single_penalty_factor) {
      34           0 :         if (first_iter || norm_e > params.θ * old_norm_e) {
      35           0 :             real_t new_Σ = std::fmin(params.Σ_max, Δ * Σ_old(0));
      36           0 :             Σ.setConstant(new_Σ);
      37           0 :         } else {
      38           0 :             Σ = Σ_old;
      39             :         }
      40           0 :     } else {
      41         136 :         for (Eigen::Index i = 0; i < e.rows(); ++i) {
      42         113 :             if (first_iter || std::abs(e(i)) > params.θ * std::abs(old_e(i))) {
      43         130 :                 Σ(i) = std::fmin(params.Σ_max,
      44         130 :                                  std::fmax(Δ * std::abs(e(i)) / norm_e, 1) *
      45          65 :                                      Σ_old(i));
      46          65 :             } else {
      47          48 :                 Σ(i) = Σ_old(i);
      48             :             }
      49         113 :         }
      50             :     }
      51          26 : }
      52             : 
      53           1 : inline void initialize_penalty(const Problem &p, const ALMParams &params,
      54             :                                crvec x0, rvec Σ) {
      55           1 :     real_t f0 = p.f(x0);
      56           1 :     vec g0(p.m);
      57           1 :     p.g(x0, g0);
      58             :     // TODO: reuse evaluations of f ang g in PANOC?
      59           2 :     real_t σ = params.σ₀ * std::max(real_t(1), std::abs(f0)) /
      60           1 :                std::max(real_t(1), 0.5 * g0.squaredNorm());
      61           1 :     σ = std::max(σ, params.Σ_min);
      62           1 :     σ = std::min(σ, params.Σ_max);
      63           1 :     Σ.fill(σ);
      64           1 : }
      65             : 
      66           0 : inline void apply_preconditioning(const Problem &problem, Problem &prec_problem,
      67             :                                   crvec x, real_t &prec_f, vec &prec_g) {
      68           0 :     vec grad_f(problem.n);
      69           0 :     vec v = vec::Zero(problem.m);
      70           0 :     vec grad_g(problem.n);
      71           0 :     problem.grad_f(x, grad_f);
      72           0 :     prec_f = 1. / std::max(grad_f.lpNorm<Eigen::Infinity>(), real_t{1});
      73           0 :     prec_g.resize(problem.m);
      74             : 
      75           0 :     for (Eigen::Index i = 0; i < problem.m; ++i) {
      76           0 :         v(i) = 1;
      77           0 :         problem.grad_g_prod(x, v, grad_g);
      78           0 :         v(i)      = 0;
      79           0 :         prec_g(i) = 1. / std::max(grad_g.lpNorm<Eigen::Infinity>(), real_t{1});
      80           0 :     }
      81             : 
      82           0 :     auto prec_f_fun      = [&](crvec x) { return problem.f(x) * prec_f; };
      83           0 :     auto prec_grad_f_fun = [&](crvec x, rvec grad_f) {
      84           0 :         problem.grad_f(x, grad_f);
      85           0 :         grad_f *= prec_f;
      86           0 :     };
      87           0 :     auto prec_g_fun = [&](crvec x, rvec g) {
      88           0 :         problem.g(x, g);
      89           0 :         g = prec_g.asDiagonal() * g;
      90           0 :     };
      91           0 :     auto prec_grad_g_fun = [&](crvec x, crvec v, rvec grad_g) {
      92           0 :         vec prec_v = prec_g.asDiagonal() * v;
      93           0 :         problem.grad_g_prod(x, prec_v, grad_g);
      94           0 :     };
      95           0 :     prec_problem = Problem{
      96           0 :         problem.n,
      97           0 :         problem.m,
      98           0 :         problem.C,
      99           0 :         Box{prec_g.asDiagonal() * problem.D.upperbound,
     100           0 :             prec_g.asDiagonal() * problem.D.lowerbound},
     101           0 :         std::move(prec_f_fun),
     102           0 :         std::move(prec_grad_f_fun),
     103           0 :         std::move(prec_g_fun),
     104           0 :         std::move(prec_grad_g_fun),
     105           0 :         [](crvec, unsigned, rvec) {
     106           0 :             throw std::logic_error("Preconditioning for second-order solvers "
     107             :                                    "has not yet been implemented");
     108           0 :         },
     109           0 :         [](crvec, crvec, crvec, rvec) {
     110           0 :             throw std::logic_error("Preconditioning for second-order solvers "
     111             :                                    "has not yet been implemented");
     112           0 :         },
     113           0 :         [](crvec, crvec, rmat) {
     114           0 :             throw std::logic_error("Preconditioning for second-order solvers "
     115             :                                    "has not yet been implemented");
     116           0 :         },
     117             :     };
     118           0 : }
     119             : 
     120             : } // namespace pa::detail