MOAB: Mesh Oriented datABase
(version 5.4.1)
|
Optimizes the objective function using the Polack-Ribiere scheme. More...
#include <ConjugateGradient.hpp>
Public Member Functions | |
MESQUITE_EXPORT | ConjugateGradient (ObjectiveFunction *objective) |
MESQUITE_EXPORT | ConjugateGradient (ObjectiveFunction *objective, MsqError &err) |
virtual MESQUITE_EXPORT | ~ConjugateGradient () |
virtual MESQUITE_EXPORT std::string | get_name () const |
Get string name for use in diagnostic and status output. | |
virtual PatchSet * | get_patch_set () |
MESQUITE_EXPORT void | set_debugging_level (int new_lev) |
Protected Member Functions | |
virtual void | initialize (PatchData &pd, MsqError &err) |
Initialize data for smoothing process. | |
virtual void | optimize_vertex_positions (PatchData &pd, MsqError &err) |
virtual void | initialize_mesh_iteration (PatchData &pd, MsqError &err) |
virtual void | terminate_mesh_iteration (PatchData &pd, MsqError &err) |
virtual void | cleanup () |
Delete arrays initially created in initialize(). | |
double | get_step (PatchData &pd, double f0, int &j, MsqError &err) |
Returns the step distance to take in the search direction. | |
Private Attributes | |
std::vector< Vector3D > | fGrad |
Culls the vertex list free_vertex_list. | |
std::vector< Vector3D > | pGrad |
std::vector< Vector3D > | fNewGrad |
PatchDataVerticesMemento * | pMemento |
int | conjGradDebug |
Optimizes the objective function using the Polack-Ribiere scheme.
Definition at line 51 of file ConjugateGradient.hpp.
Definition at line 60 of file ConjugateGradient.cpp.
: VertexMover( objective ), PatchSetUser( true ), pMemento( NULL ), conjGradDebug( 0 ) { }
MBMesquite::ConjugateGradient::ConjugateGradient | ( | ObjectiveFunction * | objective, |
MsqError & | err | ||
) |
Definition at line 65 of file ConjugateGradient.cpp.
References MBMesquite::TerminationCriterion::add_iteration_limit(), MBMesquite::QualityImprover::get_inner_termination_criterion(), MBMesquite::MsqError::INVALID_STATE, MSQ_ERRRTN, MSQ_SETERR, and set_debugging_level().
: VertexMover( objective ), PatchSetUser( true ), pMemento( NULL ), conjGradDebug( 0 ) { // Michael:: default to global? set_debugging_level( 0 ); // set the default inner termination criterion TerminationCriterion* default_crit = get_inner_termination_criterion(); if( default_crit == NULL ) { MSQ_SETERR( err ) ( "QualityImprover did not create a default inner " "termination criterion.", MsqError::INVALID_STATE ); return; } else { default_crit->add_iteration_limit( 5 );MSQ_ERRRTN( err ); } }
MBMesquite::ConjugateGradient::~ConjugateGradient | ( | ) | [virtual] |
Definition at line 86 of file ConjugateGradient.cpp.
References pMemento.
{ // Checks that cleanup() has been called. assert( pMemento == NULL ); }
void MBMesquite::ConjugateGradient::cleanup | ( | ) | [protected, virtual] |
Delete arrays initially created in initialize().
Implements MBMesquite::VertexMover.
Definition at line 284 of file ConjugateGradient.cpp.
std::string MBMesquite::ConjugateGradient::get_name | ( | ) | const [virtual] |
Get string name for use in diagnostic and status output.
Implements MBMesquite::Instruction.
Definition at line 50 of file ConjugateGradient.cpp.
{ return "ConjugateGradient"; }
PatchSet * MBMesquite::ConjugateGradient::get_patch_set | ( | ) | [virtual] |
Reimplemented from MBMesquite::PatchSetUser.
Definition at line 55 of file ConjugateGradient.cpp.
{ return PatchSetUser::get_patch_set(); }
double MBMesquite::ConjugateGradient::get_step | ( | PatchData & | pd, |
double | f0, | ||
int & | j, | ||
MsqError & | err | ||
) | [protected] |
Returns the step distance to take in the search direction.
Computes a distance to move vertices given an initial position and search direction (stored in data member pGrad).
Returns alp, the double which scales the search direction vector which when added to the old nodal positions yields the new nodal positions.
Definition at line 304 of file ConjugateGradient.cpp.
References MBMesquite::arrptr(), MBMesquite::MsqError::BARRIER_VIOLATED, MBMesquite::MsqError::clear(), MBMesquite::MsqError::error_code(), MBMesquite::OFEvaluator::evaluate(), MBMesquite::VertexMover::get_objective_function_evaluator(), MBMesquite::MsqError::INVALID_MESH, MSQ_ERRZERO, MBMesquite::MSQ_MIN, MSQ_PRINT, MSQ_SETERR, MBMesquite::PatchData::num_free_vertices(), pGrad, pMemento, MBMesquite::PatchData::recreate_vertices_memento(), MBMesquite::PatchData::set_free_vertices_constrained(), and MBMesquite::PatchData::set_to_vertices_memento().
Referenced by optimize_vertex_positions().
{ // get OF evaluator OFEvaluator& objFunc = get_objective_function_evaluator(); size_t num_vertices = pd.num_free_vertices(); // initial guess for alp double alp = 1.0; int jmax = 100; double rho = 0.5; // feasible=false implies the mesh is not in the feasible region bool feasible = false; int found = 0; // f and fnew hold the objective function value double f = 0; double fnew = 0; // Counter to avoid infinitly scaling alp j = 0; // save memento pd.recreate_vertices_memento( pMemento, err ); // if we must check feasiblility // while step takes mesh into infeasible region and ... while( j < jmax && !feasible && alp > MSQ_MIN ) { ++j; pd.set_free_vertices_constrained( pMemento, arrptr( pGrad ), num_vertices, alp, err ); feasible = objFunc.evaluate( pd, f, err ); if( err.error_code() == err.BARRIER_VIOLATED ) err.clear(); // barrier violation does not represent an actual error here MSQ_ERRZERO( err ); // if not feasible, try a smaller alp (take smaller step) if( !feasible ) { alp *= rho; } } // end while ... // if above while ended due to j>=jmax, no valid step was found. if( j >= jmax ) { MSQ_PRINT( 2 )( "\nFeasible Point Not Found" ); return 0.0; } // Message::print_info("\nOriginal f %f, first new f = %f, alp = %f",f0,f,alp); // if new f is larger than original, our step was too large if( f >= f0 ) { j = 0; while( j < jmax && found == 0 ) { ++j; alp *= rho; pd.set_free_vertices_constrained( pMemento, arrptr( pGrad ), num_vertices, alp, err ); // Get new obj value // if patch is now invalid, then the feasible region is convex or // we have an error. For now, we assume an error. if( !objFunc.evaluate( pd, f, err ) ) { MSQ_SETERR( err ) ( "Non-convex feasiblility region found.", MsqError::INVALID_MESH ); } pd.set_to_vertices_memento( pMemento, err ); MSQ_ERRZERO( err ); // if our step has now improved the objective function value if( f < f0 ) { found = 1; } } // end while j less than jmax // Message::print_info("\nj = %d found = %d f = %20.18f f0 = %20.18f\n",j,found,f,f0); // if above ended because of j>=jmax, take no step if( found == 0 ) { // Message::print_info("alp = %10.8f, but returning zero\n",alp); alp = 0.0; return alp; } j = 0; // while shrinking the step improves the objFunc value further, // scale alp down. Return alp, when scaling once more would // no longer improve the objFunc value. while( j < jmax ) { ++j; alp *= rho; // step alp in search direction from original positions pd.set_free_vertices_constrained( pMemento, arrptr( pGrad ), num_vertices, alp, err ); MSQ_ERRZERO( err ); // get new objective function value if( !objFunc.evaluate( pd, fnew, err ) ) MSQ_SETERR( err ) ( "Non-convex feasiblility region found while " "computing new f.", MsqError::INVALID_MESH ); if( fnew < f ) { f = fnew; } else { // Reset the vertices to original position pd.set_to_vertices_memento( pMemento, err ); MSQ_ERRZERO( err ); alp /= rho; return alp; } } // Reset the vertices to original position and return alp pd.set_to_vertices_memento( pMemento, err ); MSQ_ERRZERO( err ); return alp; } // else our new f was already smaller than our original else { j = 0; // check to see how large of step we can take while( j < jmax && found == 0 ) { ++j; // scale alp up (rho must be less than 1) alp /= rho; // step alp in search direction from original positions pd.set_free_vertices_constrained( pMemento, arrptr( pGrad ), num_vertices, alp, err ); MSQ_ERRZERO( err ); feasible = objFunc.evaluate( pd, fnew, err ); if( err.error_code() == err.BARRIER_VIOLATED ) err.clear(); // evaluate() error does not represent an actual problem here MSQ_ERRZERO( err ); if( !feasible ) { alp *= rho; // Reset the vertices to original position and return alp pd.set_to_vertices_memento( pMemento, err ); MSQ_ERRZERO( err ); return alp; } if( fnew < f ) { f = fnew; } else { found = 1; alp *= rho; } } // Reset the vertices to original position and return alp pd.set_to_vertices_memento( pMemento, err ); MSQ_ERRZERO( err ); return alp; } }
void MBMesquite::ConjugateGradient::initialize | ( | PatchData & | pd, |
MsqError & | err | ||
) | [protected, virtual] |
Initialize data for smoothing process.
Implements MBMesquite::VertexMover.
Definition at line 92 of file ConjugateGradient.cpp.
References MBMesquite::PatchData::create_vertices_memento(), MBMesquite::get_parallel_rank(), MBMesquite::get_parallel_size(), MSQ_DBGOUT, and pMemento.
{ if( get_parallel_size() ) { MSQ_DBGOUT( 2 ) << "\nP[" << get_parallel_rank() << "] " << "o Performing Conjugate Gradient optimization.\n"; } else { MSQ_DBGOUT( 2 ) << "\no Performing Conjugate Gradient optimization.\n"; } pMemento = pd.create_vertices_memento( err ); }
void MBMesquite::ConjugateGradient::initialize_mesh_iteration | ( | PatchData & | pd, |
MsqError & | err | ||
) | [protected, virtual] |
void MBMesquite::ConjugateGradient::optimize_vertex_positions | ( | PatchData & | pd, |
MsqError & | err | ||
) | [protected, virtual] |
Performs Conjugate gradient minimization on the PatchData, pd.
Implements MBMesquite::VertexMover.
Definition at line 109 of file ConjugateGradient.cpp.
References MBMesquite::TerminationCriterion::accumulate_inner(), MBMesquite::TerminationCriterion::accumulate_patch(), MBMesquite::arrptr(), conjGradDebug, MBMesquite::divide(), fGrad, fNewGrad, MBMesquite::QualityImprover::get_inner_termination_criterion(), MBMesquite::VertexMover::get_objective_function_evaluator(), get_step(), MBMesquite::MsqError::INVALID_MESH, MBMesquite::Linf(), MBMesquite::PatchData::move_free_vertices_constrained(), MSQ_CHKERR, MSQ_DBGOUT, MSQ_ERRRTN, MSQ_FUNCTION_TIMER, MBMesquite::MSQ_MAX_CAP, MSQ_PRINT, MSQ_SETERR, MBMesquite::PatchData::num_free_vertices(), pGrad, MBMesquite::Timer::since_birth(), MBMesquite::TerminationCriterion::terminate(), and MBMesquite::OFEvaluator::update().
{ // pd.reorder(); MSQ_FUNCTION_TIMER( "ConjugateGradient::optimize_vertex_positions" ); Timer c_timer; size_t num_vert = pd.num_free_vertices(); if( num_vert < 1 ) { MSQ_DBGOUT( 1 ) << "\nEmpty free vertex list in ConjugateGradient\n"; return; } /* //zero out arrays int zero_loop=0; while(zero_loop<arraySize){ fGrad[zero_loop].set(0,0,0); pGrad[zero_loop].set(0,0,0); fNewGrad[zero_loop].set(0,0,0); ++zero_loop; } */ // get OF evaluator OFEvaluator& objFunc = get_objective_function_evaluator(); size_t ind; // Michael cull list: possibly set soft_fixed flags here // MsqFreeVertexIndexIterator free_iter(pd, err); MSQ_ERRRTN(err); double f = 0; // Michael, this isn't equivalent to CUBIT because we only want to check // the objective function value of the 'bad' elements // if invalid initial patch set an error. bool temp_bool = objFunc.update( pd, f, fGrad, err ); assert( fGrad.size() == num_vert ); if( MSQ_CHKERR( err ) ) return; if( !temp_bool ) { MSQ_SETERR( err ) ( "Conjugate Gradient not able to get valid gradient " "and function values on intial patch.", MsqError::INVALID_MESH ); return; } double grad_norm = MSQ_MAX_CAP; if( conjGradDebug > 0 ) { MSQ_PRINT( 2 )( "\nCG's DEGUB LEVEL = %i \n", conjGradDebug ); grad_norm = Linf( arrptr( fGrad ), fGrad.size() ); MSQ_PRINT( 2 )( "\nCG's FIRST VALUE = %f,grad_norm = %f", f, grad_norm ); MSQ_PRINT( 2 )( "\n TIME %f", c_timer.since_birth() ); grad_norm = MSQ_MAX_CAP; } // Initializing pGrad (search direction). pGrad.resize( fGrad.size() ); for( ind = 0; ind < num_vert; ++ind ) pGrad[ind] = ( -fGrad[ind] ); int j = 0; // total nb of step size changes ... not used much int i = 0; // iteration counter unsigned m = 0; // double alp = MSQ_MAX_CAP; // alp: scale factor of search direction // we know inner_criterion is false because it was checked in // loop_over_mesh before being sent here. TerminationCriterion* term_crit = get_inner_termination_criterion(); // while ((i<maxIteration && alp>stepBound && grad_norm>normGradientBound) // && !inner_criterion){ while( !term_crit->terminate() ) { ++i; // std::cout<<"\Michael delete i = "<<i; int k = 0; alp = get_step( pd, f, k, err ); j += k; if( conjGradDebug > 2 ) { MSQ_PRINT( 2 )( "\n Alp initial, alp = %20.18f", alp ); } // if alp == 0, revert to steepest descent search direction if( alp == 0 ) { for( m = 0; m < num_vert; ++m ) { pGrad[m] = ( -fGrad[m] ); } alp = get_step( pd, f, k, err ); j += k; if( conjGradDebug > 1 ) { MSQ_PRINT( 2 )( "\n CG's search direction reset." ); if( conjGradDebug > 2 ) MSQ_PRINT( 2 )( "\n Alp was zero, alp = %20.18f", alp ); } } if( alp != 0 ) { pd.move_free_vertices_constrained( arrptr( pGrad ), num_vert, alp, err );MSQ_ERRRTN( err ); if( !objFunc.update( pd, f, fNewGrad, err ) ) { MSQ_SETERR( err ) ( "Error inside Conjugate Gradient, vertices moved " "making function value invalid.", MsqError::INVALID_MESH ); return; } assert( fNewGrad.size() == (unsigned)num_vert ); if( conjGradDebug > 0 ) { grad_norm = Linf( arrptr( fNewGrad ), num_vert ); MSQ_PRINT( 2 ) ( "\nCG's VALUE = %f, iter. = %i, grad_norm = %f, alp = %f", f, i, grad_norm, alp ); MSQ_PRINT( 2 )( "\n TIME %f", c_timer.since_birth() ); } double s11 = 0; double s12 = 0; double s22 = 0; // free_iter.reset(); // while (free_iter.next()) { // m=free_iter.value(); for( m = 0; m < num_vert; ++m ) { s11 += fGrad[m] % fGrad[m]; s12 += fGrad[m] % fNewGrad[m]; s22 += fNewGrad[m] % fNewGrad[m]; } // Steepest Descent (takes 2-3 times as long as P-R) // double bet=0; // Fletcher-Reeves (takes twice as long as P-R) // double bet = s22/s11; // Polack-Ribiere double bet; if( !divide( s22 - s12, s11, bet ) ) return; // gradient is zero // free_iter.reset(); // while (free_iter.next()) { // m=free_iter.value(); for( m = 0; m < num_vert; ++m ) { pGrad[m] = ( -fNewGrad[m] + ( bet * pGrad[m] ) ); fGrad[m] = fNewGrad[m]; } if( conjGradDebug > 2 ) { MSQ_PRINT( 2 ) ( " \nSEARCH DIRECTION INFINITY NORM = %e", Linf( arrptr( fNewGrad ), num_vert ) ); } } // end if on alp == 0 term_crit->accumulate_patch( pd, err );MSQ_ERRRTN( err ); term_crit->accumulate_inner( pd, f, arrptr( fGrad ), err );MSQ_ERRRTN( err ); } // end while if( conjGradDebug > 0 ) { MSQ_PRINT( 2 )( "\nConjugate Gradient complete i=%i ", i ); MSQ_PRINT( 2 )( "\n- FINAL value = %f, alp=%4.2e grad_norm=%4.2e", f, alp, grad_norm ); MSQ_PRINT( 2 )( "\n FINAL TIME %f", c_timer.since_birth() ); } }
MESQUITE_EXPORT void MBMesquite::ConjugateGradient::set_debugging_level | ( | int | new_lev | ) | [inline] |
Just for debugging purposes or for obtaining more data during the optimization process.
Definition at line 65 of file ConjugateGradient.hpp.
References conjGradDebug.
Referenced by ConjugateGradient(), run_solution_mesh_optimizer(), SphericalGeometryTest::test_cg_mesh_cond_sphere(), and PlanarGeometryTest::test_plane_tri_xz().
{ conjGradDebug = new_lev; }
void MBMesquite::ConjugateGradient::terminate_mesh_iteration | ( | PatchData & | pd, |
MsqError & | err | ||
) | [protected, virtual] |
Implements MBMesquite::VertexMover.
Definition at line 279 of file ConjugateGradient.cpp.
{
// cout << "- Executing ConjugateGradient::iteration_complete()\n";
}
int MBMesquite::ConjugateGradient::conjGradDebug [private] |
Definition at line 93 of file ConjugateGradient.hpp.
Referenced by optimize_vertex_positions(), and set_debugging_level().
std::vector< Vector3D > MBMesquite::ConjugateGradient::fGrad [private] |
Culls the vertex list free_vertex_list.
Definition at line 90 of file ConjugateGradient.hpp.
Referenced by cleanup(), and optimize_vertex_positions().
std::vector< Vector3D > MBMesquite::ConjugateGradient::fNewGrad [private] |
Definition at line 90 of file ConjugateGradient.hpp.
Referenced by cleanup(), and optimize_vertex_positions().
std::vector< Vector3D > MBMesquite::ConjugateGradient::pGrad [private] |
Definition at line 90 of file ConjugateGradient.hpp.
Referenced by cleanup(), get_step(), and optimize_vertex_positions().
Definition at line 91 of file ConjugateGradient.hpp.
Referenced by cleanup(), get_step(), initialize(), and ~ConjugateGradient().