MOAB: Mesh Oriented datABase  (version 5.2.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         { std::cout << "minimum edge length " << minEdgeLen << " maximum edge length " << maxEdgeLen << std::endl; }
00313         MSQ_PRINT( 3 )
00314         ( "minimum edge length %e    maximum edge length %e\n", minEdgeLen, maxEdgeLen );
00315     }
00316     //  setTolerance(ftol);
00317     unsigned numRow = dimension;
00318     unsigned numCol = numRow + 1;
00319     if( numRow * numCol <= simplex.max_size() )
00320     {
00321         simplex.assign( numRow * numCol, 0. );  // guard against previous simplex value
00322         double scale = minEdgeLen * mScaleDiameter;
00323         ;
00324         const MsqVertex* coord = pd.get_vertex_array( err );
00325         if( pd.num_free_vertices() > 1 )
00326         {
00327             MSQ_SETERR( err )
00328             ( "Only one free vertex per patch implemented", MsqError::NOT_IMPLEMENTED );
00329         }
00330         size_t index = 0;
00331         for( unsigned col = 0; col < numCol; col++ )
00332         {
00333             for( unsigned row = 0; row < numRow; row++ )
00334             {
00335                 simplex[row + col * numRow] = coord[index][row];
00336                 if( row == col - 1 ) { simplex[row + col * numRow] += scale / static_cast< double >( numCol ); }
00337             }
00338         }
00339     }
00340     else
00341     {
00342         MSQ_SETERR( err )( "Use patch with fewer free vertices", MsqError::OUT_OF_MEMORY );
00343         if( mNonGradDebug >= 1 ) { std::cout << "ERROR: Too many free vertices in patch" << std::endl; }
00344         // MSQ_PRINT(1)("ERROR: Too many free vertices in patch\n");
00345     }
00346 }
00347 
00348 void NonGradient::optimize_vertex_positions( PatchData& pd, MsqError& err )
00349 {
00350     MSQ_FUNCTION_TIMER( "NonGradient::optimize_vertex_positions" );
00351     int numRow = getDimension();
00352     int numCol = numRow + 1;
00353     std::vector< double > p_height( numCol );
00354 
00355     for( int col = 0; col < numCol; col++ )
00356     {
00357         p_height[col] = evaluate( &simplex[col * numRow], pd, err );  //  eval patch stuff
00358     }
00359     if( mNonGradDebug > 0 ) { printSimplex( simplex, p_height ); }
00360 
00361     // standardization
00362     TerminationCriterion* term_crit = get_inner_termination_criterion();
00363     // int maxNumEval = getMaxNumEval();
00364     // double threshold = getThreshold();
00365 
00366     //  double ftol = getTolerance();
00367     int ilo  = 0;  // height[ilo]<=...
00368     int inhi = 0;  //...<=height[inhi]<=
00369     int ihi  = 0;  //<=height[ihi]
00370     //  double rtol = 2.*ftol;
00371     double ysave;
00372     double ytry;
00373     bool afterEvaluation = false;
00374     std::vector< double > rowSum( numRow );
00375     getRowSum( numRow, numCol, simplex, rowSum );
00376     while( !( term_crit->terminate() ) )
00377     {
00378 
00379         if( mNonGradDebug > 0 ) { printSimplex( simplex, p_height ); }
00380         // std::cout << "rtol " << rtol << " ftol " << ftol << " MesquiteIter " <<
00381         // term_crit->get_iteration_count() << " Done " << term_crit->terminate() << std::endl;
00382 
00383         if( afterEvaluation )
00384         {
00385             // reflect highPt through opposite face
00386             // p_height[0] may vanish
00387             /*
00388                   if( !testRowSum( numRow, numCol, &simplex[0], &rowSum[0]) )
00389                   {
00390                     // Before uncommenting here and ...
00391                     //MSQ_SETERR(err)("Internal check sum test A failed", MsqError::INTERNAL_ERROR);
00392                     //MSQ_ERRRTN(err);
00393                   }
00394             */
00395             ytry = amotry( simplex, p_height, &rowSum[0], ihi, -1.0, pd, err );
00396             /*
00397                   if( !testRowSum( numRow, numCol, &simplex[0], &rowSum[0]) )
00398                   {
00399                      // ... here, determine a maxVal majorizing the previous as well as the current
00400                simplex.
00401                      //MSQ_SETERR(err)("Internal check sum test B failed",
00402                MsqError::INTERNAL_ERROR);
00403                      //MSQ_ERRRTN(err);
00404                   }
00405             */
00406 
00407             /*
00408                   if( p_height[0] == 0.)
00409                   {
00410                      MSQ_SETERR(err)("(B) Zero objective function value", MsqError::INTERNAL_ERROR);
00411                      exit(-1);
00412                   }
00413             */
00414             if( ytry <= p_height[ilo] )
00415             {
00416                 ytry = amotry( simplex, p_height, &rowSum[0], ihi, -2.0, pd, err );
00417                 if( mNonGradDebug >= 3 ) { std::cout << "Reflect and Expand from highPt " << ytry << std::endl; }
00418                 // MSQ_PRINT(3)("Reflect and Expand from highPt : %e\n",ytry);
00419             }
00420             else
00421             {
00422                 if( ytry >= p_height[inhi] )
00423                 {
00424                     ysave = p_height[ihi];  // Contract along highPt
00425                     ytry  = amotry( simplex, p_height, &rowSum[0], ihi, 0.5, pd, err );
00426                     if( ytry >= ysave )
00427                     {  // contract all directions toward lowPt
00428                         for( int col = 0; col < numCol; col++ )
00429                         {
00430                             if( col != ilo )
00431                             {
00432                                 for( int row = 0; row < numRow; row++ )
00433                                 {
00434                                     rowSum[row] = 0.5 * ( simplex[row + col * numRow] + simplex[row + ilo * numRow] );
00435                                     simplex[row + col * numRow] = rowSum[row];
00436                                 }
00437                                 p_height[col] = evaluate( &rowSum[0], pd, err );
00438                                 if( mNonGradDebug >= 3 )
00439                                 {
00440                                     std::cout << "Contract all directions toward lowPt value( " << col
00441                                               << " ) = " << p_height[col] << " ilo = " << ilo << std::endl;
00442                                 }
00443                                 // MSQ_PRINT(3)("Contract all directions toward lowPt value( %d ) =
00444                                 // %e    ilo = %d\n", col, p_height[col], ilo);
00445                             }
00446                         }
00447                     }
00448                 }
00449             }     // ytri > h(ilo)
00450         }         // after evaluation
00451         ilo = 1;  // conditional operator or inline if
00452         ihi = p_height[0] > p_height[1] ? ( inhi = 1, 0 ) : ( inhi = 0, 1 );
00453         for( int col = 0; col < numCol; col++ )
00454         {
00455             if( p_height[col] <= p_height[ilo] )
00456             {
00457                 ilo = col;  // ilo := argmin height
00458             }
00459             if( p_height[col] > p_height[ihi] )
00460             {
00461                 inhi = ihi;
00462                 ihi  = col;
00463             }
00464             else  // p_height[ihi] >= p_height[col]
00465                 if( col != ihi && p_height[col] > p_height[inhi] )
00466                 inhi = col;
00467         }
00468         //    rtol=2.0*fabs( p_height[ihi]-p_height[ilo] )/
00469         //         ( fabs(p_height[ihi])+fabs(p_height[ilo])+threshold );
00470         afterEvaluation = true;
00471     }  //  while not converged
00472 
00473     // Always set to current best mesh.
00474     {
00475         if( ilo != 0 )
00476         {
00477             double yTemp  = p_height[0];
00478             p_height[0]   = p_height[ilo];  // height dimension numCol
00479             p_height[ilo] = yTemp;
00480             for( int row = 1; row < numRow; row++ )
00481             {
00482                 yTemp                       = simplex[row];
00483                 simplex[row]                = simplex[row + ilo * numRow];
00484                 simplex[row + ilo * numRow] = yTemp;
00485             }
00486         }
00487     }
00488     if( pd.num_free_vertices() > 1 )
00489     {
00490         MSQ_SETERR( err )
00491         ( "Only one free vertex per patch implemented", MsqError::NOT_IMPLEMENTED );
00492     }
00493 
00494     Vector3D newPoint( &simplex[0] );
00495     size_t vertexIndex = 0;  // fix c.f. freeVertexIndex
00496     pd.set_vertex_coordinates( newPoint, vertexIndex, err );
00497     pd.snap_vertex_to_domain( vertexIndex, err );
00498     if( term_crit->terminate() )
00499     {
00500         if( mNonGradDebug >= 1 ) { std::cout << "Optimization Termination OptStatus: Max Iter Exceeded" << std::endl; }
00501         // MSQ_PRINT(1)("Optimization Termination OptStatus: Max Iter Exceeded\n");
00502     }
00503 }
00504 
00505 void NonGradient::terminate_mesh_iteration( PatchData& /*pd*/, MsqError& /*err*/ )
00506 {
00507     if( mNonGradDebug >= 2 ) { std::cout << "- Executing NonGradient::iteration_complete()" << std::endl; }
00508     // MSQ_PRINT(2)("\n - Executing NonGradient::iteration_complete() \n");
00509 }
00510 
00511 void NonGradient::cleanup()
00512 {
00513     if( mNonGradDebug >= 2 ) { std::cout << " - Executing NonGradient::iteration_end()" << std::endl; }
00514     // MSQ_PRINT(2)("\n - Executing NonGradient::iteration_end() \n");
00515 }
00516 
00517 }  // namespace MBMesquite
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines