Prev Next min_nso_linear.hpp Headings

@(@\newcommand{\W}[1]{ \; #1 \; } \newcommand{\R}[1]{ {\rm #1} } \newcommand{\B}[1]{ {\bf #1} } \newcommand{\D}[2]{ \frac{\partial #1}{\partial #2} } \newcommand{\DD}[3]{ \frac{\partial^2 #1}{\partial #2 \partial #3} } \newcommand{\Dpow}[2]{ \frac{\partial^{#1}}{\partial {#2}^{#1}} } \newcommand{\dpow}[2]{ \frac{ {\rm d}^{#1}}{{\rm d}\, {#2}^{#1}} }@)@
min_nso_linear Source Code
namespace {
     CPPAD_TESTVECTOR(double) min_nso_linear_join(
          const CPPAD_TESTVECTOR(double)& x ,
          const CPPAD_TESTVECTOR(double)& u )
     {     size_t n = x.size();
          size_t s = u.size();
          CPPAD_TESTVECTOR(double) xu(n + s);
          for(size_t j = 0; j < n; j++)
               xu[j] = x[j];
          for(size_t j = 0; j < s; j++)
               xu[n + j] = u[j];
          return xu;
     }
}

// BEGIN C++
namespace CppAD { // BEGIN_CPPAD_NAMESPACE

// BEGIN PROTOTYPE
template <class DblVector, class SizeVector>
bool min_nso_linear(
     size_t           level     ,
     ADFun<double>&   g         ,
     ADFun<double>&   a         ,
     const DblVector& epsilon   ,
     SizeVector       maxitr    ,
     double           b_in      ,
     const DblVector& x_in      ,
     DblVector&       x_out     )
// END PROTOTYPE
{
     using std::fabs;
     //
     // number of absolute value terms
     size_t s  = a.Range();
     //
     // size of domain for f
     size_t n  = g.Domain() - s;
     //
     // size of range space for f
     size_t m = g.Range() - s;
     //
     CPPAD_ASSERT_KNOWN(
          level <= 5,
          "min_nso_linear: level is not less that or equal 5"
     );
     CPPAD_ASSERT_KNOWN(
          size_t(epsilon.size()) == 2,
          "min_nso_linear: size of epsilon not equal to 2"
     );
     CPPAD_ASSERT_KNOWN(
          size_t(maxitr.size()) == 3,
          "min_nso_linear: size of maxitr not equal to 3"
     );
     CPPAD_ASSERT_KNOWN(
          g.Domain() > s && g.Range() > s,
          "min_nso_linear: g, a is not an abs-normal representation"
     );
     CPPAD_ASSERT_KNOWN(
          m == 1,
          "min_nso_linear: m is not equal to 1"
     );
     CPPAD_ASSERT_KNOWN(
          size_t(x_in.size()) == n,
          "min_nso_linear: size of x_in not equal to n"
     );
     CPPAD_ASSERT_KNOWN(
          size_t(x_out.size()) == n,
          "min_nso_linear: size of x_out not equal to n"
     );
     CPPAD_ASSERT_KNOWN(
          epsilon[0] < b_in,
          "min_nso_linear: b_in <= epsilon[0]"
     );
     if( level > 0 )
     {     std::cout << "start min_nso_linear\n";
          std::cout << "b_in = " << b_in << "\n";
          CppAD::abs_print_mat("x_in", n, 1, x_in);
     }
     // level in abs_min_linear sub-problem
     size_t level_tilde = 0;
     if( level > 0 )
          level_tilde = level - 1;
     //
     // maxitr in abs_min_linear sub-problem
     SizeVector maxitr_tilde(2);
     maxitr_tilde[0] = maxitr[1];
     maxitr_tilde[1] = maxitr[2];
     //
     // epsilon in abs_min_linear sub-problem
     DblVector eps_tilde(2);
     eps_tilde[0] = epsilon[0] / 10.;
     eps_tilde[1] = epsilon[1] / 10.;
     //
     // current bound
     double b_cur = b_in;
     //
     // initilaize the current x
     x_out = x_in;
     //
     // value of a(x) at current x
     DblVector a_cur = a.Forward(0, x_out);
     //
     // (x_out, a_cur)
     DblVector xu_cur = min_nso_linear_join(x_out, a_cur);
     //
     // value of g[ x_cur, a_cur ]
     DblVector g_cur = g.Forward(0, xu_cur);
     //
     for(size_t itr = 0; itr < maxitr[0]; itr++)
     {
          // Jacobian of g[ x_cur, a_cur ]
          DblVector g_jac = g.Jacobian(xu_cur);
          //
          // bound in abs_min_linear sub-problem
          DblVector bound_tilde(n);
          for(size_t j = 0; j < n; j++)
               bound_tilde[j] = b_cur;
          //
          DblVector delta_x(n);
          bool ok = abs_min_linear(
               level_tilde, n, m, s,
               g_cur, g_jac, bound_tilde, eps_tilde, maxitr_tilde, delta_x
          );
          if( ! ok )
          {     if( level > 0 )
                    std::cout << "end min_nso_linear: abs_min_linear failed\n";
               return false;
          }
          //
          // new candidate value for x
          DblVector x_new(n);
          double max_delta_x = 0.0;
          for(size_t j = 0; j < n; j++)
          {     x_new[j] = x_out[j] + delta_x[j];
               max_delta_x = std::max(max_delta_x, std::fabs( delta_x[j] ) );
          }
          //
          if( max_delta_x < b_cur && max_delta_x < epsilon[0] )
          {     if( level > 0 )
                    std::cout << "end min_nso_linear: delta_x is near zero\n";
               return true;
          }
          // value of abs-normal approximation at minimizer
          DblVector g_tilde = CppAD::abs_eval(n, m, s, g_cur, g_jac, delta_x);
          //
          double derivative = (g_tilde[0] - g_cur[0]) / max_delta_x;
          CPPAD_ASSERT_UNKNOWN( derivative <= 0.0 )
          if( - epsilon[1] < derivative )
          {     if( level > 0 )
                    std::cout << "end min_nso_linear: derivative near zero\n";
               return true;
          }
          //
          // value of a(x) at new x
          DblVector a_new = a.Forward(0, x_new);
          //
          // (x_new, a_new)
          DblVector xu_new = min_nso_linear_join(x_new, a_new);
          //
          // value of g[ x_new, a_new ]
          DblVector g_new = g.Forward(0, xu_new);
          //
          //
          // check for descent of objective
          double rate_new = (g_new[0] - g_cur[0]) / max_delta_x;
          if( - epsilon[1] < rate_new )
          {     // did not get sufficient descent
               b_cur /= 2.0;
               if( level > 0 )
                    std::cout << "itr = " << itr
                    << ", rate_new = " << rate_new
                    << ", b_cur = " << b_cur << "\n";
               //
          }
          else
          {     // got sufficient descent so accept candidate for x
               x_out  = x_new;
               a_cur  = a_new;
               g_cur  = g_new;
               xu_cur = xu_new;
               //
               if( level >  0 )
               {     std::cout << "itr = " << itr
                    << ", derivative = "<< derivative
                    << ", max_delta_x = "<< max_delta_x
                    << ", objective = " << g_cur[0] << "\n";
                    abs_print_mat("x_out", n, 1, x_out);
               }
          }
     }
     if( level > 0 )
          std::cout << "end min_nso_linear: maximum number of iterations exceeded\n";
     return false;
}
} // END_CPPAD_NAMESPACE

Input File: example/abs_normal/min_nso_linear.omh