1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
// IASolverInt.cpp
// Interval Assignment for Meshkit
//
#include "meshkit/IASolverInt.hpp"
#include "meshkit/IAData.hpp"
#include "meshkit/IPData.hpp"
#include "meshkit/IASolution.hpp"
#include "meshkit/IARoundingNlp.hpp"
#include "meshkit/IASolverRelaxed.hpp"
#include "meshkit/IASolverBend.hpp"
// #include "meshkit/IAMINlp.hpp"
#include "meshkit/IAIntCosNlp.hpp"
#include "meshkit/IAIntParabolaNlp.hpp"

// #include "meshkit/IAMilp.hpp" // glpk based solution too slow

#include <stdio.h>
#include <math.h>
#include <limits.h>

#include "IpIpoptApplication.hpp"

namespace MeshKit 
{
    
IASolverInt::IASolverInt(const IAData * ia_data_ptr, IASolution *relaxed_solution_ptr, 
  const bool set_silent) 
  : IASolverToolInt(ia_data_ptr, relaxed_solution_ptr, true), 
  silent(false),
//  silent(set_silent), 
  debugging(true)
{ 
  ip_data(new IPData);
  // initialize copies relaxed solution, then we can overwrite relaxed_solution_pointer
  ip_data()->initialize(relaxed_solution_ptr->x_solution); 
}

/** default destructor */
IASolverInt::~IASolverInt() 
{
  delete ip_data();
}
    
bool IASolverInt::solve_wave_workhorse(IAIntWaveNlp *mynlp)
{
  if (debugging)
  {
    printf("IASolverInt::solve_wave() - ");        
    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");
  }
  
  // solver setup  
  Ipopt::SmartPtr<Ipopt::IpoptApplication> app = IpoptApplicationFactory();
/* try leaving defaults
  // convergence parameters
  // see $IPOPTDIR/Ipopt/src/Interfaces/IpIpoptApplication.cpp
  // our real criteria are: all integer, constraints satisfied. How to test the "all_integer" part?
  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
  app->Options()->SetNumericValue("max_cpu_time", sqrt( iaData->num_variables() ) ); // max time allowed in seconds
  app->Options()->SetIntegerValue("max_iter", 3 * (10 + iaData->num_variables() ) ); // max number of iterations
  // app->Options()->SetNumericValue("primal_inf_tol", 1e-2 ); 
  app->Options()->SetNumericValue("dual_inf_tol", 1e-2 ); // how close to infeasibility? // was 1e-2
  app->Options()->SetNumericValue("constr_viol_tol", 1e-2 ); // by how much can constraints be violated?
  app->Options()->SetNumericValue("compl_inf_tol", 1e-6 ); // max norm of complementary condition // was 1e-2
  
  // second criteria convergence parameters: quit if within this tol for many iterations
  // was  app->Options()->SetIntegerValue("acceptable_iter", 4 + sqrt( iaData->num_variables() ) ); //as "tol"
  app->Options()->SetNumericValue("acceptable_tol", 1e-6 ); //as "tol" was 1e-1
  
  app->Options()->SetStringValue("mu_strategy", "adaptive");
  // print level 0 to 12, Ipopt default is 5
  const int print_level = (silent) ? 0 : 1;  // simple info is 1, debug at other values
  app->Options()->SetIntegerValue("print_level", print_level);  
  // uncomment next line to write the solution to an output file
  // app->Options()->SetStringValue("output_file", "IA.out");  
  // The following overwrites the default name (ipopt.opt) of the options file
  // app->Options()->SetStringValue("option_file_name", "IA.opt");
  
  */
  
  // Intialize the IpoptApplication and process the options
  Ipopt::ApplicationReturnStatus status;
  status = app->Initialize();
  if (status != Ipopt::Solve_Succeeded) {
    if (!silent)
      printf("\n\n*** Error during initialization!\n");
    return (int) status;
  }
  
  bool try_again = true;
  int iter = 0;
  
  // print();
  bool solution_ok = false;
  
  do {
    if (debugging)
    {
      print();
      printf("%d IntWave iteration\n", iter );
      // build the hessian, evaluate it and f at the current solution?
    }
      
    // Ask Ipopt to solve the problem
    status = app->OptimizeTNLP(mynlp); // the inherited IANlp
    
    // see /CoinIpopt/build/include/coin/IpReturnCodes_inc.h
    /*
    Solve_Succeeded=0,
    Solved_To_Acceptable_Level=1,
    Infeasible_Problem_Detected=2,
    Search_Direction_Becomes_Too_Small=3,
    Diverging_Iterates=4,
    User_Requested_Stop=5,
    Feasible_Point_Found=6,
    
    Maximum_Iterations_Exceeded=-1,
    Restoration_Failed=-2,
    Error_In_Step_Computation=-3,
    Maximum_CpuTime_Exceeded=-4,
    Not_Enough_Degrees_Of_Freedom=-10,
    Invalid_Problem_Definition=-11,
    Invalid_Option=-12,
    Invalid_Number_Detected=-13,
    
    Unrecoverable_Exception=-100,
    NonIpopt_Exception_Thrown=-101,
    Insufficient_Memory=-102,
    Internal_Error=-199
     */

    bool solved_full = false;
    bool solved_partial = false;
    bool solver_failed = false;
    bool bad_problem = false;
    // setting void just to avoid compiler warning: variable 'solver_failed' set but not used [-Wunused-but-set-variable]
    (void) solver_failed;
    (void) bad_problem;
    switch (status) {
      case Ipopt::Solve_Succeeded:
      case Ipopt::Solved_To_Acceptable_Level:
      case Ipopt::Feasible_Point_Found:
        solved_full = true;
        break;
      case Ipopt::Maximum_Iterations_Exceeded:
      case Ipopt::User_Requested_Stop:
      case Ipopt::Maximum_CpuTime_Exceeded:
        solved_partial = true;
        break;
      case Ipopt::Infeasible_Problem_Detected:
      case Ipopt::Not_Enough_Degrees_Of_Freedom:
      case Ipopt::Invalid_Problem_Definition:
      case Ipopt::Invalid_Option:
      case Ipopt::Invalid_Number_Detected:
        bad_problem = true;
        break;
      case Ipopt::Search_Direction_Becomes_Too_Small:
      case Ipopt::Restoration_Failed:
      case Ipopt::Diverging_Iterates:
      case Ipopt::Error_In_Step_Computation:
      case Ipopt::Unrecoverable_Exception:
      case Ipopt::NonIpopt_Exception_Thrown:
      case Ipopt::Insufficient_Memory:
      case Ipopt::Internal_Error:        
        solver_failed = true;
        break;
        
      default:
        break;
    }
  
    if (!silent)
    {
      if (solved_full) {
        printf("\n\n*** IntWave solved!\n");
      }
      else {
        printf("\n\n*** IntWave FAILED!\n");
      }
    }
    
    if (debugging)
    {
      printf("\nChecking solution.\n");
      bool integer_sat = solution_is_integer(true);
      bool even_sat = even_constraints( false, true);
      bool equal_sat = equal_constraints( false, true );
      printf("IntWave solution summary, %s, equal-constraints %s, even-constraints %s.\n", 
             integer_sat ? "integer" : "NON-INTEGER",
             equal_sat ? "satisfied" : "VIOLATED", 
             even_sat ? "satisfied" : "VIOLATED" );
      if (!integer_sat)
        printf("investigate integer neighborhood\n");
      if (!even_sat)
        printf("investigate even neighborhood\n");
      if (!equal_sat)
        printf("investigate equal neighborhood\n");
    }
    
    
    IASolution nlp_solution;
    nlp_solution.x_solution = ia_solution()->x_solution; // vector copy
    IASolverToolInt sti( ia_data(), &nlp_solution );
    sti.round_solution();
    if (debugging)
      printf("Checking rounded solution, ignoring even constraints.\n");
    if (sti.equal_constraints(false, debugging))
    {
      // also even constraints
      if (debugging)
        printf("Rounding worked.\n");
      
      // rounding was a valid integer solution
      ia_solution()->x_solution.swap( nlp_solution.x_solution );
      // ia_solution()->obj_value is no longer accurate, as it was for the non-rounded solution
      return true;
    }
    
    // todo: detect and act
    // may have converged to a locally optimal, but non-feasible solution
    // if so, try a new starting point
    
    // check solution feasibility, even when not debugging
    
    if ( solved_full || solved_partial )
    {
      bool integer_sat = solution_is_integer(false);
      bool even_sat = even_constraints( false, false);
      bool equal_sat = equal_constraints( false, false );
      if ( integer_sat && even_sat && equal_sat )
        return true;
    }

    // find out which vars were not integer, 
    // try moving to a farther starting point resolving
    
 
    try_again = false; 
  } while (try_again);
  
  
  // now 
  // As the SmartPtrs go out of scope, the reference count
  // will be decremented and the objects will automatically
  // be deleted.  
  return solution_ok;
  
}


bool IASolverInt::solve_round()
{
  // set up and call the separate IARoundingNlp, which has a linear objective to get a natural integer solution 
  // the intuition is this will solve integrality for  most variables all at once

  if (debugging)
  {
    printf("IASolverInt::solve_bend_workhorse() - ");        
  }

  
  // solver setup  
  Ipopt::SmartPtr<Ipopt::IpoptApplication> app = IpoptApplicationFactory();
  
  // convergence parameters
  // see $IPOPTDIR/Ipopt/src/Interfaces/IpIpoptApplication.cpp
  // our real criteria are: all integer, constraints satisfied. How to test the "all_integer" part?
  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
  app->Options()->SetNumericValue("max_cpu_time", sqrt( iaData->num_variables() ) ); // max time allowed in seconds
  app->Options()->SetIntegerValue("max_iter", 3 * iaData->num_variables() ); // max number of iterations
  // app->Options()->SetNumericValue("primal_inf_tol", 1e-2 ); 
  app->Options()->SetNumericValue("dual_inf_tol", 1e-6 ); // how close to infeasibility? // was 1e-2
  app->Options()->SetNumericValue("constr_viol_tol", 1e-6 ); // by how much can constraints be violated?
  app->Options()->SetNumericValue("compl_inf_tol", 1e-6 ); // max norm of complementary condition // was 1e-2

  // second criteria convergence parameters: quit if within this tol for many iterations
// was  app->Options()->SetIntegerValue("acceptable_iter", 4 + sqrt( iaData->num_variables() ) ); //as "tol"
  app->Options()->SetNumericValue("acceptable_tol", 1e-6 ); //as "tol" was 1e-1

  app->Options()->SetStringValue("mu_strategy", "adaptive");
  // print level 0 to 12, Ipopt default is 5
  const int print_level = (silent) ? 0 : 1; 
  app->Options()->SetIntegerValue("print_level", print_level);  
  // uncomment next line to write the solution to an output file
  // app->Options()->SetStringValue("output_file", "IA.out");  
  // The following overwrites the default name (ipopt.opt) of the options file
  // app->Options()->SetStringValue("option_file_name", "IA.opt");
  
  // Intialize the IpoptApplication and process the options
  Ipopt::ApplicationReturnStatus status;
  status = app->Initialize();
  if (status != Ipopt::Solve_Succeeded) {
    if (!silent)
      printf("\n\n*** Error during initialization!\n");
    return (int) status;
  }
  
  
  Ipopt::TNLP *tnlp = NULL;

  IARoundingNlp *myianlp = new IARoundingNlp(iaData, ipData, iaSolution, silent);
  if (debugging) 
  {          
    printf("ROUNDING problem formulation\n");
    printf("Attempting to find a naturally-integer solution by linearizing the objective function.\n");
    printf("Variables are constrained within [floor,ceil] of relaxed solution.\n");
  }
  
  // problem setup
  // a couple of different models, simplest to more complex
  // IARoundingFarNlp *myianlp = new IARoundingFarNlp(iaData, ipData, this);
  // IARoundingFar3StepNlp *myianlp = new IARoundingFar3StepNlp(iaData, ipData, this); // haven't tested this. It compiles and runs but perhaps isn't correct
  // IAIntWaveNlp *myianlp = new IAIntWaveNlp(iaData, ipData, this); // haven't tested this. It compiles and runs but perhaps isn't correct

  tnlp = myianlp;
  Ipopt::SmartPtr<Ipopt::TNLP> mynlp = tnlp; // Ipopt requires the use of smartptrs!

  bool try_again = true;
  int iter = 0;
  do {
    printf("%d rounding iteration\n", iter );
    
    // Ask Ipopt to solve the problem
    status = app->OptimizeTNLP(mynlp); // the inherited IANlp
    
    if (!silent)
    {
      if (status == Ipopt::Solve_Succeeded) {
        printf("\n\n*** The problem solved!\n");
      }
      else {
        printf("\n\n*** The problem FAILED!\n");
      }
    }
    
    // The problem should have been feasible, but it is possible that it had no integer solution.
    // figure out which variables are still integer
    
    // check solution for integrality and constraint satified    
    if (debugging)
    {
      printf("\nChecking Natural (non-rounded) solution.\n");
      bool integer_sat = solution_is_integer(true);
      bool even_sat = even_constraints( false, true);
            bool equal_sat = equal_constraints( false, true );
      printf("Natural solution summary, %s, equal-constraints %s, even-constraints %s.\n", 
             integer_sat ? "integer" : "NON-INTEGER",
             equal_sat ? "satisfied" : "VIOLATED", 
             even_sat ? "satisfied" : "VIOLATED" );
      if (!integer_sat)
        printf("investigate this\n");
    }
    
    IASolution nlp_solution;
    nlp_solution.x_solution = ia_solution()->x_solution; // vector copy
    IASolverToolInt sti( ia_data(), &nlp_solution );
    sti.round_solution();
    if (debugging)
      printf("Checking rounded solution, ignoring even constraints.\n");
    if (sti.equal_constraints(false, debugging))
    {
      // also even constraints
      if (debugging)
        printf("Rounding worked.\n");

      // rounding was a valid integer solution
      ia_solution()->x_solution.swap( nlp_solution.x_solution );
      // ia_solution()->obj_value is no longer accurate, as it was for the non-rounded solution
      return true;
    }

    // find out which vars were not integer, 
    // try rounding their weights and resolving
    // bool int_sat = solution_is_integer();
    ++iter;
    try_again = iter < 4 + sqrt(iaData->num_variables());
    if (try_again)
    {
      if (debugging)
        printf("rounding failed, randomizing weights\n");
    
      myianlp->randomize_weights_of_non_int(); // try again? debug
    }
    else if (debugging)
      printf("giving up on rounding to non-integer solution\n");

    // try_again = false; // debug
  } while (try_again);

  
  // todo: update partially-integer solution, perhaps using ipData - figure out how we're going to use it, first, for what structure makes sense.
  
  // As the SmartPtrs go out of scope, the reference count
  // will be decremented and the objects will automatically
  // be deleted.  
  return status == Ipopt::Solve_Succeeded;
  
}

void IASolverInt::cleanup()
{
  ;
}

bool IASolverInt::solve_wave(const SolverType solver_type)
{
  IAIntWaveNlp *myianlp = NULL;
  if (solver_type == COS) 
  {
    if (debugging) printf("Cosine wave.\n");
    myianlp= new IAIntCosNlp(iaData, ipData, iaSolution, silent);
  }
  else if (solver_type == PARABOLA)
  {
    if (debugging) printf("Parabola wave.\n");
    myianlp = new IAIntParabolaNlp(iaData, ipData, iaSolution, silent);
  }
  else
  {
    if (debugging) printf("Invalid wave type.\n");
    return false;
  }
    
  Ipopt::SmartPtr<Ipopt::TNLP> mynlp = myianlp; // Ipopt requires the use of smartptrs!
  return solve_wave_workhorse( myianlp );
}

bool IASolverInt::solve()
{
  
  SolverType solver_type = BEND; // BEND;
  
  // debug, try solve_intwave instead 
  // longer term, use intwave as a backup when the faster and simpler milp doesn't work.
  // unfortunately, it appears to find local minima that are far from optimal, even when starting in a well
  
  bool solved = false;
  (void) solved;
  switch (solver_type) {
    case COS:
    case PARABOLA:
      solved = solve_wave(solver_type);
      break;
    case ROUNDING:
      solved = solve_round();
      break;
    case BEND:
      {
        IASolverBend sb(iaData, iaSolution, silent);
        solved = sb.solve();
      }
      break;
    default:
      solved = false;
      break;<--- Variable 'solved' is assigned a value that is never used.
  }

  bool success = solution_is_integer();
  // also check constraints and evenality
  if (success)
  {
    if (!silent)
      printf("IASolverInt produced integer solution\n");
  }
  else
  {
    // 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)
    return false; // debug;
    // solve_rounding_heuristic();
    // success = solution_is_integer();
  }
  return success;
}

} // namespace MeshKit