MeshKit  1.0
IAInterface.cpp
Go to the documentation of this file.
00001 // IAInterface.cpp MeshKit version
00002 // Interval Assignment Data for Meshkit
00003 // Interface to the rest of Meshkit
00004 
00005 #include "meshkit/IAInterface.hpp"
00006 #include "meshkit/IASolver.hpp"
00007 #include "meshkit/IADataBuilder.hpp"
00008 #include "meshkit/IASolution.hpp"
00009 
00010 #include "meshkit/ModelEnt.hpp" // from MeshKit
00011 
00012 #include <vector>
00013 #include <iterator>
00014 #include <set>
00015 #include <cstdio>
00016 #include <assert.h>
00017 #include <math.h>
00018 
00019 namespace MeshKit 
00020 {
00021   
00022 const bool IAInterface::debugging = false;
00023   
00024 //static IAVariable counter
00025 unsigned int IAVariable::numVariables(0);
00026 
00027 moab::EntityType IAInterface_types[] = {moab::MBMAXTYPE};
00028 const moab::EntityType* IAInterface::output_types() 
00029   { return IAInterface_types; }
00030     
00031 void IAInterface::setup_this()
00032 {
00033   ; //nothing for now
00034 }
00035     
00036 IAVariable *IAInterface::get_variable( ModelEnt* model_entity, bool create_if_missing )
00037 {
00038   // already exists?
00039   if (model_entity && model_entity->ia_var())
00040   {
00041     // return without setting values
00042     return model_entity->ia_var();
00043   }
00044   
00045   if (model_entity && create_if_missing) 
00046   {
00047     // default values
00048     double goal(1.);
00049     Firmness firm(DEFAULT);
00050     
00051     const Firmness me_firm( model_entity->interval_firmness() );
00052     const double me_edge_size( model_entity->mesh_interval_size() );
00053     if ( me_firm == DEFAULT || me_edge_size <= 0.) 
00054     {
00055       // ModelEnt was unset, so use default values
00056       ;
00057     }
00058     else
00059     {
00060       const double me_size = model_entity->measure();
00061       assert(me_edge_size > 0.);
00062       goal = me_size / me_edge_size;
00063       firm = model_entity->interval_firmness(); 
00064     }
00065     return create_variable( model_entity, firm, goal);  
00066   }
00067   return NULL;
00068 }
00069 
00070 
00071 IAVariable *IAInterface::create_variable( ModelEnt* model_entity, IAVariable::Firmness set_firm, double goal_value)
00072 {
00073   IAVariable *ia_var = NULL;
00074 
00075   // already exists?
00076   if (model_entity && model_entity->ia_var())
00077   {
00078     // set new values and return
00079     ia_var = model_entity->ia_var();
00080     ia_var->set_firmness(set_firm);
00081     ia_var->set_goal(goal_value);
00082     return ia_var;
00083   }
00084       
00085   ia_var = new IAVariable(model_entity, set_firm, goal_value);
00086   assert( ia_var->uniqueId + 1 == IAVariable::numVariables );
00087   variables.push_back(ia_var);
00088   assert( IAVariable::numVariables == variables.size() );
00089 
00090   return ia_var;
00091 }
00092 
00093 IAInterface::IAVariableVec IAInterface::make_constraint_group( const MEntVector &model_entity_vec )
00094 {
00095   IAVariableVec result( model_entity_vec.size() );
00096   for (unsigned int i = 0; i < model_entity_vec.size(); ++i)
00097   {
00098     assert(model_entity_vec[i]);
00099     IAVariable *v = get_variable( model_entity_vec[i], true );
00100     result[i] = v;
00101   }
00102   return result;
00103 }
00104 
00105 void IAInterface::destroy_variable( IAVariable* ia_variable )
00106 {
00107   if (!ia_variable)
00108     return;
00109   // model_entity shouldn't point to ia_variable anymore
00110   ModelEnt *me = ia_variable->get_model_ent();
00111   if (me && me->ia_var() == ia_variable)
00112     me->ia_var(NULL);
00113   
00114   variables[ ia_variable->uniqueId ] = NULL;
00115   delete ia_variable;    
00116 }
00117 
00118 
00119 void IAInterface::destroy_data()
00120 {
00121   // destroy remaining variables
00122   for (unsigned int i = 0; i < variables.size(); ++i)
00123   {
00124     if ( variables[i] )
00125     {
00126       assert( variables[i]->uniqueId == i);
00127       destroy_variable( variables[i] );
00128     }
00129   }
00130   variables.clear();
00131 
00132   // reset global variable id? Yes if this is a singleton class, which it is right now.
00133   IAVariable::numVariables = 0;
00134 
00135   // clear constraints
00136   sumEqualConstraints1.clear();
00137   sumEqualConstraints2.clear();
00138   sumEvenConstraints.clear();
00139   // ... additional types of constraints ...
00140 }
00141 
00142 IAInterface::~IAInterface()
00143 {
00144   destroy_data();
00145 }
00146 
00147 void IAInterface::constrain_sum_even( const IAVariableVec &sum_even_vars )
00148 {
00149   sumEvenConstraints.push_back( sum_even_vars ); // vector copy
00150 }
00151 
00152 void IAInterface::constrain_sum_equal( const IAVariableVec &side_one, const IAVariableVec &side_two )
00153 {
00154   sumEqualConstraints1.push_back(side_one); //vector copy of side_one
00155   sumEqualConstraints2.push_back(side_two); //vector copy of side_two
00156 }
00157 
00158 
00159 // for global vector of variables
00160 int IAInterface::variable_to_index(const IAVariable* var) const
00161 {
00162   return var->uniqueId;
00163 }
00164 
00165 IAVariable *IAInterface::index_to_variable(int ind) const
00166 {
00167   assert(ind < (int) variables.size());
00168   return variables[ind];
00169 }
00170 
00171 void IAInterface::make_set_0_to_nm1( IndexSet &index_set, const size_t k )
00172 {
00173   std::pair< IndexSet::iterator, bool> prior =  index_set.insert( 0 );
00174   assert(prior.second);
00175   for ( size_t i = 0; i < k; ++i)
00176   {
00177     prior.first = index_set.insert( prior.first, (int) i );
00178     // assert(prior.second); // second = insert-success is not assigned by this version of insert!
00179   }
00180 }
00181 
00182 void IAInterface::make_vec_0_to_nm1( IndexVec &index_vec, const size_t k )
00183 {
00184   index_vec.clear();
00185   index_vec.reserve(k);
00186   for (size_t i = 0; i < k; ++i)
00187     index_vec.push_back( (int) i);
00188 }
00189 
00190 void IAInterface::make_vec_unset( IndexVec &index_vec, const size_t k )
00191 {
00192   index_vec.clear();
00193   index_vec.reserve(k);
00194   for (int i = 0; i < (int) k; ++i)
00195     index_vec.push_back(-1);
00196 }
00197 
00198 void IAInterface::set_variable_constraint_indices(
00199             VariableConstraintDependencies &var_con_dep, 
00200             const int i_start, 
00201             const VariableVecVec &variable_vec_vec )
00202 {
00203   unsigned int j;
00204   int i;
00205   for (j = 0, i = i_start; j < variable_vec_vec.size(); ++j, ++i)
00206   {
00207     for (unsigned int k = 0; k < variable_vec_vec[j].size(); ++k )
00208     {
00209       const IAVariable *var = variable_vec_vec[j][k];
00210       const int v = variable_to_index( var ); 
00211       var_con_dep.constraintVariables[i].push_back( v );        
00212       assert( v < (int) var_con_dep.variableConstraints.size() );
00213       var_con_dep.variableConstraints[v].push_back( i );
00214     }
00215   }
00216 }
00217    
00218 void IAInterface::set_variable_constraint_indices(VariableConstraintDependencies &var_con_dep)
00219 {
00220   // create vectors of the length of constraints and variables, but containing of empty vectors
00221   var_con_dep.constraintVariables.resize( sumEqualConstraints1.size() + sumEvenConstraints.size() );
00222   var_con_dep.variableConstraints.resize( variables.size() );
00223 
00224   assert( sumEqualConstraints1.size() == sumEqualConstraints2.size() );
00225   const int sum_even_start = (int) sumEqualConstraints1.size();
00226   //zzyk hardsets
00227   // build vector map of variable-to-constraints
00228   // build vector map of constraint-to-variables
00229   set_variable_constraint_indices( var_con_dep, 0, sumEqualConstraints1 );
00230   set_variable_constraint_indices( var_con_dep, 0, sumEqualConstraints2 );
00231   set_variable_constraint_indices( var_con_dep, sum_even_start, sumEvenConstraints );
00232 }
00233 
00234 void IAInterface::ProblemSets::find_variable_dependent_set( 
00235          const int variable_j, 
00236          const VariableConstraintDependencies &var_con_dep,
00237          ProblemSets &subsets )
00238 {
00239   // return if we've already added this variable to the subset = removed it from this set
00240   IndexSet::iterator j = variableSet.find(variable_j);
00241   if ( j == variableSet.end() )
00242     return;
00243   
00244   // add the variable to the sub-problem
00245   subsets.variableSet.insert(*j);
00246   
00247   // remove the variable from the big problem
00248   variableSet.erase(j);
00249   
00250   // recursively find the dependent constraints
00251   IndexVec::const_iterator i;
00252   assert( variable_j < (int) var_con_dep.variableConstraints.size() ); 
00253   const IndexVec &v = var_con_dep.variableConstraints[variable_j];
00254   for ( i = v.begin(); i != v.end(); ++i )
00255     find_constraint_dependent_set( *i, var_con_dep, subsets );
00256 }
00257 
00258 
00259 void IAInterface::ProblemSets::find_constraint_dependent_set( 
00260          const int constraint_i, 
00261          const VariableConstraintDependencies &var_con_dep,
00262          ProblemSets &subsets )
00263 {
00264   // return if we've already added this constraint to the subset = removed it from this set
00265   // todo, could speed this up using some sort of mark on the entity
00266   IndexSet::iterator i = constraintSet.find(constraint_i);
00267   if ( i == constraintSet.end() )
00268     return;
00269   
00270   // add the constraint to the sub-problem
00271   subsets.constraintSet.insert(*i);
00272 
00273   // remove the constraint from the global problem
00274   constraintSet.erase(i);
00275 
00276   // recursively find the dependent variables
00277   IndexVec::const_iterator j;
00278   assert( constraint_i < (int) var_con_dep.constraintVariables.size() );
00279   const IndexVec &v = var_con_dep.constraintVariables[constraint_i];
00280   for ( j = v.begin(); j != v.end(); ++j )
00281     find_variable_dependent_set( *j, var_con_dep, subsets );
00282 }
00283 
00284 void IAInterface::ProblemSets::find_dependent_set( 
00285   const VariableConstraintDependencies &var_con_dep,
00286   ProblemSets &subsets )
00287 {
00288   if ( !constraintSet.empty() )
00289   {
00290     find_constraint_dependent_set( *constraintSet.begin(), var_con_dep, subsets );
00291   }
00292   else if ( !variableSet.empty() )
00293   {
00294     find_variable_dependent_set( *variableSet.begin(), var_con_dep, subsets );
00295   }
00296 }
00297    
00298 void IAInterface::make_global_set(ProblemSets &problem_sets)
00299 {
00300   make_set_0_to_nm1( problem_sets.constraintSet, sumEqualConstraints1.size() + sumEvenConstraints.size() );
00301   make_set_0_to_nm1( problem_sets.variableSet, variables.size() ); 
00302   make_vec_unset( problem_sets.varMap, variables.size() );
00303   // todo hardsets, 
00304   // move variables from variableSet to hard_variableSet if the variable is hard-set
00305 }
00306   
00307 void IAInterface::global_to_sub_side( 
00308   const VariableVec &global_constraint, 
00309   IndexVec &global_var_map,
00310   IndexVec &local_constraint, int &rhs ) const
00311 {
00312   for (unsigned int j = 0; j < global_constraint.size(); ++j )
00313   {
00314     IAVariable *v = global_constraint[j];
00315     if (v->get_firmness() == MeshKit::HARD )
00316       rhs -= floor( v->get_goal() + 0.5 ); // goal should be integer already for hardsets
00317     else
00318     {
00319       int local_index = global_var_map[ v->uniqueId ];
00320       assert( local_index >= 0 );
00321       // assert( local_index < variableSet.size() );
00322       local_constraint.push_back( local_index ); //variable_to_index( v, sub_variables ) );
00323     }
00324   }
00325 }
00326   
00327 void IAInterface::ProblemSets::print() const
00328 {
00329   printf("ProblemSets %p:\n", this);
00330   if (!empty())
00331     printf("non-");
00332   printf("empty. %lu variables, %lu hardsets, %lu constraints\n", variableSet.size(), hardVariableSet.size(), constraintSet.size());
00333   printf("Variables: ");
00334   for ( IndexSet::const_iterator i = variableSet.begin(); i != variableSet.end(); ++i )
00335   {
00336     printf("%d ", *i);
00337   }
00338   printf("\n");
00339   printf("HardSet Variables: ");
00340   for ( IndexSet::const_iterator i = hardVariableSet.begin(); i != hardVariableSet.end(); ++i )
00341   {
00342     printf("%d ", *i);
00343   }
00344   printf("\n");  
00345   printf("Constraints: ");
00346   for ( IndexSet::const_iterator i = constraintSet.begin(); i != constraintSet.end(); ++i )
00347   {
00348     printf("%d ", *i);
00349   }
00350   printf("\n\n");
00351 }
00352 
00353 void IAInterface::VariableConstraintDependencies::print() const
00354 {
00355   printf("VariableConstraintDependencies %p:\n",this);
00356   printf("%lu variables, %lu constraints\n", variableConstraints.size(), constraintVariables.size() );
00357   for (size_t i = 0; i < variableConstraints.size(); ++i)
00358   {
00359     printf("Variable %lu depends on constraints ", i);
00360     for (size_t j = 0; j < variableConstraints[i].size(); ++j)
00361     {
00362       printf("%d ", variableConstraints[i][j]);
00363     }
00364     printf("\n");
00365   }
00366   
00367   for (size_t i = 0; i < constraintVariables.size(); ++i)
00368   {
00369     printf("Constraint %lu depends on variables ", i);
00370     for (size_t j = 0; j < constraintVariables[i].size(); ++j)
00371     {
00372       printf("%d ", constraintVariables[i][j]);
00373     }
00374     printf("\n");
00375   }
00376   printf("\n");
00377 }
00378 
00379   
00380 void IAInterface::fill_problem( ProblemSets &sub_sets, 
00381   SubProblem *sub_problem, IndexVec &global_var_map ) const
00382 {
00383   IADataBuilder data_builder( sub_problem->data ); 
00384   // number the variables from 0..n, by setting ia_index
00385   sub_problem->data->I.reserve(sub_sets.variableSet.size());
00386   sub_sets.varMap.reserve(sub_sets.variableSet.size());
00387   for (IndexSet::const_iterator i = sub_sets.variableSet.begin(); i != sub_sets.variableSet.end(); ++i)
00388   {
00389     size_t index = *i; // global index
00390     IAVariable *v = index_to_variable( (int) index );
00391     if ( v->get_firmness() != MeshKit::HARD )
00392     {
00393       // index maps
00394       global_var_map[index] = (int) sub_sets.varMap.size();
00395       sub_sets.varMap.push_back( (int) index );
00396 
00397       // add the goals
00398       data_builder.add_variable( v->goal );
00399     }
00400     else
00401     {
00402      // nothing to do, it is a constant, not a variable in the sub problem
00403     }
00404   }
00405 
00406   // convert constraints to IASolver format
00407   for (IndexSet::const_iterator i = sub_sets.constraintSet.begin(); i != sub_sets.constraintSet.end(); ++i )
00408   {
00409     unsigned int c = *i;
00410     if (c < sumEqualConstraints1.size()) // equal
00411     {
00412       IndexVec side_1, side_2;
00413       int rhs_1(0), rhs_2(0);
00414         
00415       global_to_sub_side( sumEqualConstraints1[c], global_var_map, side_1, rhs_1 );
00416       global_to_sub_side( sumEqualConstraints2[c], global_var_map, side_2, rhs_2 );
00417       
00418       data_builder.constrain_opposite_side_equal(side_1, side_2, rhs_1 - rhs_2 );   
00419       // todo check should be rhs_2 - rhs_1 ?
00420     }
00421     else // even 
00422     {
00423       c -= sumEqualConstraints1.size();
00424       assert( c >= 0 ); // not useful because it is unsigned
00425       assert( c < sumEvenConstraints.size() );
00426 
00427       IndexVec side;
00428       int rhs(0);
00429       global_to_sub_side( sumEvenConstraints[c], global_var_map, side, rhs ); // todo check if should be -rhs?
00430       data_builder.constrain_sum_even(side, rhs );   
00431     }
00432   }
00433 }
00434 
00435 void IAInterface::subdivide_problem(SubProblemVec &subproblems)
00436 {
00437   // find independent subproblems, by index of constraint and variable
00438   ProblemSets global_prob_set;
00439   make_global_set(global_prob_set);
00440   if (debugging)
00441     global_prob_set.print(); 
00442   VariableConstraintDependencies var_con_dep;
00443   set_variable_constraint_indices(var_con_dep);
00444   if (debugging)
00445     var_con_dep.print(); 
00446   
00447   while ( !global_prob_set.empty() )
00448   {
00449     // define a sub-problem by indices
00450     SubProblem *subproblem = new SubProblem;
00451     subproblems.push_back( subproblem );
00452     global_prob_set.find_dependent_set(var_con_dep, subproblem->problemSets);
00453     
00454     if (debugging)
00455     {
00456       printf("subproblem peeled out:\n");
00457       subproblem->print();
00458       printf("\nremaining global problem:\n");
00459       global_prob_set.print();
00460     }
00461     
00462     // Fill in solver with data from those indices    
00463     // convert the set into an IASolver, add data
00464     fill_problem( subproblem->problemSets, subproblem, global_prob_set.varMap );
00465     if (debugging)
00466       subproblem->print();
00467   }    
00468 }
00469 
00470 void IAInterface::subdivide_problem_one(std::vector<IAData*> &subproblems)
00471 {
00472    // placeholder - just make one subproblem
00473   IAData *sub_problem = new IAData();
00474   subproblems.push_back(sub_problem);
00475   IADataBuilder data_builder( sub_problem ); // data  
00476 
00477   // this grabs the whole problem and converts it
00478   // variables numbered 0..n
00479   for(VariableVec::const_iterator i = variables.begin(); i != variables.end(); ++i)
00480   {
00481     IAVariable *v = *i;
00482     data_builder.add_variable( v->goal );
00483   }
00484   // equal constraints
00485   for (unsigned int i = 0; i < sumEqualConstraints1.size(); ++i)
00486   {
00487     std::vector<int> side_1, side_2;
00488     for (unsigned int j = 0; j < sumEqualConstraints1[i].size(); ++j)
00489     {
00490       IAVariable *v = sumEqualConstraints1[i][j];
00491       int index = variable_to_index(v);
00492       // sanity check that the indexing was correct
00493       assert( sub_problem->I[index] == v->goal ); 
00494       side_1.push_back(index);
00495     }
00496     for (unsigned int j = 0; j < sumEqualConstraints2[i].size(); ++j)
00497     {
00498       IAVariable *v = sumEqualConstraints2[i][j];
00499       int index = variable_to_index(v);
00500       side_2.push_back(index);
00501     }
00502     data_builder.constrain_opposite_side_equal(side_1, side_2, 0.);
00503   }
00504   // even constraints
00505   for (unsigned int i = 0; i < sumEvenConstraints.size(); ++i)
00506   {
00507     std::vector<int> side;
00508     for (unsigned int j = 0; j < sumEvenConstraints[i].size(); ++j)
00509     {
00510       IAVariable *v = sumEvenConstraints[i][j];
00511       int index = variable_to_index(v);
00512       side.push_back(index);
00513     }
00514     data_builder.constrain_sum_even(side, 0.);
00515   }
00516 
00517 }
00518 
00519 void IAInterface::print_problem() const
00520 {
00521   // describe variables
00522   // describe constraints
00523 }
00524 
00525 void IAInterface::SubProblem::print() const
00526 {
00527   printf("subproblem %d:\n", id);
00528   problemSets.print();
00529   IASolverTool solver_tool(data,solution);
00530   solver_tool.print();
00531   // describe variables
00532   // describe constraints
00533 }
00534 
00535 bool IAInterface::solve_subproblem( SubProblem *subproblem )
00536 {
00537   IASolver solver( subproblem->data, subproblem->solution );
00538   return solver.solve();
00539 }
00540 
00541 void IAInterface::assign_solution( SubProblem *subproblem )
00542 {
00543   // assign solution value from subproblem to model entities
00544   std::vector<double> &x_solution = subproblem->solution->x_solution; // shorthand
00545   for (unsigned int i = 0; i < x_solution.size(); ++i)
00546   {
00547     const double x = x_solution[i];
00548     // map index i to IAVariable, which will be different when we have non-full subproblems
00549     IAVariable *v = index_to_variable( subproblem->problemSets.varMap[i] );
00550     assert(v);
00551     assert( x > 0. );
00552     const int x_int = (int) floor( x + 0.5 );
00553     v->solution = x_int;
00554     if (v->get_model_ent())
00555     {
00556       v->get_model_ent()->mesh_intervals(x_int);
00557     }
00558   }
00559 }
00560 
00561 void IAInterface::execute_this()
00562 {
00563   SubProblemVec subproblems;
00564   subdivide_problem(subproblems);
00565   for (unsigned int i = 0; i < subproblems.size(); ++i)
00566   {
00567     SubProblem *p = subproblems[i];
00568     if (debugging)
00569       printf("Solving subproblem %d\n", p->id);
00570     if ( solve_subproblem(p) )
00571       assign_solution(p);
00572     else
00573     {
00574       ; // some error statement
00575     }
00576     delete p;
00577   }
00578   subproblems.clear();
00579   return;
00580 }
00581 
00582 
00583 int IAInterface::SubProblem::max_id = 0;
00584 
00585 IAInterface::SubProblem::SubProblem()
00586   : id( max_id++ )
00587 {
00588   data = new IAData;
00589   solution = new IASolution;
00590 }
00591   
00592 IAInterface::SubProblem::~SubProblem()
00593 {
00594   delete data;
00595   data = NULL;
00596   delete solution;
00597   solution = NULL;
00598 }
00599 
00600 } // namespace MeshKit
00601 
00602 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines