MOAB: Mesh Oriented datABase  (version 5.4.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     [email protected], [email protected], [email protected],
00024     [email protected], [email protected], [email protected]
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 )
00080     {
00081         return;
00082     }
00083     else
00084     {
00085         default_crit->add_iteration_limit( 5 );
00086     }
00087 }
00088 
00089 NonGradient::NonGradient( ObjectiveFunction* of, MsqError& err )
00090     : VertexMover( of ), PatchSetUser( true ), projectGradient( false ), mDimension( 0 ), mThreshold( 0.0 ),
00091       mTolerance( 0.0 ), mMaxNumEval( 0 ), mNonGradDebug( 0 ), mUseExactPenaltyFunction( true ), mScaleDiameter( 0.1 )
00092 {
00093     set_debugging_level( 2 );
00094     // set the default inner termination criterion
00095     TerminationCriterion* default_crit = get_inner_termination_criterion();
00096     if( default_crit == NULL )
00097     {
00098         MSQ_SETERR( err )
00099         ( "QualityImprover did not create a default inner "
00100           "termination criterion.",
00101           MsqError::INVALID_STATE );
00102         return;
00103     }
00104     else
00105     {
00106         default_crit->add_iteration_limit( 5 );MSQ_ERRRTN( err );
00107     }
00108 }
00109 
00110 bool NonGradient::testRowSum( int numRow, int numCol, double* matrix, double* oldRowSum )
00111 {
00112     bool same = true;
00113     std::vector< double > rowSum( numRow, 0. );
00114     double maxVal = 0.;
00115     for( int col = 0; col < numCol; col++ )
00116     {
00117         for( int row = 0; row < numRow; row++ )
00118         {
00119             rowSum[row] += matrix[row + col * numRow];
00120             if( fabs( matrix[row + col * numRow] ) > maxVal )
00121             {
00122                 maxVal = fabs( matrix[row + col * numRow] );
00123             }
00124         }
00125     }
00126     double machEps    = 1.e-14 * static_cast< double >( numRow );  // better to use system parameters
00127     double upperBound = machEps * maxVal + 1.e-15;
00128     for( int row = 0; row < numRow; row++ )
00129     {
00130         if( fabs( rowSum[row] - oldRowSum[row] ) > upperBound )
00131         {
00132             same = false;
00133             if( mNonGradDebug >= 2 )
00134             {
00135                 std::cout << " Warning: NonGradient Row Sum " << row << " Test: value " << rowSum[row]
00136                           << " Discrepancy " << rowSum[row] - oldRowSum[row] << " maxVal " << maxVal << " numRow "
00137                           << numRow << " numCol " << numCol << std::endl;
00138             }
00139             MSQ_PRINT( 2 )
00140             ( "NonGradient Row Sum [%d] Test failure: value %22.15e  difference %22.15e \n", row, rowSum[row],
00141               rowSum[row] - oldRowSum[row] );
00142         }
00143     }
00144     return ( same );
00145 }
00146 
00147 void NonGradient::getRowSum( int numRow, int numCol, std::vector< double >& matrix, std::vector< double >& rowSum )
00148 {
00149     for( int row = 0; row < numRow; row++ )
00150     {
00151         rowSum[row] = 0.;
00152     }
00153     for( int col = 0; col < numCol; col++ )
00154     {
00155         for( int row = 0; row < numRow; row++ )
00156         {
00157             rowSum[row] += matrix[row + col * numRow];
00158         }
00159     }
00160 }
00161 
00162 double NonGradient::evaluate( double* point, PatchData& pd, MsqError& err )
00163 {
00164     if( pd.num_free_vertices() > 1 )
00165     {
00166         MSQ_SETERR( err )
00167         ( "Only one free vertex per patch implemented", MsqError::NOT_IMPLEMENTED );
00168     }
00169     const size_t vertexIndex = 0;
00170     Vector3D originalVec     = pd.vertex_by_index( vertexIndex );
00171     Vector3D pointVec;
00172     for( int index = 0; index < 3; index++ )
00173     {
00174         pointVec[index] = point[index];
00175     }
00176     pd.set_vertex_coordinates( pointVec, vertexIndex, err );
00177     pd.snap_vertex_to_domain( vertexIndex, err );
00178 
00179     OFEvaluator& obj_func           = get_objective_function_evaluator();
00180     TerminationCriterion* term_crit = get_inner_termination_criterion();
00181 
00182     double value;
00183     bool feasible = obj_func.evaluate( pd, value, err );         // MSQ_ERRZERO(err);
00184     term_crit->accumulate_inner( pd, value, 0, err );            // MSQ_CHKERR(err);
00185     if( err.error_code() == err.BARRIER_VIOLATED ) err.clear();  // barrier violation not an error here
00186     MSQ_ERRZERO( err );
00187     pd.set_vertex_coordinates( originalVec, vertexIndex, err );
00188     if( !feasible && mUseExactPenaltyFunction )
00189     {  // "value" undefined btw
00190         double ensureFiniteRtol = .25;
00191         value                   = DBL_MAX * ensureFiniteRtol;
00192         // std::cout << " Infeasible patch: " << value << std::endl; printPatch( pd, err );
00193     }
00194     return value;
00195 }
00196 
00197 // Extrapolate by a factor of fac through the face of the simplex
00198 // opposite the high point to a new point.  If the new point is
00199 // better, swap it with the high point.
00200 double NonGradient::amotry( std::vector< double >& p_simplex,
00201                             std::vector< double >& p_height,
00202                             double psum[],
00203                             int ihi,
00204                             double fac,
00205                             PatchData& pd,
00206                             MsqError& err )
00207 {
00208     int numRow = getDimension();
00209     // int numCol = numRow + 1;
00210     std::vector< double > ptry( numRow );  // does this make sense?
00211     double fac1 = ( 1.0 - fac ) / static_cast< double >( numRow );
00212     double fac2 = fac1 - fac;
00213     for( int row = 0; row < numRow; row++ )
00214     {
00215         ptry[row] = psum[row] * fac1 - p_simplex[row + ihi * numRow] * fac2;
00216     }
00217     if( mNonGradDebug >= 3 )
00218     {
00219         std::cout << "Try ";
00220     }
00221     MSQ_PRINT( 3 )( "Try" );
00222     double ytry = evaluate( &ptry[0], pd, err );  // value at trial point
00223     if( mNonGradDebug >= 3 )
00224     {
00225         std::cout << ytry << std::endl;
00226     }
00227     MSQ_PRINT( 3 )( "yTry" );
00228     if( ytry < p_height[ihi] )  // better than highest (worst)
00229     {
00230         p_height[ihi] = ytry;  // swap ihi and ytry
00231         for( int row = 0; row < numRow; row++ )
00232         {
00233             psum[row] += ( ptry[row] - p_simplex[row + ihi * numRow] );
00234             p_simplex[row + ihi * numRow] = ptry[row];
00235         }
00236     }
00237     return ytry;
00238 }
00239 
00240 void NonGradient::printPatch( const PatchData& pd, MsqError& err )
00241 {
00242     if( mNonGradDebug == 0 )
00243     {
00244         return;
00245     }
00246     const size_t numNode = pd.num_nodes();  // 27,  27 what?
00247     MSQ_PRINT( 3 )( "Number of Vertices: %d\n", (int)pd.num_nodes() );
00248     const size_t numVert = pd.num_free_vertices();  // 1
00249     MSQ_PRINT( 3 )( "Num Free = %d\n", (int)pd.num_free_vertices() );
00250     const size_t numSlaveVert = pd.num_slave_vertices();  // 0
00251     const size_t numCoin      = pd.num_corners();         // 64
00252     const MsqVertex* coord    = pd.get_vertex_array( err );
00253     MSQ_PRINT( 3 )( "Number of Vertices: %d\n", (int)pd.num_nodes() );
00254 
00255     std::cout << "Patch " << numNode << "  " << numVert << "  " << numSlaveVert << "  " << numCoin << std::endl;
00256     MSQ_PRINT( 3 )( "\n" );
00257     std::cout << "Coordinate ";
00258     std::cout << "         " << std::endl;
00259     for( size_t index = 0; index < numVert; index++ )
00260     {
00261         std::cout << coord[index][0] << "  " << coord[index][1] << "  " << coord[index][2] << std::endl;
00262     }
00263     // const size_t numElt = pd.num_elements();
00264     if( mNonGradDebug >= 3 )
00265     {
00266         std::cout << "Number of Elements: " << pd.num_elements() << std::endl;
00267     }
00268     MSQ_PRINT( 3 )( "Number of Elements: %d\n", (int)pd.num_elements() );
00269 }
00270 
00271 void NonGradient::printSimplex( std::vector< double >& p_simplex, std::vector< double >& p_height )
00272 {
00273     int numRow = getDimension();
00274     int numCol = numRow + 1;
00275     for( int col = 0; col < numCol; col++ )
00276     {
00277         // std::cout << "simplex[ " << col << "]= " ;
00278         for( int row = 0; row < numRow; row++ )
00279         {
00280             std::cout << p_simplex[row + col * numRow] << " ";
00281         }
00282         // std::cout << "           "  << p_height[col] << std::endl;
00283     }
00284     for( int col = 0; col < numCol; col++ )
00285     {
00286         std::cout << p_height[col] << " ";
00287     }
00288     std::cout << std::endl;
00289 }
00290 
00291 int NonGradient::getPatchDimension( const PatchData& pd, MsqError& err )
00292 {
00293     const size_t numElt = pd.num_elements();
00294     unsigned edimMax    = 0;  // in case numElt == 0
00295     for( size_t elementId = 0; elementId < numElt; elementId++ )
00296     {
00297         const MsqMeshEntity& element = pd.element_by_index( elementId );
00298         EntityTopology type          = element.get_element_type();
00299         unsigned edim                = TopologyInfo::dimension( type );
00300         if( elementId == 0 )
00301         {
00302             edimMax = edim;
00303         }
00304         else
00305         {
00306             if( edimMax != edim )
00307             {
00308                 MSQ_SETERR( err )
00309                 ( "A patch has elements of different dimensions", MsqError::INVALID_MESH );
00310                 std::cout << "A patch has elements of different dimensions" << edimMax << "  " << edim << std::endl;
00311                 if( edimMax < edim )
00312                 {
00313                     edimMax = edim;
00314                 }
00315             }
00316         }
00317     }
00318     return ( edimMax );
00319 }
00320 
00321 void NonGradient::initialize( PatchData& /*pd*/, MsqError& /*err*/ ) {}
00322 
00323 void NonGradient::initialize_mesh_iteration( PatchData& pd, MsqError& err )
00324 {
00325     unsigned elementDimension = getPatchDimension( pd, err );  // to do: react to error
00326     unsigned dimension        = elementDimension * pd.num_free_vertices();
00327     // printPatch( pd, err );
00328     setDimension( dimension );
00329     int maxNumEval = 100 * dimension;  // 1. Make this a user parameter
00330     setMaxNumEval( maxNumEval );
00331     double threshold = 1.e-10;  // avoid division by zero
00332     setThreshold( threshold );
00333     double minEdgeLen = 0.0;
00334     double maxEdgeLen = 0.0;
00335     //  double ftol = 0.;
00336     if( dimension > 0 )
00337     {
00338         pd.get_minmax_edge_length( minEdgeLen, maxEdgeLen );
00339         // ftol = minEdgeLen * 1.e-4; // Turn off Amoeba convergence criterion
00340         if( mNonGradDebug >= 1 )
00341         {
00342             std::cout << "minimum edge length " << minEdgeLen << " maximum edge length " << maxEdgeLen << std::endl;
00343         }
00344         MSQ_PRINT( 3 )
00345         ( "minimum edge length %e    maximum edge length %e\n", minEdgeLen, maxEdgeLen );
00346     }
00347     //  setTolerance(ftol);
00348     unsigned numRow = dimension;
00349     unsigned numCol = numRow + 1;
00350     if( numRow * numCol <= simplex.max_size() )
00351     {
00352         simplex.assign( numRow * numCol, 0. );  // guard against previous simplex value
00353         double scale = minEdgeLen * mScaleDiameter;
00354         ;
00355         const MsqVertex* coord = pd.get_vertex_array( err );
00356         if( pd.num_free_vertices() > 1 )
00357         {
00358             MSQ_SETERR( err )
00359             ( "Only one free vertex per patch implemented", MsqError::NOT_IMPLEMENTED );
00360         }
00361         size_t index = 0;
00362         for( unsigned col = 0; col < numCol; col++ )
00363         {
00364             for( unsigned row = 0; row < numRow; row++ )
00365             {
00366                 simplex[row + col * numRow] = coord[index][row];
00367                 if( row == col - 1 )
00368                 {
00369                     simplex[row + col * numRow] += scale / static_cast< double >( numCol );
00370                 }
00371             }
00372         }
00373     }
00374     else
00375     {
00376         MSQ_SETERR( err )( "Use patch with fewer free vertices", MsqError::OUT_OF_MEMORY );
00377         if( mNonGradDebug >= 1 )
00378         {
00379             std::cout << "ERROR: Too many free vertices in patch" << std::endl;
00380         }
00381         // MSQ_PRINT(1)("ERROR: Too many free vertices in patch\n");
00382     }
00383 }
00384 
00385 void NonGradient::optimize_vertex_positions( PatchData& pd, MsqError& err )
00386 {
00387     MSQ_FUNCTION_TIMER( "NonGradient::optimize_vertex_positions" );
00388     int numRow = getDimension();
00389     int numCol = numRow + 1;
00390     std::vector< double > p_height( numCol );
00391 
00392     for( int col = 0; col < numCol; col++ )
00393     {
00394         p_height[col] = evaluate( &simplex[col * numRow], pd, err );  //  eval patch stuff
00395     }
00396     if( mNonGradDebug > 0 )
00397     {
00398         printSimplex( simplex, p_height );
00399     }
00400 
00401     // standardization
00402     TerminationCriterion* term_crit = get_inner_termination_criterion();
00403     // int maxNumEval = getMaxNumEval();
00404     // double threshold = getThreshold();
00405 
00406     //  double ftol = getTolerance();
00407     int ilo  = 0;  // height[ilo]<=...
00408     int inhi = 0;  //...<=height[inhi]<=
00409     int ihi  = 0;  //<=height[ihi]
00410     //  double rtol = 2.*ftol;
00411     double ysave;
00412     double ytry;
00413     bool afterEvaluation = false;
00414     std::vector< double > rowSum( numRow );
00415     getRowSum( numRow, numCol, simplex, rowSum );
00416     while( !( term_crit->terminate() ) )
00417     {
00418 
00419         if( mNonGradDebug > 0 )
00420         {
00421             printSimplex( simplex, p_height );
00422         }
00423         // std::cout << "rtol " << rtol << " ftol " << ftol << " MesquiteIter " <<
00424         // term_crit->get_iteration_count() << " Done " << term_crit->terminate() << std::endl;
00425 
00426         if( afterEvaluation )
00427         {
00428             // reflect highPt through opposite face
00429             // p_height[0] may vanish
00430             /*
00431                   if( !testRowSum( numRow, numCol, &simplex[0], &rowSum[0]) )
00432                   {
00433                     // Before uncommenting here and ...
00434                     //MSQ_SETERR(err)("Internal check sum test A failed", MsqError::INTERNAL_ERROR);
00435                     //MSQ_ERRRTN(err);
00436                   }
00437             */
00438             ytry = amotry( simplex, p_height, &rowSum[0], ihi, -1.0, pd, err );
00439             /*
00440                   if( !testRowSum( numRow, numCol, &simplex[0], &rowSum[0]) )
00441                   {
00442                      // ... here, determine a maxVal majorizing the previous as well as the current
00443                simplex.
00444                      //MSQ_SETERR(err)("Internal check sum test B failed",
00445                MsqError::INTERNAL_ERROR);
00446                      //MSQ_ERRRTN(err);
00447                   }
00448             */
00449 
00450             /*
00451                   if( p_height[0] == 0.)
00452                   {
00453                      MSQ_SETERR(err)("(B) Zero objective function value", MsqError::INTERNAL_ERROR);
00454                      exit(-1);
00455                   }
00456             */
00457             if( ytry <= p_height[ilo] )
00458             {
00459                 ytry = amotry( simplex, p_height, &rowSum[0], ihi, -2.0, pd, err );
00460                 if( mNonGradDebug >= 3 )
00461                 {
00462                     std::cout << "Reflect and Expand from highPt " << ytry << std::endl;
00463                 }
00464                 // MSQ_PRINT(3)("Reflect and Expand from highPt : %e\n",ytry);
00465             }
00466             else
00467             {
00468                 if( ytry >= p_height[inhi] )
00469                 {
00470                     ysave = p_height[ihi];  // Contract along highPt
00471                     ytry  = amotry( simplex, p_height, &rowSum[0], ihi, 0.5, pd, err );
00472                     if( ytry >= ysave )
00473                     {  // contract all directions toward lowPt
00474                         for( int col = 0; col < numCol; col++ )
00475                         {
00476                             if( col != ilo )
00477                             {
00478                                 for( int row = 0; row < numRow; row++ )
00479                                 {
00480                                     rowSum[row] = 0.5 * ( simplex[row + col * numRow] + simplex[row + ilo * numRow] );
00481                                     simplex[row + col * numRow] = rowSum[row];
00482                                 }
00483                                 p_height[col] = evaluate( &rowSum[0], pd, err );
00484                                 if( mNonGradDebug >= 3 )
00485                                 {
00486                                     std::cout << "Contract all directions toward lowPt value( " << col
00487                                               << " ) = " << p_height[col] << " ilo = " << ilo << std::endl;
00488                                 }
00489                                 // MSQ_PRINT(3)("Contract all directions toward lowPt value( %d ) =
00490                                 // %e    ilo = %d\n", col, p_height[col], ilo);
00491                             }
00492                         }
00493                     }
00494                 }
00495             }     // ytri > h(ilo)
00496         }         // after evaluation
00497         ilo = 1;  // conditional operator or inline if
00498         ihi = p_height[0] > p_height[1] ? ( inhi = 1, 0 ) : ( inhi = 0, 1 );
00499         for( int col = 0; col < numCol; col++ )
00500         {
00501             if( p_height[col] <= p_height[ilo] )
00502             {
00503                 ilo = col;  // ilo := argmin height
00504             }
00505             if( p_height[col] > p_height[ihi] )
00506             {
00507                 inhi = ihi;
00508                 ihi  = col;
00509             }
00510             else  // p_height[ihi] >= p_height[col]
00511                 if( col != ihi && p_height[col] > p_height[inhi] ) inhi = col;
00512         }
00513         //    rtol=2.0*fabs( p_height[ihi]-p_height[ilo] )/
00514         //         ( fabs(p_height[ihi])+fabs(p_height[ilo])+threshold );
00515         afterEvaluation = true;
00516     }  //  while not converged
00517 
00518     // Always set to current best mesh.
00519     {
00520         if( ilo != 0 )
00521         {
00522             double yTemp  = p_height[0];
00523             p_height[0]   = p_height[ilo];  // height dimension numCol
00524             p_height[ilo] = yTemp;
00525             for( int row = 1; row < numRow; row++ )
00526             {
00527                 yTemp                       = simplex[row];
00528                 simplex[row]                = simplex[row + ilo * numRow];
00529                 simplex[row + ilo * numRow] = yTemp;
00530             }
00531         }
00532     }
00533     if( pd.num_free_vertices() > 1 )
00534     {
00535         MSQ_SETERR( err )
00536         ( "Only one free vertex per patch implemented", MsqError::NOT_IMPLEMENTED );
00537     }
00538 
00539     Vector3D newPoint( &simplex[0] );
00540     size_t vertexIndex = 0;  // fix c.f. freeVertexIndex
00541     pd.set_vertex_coordinates( newPoint, vertexIndex, err );
00542     pd.snap_vertex_to_domain( vertexIndex, err );
00543     if( term_crit->terminate() )
00544     {
00545         if( mNonGradDebug >= 1 )
00546         {
00547             std::cout << "Optimization Termination OptStatus: Max Iter Exceeded" << std::endl;
00548         }
00549         // MSQ_PRINT(1)("Optimization Termination OptStatus: Max Iter Exceeded\n");
00550     }
00551 }
00552 
00553 void NonGradient::terminate_mesh_iteration( PatchData& /*pd*/, MsqError& /*err*/ )
00554 {
00555     if( mNonGradDebug >= 2 )
00556     {
00557         std::cout << "- Executing NonGradient::iteration_complete()" << std::endl;
00558     }
00559     // MSQ_PRINT(2)("\n - Executing NonGradient::iteration_complete() \n");
00560 }
00561 
00562 void NonGradient::cleanup()
00563 {
00564     if( mNonGradDebug >= 2 )
00565     {
00566         std::cout << " - Executing NonGradient::iteration_end()" << std::endl;
00567     }
00568     // MSQ_PRINT(2)("\n - Executing NonGradient::iteration_end() \n");
00569 }
00570 
00571 }  // namespace MBMesquite
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines