MeshKit  1.0
IAMINlp.cpp
Go to the documentation of this file.
00001 // IAMINlp.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/IAMINlp.hpp"
00014 #include "meshkit/IAData.hpp"
00015 #include "meshkit/IPData.hpp"
00016 #include "meshkit/IASolution.hpp"
00017 #include "meshkit/IANlp.hpp"
00018 
00019 #include <math.h>
00020 
00021 // for printf
00022 #ifdef HAVE_CSTDIO
00023 # include <cstdio>
00024 #else
00025 # ifdef HAVE_STDIO_H
00026 #  include <stdio.h>
00027 # else
00028 #  error "don't have header file for stdio"
00029 # endif
00030 #endif
00031 
00032 namespace MeshKit 
00033 {
00034     
00035 double p_norm = 3; // remove this, we use linear objectives here
00036 
00037 // constructor
00038 IAMINlp::IAMINlp(const IAData *data_ptr, const IPData *ip_data_ptr, IASolution *solution_ptr): 
00039 data(data_ptr), ip_data(ip_data_ptr), solution(solution_ptr),
00040 debugging(true), verbose(true) // true
00041 {
00042 
00043   printf("\nIAMINlp Problem size:\n");
00044   printf("  number of variables: %lu\n", data->I.size());
00045   printf("  number of constraints: %lu\n\n", data->constraints.size());
00046 }
00047 
00048 
00049 // n = number of variables
00050 // m = number of constraints (not counting variable bounds)
00051 
00052 // apparent Meshkit style conventions:
00053 // ClassFileName
00054 // ClassName
00055 // #ifndef MESHKIT_CLASSNAME_HPP
00056 // #define MESHKIT_CLASSNAME_HPP
00057 // namespace IAMeshKit
00058 // public class_member_function()
00059 // protected classMemberData
00060 
00061 IAMINlp::~IAMINlp() {data = NULL;}
00062 
00063 // returns the size of the problem
00064 bool IAMINlp::get_nlp_info(Index& n, Index& m, Index& nnz_jac_g,
00065                              Index& nnz_h_lag, IndexStyleEnum& index_style)
00066 {
00067   // printf("A ");
00068   // number of variables
00069   n = data->num_variables();
00070 
00071   // number of constraints
00072   m = (int) data->constraints.size() + (int) data->sumEvenConstraints.size();
00073 
00074   // number of non-zeroes in the Jacobian of the constraints
00075   size_t num_entries=0;
00076   for (std::vector<IAData::constraintRow>::const_iterator i=data->constraints.begin(); i != data->constraints.end(); ++i)
00077   {
00078         num_entries += i->M.size();
00079   }
00080   for (std::vector<IAData::sumEvenConstraintRow>::const_iterator i=data->sumEvenConstraints.begin(); i != data->sumEvenConstraints.end(); ++i)
00081   {
00082         num_entries += i->M.size();
00083   }
00084   nnz_jac_g = (int) num_entries;
00085 
00086   // number of non-zeroes in the Hessian 
00087   // diagonal entries for the objective function +
00088   // none for the constraints, since they are all linear
00089   nnz_h_lag = data->num_variables();
00090 
00091   // C style indexing (0-based)
00092   index_style = TNLP::C_STYLE;
00093 
00094   return true;
00095 }
00096 
00097 // returns the variable bounds
00098 bool IAMINlp::get_bounds_info(Index n, Number* x_l, Number* x_u,
00099                                 Index m, Number* g_l, Number* g_u)
00100 {
00101    //printf("B ");
00102 
00103   // The n and m we gave IPOPT in get_nlp_info are passed back to us.
00104   assert(n == data->num_variables());
00105   assert(m == (int)(data->constraints.size() + data->sumEvenConstraints.size()));
00106 
00107   // future interval upper and lower bounds:
00108   //   for midpoint subdivision, the lower bound may be 2 instead
00109   //   User may specify different bounds for some intervals
00110   //   Implement this by having another vector of lower bounds and one of upper bounds
00111 
00112   // relaxed problem
00113   if (ip_data->varIntegerBound.size() == 0)
00114   {
00115     // variables have lower bounds of 1 and no upper bounds
00116     for (Index i=0; i<n; ++i) 
00117     { 
00118       x_l[i] = 1.0; 
00119       x_u[i] = MESHKIT_IA_upperUnbound; 
00120     }
00121   }
00122   else
00123   {
00124     for (Index i=0; i<n; ++i) 
00125     {
00126       const double b = ip_data->varIntegerBound[i];
00127       if (b == 0) // unconstrained so far
00128       {
00129         x_l[i] = 1.0;
00130         x_u[i] = MESHKIT_IA_upperUnbound;
00131       }
00132       else
00133       {
00134         // b is a lower bound if it is bigger than the goal
00135         if (b > data->I[i]) 
00136         {
00137           x_l[i] = b;
00138           x_u[i] = MESHKIT_IA_upperUnbound; 
00139         }
00140         // otherwise b is an upper bound
00141         else
00142         {
00143           x_l[i] = 1.0; 
00144           x_u[i] = b; 
00145         }
00146       }
00147     }
00148   }
00149   
00150   // constraint bounds
00151   for (unsigned int i = 0; i<data->constraints.size(); ++i)
00152   {
00153     g_l[i] = data->constraints[i].lowerBound; 
00154     g_u[i] = data->constraints[i].upperBound; 
00155   }
00156   for (unsigned int i = 0; i<data->sumEvenConstraints.size(); ++i)
00157   {
00158     const int j = i + (int) data->constraints.size();
00159     g_l[j] = 4;
00160     g_u[j] = MESHKIT_IA_upperUnbound;
00161   }
00162 
00163   //printf("b ");
00164   return true; //means what?
00165 }
00166 
00167 // returns the initial point for the problem
00168 bool IAMINlp::get_starting_point(Index n, bool init_x, Number* x_init,
00169                                    bool init_z, Number* z_L, Number* z_U,
00170                                    Index m, bool init_lambda,
00171                                    Number* lambda)
00172 {
00173   //printf("C ");
00174 
00175   // Minimal info is starting values for x, x_init
00176   // Improvement: we can provide starting values for the dual variables if we wish
00177   assert(init_x == true);
00178   assert(init_z == false);
00179   assert(init_lambda == false);
00180 
00181   assert(n==(int)data->I.size());
00182   assert(data->num_variables()==(int)data->I.size());
00183   
00184   // initialize x to the goals
00185   if (ip_data->varIntegerBound.size() == 0)
00186     for (unsigned int i = 0; i<data->I.size(); ++i)
00187     {
00188       x_init[i]=data->I[i];
00189     }
00190   // initialize x to the prior solution
00191   else
00192     for (unsigned int i = 0; i<data->I.size(); ++i)
00193     {
00194       // todo: test if we need to modify this for x violating the variable bounds?
00195       x_init[i]= solution->x_solution[i]; 
00196     }
00197   
00198   return true;
00199 }
00200 
00201 
00202 // r
00203 inline Number IAMINlp::eval_r_i(const Number& I_i, const Number& x_i)
00204 {
00205   //return x_i / I_i + I_i / x_i;
00206   return x_i > I_i ? 
00207     (x_i - I_i) / I_i :
00208     (I_i - x_i) / x_i;
00209 }
00210 
00211 inline Number IAMINlp::eval_grad_r_i(const Number& I_i, const Number& x_i)
00212 {
00213   // return (1.0 / I_i) - I_i / (x_i * x_i);
00214   return x_i > I_i ? 
00215     1. / I_i :
00216     -I_i / (x_i * x_i);
00217 }
00218 
00219 inline Number IAMINlp::eval_hess_r_i(const Number& I_i, const Number& x_i)
00220 {
00221   // return 2.0 * I_i / (x_i * x_i * x_i);
00222   return x_i > I_i ? 
00223     0 : 
00224     2 * I_i / (x_i * x_i * x_i);
00225 }
00226 
00227 inline Number IAMINlp::eval_R_i(const Number& I_i, const Number& x_i)
00228 {
00229   if (x_i <= 0.)
00230     return MESHKIT_IA_upperUnbound;
00231   
00232   // old function, sum
00233   // return 2. + ( x_i / I_i ) * ( x_i / I_i ) + ( I_i / x_i ) * ( I_i / x_i );
00234   //  const register double r = (x_i - I_i) / I_i;
00235   //  const register double r2 = r*r;
00236   // return 2. + r2 + 1./r2;
00237   
00238   return pow( eval_r_i(I_i,x_i), p_norm );
00239 }
00240   
00241 
00242 inline Number IAMINlp::eval_grad_R_i(const Number& I_i, const Number& x_i)
00243 {
00244   if (x_i <= 0.)
00245     return MESHKIT_IA_lowerUnbound;
00246 
00247   return p_norm * pow( eval_r_i(I_i,x_i), p_norm-1) * eval_grad_r_i(I_i, x_i);
00248 }
00249 
00250 inline Number IAMINlp::eval_hess_R_i(const Number& I_i, const Number& x_i)
00251 {
00252   if (x_i <= 0.)
00253     return MESHKIT_IA_upperUnbound;
00254   
00255   // old function, sum
00256   // const register double I2 = I_i * I_i;
00257   // return 2. * ( 1. / I2 + 3. * I2 / (x_i * x_i * x_i * x_i) );
00258 
00259   if (p_norm == 2) 
00260     return (x_i > I_i) ? 2 / (I_i * I_i) : ( I_i / (x_i * x_i * x_i) ) * ( 6. * I_i / x_i - 4. ); 
00261   
00262   if (x_i > I_i ) 
00263     return p_norm * (p_norm - 1) * pow( (x_i - I_i) / I_i, p_norm -2) / (I_i * I_i); 
00264     
00265   // p (p-1) (I/x-1)^(p-2) I^2 / x^4 + p (I/x-1)^(p-1) 2 I / x^3  
00266   const register double f = (I_i - x_i) / x_i;
00267   const register double fp2 = pow( f, p_norm - 2);
00268   const register double x3 = x_i * x_i * x_i;
00269   return p_norm * ( (p_norm - 1) * fp2 * I_i * I_i / ( x3 * x_i ) + f * fp2 * 2. * I_i / x3 );
00270   
00271 }
00272 
00273 // s
00274 inline Number IAMINlp::eval_s_i(const Number& I_i, const Number& x_i)
00275 {
00276   return eval_r_i(x_i, I_i) * x_i;
00277 }
00278 
00279 inline Number IAMINlp::eval_grad_s_i(const Number& I_i, const Number& x_i)
00280 {
00281   return eval_r_i(x_i, I_i) + x_i * eval_grad_r_i(I_i, x_i);
00282 }
00283 
00284 inline Number IAMINlp::eval_hess_s_i(const Number& I_i, const Number& x_i)
00285 {
00286   return 2. * eval_grad_r_i(I_i, x_i) + x_i * eval_hess_r_i(I_i, x_i);
00287 }
00288 
00289 inline Number IAMINlp::eval_S_i(const Number& I_i, const Number& x_i)
00290 {
00291   if (x_i <= 0.)
00292     return MESHKIT_IA_upperUnbound;
00293   const double s = eval_s_i(I_i,x_i);
00294   
00295   if (p_norm == 1)
00296     return s;
00297     
00298   if (p_norm == 2)
00299     return s*s;
00300   
00301   assert(p_norm > 2);
00302   return pow(s, p_norm );
00303 }
00304 
00305 
00306 inline Number IAMINlp::eval_grad_S_i(const Number& I_i, const Number& x_i)
00307 {
00308   if (x_i <= 0.)
00309     return MESHKIT_IA_lowerUnbound;
00310   
00311   const double s = eval_s_i(I_i,x_i);
00312   const double sp = eval_grad_s_i(I_i, x_i);
00313 
00314   if (p_norm == 1)
00315     return s;
00316   
00317   assert(p_norm>1);
00318   if (p_norm == 2)
00319     return 2. * s * sp; 
00320 
00321   assert(p_norm > 2);
00322   return p_norm * pow( s, p_norm-1) * sp;
00323 }
00324 
00325 inline Number IAMINlp::eval_hess_S_i(const Number& I_i, const Number& x_i)
00326 {
00327   if (x_i <= 0.)
00328     return MESHKIT_IA_upperUnbound;
00329 
00330   const double spp = eval_hess_s_i(I_i, x_i);
00331   
00332   if (p_norm == 1)
00333     return spp;
00334   
00335   const double s = eval_s_i(I_i, x_i);
00336   const double sp = eval_grad_s_i(I_i, x_i);
00337 
00338   if (p_norm==2)
00339     return 2. * ( sp * sp + s * spp );
00340   
00341   assert(p_norm>2);
00342   return p_norm * pow(s, p_norm - 2) * ( (p_norm - 1) * sp * sp + s * spp );
00343 }
00344 
00345 
00346 // experimental results: 
00347 // best so far is f=R and p_norm = 3
00348 // f=s doesn't work well for x < I, as it tends to push the x to 1 so the weight is small !
00349 
00350 // returns the value of the objective function
00351 bool IAMINlp::eval_f(Index n, const Number* x, bool new_x, Number& obj_value)
00352 {
00353   //printf("D ");
00354 
00355   assert(n == data->num_variables());
00356 
00357   // future: perhaps try max of the f, infinity norm
00358   
00359   double obj = 0.;
00360   for (Index i = 0; i<n; ++i)
00361   {
00362     obj += eval_R_i( data->I[i], x[i] );
00363   }
00364 
00365   obj_value = obj;
00366 
00367   return true;
00368 }
00369 
00370 // return the gradient of the objective function grad_{x} f(x)
00371 bool IAMINlp::eval_grad_f(Index n, const Number* x, bool new_x, Number* grad_f)
00372 {
00373   //printf("E ");
00374 
00375   assert(n == data->num_variables());
00376 
00377   for (Index i = 0; i<n; ++i)
00378   {
00379     grad_f[i] = eval_grad_R_i( data->I[i], x[i] ); 
00380   }
00381 
00382   return true;
00383 }
00384 
00385 // return the value of the constraints: g(x)
00386 bool IAMINlp::eval_g(Index n, const Number* x, bool new_x, Index m, Number* g)
00387 {
00388   //printf("F ");
00389 
00390   assert(n == data->num_variables());
00391   assert(m == (int)(data->constraints.size() + data->sumEvenConstraints.size()));
00392   
00393   
00394   for (Index i = 0; i<(int)data->constraints.size(); ++i)
00395   {
00396     //double g_i = constraints[i].rhs;
00397     double g_i = 0.;
00398     for (Index j = 0; j < (int)data->constraints[i].M.size(); ++j)
00399     {
00400       g_i += x[ data->constraints[i].M[j].col ] * data->constraints[i].M[j].val;
00401     }
00402     g[i] = g_i;
00403   }
00404   for (Index i = 0; i<(int)data->sumEvenConstraints.size(); ++i)
00405   {
00406     //double g_i = constraints[i].rhs;
00407     const int k = (int) data->constraints.size() + i;
00408     double g_k = 0.;
00409     for (Index j = 0; j < (int)data->sumEvenConstraints[i].M.size(); ++j)
00410     {
00411       g_k += x[ data->sumEvenConstraints[i].M[j].col ] * data->sumEvenConstraints[i].M[j].val;
00412     }
00413     g[k] = g_k;
00414   }
00415 
00416   return true;
00417 }
00418 
00419 // return the structure or values of the jacobian
00420 bool IAMINlp::eval_jac_g(Index n, const Number* x, bool new_x,
00421                            Index m, Index nele_jac, Index* iRow, Index *jCol,
00422                            Number* values)
00423 {
00424   //printf("G ");
00425 
00426   assert(m == (int)(data->constraints.size() + data->sumEvenConstraints.size()));
00427   
00428   if (values == NULL) {
00429     // return the structure of the jacobian
00430     Index k=0;
00431     for (Index i = 0; i<(int)data->constraints.size(); ++i)
00432     {
00433       for (Index j = 0; j < (int)data->constraints[i].M.size(); ++j)
00434       {
00435                 iRow[k] = i;
00436         jCol[k] = data->constraints[i].M[j].col;
00437         ++k;
00438       }
00439     }
00440     for (Index i = 0; i< (int)data->sumEvenConstraints.size(); ++i)
00441     {
00442       for (Index j = 0; j < (int)data->sumEvenConstraints[i].M.size(); ++j)
00443       {
00444                 iRow[k] = i + (int) data->constraints.size();
00445         jCol[k] = data->sumEvenConstraints[i].M[j].col;
00446         ++k;
00447       }
00448     }
00449 
00450   }
00451   else {
00452     // return the values of the jacobian of the constraints
00453     Index k=0;
00454     for (Index i = 0; i < (int)data->constraints.size(); ++i)
00455     {
00456       for (Index j = 0; j < (int)data->constraints[i].M.size(); ++j)
00457       {
00458         values[k++] = data->constraints[i].M[j].val;
00459       }
00460     }
00461     for (Index i = 0; i < (int)data->sumEvenConstraints.size(); ++i)
00462     {
00463       for (Index j = 0; j < (int)data->sumEvenConstraints[i].M.size(); ++j)
00464       {
00465         values[k++] = data->sumEvenConstraints[i].M[j].val;
00466       }
00467     }
00468   }
00469 
00470   return true;
00471 }
00472 
00473 //return the structure or values of the hessian
00474 bool IAMINlp::eval_h(Index n, const Number* x, bool new_x,
00475                        Number obj_factor, Index m, const Number* lambda,
00476                        bool new_lambda, Index nele_hess, Index* iRow,
00477                        Index* jCol, Number* values)
00478 {
00479   //printf("H ");
00480 
00481   // get_nlp_info specified the number of non-zeroes, should be
00482   assert(nele_hess == data->num_variables()); 
00483   
00484   if (values == NULL) {
00485     // return the structure. This is a symmetric matrix, fill the lower left
00486     // triangle only.
00487 
00488         // because the constraints are linear
00489     // the hessian for this constraints are actually empty
00490     
00491     // Since the objective function is separable, only the diagonal is non-zero
00492     for (Index idx=0; idx<data->num_variables(); ++idx)
00493     {
00494       iRow[idx] = idx;
00495       jCol[idx] = idx;
00496     }
00497 
00498   }
00499   else {
00500     // return the values. This is a symmetric matrix, fill the lower left
00501     // triangle only
00502 
00503         // because the constraints are linear
00504     // the hessian for this problem is actually empty
00505 
00506     // fill the objective portion
00507     for (Index idx=0; idx<data->num_variables(); ++idx)
00508     {
00509       values[idx] = obj_factor * eval_hess_R_i(data->I[idx], x[idx]);
00510     }
00511   }
00512 
00513   return true;
00514 }
00515 
00516 void IAMINlp::finalize_solution(SolverReturn status,
00517                                   Index n, const Number* x, const Number* z_L, const Number* z_U,
00518                                   Index m, const Number* g, const Number* lambda,
00519                                   Number obj_value,
00520                                   const IpoptData* ip_data,
00521                                   IpoptCalculatedQuantities* ip_cq)
00522 {
00523   // here is where we would store the solution to variables, or write to a file, etc
00524   // so we could use the solution.
00525   //printf("I ");
00526 
00527   // copy solution into local storage x_solution
00528   solution->x_solution.clear(); // clear contents
00529   std::vector<double>(solution->x_solution).swap(solution->x_solution); // zero capacity
00530   solution->x_solution.reserve(n); // space for new solution
00531   for (Index i=0; i<n; i++) 
00532   {
00533     solution->x_solution.push_back( x[i] );  // values of new solution
00534   }
00535   assert( (int)solution->x_solution.size() == n );
00536   solution->obj_value = obj_value;
00537   assert(obj_value >= 0.);
00538   
00539   if (debugging)
00540   {
00541     printf("NLP solution:\n");       
00542     printf("x[%d] = %e\n", 0, x[0]);
00543     double rhs = 0., lhs = 0.;
00544     if (verbose)
00545     {
00546       printf("legend: coeff x_i (solution, goal, ratio; f(x), f'(x); F(x), F'(x) )\n");
00547       for (unsigned int j = 0; j < data->constraints.size(); ++j)
00548       {
00549         printf("constraint %d: ", j);
00550         const IAData::constraintRow & c = data->constraints[j];
00551         for (std::vector<IAData::sparseEntry>::const_iterator i = c.M.begin(); i < c.M.end(); ++i)
00552         {
00553           const double xv = x[i->col];
00554           const double gv = data->I[i->col];
00555           const double r = xv > gv ? (xv-gv) / gv : (gv - xv) / xv;
00556           printf(" %1.0f x_%d (%1.3f, %1.1f, %1.1f; %2.2g, %2.2g; %2.2g, %2.2g) ", 
00557                  i->val, i->col, xv, gv, 
00558                  r, eval_r_i(gv, xv), eval_grad_r_i(gv, xv), 
00559                  eval_R_i(gv,xv), eval_grad_R_i(gv,xv) );
00560           if (i->val > 0) 
00561             lhs += i->val * x[i->col];
00562           else
00563             rhs += i->val * x[i->col];
00564         }
00565         if (data->constraints.front().upperBound == data->constraints.front().lowerBound)
00566           printf(" = %1.1f", data->constraints.front().upperBound);  
00567         else
00568           printf(" in [%1.1f,%1.1f]", data->constraints.front().upperBound, data->constraints.front().lowerBound );
00569         printf(" <=> %f %f in ...\n", lhs, rhs );
00570       }
00571     
00572       printf("\n\nSolution of the primal variables, x\n");
00573       for (Index i=0; i<n; i++) {
00574         printf("x[%d] = %e\n", i, x[i]);
00575       }
00576       
00577       printf("\n\nSolution of the bound multipliers, z_L and z_U\n");
00578       for (Index i=0; i<n; i++) {
00579         printf("z_L[%d] = %e\n", i, z_L[i]);
00580       }
00581       for (Index i=0; i<n; i++) {
00582         printf("z_U[%d] = %e\n", i, z_U[i]);
00583       }
00584       printf("\nFinal value of the constraints:\n");
00585       for (Index i=0; i<m ;i++) {
00586         printf("g(%d) = %e\n", i, g[i]);
00587       }
00588     }
00589     
00590     printf("\n\nObjective value\n");
00591     printf("f(x*) = %e\n", obj_value);
00592     
00593   }
00594 }
00595 
00596 } // namespace MeshKit
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines