MeshKit  1.0
IASolverEven.cpp
Go to the documentation of this file.
00001 // IASolverEven.cpp
00002 // Interval Assignment for Meshkit
00003 //
00004 #include "IASolverEven.hpp"
00005 #include "IAMINlp.hpp"
00006 #include "IPData.hpp"
00007 // #include "IAMilp.hpp"
00008 
00009 #include <stdio.h>
00010 #include <math.h>
00011 #include <limits.h>
00012 
00013 #include "IpIpoptApplication.hpp"
00014 
00015 namespace MeshKit 
00016 {
00017     
00018 IASolverEven::IASolverEven(const IAData *ia_data, const IASolution *int_solution) :
00019 iaData(ia_data), debugging(true) 
00020 { 
00021   ipData.initialize(int_solution->x_solution);
00022 }
00023 
00025 IASolverEven::~IASolverEven() { iaData=NULL; }
00026 
00027 
00028 double IASolverEven::distance_to_even(int i)
00029 {
00030   double sum = sum_even_value(i);
00031   if (is_even(sum))
00032     return 0.;
00033   return sum;
00034   /*
00035   double lower_even = 2. * floor( sum / 2.);
00036   if (lower_even < 4.) 
00037     lower_even = 4.;
00038   const double upper_even = 2. * ceil( sum / 2. );
00039   const double u = upper_even - sum;
00040   const double l = sum - lower_even;
00041   if ( u > l )
00042     return l;
00043   else
00044     return u;  
00045    */
00046 }
00047 
00048 bool IASolverEven::calculate_sumeven_value(const int i, double &obj_increase, int &y_int )
00049 {
00050   // to do: make a version of this where x can flip which side of the goal it is on.
00051   // the issue is e.g. goal = 1.1, x gets rounded to 1, then x can never flip to the other side of 1.1! but if goal was 1, then it could
00052   const double x = x_solution[i];
00053   y_int = floor(x + 0.5); // round to nearest integer - this version assumes we are already integer within roundoff.
00054   // round x away from I
00055   const double goal = iaData->I[i]; // shorthand reminder
00056   if (x>goal)
00057     ++y_int;
00058   else 
00059     --y_int;
00060   if (y_int >= 1)
00061   {
00062     // compare to relaxed solution, not the latest pseudo-integer solution
00063     double f = myianlp->eval_R_i( goal, ipData.relaxedSolution[i] ); 
00064     double g = myianlp->eval_R_i( goal, y_int );
00065     obj_increase = g - f;           
00066     // assert( obj_increase > 0. ); this is only invariant if we use f = R(x), not f = R(relaxed).
00067     return true;
00068   }
00069   //else
00070   obj_increase = 0.;
00071   return false;
00072 
00073 }
00074 
00075 bool IASolverEven::calculate_rounding_value(const int i, double &obj_increase, int &y_int )
00076 {
00077   const double x = x_solution[i];
00078   y_int = floor(x + 0.5); // round to nearest integer
00079   double y = y_int;
00080   if ( fabs(x-y) > 1e-4) // not already integer within tolerance 
00081   {
00082     // non-integer found
00083     // round x away from I
00084     const double goal = iaData->I[i];
00085     // it was forced away from its goal, round in the forced direction further
00086     // else leave it rounded to the nearest integer
00087     if (fabs(x - goal) > 1.e-2) 
00088     {
00089       y_int = (x>goal) ? ceil(x) : floor(x);
00090       y = y_int;
00091     }
00092     assert(y >= 1.);
00093     
00094     // if it is already constrained to y_int, then consider it as already satisfying it.
00095     if ( ipData.varIntegerBound[i] == y )
00096     {
00097       obj_increase = 0.;
00098       return false;
00099     }
00100 
00101     // todo: if x is very close to the goal, maybe we should retain freedom to round it the other direction, try both, etc.
00102  
00103     // is this the best non-integer we've found so far?
00104     // best means the one that changing it to integer increases the obj function the least
00105     // double f = myianlp->eval_R_i( goal, x ); 
00106     // compare to relaxed solution, not the latest pseudo-integer solution
00107     double f = myianlp->eval_R_i( goal, ipData.relaxedSolution[i] ); 
00108     double g = myianlp->eval_R_i( goal, y );
00109     obj_increase = g - f;           
00110     // assert( obj_increase > 0. ); this is only invariant if we use f = R(x), not f = R(relaxed).
00111     return true;
00112   }
00113   //else
00114   obj_increase = 0.;
00115   return false;
00116 }
00117 
00118 bool IASolverEven::find_one_non_integer(int &i_nonint, int &x_bound)
00119 {
00120   i_nonint = -1; // none found
00121   double obj_increase_nonint( std::numeric_limits<double>::max() );
00122   // find the one variable whose rounding causes the least change to the objective function
00123   // to do: alt is to pick one at random, with some chance
00124   for (int i = 0; i<iaData->num_variables(); ++i)
00125   {
00126     double obj_increase;
00127     int y_int;
00128     if (calculate_rounding_value(i, obj_increase, y_int))
00129     {
00130       if (obj_increase < obj_increase_nonint)
00131       {
00132         i_nonint = i;
00133         x_bound = y_int;
00134         obj_increase_nonint = obj_increase;
00135       }
00136     }
00137   }
00138   if (i_nonint > -1)
00139     return true;
00140   // everybody is already integer, enforce some sum-even constraint
00141   // pick an un-satisfied sum-even constraint
00142   int num_noneven_found = 0;
00143   // algorithm parameter
00144   // sum_strategy = 0 pick biggest, 1: uniformly at random 
00145   // to do: test these further. Are there bad cases where randomness helps? Some random sequences give poor-quality solutions, rounding a 1.5 up to 3 or 4.
00146   const int sum_strategy = 1; // rand() % 2; // 0
00147   printf("sum_strategy=%d ", sum_strategy);
00148   double biggest_non_even_value = 0.;
00149   int biggest_non_even_i = -1;
00150   for (unsigned int i = 0; i < iaData->sumEvenConstraints.size(); ++i)
00151   {
00152     const double d = distance_to_even(i);
00153     if (d > 0.)
00154     {
00155       ++num_noneven_found;
00156       if ((sum_strategy == 0 && d > biggest_non_even_value) || 
00157           (sum_strategy == 1 && (rand() % num_noneven_found == 0)))
00158       {
00159         biggest_non_even_value = d;
00160         biggest_non_even_i = i;
00161       }
00162     }
00163   }
00164   // now pick an curve of the sum-even-constraint to round
00165   if (biggest_non_even_i > -1)
00166   {
00167     const int c = biggest_non_even_i;
00168     // algorithm parameter
00169     // 0: pick biggest non-even , 1: uniformly at random 
00170     // to do: test these further
00171     const int curve_strategy = 1; // rand() % 2; // 0
00172     printf("curve_strategy=%d\n", curve_strategy);
00173     int num_found = 0;
00174     double smallest_increase = std::numeric_limits<double>::max();
00175     for (unsigned int i = 0; i < iaData->sumEvenConstraints[c].M.size(); ++i)
00176     {
00177       int j = iaData->sumEvenConstraints[c].M[i].col;
00178       assert( iaData->sumEvenConstraints[c].M[i].val == 1. );
00179       const double x = x_solution[j];
00180       if (!is_even(x))
00181       {
00182         ++num_found;
00183         int y_int;
00184         double obj_increase;
00185         if (calculate_sumeven_value(j, obj_increase, y_int))
00186         {
00187           if ((curve_strategy == 0 && obj_increase < smallest_increase) || 
00188               (curve_strategy == 1 && (rand() % num_found == 0)))
00189           {
00190             i_nonint = j;
00191             x_bound = y_int;
00192             smallest_increase = obj_increase;
00193           }
00194         }
00195       }
00196     }
00197     // if no non-even curves were found, pick an odd one
00198     if (i_nonint == -1)
00199     {
00200       for (unsigned int i = 0; i < iaData->sumEvenConstraints[c].M.size(); ++i)
00201       {
00202         int j = iaData->sumEvenConstraints[c].M[i].col;
00203         assert( iaData->sumEvenConstraints[c].M[i].val == 1. );
00204         ++num_found;
00205         int y_int;
00206         double obj_increase;
00207         if (calculate_sumeven_value(j, obj_increase, y_int))
00208         {
00209           if (curve_strategy == 0 && obj_increase < smallest_increase || 
00210               curve_strategy == 1 && (rand() % num_found == 0))
00211           {
00212             i_nonint = j;
00213             x_bound = y_int;
00214             smallest_increase = obj_increase;
00215           }
00216         }
00217       }
00218     }
00219     if (i_nonint == -1)
00220       // failed to achieve sum-even solution
00221       return false;
00222     return true;
00223   }
00224   return false;
00225 }
00226  
00227 void IASolverEven::constrain_integer(const int i_nonint, const int x_bound)
00228 {
00229   ipData.oldBound[i_nonint] = ipData.varIntegerBound[i_nonint];;
00230   ipData.varIntegerBound[i_nonint] = (double) x_bound;
00231 }
00232 
00233 // find the one variable whose rounding causes the least change to the objective function
00234 // unused at the moment
00235 bool IASolverEven::constrain_one_non_integer( int &i, int &b)
00236 {
00237   if (find_one_non_integer(i, b))
00238   {
00239     constrain_integer(i, b);
00240     if (debugging)
00241     {
00242       report_one_constraint(i, b);
00243     }
00244     return true;
00245   }
00246   return false;
00247 }
00248 
00249 /*
00250 if (i>num_variables())
00251 {
00252   assert(phase >= 1);
00253   const int j = i - num_variables();
00254   printf("  constraining sum[%d] from %e to %c %d\n", j, sum_even_value(j), '>', b);
00255 }
00256 else
00257 
00258  
00259  if (i>num_variables())
00260  {
00261  assert(phase >= 1);
00262  const int j = i - num_variables();
00263  printf("  new sum[%d] value is %e %c %d\n", j, sum_even_value(j), '>', b);
00264  }
00265  else
00266  {
00267 */
00268 
00269 void IASolverEven::report_one_constraint(const int i, const int b)
00270 {
00271   // assumes x_solution[i] is current non-integer solution
00272   printf("  constraining x[%d] from %e (%e) to %c %d\n", i, x_solution[i], iaData->I[i], b > iaData->I[i] ? '>' : '<', b);
00273 }
00274 
00275 void IASolverEven::report_one_non_integer(const int i, const int b)
00276 {
00277   printf("  new x[%d] value is ", i );
00278   if (fabs( floor(x_solution[i]+0.5) - x_solution[i] ) < 1.0e-2)
00279     printf("%d", (int) floor(x_solution[i]+0.5));
00280   else
00281     printf("%e", x_solution[i]);
00282   printf(" (%e) %c %d\n", iaData->I[i], b > iaData->I[i] ? '>' : '<', b);
00283 }
00284 
00285 void IASolverEven::back_off(RoundingMap &rounding_map)
00286 {
00287   // undo half of the constraints
00288   // todo: 
00289   RoundingMap::iterator undo_start = rounding_map.begin(); 
00290   int num_to_keep = (int) rounding_map.size() / 2; // roundoff, rounding to zero is OK
00291   std::advance( undo_start, num_to_keep );
00292   for(RoundingMap::const_iterator i = undo_start; i != rounding_map.end(); ++i)
00293   {
00294     const int x_i = i->second.first;
00295     if (debugging)
00296     {
00297       printf(" backing off:");
00298       report_one_constraint(x_i, ipData.oldBound[x_i]);
00299     }
00300     ipData.varIntegerBound[x_i] =  ipData.oldBound[x_i]; 
00301   }
00302   rounding_map.erase(undo_start, rounding_map.end()); 
00303 }
00304 
00305 bool IASolverEven::find_many_non_integer(RoundingMap &rounding_map)
00306 {
00307   // to do: this heuristic is bad because it easily paints us into a corner 
00308   // to do: we should instead pick some constant fraction of the variable per *side* of a constraint, 
00309   // provided there is at least one other unconstrained variable on that side. And never constrain the last k-sides (k=2, or 3? or 4?) until the end, when we just constrain one at a time.
00310   
00311   // since we don't know how many are non-integer, put them all in the map, then discard the ones we won't constrain
00312   rounding_map.clear();
00313   for (int i = 0; i<iaData->num_variables(); ++i)
00314   {
00315     double obj_increase;
00316     int y_int;
00317     if (calculate_rounding_value(i, obj_increase, y_int))
00318     {
00319       // add to map
00320       // printf("inserting obj %e, i %d, x %e, y %d\n", obj_increase, i, x_solution[i], y_int);
00321       rounding_map.insert(std::make_pair(obj_increase,std::make_pair(i,y_int)));  
00322     }
00323   }
00324   if (rounding_map.empty())    
00325     return false;
00326   //const unsigned long num_to_round = 1;
00327   const unsigned long num_to_round = 1; // 1 + rounding_map.size() / 8; // todo: experiment with the 16 part
00328   if (debugging)
00329   {
00330     printf("Found %ld of %d non-integers, rounding %ld of them.",
00331            rounding_map.size(), iaData->num_variables(), num_to_round);
00332   }
00333   // remove everything after the "num_to_round" element
00334   RoundingMap::iterator r = rounding_map.begin();
00335   advance(r, num_to_round);
00336   rounding_map.erase(r, rounding_map.end()); 
00337   // todo: if erasing (or subsequent iterations) is slow due to tree rebalancing, then instead we should transfer the first num_to_round elements to a std::vector and return that instead, and discard the map
00338 
00339   if (debugging)
00340   {
00341     if (rounding_map.size() > 1)
00342       printf(" obj_increase range %e to %e\n",
00343              rounding_map.begin()->first, rounding_map.rbegin()->first);
00344     else if (rounding_map.size() == 1)
00345       printf( " obj_increase expected %e\n",
00346              rounding_map.begin()->first);
00347     else
00348       printf( " rounding nothing\n");
00349   }
00350   return true;
00351 }
00352 
00353 bool IASolverEven::constrain_many_non_integer(RoundingMap &rounding_map)
00354 {
00355   if (find_many_non_integer(rounding_map))
00356   {
00357     for(RoundingMap::const_iterator i = rounding_map.begin(); i != rounding_map.end(); ++i)
00358     {
00359       const int x_i = i->second.first;
00360       const int b = i->second.second;
00361       constrain_integer(x_i, b); // index, new integer bound
00362       if (debugging)
00363       {
00364         report_one_constraint(x_i, b);
00365       }
00366     }
00367     // todo: heuristic to decrease the number of constrained variables if the solver returns infeasible, or it decreased these variables beyond the new upper bounds . save the old bounds to restore them.
00368     return true;
00369   }
00370   return false; 
00371 }
00372 
00373 void IASolverEven::report_many_non_integer(RoundingMap &rounding_map)
00374 {
00375   for(RoundingMap::const_iterator i = rounding_map.begin(); i != rounding_map.end(); ++i)
00376   {
00377     report_one_non_integer(i->second.first, i->second.second);
00378   }
00379 }
00380 
00381 void IASolverEven::cleanup()
00382 {
00383   myianlp = NULL; // smart pointers should take care of deleting this.
00384 }
00385 
00386 bool IASolverEven::solve()
00387 {
00388   
00389   // solve the nlp to get a non-integral solution, which we hope is close to a good integer solution    
00390   // =======================================
00391   // from HS071 ipopt example
00392 
00393   // smaller p has the effect of valuing the fidelity of shorter curves over longer curves more
00394   // larger p approaches min max
00395   myianlp = new IAMINlp(iaData, &ipData, this);
00396   SmartPtr<TNLP> mynlp = myianlp; // Ipopt requires the use of smartptrs!
00397   
00398   SmartPtr<IpoptApplication> app = IpoptApplicationFactory();
00399   app->Options()->SetNumericValue("tol", 1e-7); // 2 seems close enough, could do less, say .1
00400   app->Options()->SetStringValue("mu_strategy", "adaptive");
00401   app->Options()->SetIntegerValue("print_level", 1);  // 0 to 12, most. Default is 5
00402   // uncomment next line to write the solution to an output file
00403   // app->Options()->SetStringValue("output_file", "IA.out");  
00404   // The following overwrites the default name (ipopt.opt) of the options file
00405   // app->Options()->SetStringValue("option_file_name", "IA.opt");
00406   
00407   // Intialize the IpoptApplication and process the options
00408   Ipopt::ApplicationReturnStatus status;
00409   status = app->Initialize();
00410   if (status != Ipopt::Solve_Succeeded) {
00411     printf("\n\n*** Error during initialization!\n");
00412     cleanup();
00413     return (int) status;
00414   }
00415   
00416   // Ask Ipopt to solve the problem
00417   status = app->OptimizeTNLP(mynlp); // the inherited IANlp
00418   
00419   if (status == Solve_Succeeded) {
00420     printf("\n\n*** The relaxed problem solved!\n");
00421   }
00422   else {
00423     printf("\n\n*** The relaxed problem FAILED!\n");
00424     cleanup();
00425     return false;
00426   }
00427   
00428   // make integer
00429   // Pick the variable whose rounding to the next further integer from its goal results in the lowest individual contribution (or increase in contribution from its current value) to the objective function.
00430   // Round it that way, setting a lower bound on its value (but leaving objective, etc., the same)
00431   // Re-solve the NLP. Some old vars might have their delta directions or values flip. That is OK.
00432   // Iterate.
00433 
00434   ipData.initialize(x_solution);
00435   int iter = 0;
00436   printf("\n\n*** Solving integer solution.\n"); 
00437 
00438     
00439   /* */
00440   // many at a time speedup, could skip and still get correct output
00441   if (debugging)
00442     printf("Many-at-a-time rounding\n");
00443   RoundingMap rounding_map;
00444   // phase = -1;
00445   while (constrain_many_non_integer(rounding_map))
00446   /* */
00447   {
00448     
00449     if (debugging)
00450       printf("MINLP many-at-once iter %d\n", iter);
00451     else
00452       printf("%d ", iter);
00453     
00454     // solve for many new constraints, but reduce those constraints if the problem becomes infeasible
00455     bool try_again = false;
00456     do 
00457     {
00458       try_again = false;
00459       status = app->OptimizeTNLP(mynlp); // double-check that resolving works, w/out initialize, etc 
00460       if (status != Ipopt::Solve_Succeeded && status != Ipopt::Solved_To_Acceptable_Level)
00461       {
00462         back_off(rounding_map);
00463         try_again = rounding_map.size();
00464       }
00465     }
00466     while (try_again);
00467     if (rounding_map.size() == 0) // no progress was made
00468       break;
00469       
00470     // to do: what to do if status != Solve_Succeeded ?
00471     // assert the new constraint was satisfied, one-at-a-time
00472     // assert( (b >= I[i] && x_solution[i] >= b) || (b < I[i] && x_solution[i] <= b)); 
00473     if (debugging)
00474     {
00475       // report_one_non_integer(i, b);
00476       report_many_non_integer(rounding_map);
00477     }
00478     ++iter;
00479     // todo: quit at an iteration limit? but a single variable could be increased several times...
00480   }
00481 
00482   
00483   // one at a time
00484   if (debugging)
00485     printf("\nOne-at-a-time rounding, including sum-even constraints.\n");
00486   int i; // non-integer variable to make integer
00487   int b; // value to make it integer
00488   iter = 0;
00489 
00490   while (constrain_one_non_integer(i,b))
00491   {
00492       
00493     if (debugging)
00494       printf("MINLP one-at-a-time iter %d\n", iter);
00495     else
00496       printf("%d ", iter);
00497       
00498     status = app->OptimizeTNLP(mynlp); // double-check that resolving works, w/out initialize, etc 
00499     // to do: what to do if status != Solve_Succeeded ?
00500     const double g = iaData->I[i];
00501     const double x = x_solution[i];
00502     assert((b >= g && x + 1.e-2 >= b) || (b < g && x <= b + 1.e-2)); // assert the new constraint was satisfied, one-at-a-time
00503     if ((status != Ipopt::Solve_Succeeded && status != Ipopt::Solved_To_Acceptable_Level) || 
00504         ! ((b >= g && x + 1.e-2 >= b) || (b < g && x <= b + 1.e-2)))
00505     {
00506       // see ipopt/CoinIpopt/Ipopt/src/Interfaces/IpReturnCodes_inc.h for the enum of possible values
00507       // status = Infeasible_Problem_Detected; 
00508       if (status == Infeasible_Problem_Detected )
00509         printf("IPOPT Infeasible_Problem_Detected\n");
00510       else
00511         printf("IPOPT bad solution return code %d\n", status);
00512       break;
00513     }
00514     if (debugging)
00515     {
00516       report_one_non_integer(i, b);
00517     }
00518     ++iter;
00519     // todo: quit at an iteration limit? but a single variable could be increased several times...
00520   }
00521   
00522   if (status != Ipopt::Solve_Succeeded)
00523     printf("***Integer Fail!***\n");
00524   else
00525   {
00526     if (debugging)
00527     {
00528       // print_solution();
00529     }
00530   }
00531   
00532   // As the SmartPtrs go out of scope, the reference count
00533   // will be decremented and the objects will automatically
00534   // be deleted. 
00535   cleanup();
00536   return status == Ipopt::Solve_Succeeded;
00537   
00538 }
00539 
00540 } // namespace MeshKit
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines