MeshKit  1.0
IARoundingFar3StepNlp.cpp
Go to the documentation of this file.
00001 // IARoundingFar3StepNlp.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 "IARoundingFar3StepNlp.hpp"
00014 #include "IARoundingNlp.hpp"
00015 #include "IAData.hpp"
00016 #include "IPData.hpp"
00017 #include "IASolution.hpp"
00018 #include "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 = (1 + floor(kp)/2) * 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+ = hp is weights n..2n-1
00061  h- = hm is weights 2n..
00062  
00063  // variable layout
00064  x01 is variables  n..2n-1  : x01_start
00065  xp is variables  2n..3n-1  : xp_start
00066  xm is varaibles  3n..4n-1  : xm_start
00067  
00068  // constrataint layout
00069  sum-equal  
00070  sum-even                : sum_even_start 
00071  x = xl + x01 + xp - xm  : x_constraint_start
00072  or g(j)= x - x01 - xp + xm = xl 
00073 */
00074 
00075 
00076 
00077 // the objective function we should use for deriving the weights
00078 double IARoundingFar3StepNlp::f_x_value( double I_i, double x_i ) const
00079 {
00080   // to do revise
00081   return x_i > I_i ? 
00082   (x_i - I_i) / I_i :
00083   1.103402234045 * (I_i - x_i) / x_i; 
00084 }
00085 
00086 // destructor
00087 IARoundingFar3StepNlp::~IARoundingFar3StepNlp() 
00088 {
00089   data = NULL; 
00090   ipData = NULL;
00091   solution = NULL;
00092 }
00093 
00094 // constructor
00095 IARoundingFar3StepNlp::IARoundingFar3StepNlp(const IAData *data_ptr, const IPData *ip_data_ptr, IASolution *solution_ptr): 
00096 data(data_ptr), ipData(ip_data_ptr), solution(solution_ptr), baseNlp(data_ptr, solution_ptr),
00097 x01_start( data_ptr->num_variables() ), 
00098 xp_start( 2 * data_ptr->num_variables() ), 
00099 xm_start( 3 * data_ptr->num_variables() ),
00100 base_n ( data_ptr->num_variables() ),
00101 base_m ( (int) ( data_ptr->constraints.size() + data->sumEvenConstraints.size() ) ),
00102 problem_n( (int) 4 * data_ptr->num_variables() ),
00103 problem_m( (int) (data_ptr->constraints.size() + data_ptr->sumEvenConstraints.size() + data_ptr->num_variables())),
00104 sum_even_start( (int) data_ptr->constraints.size()), 
00105 x_constraint_start( (int) (data_ptr->constraints.size() + data_ptr->sumEvenConstraints.size())), 
00106 hess_option(ZERO),
00107 debugging(true), verbose(true) // true
00108 {
00109   printf("\nIARoundingFar3StepNlp Problem:\n");
00110   
00111   
00112   // weights for function value
00113   h0.resize(data->num_variables());
00114   hp.resize(data->num_variables());
00115   hm.resize(data->num_variables());
00116   for (int i = 0; i < data->num_variables(); ++i)
00117   {
00118     const double x = ipData->relaxedSolution[i];
00119     const double xl = floor(x);
00120     const double xh = xl + 1;
00121     const double g = data->I[i]; // goal
00122     
00123     const double fl = f_x_value(g, xl); 
00124     const double fh = f_x_value(g, xh);
00125     const double w = fh - fl;
00126     h0[i] = w;
00127 
00128     // use sqrt(3) to ensure hp and hm is larger than h0
00129     const double fp = sqrt(3) * (f_x_value(g,xh+1) - f_x_value(g,xh));
00130     double fm = 0.;
00131     
00132     if (xl > 1)
00133     {
00134       fm = sqrt(5) * (f_x_value(g,xl) - f_x_value(g,xl-1));
00135     }
00136     else // if (xl == 1)
00137     {
00138       // unused anyway, assign it something that won't screw up the distribution
00139       fm = - sqrt(5) * fabs(f_x_value(g,2) - f_x_value(g,1));
00140     }
00141 
00142     // case A
00143     if ( x > ceil(g) )
00144     {
00145       hp[i] = fp;
00146       hm[i] = -fp * sqrt(11) / 5.;      
00147       assert( fabs(hm[i]) < fabs(hp[i]) );
00148     }
00149     // case B
00150     else if ( x >= floor(g) )
00151     {
00152       hp[i] = fp;
00153       hm[i] = fm;
00154     }
00155     // case C
00156     else // x < floor(g)
00157     {
00158       hm[i] = fm;
00159       hp[i] = -fm * sqrt(7) / 3.;
00160       assert( fabs(hm[i]) > fabs(hp[i]) );
00161     }
00162         
00163     // invariants for all cases, to produce a super-linear preference for xl or xh
00164     assert( hp[i] > 0 );
00165     assert( hp[i] > h0[i] );
00166     assert( hm[i] < 0 );
00167     assert( hm[i] < h0[i] );
00168 
00169     /*
00170     //zzyk experiment with limited f' discontinuity
00171      // this makes hp < 0 and hm > 0 sometimes, so watch the other asserts
00172     if ( h0[i] > 0 )
00173     {
00174       hp[i] = h0[i]*2.;
00175       hm[i] = h0[i]*0.5;
00176     }
00177     else
00178     {
00179       hp[i] = h0[i]*0.5;
00180       hm[i] = h0[i]*2.;      
00181     }
00182      */
00183     
00184     if (0 && debugging)
00185       printf(": %f %f %f", h0[i], hp[i], hm[i]); // 0th+/i, first should be positive, second negative 
00186   }
00187   
00188   // uniquify and scale the weights = slopes
00189   const double h0_lo = 1;
00190   const double h0_hi = h0_lo * (10 + 2 * data->num_variables());
00191   // uniquify preserving h0 < hm, hp
00192   if (0)
00193   {
00194     // todo experiment with ranges
00195     // want the largest h0 to be less then the smallest hp or hm, in absolute value
00196     // so that a simple rounding from the relaxed solution is optimal, if feasible
00197     // 1.e4? 1.e3 makes the jumps in values wash out for the big chain problem, need some jumps?
00198     // idea, dynamically scale this based on the problem size, h0_hi = 2*num_variables;
00199     const double hpm_lo = h0_hi + 1.;
00200     const double hpm_hi = 1.e5;
00201     h0.uniquify(h0_lo, h0_hi);
00202     
00203     // uniquify and scale hp hm to come after h0
00204     IAWeights hpm;
00205     hpm.weightvec = hp.weightvec;                  // hpm = hp
00206     hpm.weightvec.insert( hpm.weightvec.end(), hm.weightvec.begin(), hm.weightvec.end() ); // hmp += hm
00207     hpm.uniquify_weights(hpm_lo, hpm_hi); 
00208     hp.weightvec.assign(hpm.weightvec.begin(), hpm.weightvec.begin()+hp.size()); // extract hp
00209     hm.weightvec.assign(hpm.weightvec.begin()+hp.weightvec.size(), hpm.weightvec.end());   // extract hm
00210   }
00211   // uniquify lumping h0 with hm and hp
00212   else {
00213     h0.uniquify(h0, h0_lo, h0_hi);
00214     
00215     // uniquify and scale hp hm to come after h0
00216     IAWeights h0pm( h0 );                  // h0pm = h0
00217     h0pm.insert( h0pm.end(), hp.begin(), hp.end() ); // h0pm += hp
00218     h0pm.insert( h0pm.end(), hm.begin(), hm.end() ); // h0pm += hm
00219     h0pm.uniquify_weights(h0_lo, h0_hi); 
00220     const size_t sz = data->num_variables();
00221     h0.assign(h0pm.begin(),      h0pm.begin()+sz);   // extract h0
00222     hp.assign(h0pm.begin()+sz,   h0pm.begin()+2*sz); // extract hp
00223     hm.assign(h0pm.begin()+2*sz, h0pm.end());        // extract hm
00224   
00225   }
00226   
00227   if (debugging)
00228   {
00229     printf("\n");
00230     for (int i = 0; i < data->num_variables(); ++i)
00231     {
00232       const double x = ipData->relaxedSolution[i];
00233       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]);
00234     }
00235   }
00236 }
00237 
00238 bool IARoundingFar3StepNlp::randomize_weights_of_non_int()
00239 {
00240   // either they all will be true or they all will be false, because
00241   // true means a solution value was non-integer
00242   IASolverToolInt sti(data, solution);
00243   bool changed = 
00244         randomize_weights_of_non_int( h0,  0.03124234 ) &&
00245     randomize_weights_of_non_int( hm,  0.02894834 ) &&
00246         randomize_weights_of_non_int( hp,  0.02745675 );
00247   
00248   return changed;
00249 }
00250 
00251 
00252 // returns the size of the problem
00253 bool IARoundingFar3StepNlp::get_nlp_info(Index& n, Index& m, Index& nnz_jac_g,
00254                              Index& nnz_h_lag, IndexStyleEnum& index_style)
00255 {
00256   bool base_ok = baseNlp.get_nlp_info(n, m, nnz_jac_g, nnz_h_lag, index_style);
00257   
00258   m = problem_m;
00259   n = problem_n;
00260 
00261   // nnz_jac_g == nele_jac in eval_jac_g
00262   // four coefficients for each g() = x - x01 - xp + xm = xl const constraint
00263   nnz_jac_g += 4 * base_n; 
00264   
00265   // nnz_h_lag == nele_hess in eval_h
00266   if (hess_option == ZERO)
00267     nnz_h_lag = 0.;
00268   else if (hess_option == ROUNDED)
00269     nnz_h_lag =  2 * data->num_variables(); // one for each xp, xm
00270   
00271   return true && base_ok;
00272   // need to change this if there are more variables, such as delta-minuses
00273 }
00274 
00275 // returns the variable bounds
00276 bool IARoundingFar3StepNlp::get_bounds_info(Index n, Number* x_l, Number* x_u,
00277                                 Index m, Number* g_l, Number* g_u)
00278 {
00279   const bool base_ok = baseNlp.get_bounds_info(base_n, x_l, x_u, base_m, g_l, g_u);
00280   
00281   for (Index i=0; i<data->num_variables(); ++i)
00282   {
00283     // x01
00284     const int x01i = x01_start + i;
00285     x_l[x01i] = 0;
00286     x_u[x01i] = 1;
00287     
00288     // xp
00289     const int xpi = xp_start + i;
00290     x_l[xpi] = 0;
00291     x_u[xpi] = MESHKIT_IA_upperUnbound;
00292     
00293     // xm
00294     const int xmi = xm_start + i;
00295     x_l[xmi] = 0;
00296     x_u[xmi] = ipData->get_xl(i) - 1;
00297     assert( x_u[xmi] >= 0. );
00298   }
00299 
00300   // x = xl + x01 + xp - xm  : x_constraint_start
00301   // or g(j)= x - x01 - xp + xm = xl 
00302   for (Index j=0; j<data->num_variables(); ++j)
00303   {
00304     const int xcj = x_constraint_start + j;
00305     const int xl = ipData->get_xl(j);
00306     g_l[xcj] = xl; 
00307     g_u[xcj] = xl;
00308   }
00309   
00310   return base_ok && true;
00311 }
00312 
00313 // returns the initial point for the problem
00314 bool IARoundingFar3StepNlp::get_starting_point(Index n, bool init_x, Number* x_init,
00315                                    bool init_z, Number* z_L, Number* z_U,
00316                                    Index m, bool init_lambda,
00317                                    Number* lambda)
00318 {
00319   // Minimal info is starting values for x, x_init
00320   // Improvement: we can provide starting values for the dual variables if we wish
00321   assert(init_x == true);
00322   assert(init_z == false);
00323   assert(init_lambda == false);
00324 
00325   // initialize x to the relaxed solution
00326   for (Index i=0; i<base_n; ++i) 
00327   {
00328     const double xr = ipData->relaxedSolution[i];
00329     x_init[i] = xr;
00330     x_init[x01_start + i] = xr - floor(xr);
00331     x_init[xp_start + i] = 0.;
00332     x_init[xm_start + i] = 0.;
00333   }
00334 
00335   return true;
00336 }
00337 
00338 inline 
00339 double IARoundingFar3StepNlp::get_f_xl(int i) const
00340 {
00341   return (h0[i] > 0.) ? 0. : -h0[i];
00342 }
00343 
00344 inline 
00345 double IARoundingFar3StepNlp::get_f_xh(int i) const
00346 {
00347   return (h0[i] > 0.) ? h0[i] : 0.;
00348 }
00349 
00350 
00351 // returns the value of the objective function
00352 bool IARoundingFar3StepNlp::eval_f(Index n, const Number* x, bool new_x, Number& obj_value)
00353 {
00354   assert(n == problem_n);
00355 
00356   double obj = 0.;
00357   for (Index i = 0; i<base_n; ++i)
00358   {
00359     // x itself contributes nothing, just its components x01, xp, xm
00360 
00361     // xp
00362     {
00363       const double xpv = x[xp_start + i];
00364       int kp_int;
00365       double kp_frac;
00366       IPData::get_frac( xpv, kp_int, kp_frac );
00367       const double objp = hp[i] * ( kp_int * ( 1. + (kp_int + 1.) / 4. ) + kp_frac * ( 1. + kp_int / 2. ) );
00368       assert( xpv >= 0. );
00369       assert(objp >= 0.);    
00370       obj += objp;
00371     }
00372 
00373     // xm
00374     {
00375       const double xmv = x[xm_start + i];
00376       int km_int;
00377       double km_frac;
00378       IPData::get_frac( xmv, km_int, km_frac );
00379       const double objm = -hm[i] * ( km_int * ( 1. + (km_int + 1.) / 4. ) + km_frac * ( 1. + km_int / 2. ) );
00380       assert( xmv >= 0. );
00381       assert(objm >= 0.);    
00382       obj += objm;
00383     }
00384 
00385     // x01
00386     {
00387       const double x01v = x[x01_start + i];
00388       const double obj01 = (h0[i] > 0) ? h0[i] * x01v : h0[i] * ( x01v - 1. );
00389       assert(x01v >= 0.);
00390       assert(x01v <= 1.);
00391       assert(obj01 >= 0.);
00392       obj += obj01;
00393     }
00394   }
00395   assert(obj >= 0.);
00396   obj_value = obj;
00397   printf(" %f", obj_value); // debug progress, is this going down?
00398   
00399   return true;
00400 }
00401 
00402 // return the gradient of the objective function grad_{x} f(x)
00403 bool IARoundingFar3StepNlp::eval_grad_f(Index n, const Number* x, bool new_x, Number* grad_f)
00404 {
00405   assert(n == problem_n);
00406   for (Index i = 0; i<base_n; ++i)
00407   {
00408     // x
00409     grad_f[i] = 0;
00410     
00411     // xp
00412     {
00413       const double xpv = x[xp_start + i];
00414       int kp_int;
00415       double kp_frac;
00416       IPData::get_frac( xpv, kp_int, kp_frac );    
00417       grad_f[xp_start+i] = hp[i] * ( 1. + kp_int / 2. );
00418     }
00419     
00420     // xm
00421     {
00422       const double xmv = x[xm_start + i];
00423       int km_int;
00424       double km_frac;
00425       IPData::get_frac( xmv, km_int, km_frac );
00426       grad_f[xm_start+i] = -hm[i] * ( 1. + km_int / 2. );
00427     }
00428     
00429     // x01
00430     grad_f[x01_start+i] = h0[i];
00431   }
00432   
00433   return true;
00434 }
00435 
00436 // return the value of the constraints: g(x)
00437 bool IARoundingFar3StepNlp::eval_g(Index n, const Number* x, bool new_x, Index m, Number* g)
00438 {
00439   const bool base_ok = baseNlp.eval_g(base_n, x, new_x, base_m, g);
00440   
00441   // x_constraints
00442   for (int i=0; i<base_n; ++i)
00443   {
00444     // g(k)= x - x01 - xp + xm = xl 
00445     const double xv = x[i];
00446     const double x01v = x[x01_start+i];
00447     const double xpv = x[xp_start+i];
00448     const double xmv = x[xm_start+i];
00449     const double g_k = xv - x01v - xpv + xmv;
00450     const int k = x_constraint_start + i;
00451     g[k] = g_k;    
00452   }
00453   
00454   return base_ok && true;
00455 }
00456 
00457 // return the structure or values of the jacobian
00458 bool IARoundingFar3StepNlp::eval_jac_g(Index n, const Number* x, bool new_x,
00459                            Index m, Index nele_jac, Index* iRow, Index *jCol,
00460                            Number* values)
00461 {
00462   int base_nele_jac = baseNlp.get_neleJac();
00463   bool base_ok = baseNlp.eval_jac_g(base_n, x, new_x, base_m, base_nele_jac, iRow, jCol, values);
00464   base_nele_jac = baseNlp.get_neleJac();
00465   int k = base_nele_jac;
00466 
00467   if (values == NULL) 
00468   {
00469     // return the structure of the jacobian
00470 
00471     // x_constraints
00472     for (int i=0; i<base_n; ++i)
00473     {
00474       // g(k)= x - x01 - xp + xm = xl const
00475       int xi = i + x_constraint_start;
00476       iRow[k] = xi;
00477       jCol[k] = i;             // x
00478       ++k;
00479       iRow[k] = xi;
00480       jCol[k] = x01_start + i; // x01
00481       ++k;
00482       iRow[k] = xi;
00483       jCol[k] = xp_start + i;  // xp
00484       ++k;
00485       iRow[k] = xi;
00486       jCol[k] = xm_start + i;  // xm
00487       ++k;
00488     }
00489   }
00490   else
00491   {
00492     // return the values of the jacobian of the constraints
00493 
00494     // x_constraints
00495     for (int i=0; i<base_n; ++i)
00496     {
00497       // g(k)= x - x01 - xp + xm = xl const
00498       values[k++] =  1;  // x
00499       values[k++] = -1;  // x01
00500       values[k++] = -1;  // xp
00501       values[k++] =  1;  // xm
00502     }
00503 
00504   }
00505   return base_ok && true;
00506 }
00507 
00508 //return the structure or values of the hessian
00509 bool IARoundingFar3StepNlp::eval_h(Index n, const Number* x, bool new_x,
00510                        Number obj_factor, Index m, const Number* lambda,
00511                        bool new_lambda, Index nele_hess, Index* iRow,
00512                        Index* jCol, Number* values)
00513 {
00514         // because the constraints are linear
00515   // the hessian for them is actually empty
00516 
00517   // get_nlp_info specified the number of non-zeroes, nele_hess
00518 
00519   
00520   // This is a symmetric matrix, fill the lower left triangle only.
00521   if (values == NULL) {
00522     // return the structure. 
00523     if (hess_option == ZERO)
00524     {
00525       ; // empty
00526     }
00527     else if (hess_option == ROUNDED)
00528     {
00529       // Since the objective function is separable, only the diagonal is non-zero
00530       int k = 0;
00531       for ( int i = 0; i < base_n; ++i )
00532       {
00533         iRow[k] = xp_start + i;
00534         jCol[k] = xp_start + i;
00535         ++k;
00536       }
00537       for ( int i = 0; i < base_n; ++i )
00538       {
00539         iRow[k] = xm_start + i;
00540         jCol[k] = xm_start + i;
00541         ++k;
00542       }
00543     }
00544   }
00545   else {
00546     // return the values. 
00547     int k = 0;
00548     if (hess_option == ZERO)
00549     {
00550       ; // empty
00551     }
00552     else if (hess_option == ROUNDED)
00553     {
00554       for ( int i = 0; i < base_n; ++i )
00555       {
00556         values[k] = hp[i] / 2.; 
00557         // but linearly fade to zero towards xh
00558         const double xp = x[xp_start + i];
00559         if (xp < 0.5)
00560           values[k] *= 2. * xp;
00561         
00562         ++k;
00563       }
00564       for ( int i = 0; i < base_n; ++i )
00565       {
00566         values[k] = -hm[i] / 2.; 
00567         const double xm = x[xm_start + i];
00568         if (xm < 0.5)
00569           values[k] *= 2. * xm;
00570 
00571         ++k;
00572       }
00573       
00574     }
00575     assert( k == nele_hess );
00576   }
00577 
00578   
00579   return true;
00580 }
00581 
00582 void IARoundingFar3StepNlp::finalize_solution(SolverReturn status,
00583                                       Index n, const Number* x, const Number* z_L, const Number* z_U,
00584                                       Index m, const Number* g, const Number* lambda,
00585                                       Number obj_value,
00586                                       const IpoptData* ip_data,
00587                                       IpoptCalculatedQuantities* ip_cq)
00588 {
00589   baseNlp.finalize_solution(status, base_n, x, z_L, z_U, base_m, g, lambda, obj_value, ip_data, ip_cq);
00590   // ignore the other constraints and variables, as they were just bookkeepings.
00591 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines