MeshKit  1.0
IASolverTool.cpp
Go to the documentation of this file.
00001 // IASolverTool.cpp
00002 // Interval Assignment for Meshkit
00003 //
00004 
00005 #include "meshkit/IASolverTool.hpp"
00006 #include "meshkit/IAData.hpp"
00007 #include "meshkit/IASolution.hpp"
00008 #include "meshkit/IPData.hpp"
00009 
00010 #include <cstdlib>
00011 #include <stdio.h>
00012 #include <math.h>
00013 #include <limits.h>
00014 
00015 namespace MeshKit 
00016 {
00017 
00018 IASolverTool::~IASolverTool()
00019 {
00020   iaData = NULL; 
00021   iaSolution=NULL;
00022 }
00023   
00024 bool IASolverTool::is_integer(const double x, double &x_int_double) const
00025 {
00026   x_int_double = ( floor( x + 0.5 ) );
00027   // beware cases where this epsilon is too large, e.g. sides with more than 100 curves = 1 / 1.e-2
00028   return fabs( x - x_int_double ) < 1.e-2; 
00029 }
00030 
00031 bool IASolverTool::is_integer(const double x, int &x_int) const
00032 {
00033   double x_int_double;
00034   bool is_int = is_integer( x, x_int_double );
00035   x_int = (int) ( x_int_double );
00036   return is_int; 
00037 }
00038 
00039   
00040 bool IASolverTool::is_integer(const double x) const
00041 {
00042   double x_int_double;
00043   return is_integer( x, x_int_double );
00044 }
00045 
00046   
00047 bool IASolverTool::valid_solution() const
00048 {
00049   return (iaSolution && iaData && (iaSolution->x_solution.size() >= iaData->I.size())); // could be bigger if deltas are retained
00050 }
00051   
00052     
00053 void IASolverTool::print(const bool do_print_solution, const bool do_print_nonint,
00054                          const bool do_print_equal_constraints, const bool do_print_nonequal,
00055                          const bool do_print_even_constraints, const bool do_print_noneven ) const
00056 {
00057   
00058   int info_case;
00059   /*  
00060    four cases
00061    0 nothing
00062    1 data_and_solution
00063    2 data_only
00064    3 solution_only 
00065    */
00066   if (iaData)
00067   {
00068     if (valid_solution())
00069       info_case = 1; //data_and_solution = true;
00070     else
00071       info_case = 2; // data_only = true;    
00072   }
00073   else
00074   {
00075     if (iaSolution && iaSolution->x_solution.size())
00076       info_case = 3; // solution_only = true;
00077     else
00078       info_case = 0; // nothing = true;
00079   }
00080   
00081   // header
00082   switch (info_case) {
00083     case 0:
00084       printf("no data, no solution\n");
00085       return;
00086       break;
00087 
00088     case 1:
00089       printf("\nIA data and solution:\n");
00090       break;
00091       
00092     case 2:
00093       printf("\nIA data:\n");
00094       break;
00095 
00096     case 3:
00097       printf("\nIA solution:\n");
00098       break;
00099 
00100     default:
00101       break;
00102   }
00103 
00104   // variables
00105   if ( info_case == 1 || info_case == 2 ) // data exists
00106   {
00107     printf("%d vars\n", iaData->num_variables());
00108     for (int i=0; i<iaData->num_variables(); ++i)
00109     {
00110       printf("%d x (goal %e) ",i, iaData->I[i]);
00111       if (do_print_solution && info_case == 1)
00112       {
00113         //printf(" relaxed %e solution ", ip_data.relaxedSolution[i]);
00114         const double x = iaSolution->x_solution[i];
00115         int x_int;
00116         if (is_integer(x,x_int))
00117           printf("%d\n",x_int);
00118         else
00119           printf("%e  NON-INTEGER\n",x);
00120       }
00121       else
00122       {
00123         printf("\n");
00124       }
00125     }
00126   }
00127   else if ( info_case == 3 && do_print_solution)
00128   {
00129     printf("%d x solution values\n", iaData->num_variables());
00130     for (unsigned int i=0; i<iaSolution->x_solution.size(); ++i)
00131     {
00132       printf("x_%d ",i);
00133       //printf(" relaxed %e solution ", ip_data.relaxedSolution[i]);
00134       const double x = iaSolution->x_solution[i];
00135       int x_int;
00136       if (is_integer(x,x_int))
00137         printf("%d\n",x_int);
00138       else
00139         printf("%e  NON-INTEGER\n",x);
00140     }
00141   }
00142   
00143   if ( info_case == 1 || info_case == 2 ) // data exists
00144   {
00145     if (do_print_equal_constraints || do_print_nonequal)
00146     {
00147       printf("%lu equality constraints:\n", iaData->constraints.size());
00148       equal_constraints(do_print_equal_constraints, do_print_nonequal );
00149     }
00150     if (do_print_even_constraints || do_print_noneven)
00151     {
00152       printf("%lu even constraints:\n", iaData->sumEvenConstraints.size());
00153       even_constraints(do_print_even_constraints, do_print_noneven);
00154     }
00155   }
00156   
00157   if (do_print_solution && valid_solution())
00158     printf("objective function value %e\n", iaSolution->obj_value);
00159   printf("\n");
00160 }
00161   
00162 void IASolverTool::print_solution() const
00163 {
00164   print( true, false, false, false, false, false );
00165 }
00166   
00167 void IASolverTool::print_problem() const
00168 {
00169   print( false, true, false, false, false, false );
00170 }
00171   
00172 bool IASolverTool::is_even(double y, int &y_even) const
00173 {
00174   // nearest even value, including 0 and 2
00175   y_even = (int) floor( 2. * floor( y / 2. + 0.5 ) + 1.e-2 );
00176   if ( fabs( y - y_even ) < 1.0e-2 )
00177     return true;
00178   return false;
00179 }
00180   
00181 bool IASolverTool::is_even(double y) const
00182 {
00183   int y_even;
00184   return is_even(y, y_even);
00185 }
00186   
00187 double IASolverTool::equal_value(const int i) const
00188 {
00189   double g_i = 0.; // no rhs, it is placed into the bounds instead
00190   if (valid_solution())
00191   {
00192     for (unsigned int j = 0; j < iaData->constraints[i].M.size(); ++j)
00193     {
00194       g_i += iaSolution->x_solution[ iaData->constraints[i].M[j].col ] * iaData->constraints[i].M[j].val;
00195     }
00196   }
00197   return g_i;
00198 }
00199   
00200   
00201 double IASolverTool::even_value(const int i) const
00202 {
00203   double g_i = - iaData->sumEvenConstraints[i].rhs;
00204   if (valid_solution())
00205   {
00206     for (unsigned int j = 0; j < iaData->sumEvenConstraints[i].M.size(); ++j)
00207     {
00208       double x = iaSolution->x_solution[ iaData->sumEvenConstraints[i].M[j].col ];
00209       double x_int_double;
00210       if (is_integer(x, x_int_double))
00211         x = x_int_double;          
00212       g_i += x * iaData->sumEvenConstraints[i].M[j].val;
00213     }
00214   }       
00215   return g_i;
00216 }
00217 
00218 void IASolverTool::even_floor_ceil(double s, double &s_floor, double &s_ceil) const
00219 {
00220   s_floor = 2. * floor(s / 2);
00221   s_ceil = s_floor + 2.;
00222 }
00223   
00224 void IASolverTool::int_floor_ceil(double x, double &x_floor, double &x_ceil) const
00225 {
00226   x_floor = floor(x);
00227   x_ceil = x_floor + 1.;
00228 }
00229   
00230   // compare to 
00231   // bool IANlp::eval_g(Index n, const Number* x, bool new_x, Index m, Number* g)
00232   // does the same thing with the data in ipopt format
00233 
00234 std::pair< bool, double> IASolverTool::equal_constraint(const int i, const bool print_me, const bool print_unsatisfied ) const
00235 {
00236   const double g_i = equal_value(i);
00237   // epsilons values
00238   const double epsilon_lower = 1.e-2;
00239   const double epsilon_upper = 1.e-2;
00240   bool satisfied = 
00241     (g_i + epsilon_lower > iaData->constraints[i].lowerBound) || 
00242     (g_i < iaData->constraints[i].upperBound + epsilon_upper);
00243   if (print_me || (print_unsatisfied && !satisfied))
00244   {
00245     printf("equal constraint %d:", i);
00246     for (unsigned int j = 0; j < iaData->constraints[i].M.size(); ++j)
00247     {
00248       printf( " %f x_%d", iaData->constraints[i].M[j].val, iaData->constraints[i].M[j].col );
00249       if (valid_solution())
00250         printf(" (%f)", 
00251                iaSolution->x_solution[ iaData->constraints[i].M[j].col ] * iaData->constraints[i].M[j].val );
00252     }
00253     if (valid_solution())
00254       printf(" = %f", g_i);
00255     printf(" in [%f, %f]", iaData->constraints[i].lowerBound, iaData->constraints[i].upperBound);
00256     if (!satisfied)
00257       printf("  VIOLATED");
00258     printf("\n");
00259   }
00260   return std::make_pair(satisfied, g_i);
00261 }
00262   
00263 std::pair< bool, double> IASolverTool::even_constraint(const int i, const bool print_me, const bool print_unsatisfied ) const
00264 {
00265   const double g_i = even_value(i);
00266   // epsilon values
00267   int g_even;
00268   bool satisfied = is_even(g_i, g_even);
00269   if (print_me || (print_unsatisfied && !satisfied))
00270   {
00271     printf("even constraint %d:", i);
00272     for (unsigned int j = 0; j < iaData->sumEvenConstraints[i].M.size(); ++j)
00273     {
00274       printf( " %f x_%d", iaData->sumEvenConstraints[i].M[j].val, iaData->sumEvenConstraints[i].M[j].col );
00275       if (valid_solution())
00276         printf(" (%f) ", 
00277                iaSolution->x_solution[ iaData->sumEvenConstraints[i].M[j].col ] 
00278                * iaData->sumEvenConstraints[i].M[j].val );
00279     }
00280     if (valid_solution())
00281     {
00282       if (satisfied)
00283         printf(" = %f = %d", g_i, g_even);
00284       else
00285         printf(" = %f != %d NOT-EVEN", g_i, g_even);
00286     }
00287     printf("\n");
00288   }
00289   return std::make_pair(satisfied, g_i);
00290 }
00291   
00292   
00293 bool IASolverTool::equal_constraints( bool print_me, bool print_unsatisfied ) const
00294 {
00295   bool equal_satisfied = true;
00296   for (unsigned int i = 0; i<iaData->constraints.size(); ++i)
00297   {
00298     if (! equal_constraint(i, print_me, print_unsatisfied).first )
00299       equal_satisfied = false;
00300   }
00301   return equal_satisfied;
00302 }
00303   
00304 bool IASolverTool::even_constraints( bool print_me, bool print_unsatisfied ) const
00305 {
00306   bool even_satisfied = true;
00307   for (unsigned int i = 0; i<iaData->sumEvenConstraints.size(); ++i)
00308   {
00309     if (! even_constraint(i, print_me, print_unsatisfied).first )
00310       even_satisfied = false;
00311   }
00312   return even_satisfied;
00313 }
00314   
00315 bool IASolverTool::all_constraints( bool print_me, bool print_unsatisfied ) const
00316 {
00317   bool equal_satisfied = equal_constraints( print_me, print_unsatisfied );
00318   bool even_satisfied = even_constraints( print_me, print_unsatisfied );
00319   return equal_satisfied && even_satisfied;
00320 }
00321 
00322 
00323 } // namespace MeshKit
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines