MeshKit  1.0
IARoundingFarNlp.cpp
Go to the documentation of this file.
00001 // IARoundingFarNlp.cpp
00002 // Interval Assignment for Meshkit
00003 //
00004 // Adapted from
00005 // Copyright (C) 2005, 2006 International Business Machines and others.
00006 // All Rights Reserved.
00007 // This code is published under the Eclipse Public License.
00008 //
00009 // $Id: hs071_nlp.cpp 1864 2010-12-22 19:21:02Z andreasw $
00010 //
00011 // Authors:  Carl Laird, Andreas Waechter     IBM    2005-08-16
00012 
00013 #include "meshkit/IARoundingFarNlp.hpp"
00014 #include "meshkit/IARoundingNlp.hpp"
00015 #include "meshkit/IAData.hpp"
00016 #include "meshkit/IPData.hpp"
00017 #include "meshkit/IASolution.hpp"
00018 #include "meshkit/IANlp.hpp"
00019 
00020 #include <math.h>
00021 #include <algorithm>
00022 
00023 // for printf
00024 #ifdef HAVE_CSTDIO
00025 # include <cstdio>
00026 #else
00027 # ifdef HAVE_STDIO_H
00028 #  include <stdio.h>
00029 # else
00030 #  error "don't have header file for stdio"
00031 # endif
00032 #endif
00033 
00034 /* Idea: a form of IARoundingNlp with larger variable bounds, but still with a natural integer solution.
00035  x in [1..inf]
00036  xr = x optimal relaxed solution with objective function fnlp, see IANlp.xpp 
00037  f is piecewise linear, with corners at integer values. f slopes are unique (we hope)
00038  Slope definitions
00039  for x between xl = floor xr and xh = ceil xr, we use the difference in fnlp between xl and xh
00040  
00041  case A. xr > ceil g, g is goal I[i]
00042  for x above xh, 
00043  let h+ be fnlp ( xh+1 ) - fnlp ( xh )
00044  let kp be the number of intervals x is above xh 
00045  then slope = floor kp * h+
00046  for x below xl, h- = sqrt(11) / 5 h+, and slope = floor km * h-
00047  all this is weighted by some unique weight
00048  
00049  case B. xr in [ floor g, ceil g] 
00050  h+ = fnlp ( xh+1 ) - fnlp ( xh )
00051  h- = fnlp ( xl-1 ) - fnlp ( xl ), not used if xl == 1
00052  
00053  case C. xr < floor g
00054  h- = fnlp ( xl-1 ) - fnlp ( xl )
00055  h+ = sqrt(10)/5 h-
00056  If g < 2, then h- is unused, and h+ = as in case B
00057 
00058  // representation:
00059  h0 is weights 0..n-1
00060  h+ is weights n..2n-1
00061  h- is weights 2n..
00062  
00063 */
00064 
00065 
00066 // to do
00067 // set h
00068 // uniquify h initially
00069 // randomize h if solution is non-integer
00070 // compute f
00071 // compute f'
00072 // compute f'' == 0
00073 
00074 // clean up
00075 // move the re-used part of IANlp, dealing with constraints but not obj, to some underlying class, without it being a TNLP
00076 
00077 
00078 // the objective function we should use for deriving the weights
00079 double IARoundingFarNlp::f_x_value( double I_i, double x_i ) const
00080 {
00081   // to do revise
00082   return x_i > I_i ? 
00083   (x_i - I_i) / I_i :
00084   1.103402234045 * (I_i - x_i) / x_i; 
00085 }
00086 
00087 // destructor
00088 IARoundingFarNlp::~IARoundingFarNlp() 
00089 {
00090   data = NULL; 
00091   ipData = NULL;
00092   solution = NULL;
00093 }
00094 
00095 // constructor
00096 IARoundingFarNlp::IARoundingFarNlp(const IAData *data_ptr, const IPData *ip_data_ptr, IASolution *solution_ptr): 
00097 data(data_ptr), ipData(ip_data_ptr), solution(solution_ptr), baseNlp(data_ptr, solution_ptr),
00098 debugging(true), verbose(true) // true
00099 {
00100   printf("\nIARoundingFarNlp Problem:\n");
00101   
00102   // weights for function value
00103   h0.resize(data->num_variables());
00104   hp.resize(data->num_variables());
00105   hm.resize(data->num_variables());
00106   for (int i = 0; i < data->num_variables(); ++i)
00107   {
00108     const double x = ipData->relaxedSolution[i];
00109     const double xl = floor(x);
00110     const double xh = xl + 1;
00111     const double g = data->I[i]; // goal
00112     
00113     const double fl = f_x_value(g, xl); 
00114     const double fh = f_x_value(g, xh);
00115     const double w = fh - fl;
00116     h0[i] = w;
00117 
00118     // use sqrt(3) to ensure hp and hm is larger than h0
00119     const double fp = sqrt(3) * (f_x_value(g,xh+1) - f_x_value(g,xh));
00120     double fm = 0.;
00121     
00122     if (xl > 1)
00123     {
00124       fm = sqrt(5) * (f_x_value(g,xl) - f_x_value(g,xl-1));
00125     }
00126     else // if (xl == 1)
00127     {
00128       // unused anyway, assign it something that won't screw up the distribution
00129       fm = - sqrt(5) * fabs(f_x_value(g,2) - f_x_value(g,1));
00130     }
00131 
00132     // case A
00133     if ( x > ceil(g) )
00134     {
00135       hp[i] = fp;
00136       hm[i] = -fp * sqrt(11) / 5.;      
00137       assert( fabs(hm[i]) < fabs(hp[i]) );
00138     }
00139     // case B
00140     else if ( x >= floor(g) )
00141     {
00142       hp[i] = fp;
00143       hm[i] = fm;
00144     }
00145     // case C
00146     else // x < floor(g)
00147     {
00148       hm[i] = fm;
00149       hp[i] = -fm * sqrt(7) / 3.;
00150       assert( fabs(hm[i]) > fabs(hp[i]) );
00151     }
00152         
00153     // invariants for all cases, to produce a super-linear preference for xl or xh
00154     assert( hp[i] > 0 );
00155     assert( hp[i] > h0[i] );
00156     assert( hm[i] < 0 );
00157     assert( hm[i] < h0[i] );
00158 
00159     /*
00160     //zzyk experiment with limited f' discontinuity
00161      // this makes hp < 0 and hm > 0 sometimes, so watch the other asserts
00162     if ( h0[i] > 0 )
00163     {
00164       hp[i] = h0[i]*2.;
00165       hm[i] = h0[i]*0.5;
00166     }
00167     else
00168     {
00169       hp[i] = h0[i]*0.5;
00170       hm[i] = h0[i]*2.;      
00171     }
00172      */
00173     
00174     if (0 && debugging)
00175       printf(": %f %f %f", h0[i], hp[i], hm[i]); // 0th+/i, first should be positive, second negative 
00176   }
00177   
00178   // uniquify and scale the weights = slopes
00179   const double h0_lo = 1;
00180   const double h0_hi = h0_lo * (10 + 2 * data->num_variables());
00181   // uniquify preserving h0 < hm, hp
00182   if (0)
00183   {
00184     // todo experiment with ranges
00185     // want the largest h0 to be less then the smallest hp or hm, in absolute value
00186     // so that a simple rounding from the relaxed solution is optimal, if feasible
00187     // 1.e4? 1.e3 makes the jumps in values wash out for the big chain problem, need some jumps?
00188     // idea, dynamically scale this based on the problem size, h0_hi = 2*num_variables;
00189     const double hpm_lo = h0_hi + 1.;
00190     const double hpm_hi = 1.e5;
00191     h0.uniquify(h0_lo, h0_hi);
00192     
00193     // uniquify and scale hp hm to come after h0
00194     IAWeights hpm();
00195     hpm = hp;                 // hpm = hp
00196     hpm.insert( hpm.end(), hm.begin(), hm.end() ); // hmp += hm
00197     hpm.uniquify(hpm_lo, hpm_hi); 
00198     hp.assign(hpm.begin(), hpm.begin()+hp.size()); // extract hp
00199     hm.assign(hpm.begin()+hp.size(), hpm.end());   // extract hm
00200   }
00201   // uniquify lumping h0 with hm and hp
00202   else {
00203     h0.uniquify(h0, h0_lo, h0_hi);
00204     
00205     // uniquify and scale hp hm to come after h0
00206     IAWeights h0pm;
00207     hpm = h0;                  // h0pm = h0
00208     h0pm.insert( h0pm.end(), hp.begin(), hp.end() ); // h0pm += hp
00209     h0pm.insert( h0pm.end(), hm.begin(), hm.end() ); // h0pm += hm
00210     h0pm.uniquify(h0_lo, h0_hi); 
00211     const size_t sz = data->num_variables();
00212     h0.assign(h0pm.begin(),      h0pm.begin()+sz);   // extract h0
00213     hp.assign(h0pm.begin()+sz,   h0pm.begin()+2*sz); // extract hp
00214     hm.assign(h0pm.begin()+2*sz, h0pm.end());        // extract hm
00215   
00216   }
00217   
00218   if (debugging)
00219   {
00220     printf("\n");
00221     for (int i = 0; i < data->num_variables(); ++i)
00222     {
00223       const double x = ipData->relaxedSolution[i];
00224       printf("x %d (relaxed %f, goal %e), h0 %10.4f hp %10.4f hm %10.4f\n", i, x, data->I[i], h0[i], hp[i], hm[i]);
00225     }
00226   }
00227 }
00228 
00229 
00230 // returns the size of the problem
00231 bool IARoundingFarNlp::get_nlp_info(Index& n, Index& m, Index& nnz_jac_g,
00232                              Index& nnz_h_lag, IndexStyleEnum& index_style)
00233 {
00234   bool base_ok = baseNlp.get_nlp_info(n, m, nnz_jac_g, nnz_h_lag, index_style);
00235   nnz_h_lag = 0; // obj is piecewise linear. We could define h near the corners, but lets try zero first. nele_hess
00236   return true && base_ok;
00237   // need to change this if there are more variables, such as delta-minuses
00238 }
00239 
00240 // returns the variable bounds
00241 bool IARoundingFarNlp::get_bounds_info(Index n, Number* x_l, Number* x_u,
00242                                 Index m, Number* g_l, Number* g_u)
00243 {
00244   const bool base_ok = baseNlp.get_bounds_info(n, x_l, x_u, m, g_l, g_u);
00245 
00246   // test various ranges  
00247   /*
00248   for (Index i=0; i<n; ++i) 
00249   {
00250   */
00251   
00252     /*
00253     // debug: x in [xl-1, xh+1]
00254     const double x = ipData->relaxedSolution[i];
00255     const double xl = floor(x);
00256     assert(xl >= 1.);
00257     if (xl > 1)
00258       x_l[i] = xl-1;
00259     else
00260       x_l[i] = 1;
00261     x_u[i] = xl + 2;
00262     */
00263     
00264     /*
00265     // debug: x in [xl,xh]
00266     x_l[i] = xl;
00267     x_u[i] = xl + 1.;
00268     */
00269     
00270     //if (debugging)
00271     //    printf("x %d (%f) in [%d,%d]\n", i, x, (int) x_l[i], (int) x_u[i]);
00272     
00273   /*
00274   }
00275   */
00276   
00277   return base_ok && true;
00278 }
00279 
00280 // returns the initial point for the problem
00281 bool IARoundingFarNlp::get_starting_point(Index n, bool init_x, Number* x_init,
00282                                    bool init_z, Number* z_L, Number* z_U,
00283                                    Index m, bool init_lambda,
00284                                    Number* lambda)
00285 {
00286   // Minimal info is starting values for x, x_init
00287   // Improvement: we can provide starting values for the dual variables if we wish
00288   assert(init_x == true);
00289   assert(init_z == false);
00290   assert(init_lambda == false);
00291 
00292   // initialize x to the relaxed solution
00293   for (Index i=0; i<n; ++i) 
00294   {
00295     x_init[i] = ipData->relaxedSolution[i];
00296   }
00297 
00298   return true;
00299 }
00300 
00301 inline 
00302 double IARoundingFarNlp::get_f_xl(int i) const
00303 {
00304   return (h0[i] > 0.) ? 0. : -h0[i];
00305 }
00306 
00307 inline 
00308 double IARoundingFarNlp::get_f_xh(int i) const
00309 {
00310   return (h0[i] > 0.) ? h0[i] : 0.;
00311 }
00312 
00313 
00314 // returns the value of the objective function
00315 bool IARoundingFarNlp::eval_f(Index n, const Number* x, bool new_x, Number& obj_value)
00316 {
00317   assert(n == data->num_variables());
00318 
00319   double obj = 0.;
00320   for (Index i = 0; i<n; ++i)
00321   {
00322     // find which interval we're in
00323     const double xl = ipData->get_xl(i);
00324     const double xh = xl + 1;
00325     const double x_i = x[i];  
00326     
00327     double obj_i(-1.); // bad value
00328     if (x_i >xh)  // above xh
00329     {
00330       // the slope of the kth interval is k hp[i]
00331       // f = f(xh) + sum_{k=1:kp} k hp + kp_frac (kp+1) hp
00332       //   = f(xh) + kp (kp+1) / 2  hp + kp_frac (kp+1) hp
00333       //   = f(xh) + hp (kp+1) (kp_frac + kp / 2 )
00334       int kp_int;
00335       double kp_frac;
00336       ipData->get_kp(i, x_i, kp_int, kp_frac);
00337       obj_i = get_f_xh(i) + hp[i] * (kp_int + 1) * (kp_frac + ((double) kp_int) / 2.);
00338       assert(obj_i >= 0.);
00339     }
00340     else if (x_i < xl) // below xl
00341     {
00342       // the slope of the kth interval is k hm[i]
00343       // f = f(xl) + sum_{k=1:km} k hm + km_frac (km+1) hm
00344       //   = f(xl) + km (km+1) / 2  hm + km_frac (km+1) hm
00345       //   = f(xl) + hm (km+1) (k_frac + km / 2 )
00346       int km_int;
00347       double km_frac;
00348       ipData->get_km(i, x_i, km_int, km_frac);
00349       obj_i = get_f_xl(i) + fabs(hm[i]) * (km_int + 1) * (km_frac + ((double) km_int) / 2.); 
00350       assert(obj_i >= 0.);
00351     }
00352     else // middle
00353     { 
00354     // zzyk - testing truncating x within an interval around xl, continuity
00355     /*
00356     static double max_over = 0.;
00357     static double max_under = 0.;
00358     
00359     if (x_i < xl)
00360     {
00361       if ( xl - x_i > max_under )
00362       {
00363         max_under = xl - x_i;
00364         printf("max under %g\n", max_under);
00365       }
00366       x_i = xl;
00367     }
00368     else if (x_i > xh)
00369     {
00370       if ( x_i - xh > max_over )
00371       {
00372         max_over = x_i - xh;
00373         printf("max OVER %g\n", max_over);
00374       }
00375       x_i = xh;      
00376     } */
00377     
00378       obj_i = (h0[i] > 0) ? h0[i] * ( x_i - xl ) : -h0[i] * ( xh - x_i );
00379       assert(obj_i >= 0.);
00380 
00381     } 
00382     assert(obj_i >= 0.);
00383 
00384     obj += obj_i;
00385   }
00386   assert(obj >= 0.);
00387   obj_value = obj;
00388   printf(" %f", obj_value); // debug progress, is this going anywhere?
00389   
00390   return true;
00391 }
00392 
00393 // return the gradient of the objective function grad_{x} f(x)
00394 bool IARoundingFarNlp::eval_grad_f(Index n, const Number* x, bool new_x, Number* grad_f)
00395 {
00396   assert(n == data->num_variables());
00397   for (Index i = 0; i<n; ++i)
00398   {
00399     // find which interval we're in
00400     if (x[i] > ipData->get_xh(i)) // above xh
00401     {
00402       // the slope of the kth interval is k hp[i]
00403       int kp_int;
00404       double kp_frac;
00405       ipData->get_kp(i, x[i], kp_int, kp_frac);
00406       grad_f[i] = (1+kp_int) * hp[i];
00407     }
00408     else if (x[i] < ipData->get_xl(i)) // below xl
00409     {
00410       // the slope of the kth interval is k hm[i]
00411       int km_int;
00412       double km_frac;
00413       ipData->get_km(i, x[i], km_int, km_frac);
00414       grad_f[i] = (1+km_int) * hm[i];
00415     }
00416     else // between xl and xh
00417     { 
00418       // slope is h0
00419       grad_f[i] = h0[i];
00420     }
00421   }
00422   return true;
00423 }
00424 
00425 // return the value of the constraints: g(x)
00426 bool IARoundingFarNlp::eval_g(Index n, const Number* x, bool new_x, Index m, Number* g)
00427 {
00428   return baseNlp.eval_g(n, x, new_x, m, g);
00429 }
00430 
00431 // return the structure or values of the jacobian
00432 bool IARoundingFarNlp::eval_jac_g(Index n, const Number* x, bool new_x,
00433                            Index m, Index nele_jac, Index* iRow, Index *jCol,
00434                            Number* values)
00435 {
00436   return baseNlp.eval_jac_g(n, x, new_x, m, nele_jac, iRow, jCol, values);
00437 }
00438 
00439 //return the structure or values of the hessian
00440 bool IARoundingFarNlp::eval_h(Index n, const Number* x, bool new_x,
00441                        Number obj_factor, Index m, const Number* lambda,
00442                        bool new_lambda, Index nele_hess, Index* iRow,
00443                        Index* jCol, Number* values)
00444 {
00445         // because the constraints are linear
00446   // and the objective function is piecewise linear
00447   // the hessian for this problem is actually empty
00448 
00449   // to do: may need to fill with zeros so ipopt doesn't think it is globally linear
00450   assert(nele_hess == 0); 
00451   
00452   // This is a symmetric matrix, fill the lower left triangle only.
00453   if (values == NULL) {
00454     // return the structure. 
00455     ;
00456   }
00457   else {
00458     // return the values. 
00459     ;
00460   }
00461 
00462   return true;
00463 }
00464 
00465 void IARoundingFarNlp::finalize_solution(SolverReturn status,
00466                                       Index n, const Number* x, const Number* z_L, const Number* z_U,
00467                                       Index m, const Number* g, const Number* lambda,
00468                                       Number obj_value,
00469                                       const IpoptData* ip_data,
00470                                       IpoptCalculatedQuantities* ip_cq)
00471 {
00472   baseNlp.finalize_solution(status, n, x, z_L, z_U, m, g, lambda, obj_value, ip_data, ip_cq);
00473 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines