MOAB: Mesh Oriented datABase  (version 5.3.1)
NonGradient.cpp
Go to the documentation of this file.
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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines