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