MeshKit  1.0
IASolverInt.cpp
Go to the documentation of this file.
00001 
00002 // IASolverInt.cpp
00003 // Interval Assignment for Meshkit
00004 //
00005 #include "IASolverInt.hpp"
00006 #include "IAData.hpp"
00007 #include "IPData.hpp"
00008 #include "IASolution.hpp"
00009 #include "IARoundingNlp.hpp"
00010 #include "IASolverRelaxed.hpp"
00011 #include "IASolverBend.hpp"
00012 // #include "IAMINlp.hpp"
00013 #include "IAIntCosNlp.hpp"
00014 #include "IAIntParabolaNlp.hpp"
00015 
00016 // #include "IAMilp.hpp" // glpk based solution too slow
00017 
00018 #include <stdio.h>
00019 #include <math.h>
00020 #include <limits.h>
00021 
00022 #include "IpIpoptApplication.hpp"
00023 
00024 namespace MeshKit 
00025 {
00026     
00027 IASolverInt::IASolverInt(const IAData * ia_data_ptr, IASolution *relaxed_solution_ptr, 
00028   const bool set_silent) 
00029   : IASolverToolInt(ia_data_ptr, relaxed_solution_ptr, true), 
00030   silent(false),
00031 //  silent(set_silent), 
00032   debugging(true)
00033 { 
00034   ip_data(new IPData);
00035   // initialize copies relaxed solution, then we can overwrite relaxed_solution_pointer
00036   ip_data()->initialize(relaxed_solution_ptr->x_solution); 
00037 }
00038 
00040 IASolverInt::~IASolverInt() 
00041 {
00042   delete ip_data();
00043 }
00044     
00045 bool IASolverInt::solve_wave_workhorse(IAIntWaveNlp *mynlp)
00046 {
00047   if (debugging)
00048   {
00049     printf("IASolverInt::solve_wave() - ");        
00050     printf("Attempting to enforce an integer and even solution to the relaxed NLP by adding constraints that repeat wave-like at each integer lattice point.\n");
00051   }
00052   
00053   // solver setup  
00054   Ipopt::SmartPtr<Ipopt::IpoptApplication> app = IpoptApplicationFactory();
00055 /* try leaving defaults
00056   // convergence parameters
00057   // see $IPOPTDIR/Ipopt/src/Interfaces/IpIpoptApplication.cpp
00058   // our real criteria are: all integer, constraints satisfied. How to test the "all_integer" part?
00059   app->Options()->SetNumericValue("tol", 1e-6); //"converged" if NLP error<this, default is 1e-7. Obj are scaled to be >1, so e-2 is plenty // was 1e-2
00060   app->Options()->SetNumericValue("max_cpu_time", sqrt( iaData->num_variables() ) ); // max time allowed in seconds
00061   app->Options()->SetIntegerValue("max_iter", 3 * (10 + iaData->num_variables() ) ); // max number of iterations
00062   // app->Options()->SetNumericValue("primal_inf_tol", 1e-2 ); 
00063   app->Options()->SetNumericValue("dual_inf_tol", 1e-2 ); // how close to infeasibility? // was 1e-2
00064   app->Options()->SetNumericValue("constr_viol_tol", 1e-2 ); // by how much can constraints be violated?
00065   app->Options()->SetNumericValue("compl_inf_tol", 1e-6 ); // max norm of complementary condition // was 1e-2
00066   
00067   // second criteria convergence parameters: quit if within this tol for many iterations
00068   // was  app->Options()->SetIntegerValue("acceptable_iter", 4 + sqrt( iaData->num_variables() ) ); //as "tol"
00069   app->Options()->SetNumericValue("acceptable_tol", 1e-6 ); //as "tol" was 1e-1
00070   
00071   app->Options()->SetStringValue("mu_strategy", "adaptive");
00072   // print level 0 to 12, Ipopt default is 5
00073   const int print_level = (silent) ? 0 : 1;  // simple info is 1, debug at other values
00074   app->Options()->SetIntegerValue("print_level", print_level);  
00075   // uncomment next line to write the solution to an output file
00076   // app->Options()->SetStringValue("output_file", "IA.out");  
00077   // The following overwrites the default name (ipopt.opt) of the options file
00078   // app->Options()->SetStringValue("option_file_name", "IA.opt");
00079   
00080   */
00081   
00082   // Intialize the IpoptApplication and process the options
00083   Ipopt::ApplicationReturnStatus status;
00084   status = app->Initialize();
00085   if (status != Ipopt::Solve_Succeeded) {
00086     if (!silent)
00087       printf("\n\n*** Error during initialization!\n");
00088     return (int) status;
00089   }
00090   
00091   bool try_again = true;
00092   int iter = 0;
00093   
00094   // print();
00095   bool solution_ok = false;
00096   
00097   do {
00098     if (debugging)
00099     {
00100       print();
00101       printf("%d IntWave iteration\n", iter );
00102       // build the hessian, evaluate it and f at the current solution?
00103     }
00104       
00105     // Ask Ipopt to solve the problem
00106     status = app->OptimizeTNLP(mynlp); // the inherited IANlp
00107     
00108     // see /CoinIpopt/build/include/coin/IpReturnCodes_inc.h
00109     /*
00110     Solve_Succeeded=0,
00111     Solved_To_Acceptable_Level=1,
00112     Infeasible_Problem_Detected=2,
00113     Search_Direction_Becomes_Too_Small=3,
00114     Diverging_Iterates=4,
00115     User_Requested_Stop=5,
00116     Feasible_Point_Found=6,
00117     
00118     Maximum_Iterations_Exceeded=-1,
00119     Restoration_Failed=-2,
00120     Error_In_Step_Computation=-3,
00121     Maximum_CpuTime_Exceeded=-4,
00122     Not_Enough_Degrees_Of_Freedom=-10,
00123     Invalid_Problem_Definition=-11,
00124     Invalid_Option=-12,
00125     Invalid_Number_Detected=-13,
00126     
00127     Unrecoverable_Exception=-100,
00128     NonIpopt_Exception_Thrown=-101,
00129     Insufficient_Memory=-102,
00130     Internal_Error=-199
00131      */
00132 
00133     bool solved_full = false;
00134     bool solved_partial = false;
00135     bool solver_failed = false;
00136     bool bad_problem = false;
00137 
00138     switch (status) {
00139       case Ipopt::Solve_Succeeded:
00140       case Ipopt::Solved_To_Acceptable_Level:
00141       case Ipopt::Feasible_Point_Found:
00142         solved_full = true;
00143         break;
00144       case Ipopt::Maximum_Iterations_Exceeded:
00145       case Ipopt::User_Requested_Stop:
00146       case Ipopt::Maximum_CpuTime_Exceeded:
00147         solved_partial = true;
00148         break;
00149       case Ipopt::Infeasible_Problem_Detected:
00150       case Ipopt::Not_Enough_Degrees_Of_Freedom:
00151       case Ipopt::Invalid_Problem_Definition:
00152       case Ipopt::Invalid_Option:
00153       case Ipopt::Invalid_Number_Detected:
00154         bad_problem = true;
00155         break;
00156       case Ipopt::Search_Direction_Becomes_Too_Small:
00157       case Ipopt::Restoration_Failed:
00158       case Ipopt::Diverging_Iterates:
00159       case Ipopt::Error_In_Step_Computation:
00160       case Ipopt::Unrecoverable_Exception:
00161       case Ipopt::NonIpopt_Exception_Thrown:
00162       case Ipopt::Insufficient_Memory:
00163       case Ipopt::Internal_Error:        
00164         solver_failed = true;
00165         break;
00166         
00167       default:
00168         break;
00169     }
00170   
00171     if (!silent)
00172     {
00173       if (solved_full) {
00174         printf("\n\n*** IntWave solved!\n");
00175       }
00176       else {
00177         printf("\n\n*** IntWave FAILED!\n");
00178       }
00179     }
00180     
00181     if (debugging)
00182     {
00183       printf("\nChecking solution.\n");
00184       bool integer_sat = solution_is_integer(true);
00185       bool even_sat = even_constraints( false, true);
00186       bool equal_sat = equal_constraints( false, true );
00187       printf("IntWave solution summary, %s, equal-constraints %s, even-constraints %s.\n", 
00188              integer_sat ? "integer" : "NON-INTEGER",
00189              equal_sat ? "satisfied" : "VIOLATED", 
00190              even_sat ? "satisfied" : "VIOLATED" );
00191       if (!integer_sat)
00192         printf("investigate integer neighborhood\n");
00193       if (!even_sat)
00194         printf("investigate even neighborhood\n");
00195       if (!equal_sat)
00196         printf("investigate equal neighborhood\n");
00197     }
00198     
00199     
00200     IASolution nlp_solution;
00201     nlp_solution.x_solution = ia_solution()->x_solution; // vector copy
00202     IASolverToolInt sti( ia_data(), &nlp_solution );
00203     sti.round_solution();
00204     if (debugging)
00205       printf("Checking rounded solution, ignoring even constraints.\n");
00206     if (sti.equal_constraints(false, debugging))
00207     {
00208       // also even constraints
00209       if (debugging)
00210         printf("Rounding worked.\n");
00211       
00212       // rounding was a valid integer solution
00213       ia_solution()->x_solution.swap( nlp_solution.x_solution );
00214       // ia_solution()->obj_value is no longer accurate, as it was for the non-rounded solution
00215       return true;
00216     }
00217     
00218     // todo: detect and act
00219     // may have converged to a locally optimal, but non-feasible solution
00220     // if so, try a new starting point
00221     
00222     // check solution feasibility, even when not debugging
00223     
00224     if ( solved_full || solved_partial )
00225     {
00226       bool integer_sat = solution_is_integer(false);
00227       bool even_sat = even_constraints( false, false);
00228       bool equal_sat = equal_constraints( false, false );
00229       if ( integer_sat && even_sat && equal_sat )
00230         return true;
00231     }
00232 
00233     // find out which vars were not integer, 
00234     // try moving to a farther starting point resolving
00235     
00236  
00237     try_again = false; 
00238   } while (try_again);
00239   
00240   
00241   // now 
00242   // As the SmartPtrs go out of scope, the reference count
00243   // will be decremented and the objects will automatically
00244   // be deleted.  
00245   return solution_ok;
00246   
00247 }
00248 
00249 
00250 bool IASolverInt::solve_round()
00251 {
00252   // set up and call the separate IARoundingNlp, which has a linear objective to get a natural integer solution 
00253   // the intuition is this will solve integrality for  most variables all at once
00254 
00255   if (debugging)
00256   {
00257     printf("IASolverInt::solve_bend_workhorse() - ");        
00258   }
00259 
00260   
00261   // solver setup  
00262   Ipopt::SmartPtr<Ipopt::IpoptApplication> app = IpoptApplicationFactory();
00263   
00264   // convergence parameters
00265   // see $IPOPTDIR/Ipopt/src/Interfaces/IpIpoptApplication.cpp
00266   // our real criteria are: all integer, constraints satisfied. How to test the "all_integer" part?
00267   app->Options()->SetNumericValue("tol", 1e-6); //"converged" if NLP error<this, default is 1e-7. Obj are scaled to be >1, so e-2 is plenty // was 1e-2
00268   app->Options()->SetNumericValue("max_cpu_time", sqrt( iaData->num_variables() ) ); // max time allowed in seconds
00269   app->Options()->SetIntegerValue("max_iter", 3 * iaData->num_variables() ); // max number of iterations
00270   // app->Options()->SetNumericValue("primal_inf_tol", 1e-2 ); 
00271   app->Options()->SetNumericValue("dual_inf_tol", 1e-6 ); // how close to infeasibility? // was 1e-2
00272   app->Options()->SetNumericValue("constr_viol_tol", 1e-6 ); // by how much can constraints be violated?
00273   app->Options()->SetNumericValue("compl_inf_tol", 1e-6 ); // max norm of complementary condition // was 1e-2
00274 
00275   // second criteria convergence parameters: quit if within this tol for many iterations
00276 // was  app->Options()->SetIntegerValue("acceptable_iter", 4 + sqrt( iaData->num_variables() ) ); //as "tol"
00277   app->Options()->SetNumericValue("acceptable_tol", 1e-6 ); //as "tol" was 1e-1
00278 
00279   app->Options()->SetStringValue("mu_strategy", "adaptive");
00280   // print level 0 to 12, Ipopt default is 5
00281   const int print_level = (silent) ? 0 : 1; 
00282   app->Options()->SetIntegerValue("print_level", print_level);  
00283   // uncomment next line to write the solution to an output file
00284   // app->Options()->SetStringValue("output_file", "IA.out");  
00285   // The following overwrites the default name (ipopt.opt) of the options file
00286   // app->Options()->SetStringValue("option_file_name", "IA.opt");
00287   
00288   // Intialize the IpoptApplication and process the options
00289   Ipopt::ApplicationReturnStatus status;
00290   status = app->Initialize();
00291   if (status != Ipopt::Solve_Succeeded) {
00292     if (!silent)
00293       printf("\n\n*** Error during initialization!\n");
00294     return (int) status;
00295   }
00296   
00297   
00298   Ipopt::TNLP *tnlp = NULL;
00299 
00300   IARoundingNlp *myianlp = new IARoundingNlp(iaData, ipData, iaSolution, silent);
00301   if (debugging) 
00302   {          
00303     printf("ROUNDING problem formulation\n");
00304     printf("Attempting to find a naturally-integer solution by linearizing the objective function.\n");
00305     printf("Variables are constrained within [floor,ceil] of relaxed solution.\n");
00306   }
00307   
00308   // problem setup
00309   // a couple of different models, simplest to more complex
00310   // IARoundingFarNlp *myianlp = new IARoundingFarNlp(iaData, ipData, this);
00311   // IARoundingFar3StepNlp *myianlp = new IARoundingFar3StepNlp(iaData, ipData, this); // haven't tested this. It compiles and runs but perhaps isn't correct
00312   // IAIntWaveNlp *myianlp = new IAIntWaveNlp(iaData, ipData, this); // haven't tested this. It compiles and runs but perhaps isn't correct
00313 
00314   tnlp = myianlp;
00315   Ipopt::SmartPtr<Ipopt::TNLP> mynlp = tnlp; // Ipopt requires the use of smartptrs!
00316 
00317   bool try_again = true;
00318   int iter = 0;
00319   do {
00320     printf("%d rounding iteration\n", iter );
00321     
00322     // Ask Ipopt to solve the problem
00323     status = app->OptimizeTNLP(mynlp); // the inherited IANlp
00324     
00325     if (!silent)
00326     {
00327       if (status == Ipopt::Solve_Succeeded) {
00328         printf("\n\n*** The problem solved!\n");
00329       }
00330       else {
00331         printf("\n\n*** The problem FAILED!\n");
00332       }
00333     }
00334     
00335     // The problem should have been feasible, but it is possible that it had no integer solution.
00336     // figure out which variables are still integer
00337     
00338     // check solution for integrality and constraint satified    
00339     if (debugging)
00340     {
00341       printf("\nChecking Natural (non-rounded) solution.\n");
00342       bool integer_sat = solution_is_integer(true);
00343       bool even_sat = even_constraints( false, true);
00344             bool equal_sat = equal_constraints( false, true );
00345       printf("Natural solution summary, %s, equal-constraints %s, even-constraints %s.\n", 
00346              integer_sat ? "integer" : "NON-INTEGER",
00347              equal_sat ? "satisfied" : "VIOLATED", 
00348              even_sat ? "satisfied" : "VIOLATED" );
00349       if (!integer_sat)
00350         printf("investigate this\n");
00351     }
00352     
00353     IASolution nlp_solution;
00354     nlp_solution.x_solution = ia_solution()->x_solution; // vector copy
00355     IASolverToolInt sti( ia_data(), &nlp_solution );
00356     sti.round_solution();
00357     if (debugging)
00358       printf("Checking rounded solution, ignoring even constraints.\n");
00359     if (sti.equal_constraints(false, debugging))
00360     {
00361       // also even constraints
00362       if (debugging)
00363         printf("Rounding worked.\n");
00364 
00365       // rounding was a valid integer solution
00366       ia_solution()->x_solution.swap( nlp_solution.x_solution );
00367       // ia_solution()->obj_value is no longer accurate, as it was for the non-rounded solution
00368       return true;
00369     }
00370 
00371     // find out which vars were not integer, 
00372     // try rounding their weights and resolving
00373     // bool int_sat = solution_is_integer();
00374     ++iter;
00375     try_again = iter < 4 + sqrt(iaData->num_variables());
00376     if (try_again)
00377     {
00378       if (debugging)
00379         printf("rounding failed, randomizing weights\n");
00380     
00381       myianlp->randomize_weights_of_non_int(); // try again? debug
00382     }
00383     else if (debugging)
00384       printf("giving up on rounding to non-integer solution\n");
00385 
00386     // try_again = false; // debug
00387   } while (try_again);
00388 
00389   
00390   // todo: update partially-integer solution, perhaps using ipData - figure out how we're going to use it, first, for what structure makes sense.
00391   
00392   // As the SmartPtrs go out of scope, the reference count
00393   // will be decremented and the objects will automatically
00394   // be deleted.  
00395   return status == Ipopt::Solve_Succeeded;
00396   
00397 }
00398 
00399 void IASolverInt::cleanup()
00400 {
00401   ;
00402 }
00403 
00404 bool IASolverInt::solve_wave(const SolverType solver_type)
00405 {
00406   IAIntWaveNlp *myianlp = NULL;
00407   if (solver_type == COS) 
00408   {
00409     if (debugging) printf("Cosine wave.\n");
00410     myianlp= new IAIntCosNlp(iaData, ipData, iaSolution, silent);
00411   }
00412   else if (solver_type == PARABOLA)
00413   {
00414     if (debugging) printf("Parabola wave.\n");
00415     myianlp = new IAIntParabolaNlp(iaData, ipData, iaSolution, silent);
00416   }
00417   else
00418   {
00419     if (debugging) printf("Invalid wave type.\n");
00420     return false;
00421   }
00422     
00423   Ipopt::SmartPtr<Ipopt::TNLP> mynlp = myianlp; // Ipopt requires the use of smartptrs!
00424   return solve_wave_workhorse( myianlp );
00425 }
00426 
00427 bool IASolverInt::solve()
00428 {
00429   
00430   SolverType solver_type = BEND; // BEND;
00431   
00432   // debug, try solve_intwave instead 
00433   // longer term, use intwave as a backup when the faster and simpler milp doesn't work.
00434   // unfortunately, it appears to find local minima that are far from optimal, even when starting in a well
00435   
00436   bool solved = false;
00437   switch (solver_type) {
00438     case COS:
00439     case PARABOLA:
00440       solved = solve_wave(solver_type);
00441       break;
00442     case ROUNDING:
00443       solved = solve_round();
00444       break;
00445     case BEND:
00446       {
00447         IASolverBend sb(iaData, iaSolution, silent);
00448         solved = sb.solve();
00449       }
00450       break;
00451     default:
00452       solved = false;
00453       break;
00454   }
00455 
00456   bool success = solution_is_integer();
00457   // also check constraints and evenality
00458   if (success)
00459   {
00460     if (!silent)
00461       printf("IASolverInt produced integer solution\n");
00462   }
00463   else
00464   {
00465     // todo: rather than applying the rounding heuristic, implement a form of IARoundingNlp with larger variable bounds, but still with a natural integer solution, by extending x to also depend on a delta_plus and delta_minus extending x beyond xl and xl+1, i.e. x = xl (const) + xh (0-1 var) + delta_plus - delta_minus. With linear objective with weight for xh as before, but weight for delta_plus to be f( xl + 2 ) - f (xl + 1), delta_minus f( xl - 2) - f(xl -1)
00466     return false; // debug;
00467     // solve_rounding_heuristic();
00468     // success = solution_is_integer();
00469   }
00470   return success;
00471 }
00472 
00473 } // namespace MeshKit
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines