MeshKit  1.0
IASolverBend.cpp
Go to the documentation of this file.
00001 // IASolverBend.cpp
00002 // Interval Assignment for Meshkit
00003 //
00004 #include "meshkit/IASolverBend.hpp"
00005 #include "meshkit/IAData.hpp"
00006 #include "meshkit/IPData.hpp"
00007 #include "meshkit/IPBend.hpp"
00008 #include "meshkit/IASolution.hpp"
00009 #include "meshkit/IABendNlp.hpp"
00010 
00011 #include <stdio.h>
00012 #include <math.h>
00013 #include <limits.h>
00014 
00015 #include "IpIpoptApplication.hpp"
00016 
00017 namespace MeshKit 
00018 {
00019     
00020 IASolverBend::IASolverBend(const IAData * ia_data_ptr, IASolution *relaxed_solution_ptr, 
00021   const bool set_silent) 
00022   : IASolverToolInt(ia_data_ptr, relaxed_solution_ptr, true), evenConstraintsActive(false),
00023     silent(set_silent), debugging(true)
00024 //    silent(set_silent), debugging(false)
00025 { 
00026   ip_data(new IPData);
00027   // initialize copies relaxed solution, then we can overwrite relaxed_solution_pointer with our integer solution 
00028   ip_data()->initialize(relaxed_solution_ptr->x_solution); 
00029 }
00030 
00032 IASolverBend::~IASolverBend() 
00033 {
00034   delete ip_data();
00035 }
00036     
00037 bool IASolverBend::solve_nlp() // IABendNlp *mynlp
00038 {
00039   if (debugging)
00040   {
00041     printf("IASolverBend::solve_nlp() == ");        
00042     printf("BEND problem formulation\n");
00043     printf("Attempting to find a naturally-integer solution by linearizing and bending the objective function at integer values.\n");
00044     printf("x = sum of positive and negative deltas around the floor of the relaxed solution. Deltas are within [0,1]. Deltas are dynamically added. Objective is linear function of weighted deltas, randomized to break ties.\n");
00045   }
00046   
00047   // solver setup  
00048   Ipopt::SmartPtr<Ipopt::IpoptApplication> app = IpoptApplicationFactory();
00049 
00050   //  jac_d_constant
00051   if (evenConstraintsActive && !ia_data()->sumEvenConstraints.empty() )
00052   {
00053     app->Options()->SetStringValue("jac_d_constant", "no"); // default
00054   }
00055   else
00056   {
00057     app->Options()->SetStringValue("jac_d_constant", "yes");
00058   }
00059   // even with the even-constraints, the hessian is constant
00060   app->Options()->SetStringValue("hessian_constant", "yes"); // no by default
00061 
00062 /* try leaving defaults
00063   // convergence parameters
00064   // see $IPOPTDIR/Ipopt/src/Interfaces/IpIpoptApplication.cpp
00065   // our real criteria are: all integer, constraints satisfied. How to test the "all_integer" part?
00066   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
00067   app->Options()->SetNumericValue("max_cpu_time", sqrt( iaData->num_variables() ) ); // max time allowed in seconds
00068   app->Options()->SetIntegerValue("max_iter", 3 * (10 + iaData->num_variables() ) ); // max number of iterations
00069   // app->Options()->SetNumericValue("primal_inf_tol", 1e-2 ); 
00070   app->Options()->SetNumericValue("dual_inf_tol", 1e-2 ); // how close to infeasibility? // was 1e-2
00071   app->Options()->SetNumericValue("constr_viol_tol", 1e-2 ); // by how much can constraints be violated?
00072   app->Options()->SetNumericValue("compl_inf_tol", 1e-6 ); // max norm of complementary condition // was 1e-2
00073   
00074   // second criteria convergence parameters: quit if within this tol for many iterations
00075   // was  app->Options()->SetIntegerValue("acceptable_iter", 4 + sqrt( iaData->num_variables() ) ); //as "tol"
00076   app->Options()->SetNumericValue("acceptable_tol", 1e-6 ); //as "tol" was 1e-1
00077   
00078   app->Options()->SetStringValue("mu_strategy", "adaptive");
00079   // print level 0 to 12, Ipopt default is 5
00080   const int print_level = (silent) ? 0 : 1;  // simple info is 1, debug at other values
00081   app->Options()->SetIntegerValue("print_level", print_level);  
00082   // uncomment next line to write the solution to an output file
00083   // app->Options()->SetStringValue("output_file", "IA.out");  
00084   // The following overwrites the default name (ipopt.opt) of the options file
00085   // app->Options()->SetStringValue("option_file_name", "IA.opt");
00086   
00087   */
00088   const int print_level = 2;  // simple info is 1, debug at other values
00089   app->Options()->SetIntegerValue("print_level", print_level);  
00090   
00091   // Intialize the IpoptApplication and process the options
00092   Ipopt::ApplicationReturnStatus status;
00093   status = app->Initialize();
00094   if (status != Ipopt::Solve_Succeeded) {
00095     if (!silent)
00096       printf("\n\n*** Error during initialization!\n");
00097     return (int) status;
00098   }
00099   
00100   bool try_again = false; // only try again if we didn't converge and want to spend more time
00101   int iter = 0;
00102   
00103   // print();
00104   bool solution_ok = false;
00105   
00106   myianlp->even_constraints_active( evenConstraintsActive );
00107   
00108   do {
00109     if (debugging)
00110     {
00111       print();
00112       printf("%d Bend iteration\n", iter );
00113       // build the hessian, evaluate it and f at the current solution?
00114     }
00115       
00116     // Ask Ipopt to solve the problem
00117     status = app->OptimizeTNLP(myianlp); // the inherited IANlp
00118     
00119     // see /CoinIpopt/build/include/coin/IpReturnCodes_inc.h
00120     /*
00121     Solve_Succeeded=0,
00122     Solved_To_Acceptable_Level=1,
00123     Infeasible_Problem_Detected=2,
00124     Search_Direction_Becomes_Too_Small=3,
00125     Diverging_Iterates=4,
00126     User_Requested_Stop=5,
00127     Feasible_Point_Found=6,
00128     
00129     Maximum_Iterations_Exceeded=-1,
00130     Restoration_Failed=-2,
00131     Error_In_Step_Computation=-3,
00132     Maximum_CpuTime_Exceeded=-4,
00133     Not_Enough_Degrees_Of_Freedom=-10,
00134     Invalid_Problem_Definition=-11,
00135     Invalid_Option=-12,
00136     Invalid_Number_Detected=-13,
00137     
00138     Unrecoverable_Exception=-100,
00139     NonIpopt_Exception_Thrown=-101,
00140     Insufficient_Memory=-102,
00141     Internal_Error=-199
00142      */
00143 
00144     bool solved_full = false;
00145     bool solved_partial = false;
00146     bool solved_failure = false;
00147     bool problem_bad = false;
00148     bool problem_unbounded = false;
00149     (void) solved_failure;
00150     (void) problem_bad;
00151     (void) problem_unbounded;
00152 
00153     switch (status) {
00154       case Ipopt::Solve_Succeeded:
00155       case Ipopt::Solved_To_Acceptable_Level:
00156       case Ipopt::Feasible_Point_Found:
00157         solved_full = true;
00158         break;
00159       case Ipopt::Maximum_Iterations_Exceeded:
00160       case Ipopt::User_Requested_Stop:
00161       case Ipopt::Maximum_CpuTime_Exceeded:
00162       case Ipopt::Search_Direction_Becomes_Too_Small:
00163         solved_partial = true;
00164         break;
00165       case Ipopt::Infeasible_Problem_Detected:
00166       case Ipopt::Not_Enough_Degrees_Of_Freedom:
00167       case Ipopt::Invalid_Problem_Definition:
00168       case Ipopt::Invalid_Option:
00169       case Ipopt::Invalid_Number_Detected:
00170         problem_bad = true;
00171         break;
00172       case Ipopt::Diverging_Iterates:
00173         problem_unbounded = true;
00174         solved_partial = true;
00175         break;
00176       case Ipopt::Restoration_Failed:
00177       case Ipopt::Error_In_Step_Computation:
00178       case Ipopt::Unrecoverable_Exception:
00179       case Ipopt::NonIpopt_Exception_Thrown:
00180       case Ipopt::Insufficient_Memory:
00181       case Ipopt::Internal_Error:        
00182         solved_failure = true;
00183         break;
00184         
00185       default:
00186         break;
00187     }
00188   
00189     if (!silent)
00190     {
00191       if (solved_full) {
00192         printf("\n\n*** BendNlp solved!\n");
00193       }
00194       else if (solved_partial) {
00195         printf("\n\n*** BendNlp partial success!\n");
00196       }
00197       else {
00198         printf("\n\n*** BendNlp FAILED!\n");
00199       }
00200     }
00201     
00202     if (debugging)
00203     {
00204       printf("\nChecking solution.\n");
00205       bool integer_sat = solution_is_integer(true);
00206       bool even_sat = even_constraints( false, true);
00207       bool equal_sat = equal_constraints( false, true );
00208       printf("Bend solution summary, %s, equal-constraints %s, even-constraints %s.\n", 
00209              integer_sat ? "integer" : "NON-INTEGER",
00210              equal_sat ? "satisfied" : "VIOLATED", 
00211              even_sat ? "satisfied" : "VIOLATED" );
00212     }
00213     try_again = false; 
00214   
00215     if ( solved_full || solved_partial )
00216     {
00217       return true;
00218     }
00219     else
00220     {
00221       // todo: tweak the problem and try again
00222       return false;
00223     }
00224     
00225   } while (try_again);
00226   
00227   return solution_ok;
00228   
00229 }
00230 
00231 bool IASolverBend::round_solution()
00232 {
00233   IASolution nlp_solution;
00234   nlp_solution.x_solution = ia_solution()->x_solution; // vector copy
00235   IASolverToolInt sti( ia_data(), &nlp_solution );
00236   sti.round_solution();
00237   if (debugging)
00238     printf("Checking rounded bend solution.\n");
00239   if (sti.equal_constraints(false, debugging) && sti.even_constraints(false, debugging) )
00240   {
00241   if (debugging)
00242     printf("Rounding worked.\n");
00243     
00244     // rounding was a valid integer solution
00245     ia_solution()->x_solution.swap( nlp_solution.x_solution );
00246     // ia_solution()->obj_value is no longer accurate, as it was for the non-rounded solution
00247     return true;
00248   }
00249   return false;
00250 }
00251 
00252 void IASolverBend::cleanup()
00253 {
00254   // the solution includes the delta values.
00255   // remove those
00256   iaSolution->x_solution.resize( (unsigned) iaData->num_variables());
00257 }
00258 
00259 
00260 // the objective function we should use for deriving the weights
00261 double IASolverBend::f_x_value( double I_i, double x_i ) const
00262 {
00263   return myianlp->f_x_value(I_i, x_i);
00264 }
00265 
00266   
00267 double IASolverBend::fpow(double f) const
00268 {
00269   return f*f*f; // f^3
00270   // todo, experiment with squaring, or less
00271 }
00272   
00273  
00274   /* Idea: a form of IARoundingNlp with larger variable bounds, but still with a natural integer solution.
00275    x in [1..inf]
00276    xr = x optimal relaxed solution with objective function fnlp, see IANlp.xpp 
00277    f is piecewise linear, with corners at integer values. f slopes are unique (we hope)
00278    Slope definitions
00279    for x between xl = floor xr and xh = ceil xr, we use the difference in fnlp between xl and xh
00280    
00281    case A. xr > ceil g, g is goal I[i]
00282    for x above xh, 
00283    let h+ be fnlp ( xh+1 ) - fnlp ( xh )
00284    let kp be the number of intervals x is above xh 
00285    then slope = floor kp * h+
00286    for x below xl, h- = sqrt(11) / 5 h+, and slope = floor km * h-
00287    all this is weighted by some unique weight
00288    
00289    case B. xr in [ floor g, ceil g] 
00290    h+ = fnlp ( xh+1 ) - fnlp ( xh )
00291    h- = fnlp ( xl-1 ) - fnlp ( xl ), not used if xl == 1
00292    
00293    case C. xr < floor g
00294    h- = fnlp ( xl-1 ) - fnlp ( xl )
00295    h+ = sqrt(10)/5 h-
00296    If g < 2, then h- is unused, and h+ = as in case B
00297    
00298    // representation:
00299    h0 is weights 0..n-1
00300    h+ is weights n..2n-1
00301    h- is weights 2n..
00302    
00303    */
00304 
00305   
00306 /* try parabolic constraints instead
00307  
00308 void IASolverBend::add_bend_sum_weights(unsigned int i, const double factor)
00309 {
00310   // scaling issue, 
00311   // in order to overcome the slopes of the integers, the minimum slope should be larger
00312   // than the largest active slope, where active means the weight of the bend delta at the current solution
00313   
00314   IPBend &bend = bendData.bendVec[i];
00315   // designed so xl is the "ideal"
00316   // const double xl = bend.xl; 
00317   
00318   const double s = even_value(i);
00319   const double g = s/2.;
00320 
00321   // pluses
00322   for (int j = 0; j < bend.numDeltaPlus; ++j)
00323   {
00324     double w = bendData.maxActiveVarWeight;
00325     // if current sum is between xl and xl + 1, underweight the first delta
00326     if ( (j == 0) && (g > bend.xl) )
00327     {
00328       double f = ceil(g) - g;
00329       assert( f >= 0.499 ); // xl was chosen to be the closest integer to s/2
00330       w *= f;
00331     }
00332     w *= factor * 2.2 * sqrt( 1.1 * j + 1.1 );    
00333     weights.push_back(w);
00334   }
00335   
00336   // minuses
00337   for (int j = 0; j < bend.numDeltaMinus; ++j)
00338   {
00339     double w = bendData.maxActiveVarWeight;
00340     if ( (j == 0) && (g < bend.xl) )
00341     {
00342       double f = g - floor(g);
00343       assert( f >= 0.499 ); // xl was chosen to be the closest integer to s/2
00344       w *= f;
00345     }
00346     w *= factor * 2.1 * sqrt( 1.02 * j + 1.02 );    
00347     weights.push_back(w);
00348   }
00349 }
00350   
00351  */
00352   
00353   /*
00354   deltapluses
00355   // now modify weights based on tilts
00356   for (std::vector<IPBend::IPTilt>::iterator t = bend.plusTilts.begin();
00357        t != bend.plusTilts.end(); ++t)
00358   {
00359     double tbig = t->first;
00360     double tlit = tbig - 1;
00361     if (tlit < g)
00362       tlit = g;
00363       const double ftlit = f_x_value(g, tlit); 
00364       const double ftbig = f_x_value(g, tbig);
00365       double slope = fpow(ftbig) - fpow(ftlit); 
00366       slope *= t->second;
00367       
00368       if (xbig >= tbig) //if +1
00369       {
00370         assert(w>0.);
00371         assert(slope>0.);
00372         w += fabs(slope);
00373       }
00374   }
00375   for (std::vector<IPBend::IPTilt>::iterator t = bend.minusTilts.begin();
00376        t != bend.minusTilts.end(); ++t)
00377   {
00378     double tbig = t->first;
00379     double tlit = tbig + 1;
00380     if (tlit > g)
00381       tlit = g;
00382       const double ftlit = f_x_value(g, tlit); 
00383       const double ftbig = f_x_value(g, tbig);
00384       double slope = fpow(ftbig) - fpow(ftlit); // always positive
00385       slope *= t->second;
00386       
00387       if (xlit <= tbig)
00388       {
00389         assert(w<0.);
00390         assert(slope>0.);
00391         w -= fabs(slope);
00392       }
00393   }
00394    deltaminuses
00395   // done with tilts
00396    // now modify weights based on tilts
00397    // this is slow and could be sped up by saving prior calculations
00398    for (std::vector<IPBend::IPTilt>::iterator t = bend.plusTilts.begin();
00399    t != bend.plusTilts.end(); ++t)
00400    {
00401    double tbig = t->first;
00402    double tlit = tbig - 1;
00403    if (tlit < g)
00404    tlit = g;
00405    const double ftlit = f_x_value(g, tlit); 
00406    const double ftbig = f_x_value(g, tbig);
00407    double slope = fpow(ftbig) - fpow(ftlit); 
00408    slope *= t->second;
00409    
00410    if (xbig >= tbig)
00411    {
00412    assert(w<0.);
00413    assert(slope>0.);
00414    w -= fabs(slope);
00415    }
00416    }
00417    for (std::vector<IPBend::IPTilt>::iterator t = bend.minusTilts.begin();
00418    t != bend.minusTilts.end(); ++t)
00419    {
00420    double tbig = t->first;
00421    double tlit = tbig + 1;
00422    if (tlit > g)
00423    tlit = g;
00424    const double ftlit = f_x_value(g, tlit); 
00425    const double ftbig = f_x_value(g, tbig);
00426    double slope = fpow(ftbig) - fpow(ftlit); // always positive
00427    slope *= t->second;
00428    
00429    if (xlit <= tbig)
00430    {
00431    assert(w>0.);
00432    assert(slope>0.);
00433    w += fabs(slope);
00434    }
00435    }
00436    // done with tilts
00437    
00438 
00439   */
00440  
00441 void IASolverBend::merge_tilts(IPBend::TiltVec &tilts)
00442 {
00443   if (tilts.size() < 2)
00444     return;
00445   
00446   std::sort( tilts.begin(), tilts.end() ); // sorts by .first
00447 
00448   // copy into a new vec
00449   IPBend::TiltVec new_tilts;
00450   new_tilts.reserve(tilts.size());
00451   
00452   // multiply tilts for the same index together, increasing them faster than additively
00453   for (unsigned int t = 0; t < tilts.size()-1; )
00454   {
00455     unsigned int u = t + 1;
00456     while (u < tilts.size() && tilts[t].first == tilts[u].first)
00457     {
00458       tilts[t].second *= tilts[u].second;
00459       u++;
00460     }
00461     t = u;
00462   }
00463   tilts.swap(new_tilts);  
00464 }
00465   
00466 void IASolverBend::tilt_weight(const IPBend::TiltVec &tilts, const int tilt_direction, const double g, const double xlit, const double xbig, const int delta_direction, double &w)
00467 {
00468   // O( num_deltas * num_tilts)
00469   // this could be sped up some by saving prior calculations or
00470   // by using the sorted order of the tilts to skip some deltas.
00471   for (IPBend::TiltVec::const_iterator t = tilts.begin();
00472        t != tilts.end(); ++t)
00473   {
00474     if (0 && debugging)
00475     {
00476       printf("  tilt (%d) %f at %d\n", tilt_direction, t->second, t->first);
00477     }
00478     double tbig = t->first;                   
00479     double tlit = tbig - tilt_direction;
00480     // if tilt_direction is negative, then tlit index is larger than tbig, 
00481     // but always tlit has smaller f value than tbig
00482     double slope = (tilt_direction > 0) ? raw_weight(g, tlit, tbig) : -raw_weight(g, tbig, tlit);
00483     assert(slope > 0. ); // always defined this way because of tilt_direction
00484     assert(t->second > 0.);
00485     
00486     slope *= t->second;
00487     
00488     // impose an absolute minimum slope of 0.1, to get out of flat regions around g
00489     
00490     if (tilt_direction > 0 && (xbig >= tbig))
00491     {
00492       if (delta_direction > 0)
00493       {
00494         assert(w>0.);
00495         w += fabs(slope);
00496       }
00497       else 
00498       {
00499         assert(w<0.);
00500         w -= fabs(slope); 
00501       }
00502     }
00503     else if (tilt_direction < 0 && (xlit <= tbig))
00504     {
00505       if (delta_direction > 0)
00506       {
00507         assert(w<0.);
00508         w -= fabs(slope);
00509       }          
00510       else {
00511         assert(w>0.);
00512         w += fabs(slope);
00513       }
00514     }
00515   }
00516 }
00517   
00518 double IASolverBend::raw_weight( const double g, double xlit, double xbig)
00519 {
00520   assert(xlit < xbig);
00521   // special handling when crossing goal to avoid zero slope
00522   if ( xlit < g && xbig > g )
00523   {
00524     if ( g - xlit < xbig - g) 
00525       xlit = g;
00526     else
00527       xbig = g;
00528   }
00529   assert( xbig >= 1.);
00530   assert( xlit >= 1.);
00531 
00532   const double flit = f_x_value(g, xlit); 
00533   const double fbig = f_x_value(g, xbig);
00534   const double w = fpow(fbig) - fpow(flit);
00535   return w;
00536 }
00537   
00538 void IASolverBend::add_bend_weights(unsigned int i)
00539 {
00540   // shorthands
00541   IPBend &bend = bendData.bendVec[i];
00542   const double xl = bend.xl;
00543   const double g = iaData->I[i]; // goal
00544   
00545   // current x solution for finding active weight
00546   const double x = iaSolution->x_solution[i]; 
00547 
00548   // pluses
00549   for (int j = 0; j < bend.numDeltaPlus; ++j)
00550   {
00551     double xlit = xl + j;
00552     double xbig = xlit + 1.;
00553     double w = raw_weight( g, xlit, xbig );
00554         
00555     if (0 && debugging)
00556     {
00557       printf("x[%u], g %f, xl %f, dplus[%u] raw_w %f\n", i, g, xl, j, w);
00558     }
00559     // now modify weights based on tilts
00560     tilt_weight(bend.plusTilts, 1, g, xlit, xbig, 1, w);
00561     tilt_weight(bend.minusTilts, -1, g, xlit, xbig, 1, w);
00562     if (0 && debugging)
00563     {
00564       if ( bend.plusTilts.size() || bend.minusTilts.size() )
00565         printf("              tilted_w %f\n", w);
00566     }
00567 
00568     weights.push_back(w);
00569 
00570     // active? or nearly active
00571     if (x <= xbig + 1. && x >= xlit)
00572     {
00573       // biggest?
00574       if (fabs(w) > bendData.maxActiveVarWeight)
00575         bendData.maxActiveVarWeight = fabs(w);
00576     }
00577   }
00578   
00579   // minuses
00580   for (int j = 0; j < bend.numDeltaMinus; ++j)
00581   {
00582     double xbig = xl - j;   
00583     double xlit = xbig - 1.;
00584     assert(xlit >= 1.); // if this fails, then the numDeltaMinus is too large
00585     
00586     double w = - raw_weight( g, xlit, xbig );
00587     
00588     if (0 && debugging)
00589     {
00590       printf("x[%u], g %f, xl %f, dminus[%u] raw_w %f\n", i, g, xl, j, w);
00591     }
00592     // now modify weights based on tilts
00593     tilt_weight(bend.plusTilts, 1, g, xlit, xbig, -1, w);
00594     tilt_weight(bend.minusTilts, -1, g, xlit, xbig, -1, w);    
00595     if (0 && debugging)
00596     {
00597       if ( bend.plusTilts.size() || bend.minusTilts.size() )
00598         printf("             tilted_w %f\n", w);
00599     }
00600     
00601     weights.push_back(w);
00602 
00603     // active? or nearly active
00604     if (x <= xbig && x >= xlit - 1.)
00605     {
00606       // biggest?
00607       if (fabs(w) > bendData.maxActiveVarWeight)
00608         bendData.maxActiveVarWeight = fabs(w);
00609     }
00610   }
00611 
00612 }
00613   
00614 
00615 void IASolverBend::initialize_ip_bends()
00616 {
00617   // problem layout:
00618   // vars:
00619   //   x variables that are supposed to be integer
00620   //   sums that are supposed to be even
00621   //   for i = 0 .. iaData->num_variables()
00622   //     x[i] delta pluses
00623   //     x[i] delta minuses
00624   // 
00625   // constraints
00626   //    base-problem
00627   //      equal
00628   //      even
00629   //    x=deltas, x = xl(const) + deltasplus - deltasminus, <==> x - deltasplus + deltasminus = -xl
00630   //    s=deltas
00631   //
00632   bendData.numSumVars = (int) iaData->sumEvenConstraints.size();
00633   bendData.sumVarStart = iaData->num_variables();
00634   
00635   // loop over the vars, 
00636   // finding upper and lower rounding values
00637   // assign weights
00638   bendData.bendVec.resize( iaData->num_variables() ); // + bendData.numSumVars if we want to do those using bends rather than waves
00639   weights.reserve(bendData.bendVec.size()*4);
00640   
00641   bendData.maxActiveVarWeight = 0.;
00642   
00643   int d_start = iaData->num_variables();
00644   for (int i = 0; i < iaData->num_variables(); ++i)
00645   {
00646     IPBend &bend = bendData.bendVec[i]; // shorthand
00647     bend.deltaIStart = d_start;
00648     
00649     double x = ipData->relaxedSolution[i];
00650     // fix any out-of-bounds roundoff issue
00651     if ( x < 1. ) 
00652       x = 1.;
00653     double xl = bend.xl = floor(x);
00654     assert(xl >= 1.);
00655     // to do, experimenting with starting with +2, -1 
00656 //    if ( x - floor(x) > 0.5 )
00657     {
00658       bend.numDeltaPlus = 2;
00659       bend.numDeltaMinus = 1; 
00660     }
00661     // just one bend, two deltas, is a bad idea, because it can lead to an unbounded objective funtion
00662 //    else
00663 //    {
00664 //      bend.numDeltaPlus = 1;
00665 //      bend.numDeltaMinus = 1;      
00666 //    }
00667 
00668     /*
00669     // debug, test negative deltas by starting with an xl that's too high
00670     x += 2.;
00671     xl = bend.xl = floor(x);
00672     bend.numDeltaPlus = 1;
00673     bend.numDeltaMinus = 1;      
00674     */
00675 
00676     // avoid minus deltas that allow the x solution to be less than 1
00677     if (bend.numDeltaMinus > xl - 1. )
00678       bend.numDeltaMinus = xl - 1;
00679     
00680     // always have at least one deltaMinus and one deltaPlus, so that the full range of x can be represented
00681     assert( ( xl == 1 ) || (bend.numDeltaMinus >= 1) );
00682     assert( bend.numDeltaPlus >= 1 );
00683     
00684     d_start += bend.num_deltas();
00685 
00686     add_bend_weights( i );
00687   }
00688   
00689   /* handle the sum-evens using parabolic constraints intead.
00690 
00691    // initialize all the sum-evens to have no bends - try to enforce those after iteration 1
00692   // as the initial int rounding could change the sum values a lot,
00693   // and unlike the x's there is no intrinsic goal for them.
00694 
00695   for (int i = 0; i < bendData.numSumVars; ++i)
00696   {
00697     IPBend &bend = bendData.bendVec[i + bendData.sumVarStart]; // shorthand
00698     // get current sum-even-value
00699     double s = even_value(i);
00700     bend.xl = floor( s / 2. + 0.5 ); // closest integer to s/2
00701     bend.numDeltaMinus = 1;
00702     bend.numDeltaPlus = 1;
00703     bend.deltaIStart = d_start;
00704 
00705     add_bend_sum_weights( i, 0. );
00706 
00707     d_start += bend.num_deltas();
00708   }
00709    */
00710   
00711   // gather the weights in a vector 
00712   
00713   // uniquify the weights
00714   weights.uniquify(1., 1.e6); // 1.e4 ??
00715   
00716   
00717   // also do that later in a lazy fasion after we see which vars are not integer
00718   
00719   
00720 }
00721   
00722  
00723 bool IASolverBend::update_ip_bends()
00724 {
00725 
00726   // ===============
00727   // check that the structure of the deltas is as expected
00728   // 1, 1, 1, 0.5, 0, 0, 0   or 
00729   // 1, 1, 1,   1, 1, 1, 3
00730   if (debugging)
00731   {
00732     for (int i = 0; i < iaData->num_variables(); ++i)
00733     {
00734       IPBend &bend = bendData.bendVec[i]; // shorthand
00735       
00736       bool plus_deltas = false;
00737       (void) plus_deltas;
00738       double xprior = 1.;
00739       for (int j = 0; j < bend.numDeltaPlus-1; ++j)
00740       {
00741         const int di = bend.deltaIStart + j;
00742         const double xp = iaSolution->x_solution[ di ];
00743         assert(xp <= xprior + 0.1 );
00744         xprior = xp;
00745         if (xprior < 0.9)
00746           xprior = 0.;
00747         if (xp > 0.1) 
00748           plus_deltas = true;
00749       }
00750       const int diplast = bend.deltaIStart + bend.numDeltaPlus-1;
00751       const double xp = iaSolution->x_solution[ diplast ];
00752       assert( xprior > 0.9 || xp < 1. );
00753       if ( xp > 0.1 )
00754         plus_deltas = true;
00755       
00756       bool minus_deltas = false;
00757       (void) minus_deltas;
00758       xprior = 1.;
00759       for (int j = 0; j < bend.numDeltaMinus-1; ++j)
00760       {
00761         const int di = bend.deltaIStart + bend.numDeltaPlus + j;
00762         const double xm = iaSolution->x_solution[ di ];
00763         assert( xm <= xprior + 0.1);
00764         xprior = xm;
00765         if (xprior < 0.9)
00766           xprior = 0.;
00767         if (xm > 0.1)
00768           minus_deltas = true;
00769       }
00770       const int dimlast = bend.deltaIStart + bend.numDeltaPlus + bend.numDeltaMinus-1;
00771       const double xm = iaSolution->x_solution[ dimlast ];
00772       assert( xprior > 0.9 || xm < 1. );
00773       if (xm > 0.1)
00774         minus_deltas = true;
00775       assert( ! (plus_deltas && minus_deltas) );
00776     }
00777   }
00778   // ===============
00779   
00780   
00781   // real algorithm
00782 
00783   // ====== new bends
00784   bool new_bend = false; // was at least one new bend added?
00785   bool randomized = false;
00786   
00787   int d_start = iaData->num_variables();
00788   for (int i = 0; i < iaData->num_variables(); ++i)
00789   {
00790     IPBend &bend = bendData.bendVec[i]; // shorthand
00791 
00792     // delta > 1? 
00793     // add more deltas
00794     const int num_delta_plus_old = bend.numDeltaPlus;
00795     const int num_delta_minus_old = bend.numDeltaMinus;
00796     if (bend.numDeltaPlus > 0)
00797     {
00798       const int di = bend.deltaIStart + num_delta_plus_old - 1;
00799       const double xp = iaSolution->x_solution[ di ]; 
00800       if (xp > 1.01)
00801       {
00802         double num_dp_added = floor(xp + 0.1);
00803 
00804         // sanity bound checks
00805         // at most double the number of delta pluses
00806         // there are two to start with, so this adds at least two
00807         double max_add = bend.numDeltaPlus;
00808         if (num_dp_added > max_add)
00809           num_dp_added = max_add;
00810         if (bend.numDeltaPlus > bend.num_deltas_max() - num_dp_added)
00811           num_dp_added = bend.num_deltas_max() - bend.numDeltaPlus;
00812         
00813         bend.numDeltaPlus += num_dp_added;
00814         
00815         if (bend.numDeltaPlus > num_delta_plus_old)
00816           new_bend = true;
00817 
00818         if (debugging)
00819         {
00820           const double xl = bend.xl; // relaxed solution
00821           const double g = iaData->I[i]; // goal
00822           printf("%d x (goal %f relaxed_floor %g) %f delta_plus[%d]=%f -> %g more delta pluses.\n", i, g, xl, iaSolution->x_solution[i], num_delta_plus_old-1, xp, num_dp_added );
00823         }
00824       }
00825     }
00826     if (bend.numDeltaMinus > 0)
00827     {
00828       const int di = bend.deltaIStart + num_delta_plus_old + num_delta_minus_old - 1;
00829       const double xm = iaSolution->x_solution[ di ]; 
00830       if (xm > 1.01 && bend.numDeltaMinus < bend.xl - 1)
00831       {
00832         double num_dm_added = floor(xm + 0.1);
00833         
00834         // sanity bound checks
00835         // at most double +1 the number of delta minuses
00836         // there is one to start with, so this adds at least two
00837         double max_add = bend.numDeltaMinus+1;
00838         if (num_dm_added > max_add)
00839           num_dm_added = max_add;
00840         if (bend.numDeltaMinus > bend.num_deltas_max() - num_dm_added)
00841           num_dm_added = bend.num_deltas_max() - bend.numDeltaMinus;
00842 
00843         // don't let x go below 1
00844         bend.numDeltaMinus += num_dm_added;
00845         if ( bend.numDeltaMinus > bend.xl - 1 )
00846         {
00847           bend.numDeltaMinus = bend.xl - 1;      
00848           num_dm_added = bend.numDeltaMinus - num_delta_minus_old;
00849         }
00850 
00851         if (bend.numDeltaMinus > num_delta_minus_old)
00852           new_bend = true;
00853 
00854         if (debugging)
00855         {
00856           const double xl = bend.xl; // relaxed solution
00857           const double g = iaData->I[i]; // goal
00858           printf("%d x (goal %f relaxed_floor %g) %f delta_minus[%d]=%f -> %g more delta minuses.\n", i, g, xl, iaSolution->x_solution[i], num_delta_minus_old-1, xm, num_dm_added );
00859         }
00860       }
00861     }
00862     bend.deltaIStart = d_start;
00863     
00864     d_start += bend.num_deltas();
00865   }
00866   
00867   // ======= new tilts
00868   // check for fractional deltas 
00869   // only do this if no new bends were added
00870   bool new_tilt = false; // was at least one new tilt added?
00871   if (!new_bend)
00872   {
00873     // the weight indices being current relies on there being no new bends added in the prior loop
00874     //IAWeights::iterator wi = weights.begin();
00875 
00876     int num_new_tilts = 0;
00877     for (int i = 0; i < iaData->num_variables(); ++i)
00878     {
00879 
00880       const double x = iaSolution->x_solution[i]; // current solution
00881       if (!is_integer(x))
00882       {
00883         IPBend &bend = bendData.bendVec[i]; // shorthand
00884 
00885         const double g = iaData->I[i]; // goal
00886         
00887         const int xf = floor(x); // int below x
00888         const int xc = xf+1.;  // int above x
00889         
00890         IPBend::IPTilt tilt;
00891 
00892         double r = ((double) rand() / RAND_MAX); // in 0,1
00893         tilt.second =  1.8 + r/2.;
00894         
00895         // if we are on a very flat part of the curve, straddling a goal,
00896         // increase the slope by a lot more
00897         if (fabs(x-g) < 1.1)
00898           tilt.second *= 3.756;
00899         
00900         // xc is farther from the goal, checks for case that the interval straddles the goal
00901         // tilt the positive branch
00902         if (xc - g > g - xf)
00903         {
00904           tilt.first = xc;
00905           tilt.second += x - xf;
00906           bend.plusTilts.push_back(tilt);
00907           
00908           merge_tilts(bend.plusTilts);
00909         }
00910         // tilt the negative branch
00911         else
00912         {
00913           tilt.first = xf;
00914           tilt.second += xc - x;
00915           
00916           bend.minusTilts.push_back(tilt);
00917           std::sort( bend.minusTilts.begin(), bend.minusTilts.end() ); // sorts by .first
00918           merge_tilts(bend.minusTilts);
00919         }
00920         ++num_new_tilts;
00921       }
00922     } // for i
00923     new_tilt = (num_new_tilts > 0);
00924     if (debugging && num_new_tilts)
00925     {
00926       printf("%d stubborn non-integer delta in solution\n", num_new_tilts);
00927     }
00928   }
00929   
00930   // if nothing changed, no need to redo weights, unless randomizing
00931   if ( !(new_bend || new_tilt || randomized) )
00932     return false;
00933       
00934   // generate new deltas and weights
00935   weights.clear();
00936   for (int i = 0; i < iaData->num_variables(); ++i)
00937   {
00938     add_bend_weights(i);
00939   }
00940   
00941     // if delta non-integer, apply tilts
00942     // randomize weight w/ excluded middle
00943 //    for (unsigned int j = 0; j < dp_non_int.size(); ++j)
00944 //    {
00945 //      int k = bend.deltaIStart - iaData->num_variables() + dp_non_int[j];
00946 //      // this might make weights tend to zero. Revisit later
00947 //      weights[ k ] *= 1. + 0.2 * IAWeights::rand_excluded_middle();
00948 //    }
00949 //    for (unsigned int j = 0; j < dm_non_int.size(); ++j)
00950 //    {
00951 //      int k = bend.deltaIStart + bend.numDeltaPlus - iaData->num_variables() + dm_non_int[j];
00952 //      // this might make weights tend to zero. Revisit later
00953 //      weights[ k ] *= 1. + 0.2 * IAWeights::rand_excluded_middle();
00954 //    }
00955 //    
00956     // todo: check that weights are still increasing with increasing k
00957   
00958    weights.uniquify(1., 1e6); // 1e4?
00959 
00960   // return if something changed
00961   return new_bend || new_tilt || randomized;
00962 
00963 }
00964   
00965 bool IASolverBend::solve()
00966 {
00967   if (debugging)
00968   {
00969   }
00970 
00971   myianlp = new IABendNlp(iaData, ipData, &bendData, iaSolution, &weights, silent);
00972   Ipopt::SmartPtr<Ipopt::TNLP> mynlp = myianlp; // Ipopt requires the use of smartptrs!
00973 
00974   // set initial ip bends from relaxed solution
00975   initialize_ip_bends();
00976 
00977   int iter = 0, bend_updates = 0;
00978   bool try_again = true;
00979   bool success = false;
00980   evenConstraintsActive = false;
00981   const int max_last_iter = 4 * ( 2 + iaData->num_variables() );
00982   const int max_first_iter = max_last_iter / 2;
00983   do 
00984   {
00985   
00986     // call the nlp solver
00987     bool solved = solve_nlp();
00988     
00989     try_again = false;
00990     if (solved)
00991     {
00992       
00993       // update bends based on solution
00994       bool changed = update_ip_bends();
00995       if (changed)
00996       {
00997         ++bend_updates;
00998         try_again = true;
00999       }
01000       // try again if sum-evens not already satisfied
01001       if (!even_constraints())
01002         try_again = true;
01003       // if nothing changed, or we've had enough initial iterations,
01004       // then activate the sum-even parabola constraints
01005       if (!changed || (iter > max_first_iter))
01006       {
01007         evenConstraintsActive = true;
01008         // zzyk switch over to the augmented problem, 
01009         // with linear constraints and standard goals for the sum-even variables.
01010         
01011       }
01012     }
01013     // avoid infinite loop if the method isn't working
01014     ++iter;
01015     if ( iter > max_last_iter )
01016       try_again = false;
01017     
01018   }
01019   while (try_again);
01020 
01021   //zzyk debug performance
01022   std::cout << iter << " bend iterations, " << bend_updates << " bend updates" << std::endl;
01023 
01024   cleanup();
01025   
01026   success = solution_is_integer() && all_constraints();
01027   if (success)
01028   {
01029     if (!silent)
01030       printf("IASolverBend produced integer and even solution\n");
01031   }
01032   else
01033   {
01034     return false; 
01035   }
01036   return success;
01037 }
01038 
01039 } // namespace MeshKit
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines