MeshKit  1.0
IABendNlp.cpp
Go to the documentation of this file.
00001 // IABendNlp.cpp
00002 // Interval Assignment for Meshkit
00003 //
00004 // Adapted from
00005 // Copyright (C) 2005, 2006 International Business Machines and others.
00006 // All Rights Reserved.
00007 // This code is published under the Eclipse Public License.
00008 //
00009 // $Id: hs071_nlp.cpp 1864 2010-12-22 19:21:02Z andreasw $
00010 //
00011 // Authors:  Carl Laird, Andreas Waechter     IBM    2005-08-16
00012 
00013 #include "meshkit/IABendNlp.hpp"
00014 #include "meshkit/IAData.hpp"
00015 #include "meshkit/IPData.hpp"
00016 #include "meshkit/IPBend.hpp"
00017 #include "meshkit/IASolution.hpp"
00018 #include "meshkit/IANlp.hpp"
00019 #include "meshkit/IASolverToolInt.hpp"
00020 
00021 #include <math.h>
00022 #include <limits>
00023 #include <algorithm>
00024 
00025 // for printf
00026 #ifdef HAVE_CSTDIO
00027 # include <cstdio>
00028 #else
00029 # ifdef HAVE_STDIO_H
00030 #  include <stdio.h>
00031 # else
00032 #  error "don't have header file for stdio"
00033 # endif
00034 #endif
00035 
00036 
00037 namespace MeshKit 
00038 {
00039   
00040 
00041 // the objective function we should use for weights
00042 double IABendNlp::f_x_value( double I_i, double x_i )
00043 {
00044   return x_i > I_i ? 
00045   (x_i - I_i) / I_i :
00046   /* 1.133402234045 * */(I_i - x_i) / x_i;
00047   // I don't think this extra constant factor is necessary to break ties.
00048   // could just call baseNlp::eval_r_i
00049 
00050   // expected is 1/100 to 4 or so
00051 
00052   //??? is the following true?
00053   // the range of magnitudes of these weights is too high, the weights are too non-linear if we take them to a power
00054   // was
00055   // const double fh = IANlp::eval_R_i(data->I[i], xl+1);
00056 }
00057 
00058 // constructor
00059 IABendNlp::IABendNlp(const IAData *data_ptr, const IPData *ip_data_ptr, const IPBendData *ip_bend_data_ptr,
00060                      IASolution *solution_ptr, IAWeights *weight_ptr, const bool set_silent): 
00061                      data(data_ptr), ipData(ip_data_ptr), ipBendData(ip_bend_data_ptr), solution(solution_ptr), 
00062                      baseNlp(data_ptr, solution_ptr, set_silent),
00063                      base_n((int)data_ptr->num_variables()),
00064                      base_m((int)(data_ptr->constraints.size() + data_ptr->sumEvenConstraints.size())),
00065                      problem_n(data_ptr->num_variables() + (int) weight_ptr->size()),
00066                      problem_m( (int)(data_ptr->constraints.size() + data_ptr->num_variables() + ( evenConstraintsActive ? 2 : 1) * data_ptr->sumEvenConstraints.size())),  
00067                      weights( weight_ptr ),
00068 //  silent(set_silent), debugging(true), verbose(true) // true
00069                     silent(set_silent), debugging(false), verbose(false) // true
00070 {
00071   if (!silent)
00072   { 
00073     printf("\nIABendNlp Problem size:\n");
00074     printf("  number of variables: %lu\n", data->I.size());
00075     printf("  number of equality constraints: %lu\n", data->constraints.size());
00076     printf("  number of even constraints: %lu\n\n", data->sumEvenConstraints.size());
00077     // report on bends?
00078   }
00079 }
00080 
00081 
00082 // n = number of variables
00083 // m = number of constraints (not counting variable bounds)
00084 
00085 
00086 IABendNlp::~IABendNlp() {data = NULL;}
00087 
00088 // returns the size of the problem
00089 bool IABendNlp::get_nlp_info(Index& n, Index& m, Index& nnz_jac_g,
00090                              Index& nnz_h_lag, IndexStyleEnum& index_style)
00091 {
00092   if (debugging) printf("IABendNlp::get_nlp_info\n");
00093   
00094   problem_m = (int) (data->constraints.size() + data->num_variables() + ( evenConstraintsActive ? 2 : 1) * data->sumEvenConstraints.size());
00095   wave_even_constraint_start = base_m + data->num_variables(); // base + deltas
00096 
00097   bool base_ok = baseNlp.get_nlp_info(n, m, nnz_jac_g, nnz_h_lag, index_style);
00098   assert(n == base_n);
00099   assert(m == base_m);
00100   
00101   // add n for delta variables
00102   n += (Index) weights->size();
00103   
00104   // add m for delta-computing constraints
00105   m += (Index) data->num_variables();
00106   
00107   // add m for even-constraints
00108   if ( evenConstraintsActive )
00109     m += (Index) data->sumEvenConstraints.size();
00110   
00111   // Jacobian entries
00112 
00113   // data assumption checks
00114   if (debugging)
00115   {
00116     int num_deltas = 0;
00117     for (IPBendVec::const_iterator i=ipBendData->bendVec.begin(); 
00118          i != ipBendData->bendVec.end(); ++i )
00119     {
00120       num_deltas += i->num_deltas();
00121       printf("num_deltas: %d = %d + %d \n", i->num_deltas(), i->numDeltaPlus, i->numDeltaMinus);
00122     }
00123     assert( num_deltas == (int) weights->size() );
00124   }
00125   
00126   // base_program adds equality and even>4 constraints
00127   
00128   // delta constraints, x = sum deltas
00129   nnz_jac_g += weights->size() + data->num_variables();
00130   
00131   // wave even constraints
00132   // const Index base_nnz_jac_g = nnz_jac_g;
00133   if (evenConstraintsActive)
00134   {
00135     int num_even_entries = 0;
00136     for (std::vector<IAData::sumEvenConstraintRow>::const_iterator i=data->sumEvenConstraints.begin(); i != data->sumEvenConstraints.end(); ++i)
00137     {
00138       num_even_entries += i->M.size();
00139     }
00140     nnz_jac_g += num_even_entries;
00141   }
00142 
00143   // hessian for even constraints
00144   // objective functions for deltas are all linearized, so second derivatives are zero
00145   nnz_h_lag = 0; // nele_h // overwrite the non-linear objective
00146   build_hessian();
00147   nnz_h_lag = (Index) hessian_vector.size();
00148   
00149   if (debugging)
00150   {
00151     printf("IABendNlp::get_nlp_info: n %d, m %d, nnz_jac_g %d, nnz_h_lag, %d\n",
00152            n, m, nnz_jac_g, nnz_h_lag);
00153   }
00154   
00155   return true && base_ok;
00156 }
00157 
00158 // returns the variable bounds
00159 bool IABendNlp::get_bounds_info(Index n, Number* x_l, Number* x_u,
00160                                 Index m, Number* g_l, Number* g_u)
00161 {
00162   if (debugging) printf("IABendNlp::get_bounds_info\n");
00163   
00164   const bool base_ok = baseNlp.get_bounds_info(base_n, x_l, x_u, base_m, g_l, g_u);
00165 
00166   // delta variable bounds
00167   for (IPBendVec::const_iterator i=ipBendData->bendVec.begin(); 
00168        i != ipBendData->bendVec.end(); ++i )
00169   {
00170     for (int j = 0; j < i->numDeltaPlus; ++j)
00171     {
00172       int k = i->deltaIStart + j;
00173       assert(k<n);
00174       x_l[k] = 0;
00175       x_u[k] = 1.;      
00176       //    if (debugging)
00177       //      printf("x %d (%f) delta+_%d in [%d,%d]\n", *i, x, j, (int) x_l[k], (int) x_u[k]);
00178     }
00179     // last delta is unlimitted in theory, 
00180     // but in practice to avoid an unbounded above solution arising from the flattening, 
00181     // bound it by some large but reasonable value compared to the number of deltas
00182     const int last_plus_k = i->deltaIStart + i->numDeltaPlus - 1;
00183     assert(last_plus_k < n);
00184     // x_u[ last_plus_k ] = MESHKIT_IA_upperUnbound; // too big
00185     x_u[ last_plus_k ] = i->numDeltaPlus + 4.5; 
00186 
00187     for (int j = 0; j < i->numDeltaMinus; ++j)
00188     {
00189       int k = i->deltaIStart + i->numDeltaPlus + j;
00190       assert(k<n);
00191       x_l[k] = 0;
00192       x_u[k] = 1.;      
00193       //    if (debugging)
00194       //      printf("x %d (%f) delta+_%d in [%d,%d]\n", *i, x, j, (int) x_l[k], (int) x_u[k]);
00195     }
00196     // last delta is only limited by x >= 1 in theory,
00197     // but in practice to avoid an unbounded below solution arising from the flattening, 
00198     // bound it by some large but reasonable value compared to the number of deltas
00199     const int last_minus_k = i->deltaIStart + i->numDeltaPlus + i->numDeltaMinus - 1;
00200     assert(last_minus_k < n);
00201     // bound by the min of the next two values
00202     const double max_dm_by_relative_change = i->numDeltaMinus + 4.5;
00203     const double max_dm_by_x_above_1 = i->xl - i->numDeltaMinus;
00204     const double max_dm = (max_dm_by_relative_change < max_dm_by_x_above_1) ? max_dm_by_relative_change : max_dm_by_x_above_1;
00205     x_u[ last_minus_k ] = max_dm;
00206    
00207   }
00208   
00209   // delta constraint bounds : they are equality constraints
00210   for (Index i = 0; i < (int) ipBendData->bendVec.size(); ++i)
00211   {
00212     const int k = i + base_m;
00213     assert(k < m);
00214     g_l[k] = 0.;
00215     g_u[k] = 0.;
00216   }
00217   
00218   // sum-even wave constraint bounds
00219   if (evenConstraintsActive)
00220   {
00221     for (unsigned int i = 0; i < data->sumEvenConstraints.size(); ++i)
00222     {
00223       // parabola(s) == 1 
00224       // equality makes ipopt think it is overdetermined, so use inequality :) 
00225       const int k = i + wave_even_constraint_start;
00226       assert(k < m);
00227       g_l[k] = 0.9999; // 1 - 1.e-4 means the sum will be within 1.e-2 ( because 2 = 4/2 ) of an integer. to do: tweak this tolerance?
00228       g_u[k] = MESHKIT_IA_upperUnbound; // 1.
00229       // the actual value of g will be below 1, so this contrains it to be 1 just as well
00230     }
00231     assert(wave_even_constraint_start + (int) data->sumEvenConstraints.size() == m);
00232   }
00233   else
00234   {
00235     assert((int) ipBendData->bendVec.size() + base_m == m);
00236   }
00237   
00238   if (debugging)
00239   {
00240     for (int i = 0; i < n; ++i)
00241     {
00242       printf("x[%d] in [%g,%g]\n", i, x_l[i], x_u[i]);
00243     }
00244     for (int j = 0; j < m; ++j)
00245     {
00246       printf("g[%d] in [%g,%g]\n", j, g_l[j], g_u[j]);
00247     }
00248   }
00249   return true && base_ok;
00250 }
00251 
00252 // returns the initial point for the problem
00253 bool IABendNlp::get_starting_point(Index n, bool init_x, Number* x_init,
00254                                    bool init_z, Number* z_L, Number* z_U,
00255                                    Index m, bool init_lambda,
00256                                    Number* lambda)
00257 {
00258   if (debugging) printf("IABendNlp::get_starting_point\n");
00259 
00260   // idea: use last solution, not necessarily relaxed one ?? would that be better or worse?
00261 
00262   // Minimal info is starting values for x, x_init
00263   // Improvement: we can provide starting values for the dual variables if we wish
00264   assert(init_x == true);
00265   assert(init_z == false);
00266   assert(init_lambda == false);
00267 
00268   // initialize x to the relaxed solution
00269   {
00270     Index i;
00271     for (i=0; i<data->num_variables(); ++i) 
00272     {
00273       assert(i<n);
00274       x_init[i] = ipData->relaxedSolution[i];
00275     }
00276   
00277     // alternative:
00278     // deltas are all zero
00279     // that point is actually not feasible because the delta constraints are not satisfied
00280     for ( ; i < n; ++i )
00281     {
00282       assert(i<n);
00283       x_init[i] = 0.;
00284     }
00285   }
00286   
00287   // delta - add non-integer part of x to the delta
00288   {
00289     assert( (int) ipBendData->bendVec.size() == data->num_variables() );
00290     for (int i = 0; i < (int) ipBendData->bendVec.size(); ++i )
00291     {
00292       double d = ipData->relaxedSolution[i] - ipBendData->bendVec[i].xl;
00293       // assert( fabs(d) < 1.01 );
00294       if ( d > 0. )
00295       {
00296         for (int di = 0; di < ipBendData->bendVec[i].numDeltaPlus; ++di)
00297         {
00298           const int k = ipBendData->bendVec[i].deltaIStart + di;
00299           assert(k<n);
00300           if (( d <= 1. ) || ( di + 1 == ipBendData->bendVec[i].numDeltaPlus))
00301           {
00302             x_init[k] = d;
00303             break;
00304           }
00305           else
00306           {
00307             x_init[k] = 1.;
00308             d -= 1.;
00309           }
00310         }
00311       }
00312       else
00313       {
00314         d = fabs(d);
00315         for (int di = 0; di < ipBendData->bendVec[i].numDeltaMinus; ++di)
00316         {
00317           const int k = ipBendData->bendVec[i].deltaIStart + ipBendData->bendVec[i].numDeltaPlus + di;
00318           assert(k<n);
00319           if (( d <= 1. ) || ( di + 1 == ipBendData->bendVec[i].numDeltaMinus))
00320           {
00321             x_init[k] = fabs(d);
00322             d = 0;
00323             break;
00324           }
00325           else
00326           {
00327             x_init[k] = 1.;
00328             d -= 1.;
00329           }
00330         }
00331         
00332       }
00333     }
00334   }
00335   return true;
00336 }
00337 
00338 
00339 // returns the value of the objective function
00340 bool IABendNlp::eval_f(Index n, const Number* x, bool new_x, Number& obj_value)
00341 {
00342   if (debugging) printf("IABendNlp::eval_f\n");
00343   assert(n == data->num_variables() + (int) weights->size());
00344 
00345   // only the deltas contribute
00346   double obj = 0.;
00347   Index i(0);
00348   Index k(data->num_variables());
00349   for ( ; k<n; ++i, ++k)
00350   {
00351     obj += (*weights)[i] * x[k];
00352   }
00353   obj_value = obj;
00354 
00355   return true;
00356 }
00357 
00358 // return the gradient of the objective function grad_{x} f(x)
00359 bool IABendNlp::eval_grad_f(Index n, const Number* x, bool new_x, Number* grad_f)
00360 {
00361   const bool debug_loc = false;
00362   if (debugging && debug_loc) printf("IABendNlp::eval_grad_f\n");
00363 
00364   assert(n == data->num_variables() + (int) weights->size());
00365 
00366   // only the deltas contribute
00367   Index i;
00368   for (i = 0; i<data->num_variables(); ++i)
00369   {
00370     grad_f[i] = 0;
00371   }
00372 
00373   // deltas contribute their weight
00374   Index k(data->num_variables());
00375   for (i = 0; i < (int) weights->size(); ++i, ++k)
00376   {
00377     grad_f[k] = (*weights)[i];
00378   }
00379 
00380   return true;
00381 }
00382 
00383 // return the value of the constraints: g(x)
00384 bool IABendNlp::eval_g(Index n, const Number* x, bool new_x, Index m, Number* g)
00385 {
00386   const bool debug_loc = false;
00387   if (debugging && debug_loc) printf("IABendNlp::eval_g\n");
00388   bool base_ok = baseNlp.eval_g(base_n, x, new_x, base_m, g);
00389 
00390   // sum x - deltas
00391   //       x = xl + delta_plusses - delta_minuses
00392   // g = (     xl + delta_plusses - delta_minuses - x  ), constrained to be zero
00393   for (Index i = 0; i < (int) ipBendData->bendVec.size(); ++i)
00394   {
00395     const int k = i + base_m;
00396     const IPBend &bend = ipBendData->bendVec[i];
00397     double gk = bend.xl - x[i];
00398     if (debugging && debug_loc)
00399       printf("constraint %d: sum of deltas for x_%d: xl (%d) - x (%f) ", k, i, bend.xl, x[i]);
00400     for (int j = 0; j < bend.numDeltaPlus; ++j )
00401     {
00402       int di = bend.deltaIStart + j;
00403       gk += x[di];
00404       if (debugging && debug_loc)
00405         printf(" + %f", x[di]);
00406     }
00407     for (int j = 0; j < bend.numDeltaMinus; ++j )
00408     {
00409       int di = bend.deltaIStart + bend.numDeltaPlus + j;
00410       gk -= x[di];
00411       if (debugging && debug_loc)
00412         printf(" - %f", x[di]);
00413     }
00414     g[k] = gk;
00415     if (debugging && debug_loc)
00416       printf(" = %f\n", gk);
00417   }
00418   
00419   // sum evens, parabola value
00420   if (evenConstraintsActive)
00421   {
00422     for (unsigned int i = 0; i < data->sumEvenConstraints.size(); ++i)
00423     {
00424       const int k = i + wave_even_constraint_start;
00425       double s = baseNlp.eval_even_sum(i,x);
00426       const double gk = eval_g_int_s(s); // e.g. (s - nearest_even_number)^2;
00427       g[k] = gk;    
00428       if (debugging)
00429       {
00430         printf("constraint %d: wave even %d, sum = %f, g = %f\n",
00431                k, i, s, gk );
00432       }
00433     }
00434 
00435   }
00436   
00437   return true && base_ok;
00438 }
00439 
00440 // return the structure or values of the jacobian
00441 bool IABendNlp::eval_jac_g(Index n, const Number* x, bool new_x,
00442                            Index m, Index nele_jac, Index* iRow, Index *jCol,
00443                            Number* values)
00444 {
00445   const bool debug_loc = false;
00446   if (debugging && debug_loc) printf("IABendNlp::eval_jac_g\n");
00447   
00448   int base_nele_jac = baseNlp.get_neleJac();
00449   bool base_ok = baseNlp.eval_jac_g(base_n, x, new_x, base_m, base_nele_jac, iRow, jCol, values);
00450   base_nele_jac = baseNlp.get_neleJac();
00451   int kk = base_nele_jac;
00452   
00453   // sum x - deltas
00454   //       x = xl + delta_plusses - delta_minuses
00455   // g = (     xl + delta_plusses - delta_minuses - x  ), constrained to be zero
00456   // g' =       0 + 1 dp_indexes  - 1 dm_indexes - 1 x_i
00457   
00458 
00459   if (values == NULL) 
00460   {
00461     if (debugging && debug_loc)
00462     {
00463       printf("jacobian g structure:\n");
00464     }
00465     // deltas
00466     for (Index i = 0; i < (int) ipBendData->bendVec.size(); ++i)
00467     {
00468       const int k = i + base_m;
00469       const IPBend &bend = ipBendData->bendVec[i];
00470       if (debugging && debug_loc)
00471       {
00472         printf("constraint %d: jac for sum-deltas for x_%d, ", k, i);
00473       }
00474       iRow[kk] = k;
00475       jCol[kk] = i; // x[i]
00476       ++kk;
00477       for (int j = 0; j < bend.numDeltaPlus; ++j )
00478       {
00479         int di = bend.deltaIStart + j;
00480         iRow[kk] = k;
00481         jCol[kk] = di;
00482         if (debugging && debug_loc)
00483         {
00484           printf(" deltaplus %d (x[%d])", j, di );
00485         }
00486         ++kk;
00487       }
00488       for (int j = 0; j < bend.numDeltaMinus; ++j )
00489       {
00490         int di = bend.deltaIStart + bend.numDeltaPlus + j;
00491         iRow[kk] = k;
00492         jCol[kk] = di;
00493         if (debugging && debug_loc)
00494         {
00495           printf(" deltaminus %d (x[%d])", j, di );
00496         }
00497         ++kk;
00498       }
00499       if (debugging && debug_loc) printf("\n");
00500     }
00501     
00502     // sum evens
00503     if (evenConstraintsActive)
00504     {
00505       for (unsigned int i = 0; i< data->sumEvenConstraints.size(); ++i)
00506       {
00507         const int k = i + wave_even_constraint_start;
00508         if (debugging)
00509         {
00510           printf("wave even constraint %d: jac, for natural even constraint %d, ", k, i);
00511         }
00512         
00513         for (unsigned int j = 0; j < data->sumEvenConstraints[i].M.size(); ++j)
00514         {
00515           iRow[kk] = k;
00516           jCol[kk] = data->sumEvenConstraints[i].M[j].col;
00517           if (debugging)
00518           {
00519             printf(" kk %d (i %d, j %d)", kk, iRow[kk], jCol[kk]);
00520           }
00521           ++kk;
00522         }
00523       }
00524       if (debugging) printf("\n");
00525     }
00526     assert( kk == nele_jac );
00527     if (debugging && debug_loc)
00528     {
00529         printf("Jacobian g sparse structure:\n");
00530       for (int i = 0; i < nele_jac; ++i)
00531       {
00532         printf("entry %d: row %d col %d\n", i, iRow[i], jCol[i] );
00533               
00534       }
00535     }
00536   }
00537   else
00538   {
00539     if (debugging)
00540     {
00541       printf("jacobian g values:\n");
00542     }
00543 
00544     // deltas
00545     for (Index i = 0; i < (int) ipBendData->bendVec.size(); ++i)
00546     {
00547       if (debugging && debug_loc)
00548       {
00549         printf("constraint %d: jac for sum-deltas for x_%d, ", i + base_m, i);
00550       }
00551 
00552       const IPBend &bend = ipBendData->bendVec[i];
00553       values[kk++] = -1; // - x[i]
00554       if (debugging) printf("-1 (x_%d) ", i);
00555 
00556       for (int j = 0; j < bend.numDeltaPlus; ++j )
00557       {
00558         values[kk++] = 1; // + deltaplus_j
00559         if (debugging) printf("+1 (d+_%d) ", j);
00560       }
00561       for (int j = 0; j < bend.numDeltaMinus; ++j )
00562       {
00563         values[kk++] = -1; // - deltaplus_j
00564         if (debugging) printf("-1 (d-_%d) ", j);
00565       }
00566       if (debugging) printf("\n");
00567     }
00568     
00569     // sum evens
00570     if (evenConstraintsActive)
00571     {
00572       
00573       for (unsigned int i = 0; i< data->sumEvenConstraints.size(); ++i)
00574       {
00575         const int k = i + wave_even_constraint_start;
00576         if (debugging)
00577         {
00578           printf("wave even constraint %d: jac, for natural even constraint %d, ", k, i);
00579         }
00580         const double s = baseNlp.eval_even_sum(i, x);
00581         const double jac_gk = eval_jac_int_s(s); // e.g. -PI * cos( PI * s );
00582         for (unsigned int j = 0; j < data->sumEvenConstraints[i].M.size(); ++j)
00583         {
00584           const double coeff = data->sumEvenConstraints[i].M[j].val;
00585           values[kk++] = coeff * jac_gk;
00586           if (debugging)
00587           {
00588             printf("  %d: x_%d gradient %f * %f = %f\n", kk-1,
00589                    data->sumEvenConstraints[i].M[j].col, coeff, jac_gk, coeff * jac_gk);
00590           }
00591         }
00592         if (debugging) printf("\n");
00593       }
00594     }
00595     
00596     assert( kk == nele_jac );
00597     if (debugging && debug_loc)
00598     {
00599       printf("Jacobian g sparse values:\n");
00600       for (int i = 0; i < nele_jac; ++i)
00601       {
00602         printf("entry %d: %f\n", i, values[i]);        
00603       }
00604     }
00605 
00606   }
00607   
00608   return true && base_ok;
00609 }
00610 
00611   
00612   
00613 IABendNlp::SparseMatrixEntry::SparseMatrixEntry(const int iset, const int jset, const int kset)
00614 {
00615   if ( jset > iset )
00616   {
00617     i = jset;
00618     j = iset;
00619   }
00620   else
00621   {
00622     i = iset;
00623     j = jset;
00624   }
00625   k = kset;
00626   assert(j <= i);
00627 }
00628   
00629 void IABendNlp::add_hessian_entry( int i, int j, int &k )
00630 {
00631   SparseMatrixEntry sme(i, j, k);
00632   
00633   if (hessian_map.insert( std::make_pair(sme.key(), sme) ).second)
00634   {
00635     hessian_vector.push_back( sme );
00636     ++k;
00637     assert( (int) hessian_vector.size() == k );
00638     assert( (int) hessian_map.size() == k );
00639   }
00640   // else it was already inserted, nothing to do
00641 }
00642   
00643 int IABendNlp::SparseMatrixEntry::n(0);
00644   
00645 void IABendNlp::build_hessian() 
00646 {
00647   if (debugging) printf("IABendNlp::build_hessian\n");
00648 
00649   // only build once
00650   if (hessian_vector.size())
00651     return;
00652   
00653   hessian_vector.clear();
00654   hessian_map.clear();
00655   SparseMatrixEntry::n = data->num_variables() + (int) weights->size(); // variables + deltas
00656   int kk = 0;
00657   
00658   if (!evenConstraintsActive)
00659   {
00660     if (debugging)
00661       printf("even constraints not active => empty hessian\n");
00662     return;
00663   }
00664   
00665   // sum_evens
00666   // g = cos( pi * sum_evens) == 1
00667   // g' = -pi * sin ( pi * sum_evens ), for each of the variables contributing to the sum
00668   // g''= -pi^2 cos ( pi * sum_evens ), for each pair of variables contributing to the sum
00669   // assuming all the coefficients are 1    
00670   for (unsigned int c = 0; c< data->sumEvenConstraints.size(); ++c)
00671   {
00672     for (unsigned int i = 0; i < data->sumEvenConstraints[c].M.size(); ++i)
00673     {
00674       int ii = data->sumEvenConstraints[c].M[i].col;
00675       for (unsigned int j = 0; j <=i; ++j)
00676       {
00677         int jj = data->sumEvenConstraints[c].M[j].col;
00678         if (debugging)
00679           printf("hessian adding %d %d %d\n", ii, jj, kk);
00680         add_hessian_entry( ii, jj, kk );
00681       }
00682     }
00683   }
00684   
00685   // nele_hess = hessian_vector.size();
00686   if (debugging)
00687   {
00688     printf("==========\nBuilt Hessian\n");
00689     print_hessian();
00690     printf("==========\n");
00691   }
00692 }
00693   
00694 int IABendNlp::get_hessian_k( int i, int j )
00695 {
00696   if ( i == j )
00697     return i;
00698   SparseMatrixEntry sme(i, j, -1 );
00699   SparseMatrixEntry &entry = hessian_map[ sme.key() ]; 
00700   return entry.k;
00701 }
00702   
00703 void IABendNlp::print_hessian()
00704 {
00705   printf("Packed Hessian:\n");
00706   for (unsigned int kk = 0; kk < hessian_vector.size(); ++kk)
00707   {
00708     SparseMatrixEntry &sme = hessian_vector[kk];
00709     printf(" %d: (%d, %d)\n", sme.k, sme.i, sme.j);
00710   }
00711   printf("\n");          
00712   
00713   printf("Random Access Hessian in sequence\n");
00714   {
00715     int k = 0;
00716     SparseMatrixMap::iterator iter;
00717     for (iter = hessian_map.begin(); iter != hessian_map.end(); ++iter, ++k)
00718     {
00719       const SparseMatrixEntry &sme = iter->second;
00720       printf(" %d: (%d, %d) k:%d key:%d\n", k, sme.i, sme.j, sme.k, sme.key() );
00721     }
00722     printf("\n");          
00723   }
00724   
00725   printf("Random Access Hessian in square:\n");
00726   for (int i = 0; i < data->num_variables(); ++i)
00727   {
00728     for (int j = 0; j < data->num_variables(); ++j)
00729     {
00730       SparseMatrixEntry sme_search(i, j, -1 );
00731       SparseMatrixMap::iterator iter = hessian_map.find(sme_search.key());
00732       if (iter == hessian_map.end())
00733         printf("   . ");
00734       else
00735         printf(" %3d ", iter->second.k);
00736     }
00737     printf("\n");          
00738   }
00739   printf("\n");          
00740   
00741 }
00742   
00743 //return the structure or values of the hessian
00744 bool IABendNlp::eval_h(Index n, const Number* x, bool new_x,
00745                        Number obj_factor, Index m, const Number* lambda,
00746                        bool new_lambda, Index nele_hess, Index* iRow,
00747                        Index* jCol, Number* values)
00748 {
00749   if (debugging) printf("IABendNlp::eval_h\n");
00750 
00751   // fill values with zeros
00752   if (values)
00753   {
00754     for (unsigned int kk = 0; kk < hessian_vector.size(); ++kk)
00755     {
00756       values[kk] = 0.;
00757     }
00758     // debug, print x
00759   }
00760   
00761   // base hessian is empty, as all those constraints are linear and the obj. function is linearized 
00762   
00763   // hessian entry i,j is:
00764   // obj_factor fij + sum_k lambda[k] gkij
00765   // where fij =   d^2 f / d x_i d x_j
00766   //       gkij =  d^2 g[k] / d x_i d x_j
00767   // and d denotes partial derivative
00768   
00769   // This is a symmetric matrix, fill the lower left triangle only.
00770   if (values == NULL) {
00771     assert((int) hessian_vector.size() == nele_hess);
00772     // return the structure. 
00773     for (unsigned int kk = 0; kk < hessian_vector.size(); ++kk)
00774     {
00775       iRow[kk] = hessian_vector[kk].i;
00776       jCol[kk] = hessian_vector[kk].j;
00777     } 
00778   } // structure
00779   else {
00780     // return the values. 
00781     // g =    ( sum_evens - nearest_sum_even )^2 == 0
00782     // g' = 2 ( sum_evens - nearest_sum_even )  for each of the variables contributing to the sum
00783     // g''= 2                                   for each pair of variables contributing to the sum
00784     // assuming all the coefficients are 1
00785     if (evenConstraintsActive)
00786     {
00787       if (debugging)
00788       {
00789         printf("\nwave even non-zero hessian values:");
00790       }
00791       for (unsigned int i = 0; i< data->sumEvenConstraints.size(); ++i)
00792       {
00793         if (debugging)
00794         {
00795           printf("\n%d constraint: ", i);
00796         }
00797         const int k = i + wave_even_constraint_start; // index of the constraint in the problem, = lambda to use
00798         
00799         const double s = baseNlp.eval_even_sum(i, x); // sum of the variable values
00800         // second derivative of wave function, assuming coefficients of one
00801         const double wavepp = eval_hess_int_s(s); // e.g. ( -PI * PI * cos( PI * s ) ); 
00802         const double h_value = lambda[k] * wavepp;        
00803         // assign that value to all pairs, weighted by coeff_i * coeff_j
00804         for (unsigned int ii = 0; ii < data->sumEvenConstraints[i].M.size(); ++ii)
00805         {
00806           const int var_i_index = data->sumEvenConstraints[i].M[ii].col;
00807           const double coeff_i =  data->sumEvenConstraints[i].M[ii].val;
00808           for (unsigned int jj = 0; jj < data->sumEvenConstraints[i].M.size(); ++jj)
00809           {
00810             const int var_j_index = data->sumEvenConstraints[i].M[jj].col;
00811             const double coeff_j =  data->sumEvenConstraints[i].M[jj].val;
00812             
00813             const int kk = get_hessian_k(var_i_index, var_j_index);
00814             assert( kk >= 0 && kk < nele_hess );
00815             values[kk] += coeff_i * coeff_j * h_value;
00816             if (debugging)
00817             {
00818               printf("  lambda[%d] * coeff_%d * coeff_%d * d^2 wave / d x_%d d x_%d = %f * %f * %f * %f\n", 
00819                      k, var_i_index, var_j_index, var_i_index, var_j_index, 
00820                      lambda[k], coeff_i, coeff_j, wavepp );
00821             }
00822           }
00823         }
00824       }
00825       if (debugging)
00826       {
00827         printf("\n");
00828       }
00829       assert((int) hessian_vector.size() == nele_hess);
00830     }
00831 
00832   } // values
00833   
00834   return true;
00835 }
00836 
00837 void IABendNlp::finalize_solution(SolverReturn status,
00838                                       Index n, const Number* x, const Number* z_L, const Number* z_U,
00839                                       Index m, const Number* g, const Number* lambda,
00840                                       Number obj_value,
00841                                       const IpoptData* ip_data,
00842                                       IpoptCalculatedQuantities* ip_cq)
00843 {
00844   if (debugging) printf("IABendNlp::finalize_solution\n");
00845   // by passing n and not base_n, m and not base_m,
00846   // we all get the deltas and sum-even values
00847   baseNlp.finalize_solution(status, n, x, z_L, z_U, m, g, lambda, obj_value, ip_data, ip_cq);
00848 }
00849 
00850 } // namespace MeshKit
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines