MOAB: Mesh Oriented datABase  (version 5.3.1)
SteepestDescent.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   SteepestDescent.cpp
00029   \brief
00030 
00031   Implements the SteepestDescent class member functions.
00032 
00033   \author Thomas Leurent
00034   \date   2002-06-13
00035 */
00036 
00037 #include "SteepestDescent.hpp"
00038 #include "MsqFreeVertexIndexIterator.hpp"
00039 #include "MsqTimer.hpp"
00040 #include "MsqDebug.hpp"
00041 
00042 #include <memory>
00043 
00044 namespace MBMesquite
00045 {
00046 
00047 std::string SteepestDescent::get_name() const
00048 {
00049     return "SteepestDescent";
00050 }
00051 
00052 PatchSet* SteepestDescent::get_patch_set()
00053 {
00054     return PatchSetUser::get_patch_set();
00055 }
00056 
00057 SteepestDescent::SteepestDescent( ObjectiveFunction* of )
00058     : VertexMover( of ), PatchSetUser( true ), projectGradient( false )  //,
00059 // cosineStep(false)
00060 {
00061 }
00062 
00063 void SteepestDescent::initialize( PatchData& /*pd*/, MsqError& /*err*/ ) {}
00064 
00065 void SteepestDescent::initialize_mesh_iteration( PatchData& /*pd*/, MsqError& /*err*/ ) {}
00066 
00067 void SteepestDescent::optimize_vertex_positions( PatchData& pd, MsqError& err )
00068 {
00069     MSQ_FUNCTION_TIMER( "SteepestDescent::optimize_vertex_positions" );
00070 
00071     const int SEARCH_MAX = 100;
00072     const double c1      = 1e-4;
00073     // std::vector<Vector3D> unprojected(pd.num_free_vertices());
00074     std::vector< Vector3D > gradient( pd.num_free_vertices() );
00075     bool feasible = true;  // bool for OF values
00076     double min_edge_len, max_edge_len;
00077     double step_size = 0, original_value = 0, new_value = 0;
00078     double norm_squared = 0;
00079     PatchDataVerticesMemento* pd_previous_coords;
00080     TerminationCriterion* term_crit = get_inner_termination_criterion();
00081     OFEvaluator& obj_func           = get_objective_function_evaluator();
00082 
00083     // get vertex memento so we can restore vertex coordinates for bad steps.
00084     pd_previous_coords = pd.create_vertices_memento( err );MSQ_ERRRTN( err );
00085     // use unique_ptr to automatically delete memento when we exit this function
00086     std::unique_ptr< PatchDataVerticesMemento > memento_deleter( pd_previous_coords );
00087 
00088     // Evaluate objective function.
00089     //
00090     // Always use 'update' version when beginning optimization so that
00091     // if doing block coordinate descent the OF code knows the set of
00092     // vertices we are modifying during the optimziation (the subset
00093     // of the mesh contained in the current patch.)  This has to be
00094     // done up-front because typically an OF will just store the portion
00095     // of the OF value (e.g. the numeric contribution to the sum for an
00096     // averaging OF) for the initial patch.
00097     feasible = obj_func.update( pd, original_value, gradient, err );MSQ_ERRRTN( err );
00098     // calculate gradient dotted with itself
00099     norm_squared = length_squared( gradient );
00100 
00101     // set an error if initial patch is invalid.
00102     if( !feasible )
00103     {
00104         MSQ_SETERR( err )( "SteepestDescent passed invalid initial patch.", MsqError::INVALID_ARG );
00105         return;
00106     }
00107 
00108     // use edge length as an initial guess for for step size
00109     pd.get_minmax_edge_length( min_edge_len, max_edge_len );
00110     // step_size = max_edge_len / std::sqrt(norm_squared);
00111     // if (!finite(step_size))  // zero-length gradient
00112     //  return;
00113     //  if (norm_squared < DBL_EPSILON)
00114     //    return;
00115     if( norm_squared >= DBL_EPSILON ) step_size = max_edge_len / std::sqrt( norm_squared ) * pd.num_free_vertices();
00116 
00117     // The steepest descent loop...
00118     // We loop until the user-specified termination criteria are met.
00119     while( !term_crit->terminate() )
00120     {
00121         MSQ_DBGOUT( 3 ) << "Iteration " << term_crit->get_iteration_count() << std::endl;
00122         MSQ_DBGOUT( 3 ) << "  o  original_value: " << original_value << std::endl;
00123         MSQ_DBGOUT( 3 ) << "  o  grad norm suqared: " << norm_squared << std::endl;
00124 
00125         // Save current vertex coords so that they can be restored if
00126         // the step was bad.
00127         pd.recreate_vertices_memento( pd_previous_coords, err );MSQ_ERRRTN( err );
00128 
00129         // Reduce step size until it satisfies Armijo condition
00130         int counter = 0;
00131         for( ;; )
00132         {
00133             if( ++counter > SEARCH_MAX || step_size < DBL_EPSILON )
00134             {
00135                 MSQ_DBGOUT( 3 ) << "    o  No valid step found.  Giving Up." << std::endl;
00136                 return;
00137             }
00138 
00139             // Move vertices to new positions.
00140             // Note: step direction is -gradient so we pass +gradient and
00141             //       -step_size to achieve the same thing.
00142             pd.move_free_vertices_constrained( arrptr( gradient ), gradient.size(), -step_size, err );MSQ_ERRRTN( err );
00143             // Evaluate objective function for new vertices.  We call the
00144             // 'evaluate' form here because we aren't sure yet if we want to
00145             // keep these vertices.  Until we call 'update', we have the option
00146             // of reverting a block coordinate decent objective function's state
00147             // to that of the initial vertex coordinates.  However, for block
00148             // coordinate decent to work correctly, we will need to call an
00149             // 'update' form if we decide to keep the new vertex coordinates.
00150             feasible = obj_func.evaluate( pd, new_value, err );
00151             if( err.error_code() == err.BARRIER_VIOLATED )
00152                 err.clear();  // barrier violated does not represent an actual error here
00153             MSQ_ERRRTN( err );
00154             MSQ_DBGOUT( 3 ) << "    o  step_size: " << step_size << std::endl;
00155             MSQ_DBGOUT( 3 ) << "    o  new_value: " << new_value << std::endl;
00156 
00157             if( !feasible )
00158             {
00159                 // OF value is invalid, decrease step_size a lot
00160                 step_size *= 0.2;
00161             }
00162             else if( new_value > original_value - c1 * step_size * norm_squared )
00163             {
00164                 // Armijo condition not met.
00165                 step_size *= 0.5;
00166             }
00167             else
00168             {
00169                 // Armijo condition met, stop
00170                 break;
00171             }
00172 
00173             // undo previous step : restore vertex coordinates
00174             pd.set_to_vertices_memento( pd_previous_coords, err );MSQ_ERRRTN( err );
00175         }
00176 
00177         // Re-evaluate objective function to get gradient.
00178         // Calling the 'update' form here incorporates the new vertex
00179         // positions into the 'accumulated' value if we are doing a
00180         // block coordinate descent optimization.
00181         obj_func.update( pd, original_value, gradient, err );MSQ_ERRRTN( err );
00182         if( projectGradient )
00183         {
00184             // if (cosineStep) {
00185             //  unprojected = gradient;
00186             //  pd.project_gradient( gradient, err ); MSQ_ERRRTN(err);
00187             //  double dot = inner_product( arrptr(gradient), arrptr(unprojected), gradient.size()
00188             //  ); double lensqr1 = length_squared( gradient ); double lensqr2 = length_squared(
00189             //  unprojected ); double cossqr = dot * dot / lensqr1 / lensqr2; step_size *=
00190             //  sqrt(cossqr);
00191             //}
00192             // else {
00193             pd.project_gradient( gradient, err );MSQ_ERRRTN( err );
00194             //}
00195         }
00196 
00197         // Update terination criterion for next iteration.
00198         // This is necessary for efficiency.  Some values can be adjusted
00199         // for each iteration so we don't need to re-caculate the value
00200         // over the entire mesh.
00201         term_crit->accumulate_patch( pd, err );MSQ_ERRRTN( err );
00202         term_crit->accumulate_inner( pd, original_value, arrptr( gradient ), err );MSQ_ERRRTN( err );
00203 
00204         // Calculate initial step size for next iteration using step size
00205         // from this iteration
00206         step_size *= norm_squared;
00207         norm_squared = length_squared( gradient );
00208         //    if (norm_squared < DBL_EPSILON)
00209         //      break;
00210         if( norm_squared >= DBL_EPSILON ) step_size /= norm_squared;
00211     }
00212 }
00213 
00214 void SteepestDescent::terminate_mesh_iteration( PatchData& /*pd*/, MsqError& /*err*/ )
00215 {
00216     //  cout << "- Executing SteepestDescent::iteration_complete()\n";
00217 }
00218 
00219 void SteepestDescent::cleanup()
00220 {
00221     //  cout << "- Executing SteepestDescent::iteration_end()\n";
00222 }
00223 
00224 }  // namespace MBMesquite
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines