MOAB: Mesh Oriented datABase
(version 5.3.1)
|
00001 /****************************************************************** 00002 MESQUITE -- The Mesh Quality Improvement Toolkit 00003 00004 Copyright 2004 Sandia Corporation and Argonne National 00005 Laboratory. Under the terms of Contract DE-AC04-94AL85000 00006 with Sandia Corporation, the U.S. Government retains certain 00007 rights in this software. 00008 00009 This library is free software; you can redistribute it and/or 00010 modify it under the terms of the GNU Lesser General Public 00011 License as published by the Free Software Foundation; either 00012 version 2.1 of the License, or (at your option) any later version. 00013 00014 This library is distributed in the hope that it will be useful, 00015 but WITHOUT ANY WARRANTY; without even the implied warranty of 00016 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 00017 Lesser General Public License for more details. 00018 00019 You should have received a copy of the GNU Lesser General Public License 00020 (lgpl.txt) along with this library; if not, write to the Free Software 00021 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 00022 00023 diachin2@llnl.gov, djmelan@sandia.gov, mbrewer@sandia.gov, 00024 pknupp@sandia.gov, tleurent@mcs.anl.gov, tmunson@mcs.anl.gov 00025 00026 ******************************************************************/ 00027 /*! 00028 \file NonGradient.cpp 00029 \brief 00030 The NonGradient class is also a concrete vertex mover 00031 which performs derivative free minimization 00032 based on the Amoeba Method, as implemented in Numerical Recipes in C. 00033 00034 The other optimization methods either accept or reject a point immediately. 00035 Amoeba may not accept a point for several iterations/evaluations. 00036 Amoeba does not check for BOUNDED_VERTEX_MOVERMENT. 00037 rtol=2.0*fabs( height[ihi]-height[ilo] )/ 00038 ( fabs(height[ihi])+fabs(height[ilo])+threshold ); 00039 00040 We are aware of the problem of choosing the initial simplex, but have not 00041 sovled it. How has this been solved? 00042 There is a problem of letting Mesquite handle convergence. The idea here is to 00043 call accumulate_? only if the best point changes. But does it matter? Where 00044 are the errors determined? What does it mean to cull? 00045 What does Pat... compare and contrast the differnt vertex movement criteria 00046 with the Amoeba criterion. 00047 */ 00048 00049 #include "NonGradient.hpp" 00050 #include "MsqFreeVertexIndexIterator.hpp" 00051 #include "MsqTimer.hpp" 00052 #include "MsqDebug.hpp" 00053 #include "MsqError.hpp" 00054 #include <cmath> 00055 #include <iostream> 00056 #include "../../ObjectiveFunction/ObjectiveFunction.hpp" 00057 #include <cfloat> 00058 00059 namespace MBMesquite 00060 { 00061 00062 std::string NonGradient::get_name() const 00063 { 00064 return "NonGradient"; 00065 } 00066 00067 PatchSet* NonGradient::get_patch_set() 00068 { 00069 return PatchSetUser::get_patch_set(); 00070 } 00071 00072 NonGradient::NonGradient( ObjectiveFunction* of ) 00073 : VertexMover( of ), PatchSetUser( true ), projectGradient( false ), mDimension( 0 ), mThreshold( 0.0 ), 00074 mTolerance( 0.0 ), mMaxNumEval( 0 ), mNonGradDebug( 0 ), mUseExactPenaltyFunction( true ), mScaleDiameter( 0.1 ) 00075 { 00076 set_debugging_level( 2 ); 00077 // set the default inner termination criterion 00078 TerminationCriterion* default_crit = get_inner_termination_criterion(); 00079 if( default_crit == NULL ) { return; } 00080 else 00081 { 00082 default_crit->add_iteration_limit( 5 ); 00083 } 00084 } 00085 00086 NonGradient::NonGradient( ObjectiveFunction* of, MsqError& err ) 00087 : VertexMover( of ), PatchSetUser( true ), projectGradient( false ), mDimension( 0 ), mThreshold( 0.0 ), 00088 mTolerance( 0.0 ), mMaxNumEval( 0 ), mNonGradDebug( 0 ), mUseExactPenaltyFunction( true ), mScaleDiameter( 0.1 ) 00089 { 00090 set_debugging_level( 2 ); 00091 // set the default inner termination criterion 00092 TerminationCriterion* default_crit = get_inner_termination_criterion(); 00093 if( default_crit == NULL ) 00094 { 00095 MSQ_SETERR( err ) 00096 ( "QualityImprover did not create a default inner " 00097 "termination criterion.", 00098 MsqError::INVALID_STATE ); 00099 return; 00100 } 00101 else 00102 { 00103 default_crit->add_iteration_limit( 5 );MSQ_ERRRTN( err ); 00104 } 00105 } 00106 00107 bool NonGradient::testRowSum( int numRow, int numCol, double* matrix, double* oldRowSum ) 00108 { 00109 bool same = true; 00110 std::vector< double > rowSum( numRow, 0. ); 00111 double maxVal = 0.; 00112 for( int col = 0; col < numCol; col++ ) 00113 { 00114 for( int row = 0; row < numRow; row++ ) 00115 { 00116 rowSum[row] += matrix[row + col * numRow]; 00117 if( fabs( matrix[row + col * numRow] ) > maxVal ) { maxVal = fabs( matrix[row + col * numRow] ); } 00118 } 00119 } 00120 double machEps = 1.e-14 * static_cast< double >( numRow ); // better to use system parameters 00121 double upperBound = machEps * maxVal + 1.e-15; 00122 for( int row = 0; row < numRow; row++ ) 00123 { 00124 if( fabs( rowSum[row] - oldRowSum[row] ) > upperBound ) 00125 { 00126 same = false; 00127 if( mNonGradDebug >= 2 ) 00128 { 00129 std::cout << " Warning: NonGradient Row Sum " << row << " Test: value " << rowSum[row] 00130 << " Discrepancy " << rowSum[row] - oldRowSum[row] << " maxVal " << maxVal << " numRow " 00131 << numRow << " numCol " << numCol << std::endl; 00132 } 00133 MSQ_PRINT( 2 ) 00134 ( "NonGradient Row Sum [%d] Test failure: value %22.15e difference %22.15e \n", row, rowSum[row], 00135 rowSum[row] - oldRowSum[row] ); 00136 } 00137 } 00138 return ( same ); 00139 } 00140 00141 void NonGradient::getRowSum( int numRow, int numCol, std::vector< double >& matrix, std::vector< double >& rowSum ) 00142 { 00143 for( int row = 0; row < numRow; row++ ) 00144 { 00145 rowSum[row] = 0.; 00146 } 00147 for( int col = 0; col < numCol; col++ ) 00148 { 00149 for( int row = 0; row < numRow; row++ ) 00150 { 00151 rowSum[row] += matrix[row + col * numRow]; 00152 } 00153 } 00154 } 00155 00156 double NonGradient::evaluate( double* point, PatchData& pd, MsqError& err ) 00157 { 00158 if( pd.num_free_vertices() > 1 ) 00159 { 00160 MSQ_SETERR( err ) 00161 ( "Only one free vertex per patch implemented", MsqError::NOT_IMPLEMENTED ); 00162 } 00163 const size_t vertexIndex = 0; 00164 Vector3D originalVec = pd.vertex_by_index( vertexIndex ); 00165 Vector3D pointVec; 00166 for( int index = 0; index < 3; index++ ) 00167 { 00168 pointVec[index] = point[index]; 00169 } 00170 pd.set_vertex_coordinates( pointVec, vertexIndex, err ); 00171 pd.snap_vertex_to_domain( vertexIndex, err ); 00172 00173 OFEvaluator& obj_func = get_objective_function_evaluator(); 00174 TerminationCriterion* term_crit = get_inner_termination_criterion(); 00175 00176 double value; 00177 bool feasible = obj_func.evaluate( pd, value, err ); // MSQ_ERRZERO(err); 00178 term_crit->accumulate_inner( pd, value, 0, err ); // MSQ_CHKERR(err); 00179 if( err.error_code() == err.BARRIER_VIOLATED ) err.clear(); // barrier violation not an error here 00180 MSQ_ERRZERO( err ); 00181 pd.set_vertex_coordinates( originalVec, vertexIndex, err ); 00182 if( !feasible && mUseExactPenaltyFunction ) 00183 { // "value" undefined btw 00184 double ensureFiniteRtol = .25; 00185 value = DBL_MAX * ensureFiniteRtol; 00186 // std::cout << " Infeasible patch: " << value << std::endl; printPatch( pd, err ); 00187 } 00188 return value; 00189 } 00190 00191 // Extrapolate by a factor of fac through the face of the simplex 00192 // opposite the high point to a new point. If the new point is 00193 // better, swap it with the high point. 00194 double NonGradient::amotry( std::vector< double >& p_simplex, std::vector< double >& p_height, double psum[], int ihi, 00195 double fac, PatchData& pd, MsqError& err ) 00196 { 00197 int numRow = getDimension(); 00198 // int numCol = numRow + 1; 00199 std::vector< double > ptry( numRow ); // does this make sense? 00200 double fac1 = ( 1.0 - fac ) / static_cast< double >( numRow ); 00201 double fac2 = fac1 - fac; 00202 for( int row = 0; row < numRow; row++ ) 00203 { 00204 ptry[row] = psum[row] * fac1 - p_simplex[row + ihi * numRow] * fac2; 00205 } 00206 if( mNonGradDebug >= 3 ) { std::cout << "Try "; } 00207 MSQ_PRINT( 3 )( "Try" ); 00208 double ytry = evaluate( &ptry[0], pd, err ); // value at trial point 00209 if( mNonGradDebug >= 3 ) { std::cout << ytry << std::endl; } 00210 MSQ_PRINT( 3 )( "yTry" ); 00211 if( ytry < p_height[ihi] ) // better than highest (worst) 00212 { 00213 p_height[ihi] = ytry; // swap ihi and ytry 00214 for( int row = 0; row < numRow; row++ ) 00215 { 00216 psum[row] += ( ptry[row] - p_simplex[row + ihi * numRow] ); 00217 p_simplex[row + ihi * numRow] = ptry[row]; 00218 } 00219 } 00220 return ytry; 00221 } 00222 00223 void NonGradient::printPatch( const PatchData& pd, MsqError& err ) 00224 { 00225 if( mNonGradDebug == 0 ) { return; } 00226 const size_t numNode = pd.num_nodes(); // 27, 27 what? 00227 MSQ_PRINT( 3 )( "Number of Vertices: %d\n", (int)pd.num_nodes() ); 00228 const size_t numVert = pd.num_free_vertices(); // 1 00229 MSQ_PRINT( 3 )( "Num Free = %d\n", (int)pd.num_free_vertices() ); 00230 const size_t numSlaveVert = pd.num_slave_vertices(); // 0 00231 const size_t numCoin = pd.num_corners(); // 64 00232 const MsqVertex* coord = pd.get_vertex_array( err ); 00233 MSQ_PRINT( 3 )( "Number of Vertices: %d\n", (int)pd.num_nodes() ); 00234 00235 std::cout << "Patch " << numNode << " " << numVert << " " << numSlaveVert << " " << numCoin << std::endl; 00236 MSQ_PRINT( 3 )( "\n" ); 00237 std::cout << "Coordinate "; 00238 std::cout << " " << std::endl; 00239 for( size_t index = 0; index < numVert; index++ ) 00240 { 00241 std::cout << coord[index][0] << " " << coord[index][1] << " " << coord[index][2] << std::endl; 00242 } 00243 // const size_t numElt = pd.num_elements(); 00244 if( mNonGradDebug >= 3 ) { std::cout << "Number of Elements: " << pd.num_elements() << std::endl; } 00245 MSQ_PRINT( 3 )( "Number of Elements: %d\n", (int)pd.num_elements() ); 00246 } 00247 00248 void NonGradient::printSimplex( std::vector< double >& p_simplex, std::vector< double >& p_height ) 00249 { 00250 int numRow = getDimension(); 00251 int numCol = numRow + 1; 00252 for( int col = 0; col < numCol; col++ ) 00253 { 00254 // std::cout << "simplex[ " << col << "]= " ; 00255 for( int row = 0; row < numRow; row++ ) 00256 { 00257 std::cout << p_simplex[row + col * numRow] << " "; 00258 } 00259 // std::cout << " " << p_height[col] << std::endl; 00260 } 00261 for( int col = 0; col < numCol; col++ ) 00262 { 00263 std::cout << p_height[col] << " "; 00264 } 00265 std::cout << std::endl; 00266 } 00267 00268 int NonGradient::getPatchDimension( const PatchData& pd, MsqError& err ) 00269 { 00270 const size_t numElt = pd.num_elements(); 00271 unsigned edimMax = 0; // in case numElt == 0 00272 for( size_t elementId = 0; elementId < numElt; elementId++ ) 00273 { 00274 const MsqMeshEntity& element = pd.element_by_index( elementId ); 00275 EntityTopology type = element.get_element_type(); 00276 unsigned edim = TopologyInfo::dimension( type ); 00277 if( elementId == 0 ) { edimMax = edim; } 00278 else 00279 { 00280 if( edimMax != edim ) 00281 { 00282 MSQ_SETERR( err ) 00283 ( "A patch has elements of different dimensions", MsqError::INVALID_MESH ); 00284 std::cout << "A patch has elements of different dimensions" << edimMax << " " << edim << std::endl; 00285 if( edimMax < edim ) { edimMax = edim; } 00286 } 00287 } 00288 } 00289 return ( edimMax ); 00290 } 00291 00292 void NonGradient::initialize( PatchData& /*pd*/, MsqError& /*err*/ ) {} 00293 00294 void NonGradient::initialize_mesh_iteration( PatchData& pd, MsqError& err ) 00295 { 00296 unsigned elementDimension = getPatchDimension( pd, err ); // to do: react to error 00297 unsigned dimension = elementDimension * pd.num_free_vertices(); 00298 // printPatch( pd, err ); 00299 setDimension( dimension ); 00300 int maxNumEval = 100 * dimension; // 1. Make this a user parameter 00301 setMaxNumEval( maxNumEval ); 00302 double threshold = 1.e-10; // avoid division by zero 00303 setThreshold( threshold ); 00304 double minEdgeLen = 0.0; 00305 double maxEdgeLen = 0.0; 00306 // double ftol = 0.; 00307 if( dimension > 0 ) 00308 { 00309 pd.get_minmax_edge_length( minEdgeLen, maxEdgeLen ); 00310 // ftol = minEdgeLen * 1.e-4; // Turn off Amoeba convergence criterion 00311 if( mNonGradDebug >= 1 ) 00312 { 00313 std::cout << "minimum edge length " << minEdgeLen << " maximum edge length " << maxEdgeLen << std::endl; 00314 } 00315 MSQ_PRINT( 3 ) 00316 ( "minimum edge length %e maximum edge length %e\n", minEdgeLen, maxEdgeLen ); 00317 } 00318 // setTolerance(ftol); 00319 unsigned numRow = dimension; 00320 unsigned numCol = numRow + 1; 00321 if( numRow * numCol <= simplex.max_size() ) 00322 { 00323 simplex.assign( numRow * numCol, 0. ); // guard against previous simplex value 00324 double scale = minEdgeLen * mScaleDiameter; 00325 ; 00326 const MsqVertex* coord = pd.get_vertex_array( err ); 00327 if( pd.num_free_vertices() > 1 ) 00328 { 00329 MSQ_SETERR( err ) 00330 ( "Only one free vertex per patch implemented", MsqError::NOT_IMPLEMENTED ); 00331 } 00332 size_t index = 0; 00333 for( unsigned col = 0; col < numCol; col++ ) 00334 { 00335 for( unsigned row = 0; row < numRow; row++ ) 00336 { 00337 simplex[row + col * numRow] = coord[index][row]; 00338 if( row == col - 1 ) { simplex[row + col * numRow] += scale / static_cast< double >( numCol ); } 00339 } 00340 } 00341 } 00342 else 00343 { 00344 MSQ_SETERR( err )( "Use patch with fewer free vertices", MsqError::OUT_OF_MEMORY ); 00345 if( mNonGradDebug >= 1 ) { std::cout << "ERROR: Too many free vertices in patch" << std::endl; } 00346 // MSQ_PRINT(1)("ERROR: Too many free vertices in patch\n"); 00347 } 00348 } 00349 00350 void NonGradient::optimize_vertex_positions( PatchData& pd, MsqError& err ) 00351 { 00352 MSQ_FUNCTION_TIMER( "NonGradient::optimize_vertex_positions" ); 00353 int numRow = getDimension(); 00354 int numCol = numRow + 1; 00355 std::vector< double > p_height( numCol ); 00356 00357 for( int col = 0; col < numCol; col++ ) 00358 { 00359 p_height[col] = evaluate( &simplex[col * numRow], pd, err ); // eval patch stuff 00360 } 00361 if( mNonGradDebug > 0 ) { printSimplex( simplex, p_height ); } 00362 00363 // standardization 00364 TerminationCriterion* term_crit = get_inner_termination_criterion(); 00365 // int maxNumEval = getMaxNumEval(); 00366 // double threshold = getThreshold(); 00367 00368 // double ftol = getTolerance(); 00369 int ilo = 0; // height[ilo]<=... 00370 int inhi = 0; //...<=height[inhi]<= 00371 int ihi = 0; //<=height[ihi] 00372 // double rtol = 2.*ftol; 00373 double ysave; 00374 double ytry; 00375 bool afterEvaluation = false; 00376 std::vector< double > rowSum( numRow ); 00377 getRowSum( numRow, numCol, simplex, rowSum ); 00378 while( !( term_crit->terminate() ) ) 00379 { 00380 00381 if( mNonGradDebug > 0 ) { printSimplex( simplex, p_height ); } 00382 // std::cout << "rtol " << rtol << " ftol " << ftol << " MesquiteIter " << 00383 // term_crit->get_iteration_count() << " Done " << term_crit->terminate() << std::endl; 00384 00385 if( afterEvaluation ) 00386 { 00387 // reflect highPt through opposite face 00388 // p_height[0] may vanish 00389 /* 00390 if( !testRowSum( numRow, numCol, &simplex[0], &rowSum[0]) ) 00391 { 00392 // Before uncommenting here and ... 00393 //MSQ_SETERR(err)("Internal check sum test A failed", MsqError::INTERNAL_ERROR); 00394 //MSQ_ERRRTN(err); 00395 } 00396 */ 00397 ytry = amotry( simplex, p_height, &rowSum[0], ihi, -1.0, pd, err ); 00398 /* 00399 if( !testRowSum( numRow, numCol, &simplex[0], &rowSum[0]) ) 00400 { 00401 // ... here, determine a maxVal majorizing the previous as well as the current 00402 simplex. 00403 //MSQ_SETERR(err)("Internal check sum test B failed", 00404 MsqError::INTERNAL_ERROR); 00405 //MSQ_ERRRTN(err); 00406 } 00407 */ 00408 00409 /* 00410 if( p_height[0] == 0.) 00411 { 00412 MSQ_SETERR(err)("(B) Zero objective function value", MsqError::INTERNAL_ERROR); 00413 exit(-1); 00414 } 00415 */ 00416 if( ytry <= p_height[ilo] ) 00417 { 00418 ytry = amotry( simplex, p_height, &rowSum[0], ihi, -2.0, pd, err ); 00419 if( mNonGradDebug >= 3 ) { std::cout << "Reflect and Expand from highPt " << ytry << std::endl; } 00420 // MSQ_PRINT(3)("Reflect and Expand from highPt : %e\n",ytry); 00421 } 00422 else 00423 { 00424 if( ytry >= p_height[inhi] ) 00425 { 00426 ysave = p_height[ihi]; // Contract along highPt 00427 ytry = amotry( simplex, p_height, &rowSum[0], ihi, 0.5, pd, err ); 00428 if( ytry >= ysave ) 00429 { // contract all directions toward lowPt 00430 for( int col = 0; col < numCol; col++ ) 00431 { 00432 if( col != ilo ) 00433 { 00434 for( int row = 0; row < numRow; row++ ) 00435 { 00436 rowSum[row] = 0.5 * ( simplex[row + col * numRow] + simplex[row + ilo * numRow] ); 00437 simplex[row + col * numRow] = rowSum[row]; 00438 } 00439 p_height[col] = evaluate( &rowSum[0], pd, err ); 00440 if( mNonGradDebug >= 3 ) 00441 { 00442 std::cout << "Contract all directions toward lowPt value( " << col 00443 << " ) = " << p_height[col] << " ilo = " << ilo << std::endl; 00444 } 00445 // MSQ_PRINT(3)("Contract all directions toward lowPt value( %d ) = 00446 // %e ilo = %d\n", col, p_height[col], ilo); 00447 } 00448 } 00449 } 00450 } 00451 } // ytri > h(ilo) 00452 } // after evaluation 00453 ilo = 1; // conditional operator or inline if 00454 ihi = p_height[0] > p_height[1] ? ( inhi = 1, 0 ) : ( inhi = 0, 1 ); 00455 for( int col = 0; col < numCol; col++ ) 00456 { 00457 if( p_height[col] <= p_height[ilo] ) 00458 { 00459 ilo = col; // ilo := argmin height 00460 } 00461 if( p_height[col] > p_height[ihi] ) 00462 { 00463 inhi = ihi; 00464 ihi = col; 00465 } 00466 else // p_height[ihi] >= p_height[col] 00467 if( col != ihi && p_height[col] > p_height[inhi] ) 00468 inhi = col; 00469 } 00470 // rtol=2.0*fabs( p_height[ihi]-p_height[ilo] )/ 00471 // ( fabs(p_height[ihi])+fabs(p_height[ilo])+threshold ); 00472 afterEvaluation = true; 00473 } // while not converged 00474 00475 // Always set to current best mesh. 00476 { 00477 if( ilo != 0 ) 00478 { 00479 double yTemp = p_height[0]; 00480 p_height[0] = p_height[ilo]; // height dimension numCol 00481 p_height[ilo] = yTemp; 00482 for( int row = 1; row < numRow; row++ ) 00483 { 00484 yTemp = simplex[row]; 00485 simplex[row] = simplex[row + ilo * numRow]; 00486 simplex[row + ilo * numRow] = yTemp; 00487 } 00488 } 00489 } 00490 if( pd.num_free_vertices() > 1 ) 00491 { 00492 MSQ_SETERR( err ) 00493 ( "Only one free vertex per patch implemented", MsqError::NOT_IMPLEMENTED ); 00494 } 00495 00496 Vector3D newPoint( &simplex[0] ); 00497 size_t vertexIndex = 0; // fix c.f. freeVertexIndex 00498 pd.set_vertex_coordinates( newPoint, vertexIndex, err ); 00499 pd.snap_vertex_to_domain( vertexIndex, err ); 00500 if( term_crit->terminate() ) 00501 { 00502 if( mNonGradDebug >= 1 ) { std::cout << "Optimization Termination OptStatus: Max Iter Exceeded" << std::endl; } 00503 // MSQ_PRINT(1)("Optimization Termination OptStatus: Max Iter Exceeded\n"); 00504 } 00505 } 00506 00507 void NonGradient::terminate_mesh_iteration( PatchData& /*pd*/, MsqError& /*err*/ ) 00508 { 00509 if( mNonGradDebug >= 2 ) { std::cout << "- Executing NonGradient::iteration_complete()" << std::endl; } 00510 // MSQ_PRINT(2)("\n - Executing NonGradient::iteration_complete() \n"); 00511 } 00512 00513 void NonGradient::cleanup() 00514 { 00515 if( mNonGradDebug >= 2 ) { std::cout << " - Executing NonGradient::iteration_end()" << std::endl; } 00516 // MSQ_PRINT(2)("\n - Executing NonGradient::iteration_end() \n"); 00517 } 00518 00519 } // namespace MBMesquite