MOAB: Mesh Oriented datABase  (version 5.2.1)
ObjectiveFunction.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   ObjectiveFunction.cpp
00029   \brief
00030 
00031   \author Michael Brewer
00032   \author Thomas Leurent
00033 
00034   \date   2002-08-02
00035 */
00036 
00037 #include "ObjectiveFunction.hpp"
00038 #include "MsqVertex.hpp"
00039 #include "MsqDebug.hpp"
00040 #include "PatchData.hpp"
00041 #include "MsqError.hpp"
00042 #include "MsqHessian.hpp"
00043 #include "SymMatrix3D.hpp"
00044 #include <memory>  // unique_ptr
00045 
00046 namespace MBMesquite
00047 {
00048 
00049 ObjectiveFunction::~ObjectiveFunction() {}
00050 
00051 /*!Returns an appropiate value (eps) to use as a delta step for
00052   MsqVertex vertex in dimension k (i.e. k=0 -> x, k=1 -> y, k=2 -> z).
00053   The objective function value at the perturbed vertex position is given
00054   in local_val.
00055 */
00056 double ObjectiveFunction::get_eps( PatchData& pd, EvalType type, double& local_val, int dim, size_t vertex_index,
00057                                    MsqError& err )
00058 {
00059     double eps       = 1.e-07;
00060     const double rho = 0.5;
00061     const int imax   = 20;
00062     bool feasible    = false;
00063     double tmp_var   = 0.0;
00064     Vector3D delta( 0, 0, 0 );
00065     for( int i = 0; i < imax; ++i )
00066     {
00067         // perturb kth coord val and check feas if needed
00068         tmp_var    = pd.vertex_by_index( vertex_index )[dim];
00069         delta[dim] = eps;
00070         pd.move_vertex( delta, vertex_index, err );
00071         feasible = evaluate( type, pd, local_val, OF_FREE_EVALS_ONLY, err );
00072         MSQ_ERRZERO( err );
00073         // revert kth coord val
00074         delta      = pd.vertex_by_index( vertex_index );
00075         delta[dim] = tmp_var;
00076         pd.set_vertex_coordinates( delta, vertex_index, err );
00077         // if step was too big, shorten it and go again
00078         if( feasible )
00079             return eps;
00080         else
00081             eps *= rho;
00082     }  // end while looking for feasible eps
00083     return 0.0;
00084 }  // end function get_eps
00085 
00086 bool ObjectiveFunction::compute_subpatch_numerical_gradient( EvalType type, EvalType subtype, PatchData& pd,
00087                                                              double& flocal, Vector3D& grad, MsqError& err )
00088 {
00089     assert( pd.num_free_vertices() == 1 );
00090 
00091     double flocald = 0;
00092     double eps     = 0;
00093 
00094     bool b = evaluate( type, pd, flocal, OF_FREE_EVALS_ONLY, err );
00095     if( MSQ_CHKERR( err ) || !b ) { return false; }
00096 
00097     // loop over the three coords x,y,z
00098     for( int j = 0; j < 3; ++j )
00099     {
00100         eps = get_eps( pd, subtype, flocald, j, 0, err );
00101         MSQ_ERRZERO( err );
00102         if( eps == 0 )
00103         {
00104             MSQ_SETERR( err )
00105             ( "Dividing by zero in Objective Functions numerical grad", MsqError::INVALID_STATE );
00106             return false;
00107         }
00108         grad[j] = ( flocald - flocal ) / eps;
00109     }
00110     return true;
00111 }
00112 
00113 bool ObjectiveFunction::compute_patch_numerical_gradient( EvalType type, EvalType subtype, PatchData& pd,
00114                                                           double& flocal, std::vector< Vector3D >& grad, MsqError& err )
00115 {
00116     double flocald = 0;
00117     double eps     = 0;
00118 
00119     bool b = evaluate( type, pd, flocal, OF_FREE_EVALS_ONLY, err );
00120     if( MSQ_CHKERR( err ) || !b ) { return false; }
00121 
00122     for( size_t i = 0; i < pd.num_free_vertices(); ++i )
00123     {
00124         // loop over the three coords x,y,z
00125         for( int j = 0; j < 3; ++j )
00126         {
00127             eps = get_eps( pd, subtype, flocald, j, i, err );
00128             MSQ_ERRZERO( err );
00129             if( eps == 0 )
00130             {
00131                 MSQ_SETERR( err )
00132                 ( "Dividing by zero in Objective Functions numerical grad", MsqError::INVALID_STATE );
00133                 return false;
00134             }
00135             grad[i][j] = ( flocald - flocal ) / eps;
00136         }
00137     }
00138 
00139     return true;
00140 }
00141 
00142 /*!
00143   Numerically Calculates the gradient of the ObjectiveFunction for the
00144   free vertices in the patch.  Returns 'false' if the patch is outside
00145   of a required feasible region, returns 'ture' otherwise.
00146   The behavior of the function depends on the value of the boolean
00147   useLocalGradient.  If useLocalGradient is set to
00148   'true', compute_numerical_gradient creates a sub-patch around a free
00149   vertex, and then perturbs that vertex in one of the coordinate directions.
00150   Only the ObjectiveFunction value on the local sub-patch is used in the
00151   computation of the gradient.  Therefore, useLocalGradient should only
00152   be set to 'true' for ObjectiveFunctions which can use this method.  Unless
00153   the concrete ObjectiveFunction sets useLocalGradient to 'true' in its
00154   constructor, the value will be 'false'.  In this case, the objective
00155   function value for the entire patch is used in the calculation of the
00156   gradient.  This is computationally expensive, but it is numerically
00157   correct for all (C_1) functions.
00158   \param pd  PatchData on which the gradient is taken.
00159   \param grad  Array of Vector3D of length the number of vertices used to store gradient.
00160   \param OF_val will be set to the objective function value.
00161  */
00162 bool ObjectiveFunction::evaluate_with_gradient( EvalType eval_type, PatchData& pd, double& OF_val,
00163                                                 std::vector< Vector3D >& grad, MsqError& err )
00164 {
00165     bool b;
00166     grad.resize( pd.num_free_vertices() );
00167 
00168     // Fast path for single-free-vertex patch
00169     if( pd.num_free_vertices() == 1 )
00170     {
00171         const EvalType sub_type = ( eval_type == CALCULATE ) ? CALCULATE : TEMPORARY;
00172         b                       = compute_subpatch_numerical_gradient( eval_type, sub_type, pd, OF_val, grad[0], err );
00173         return !MSQ_CHKERR( err ) && b;
00174     }
00175 
00176     ObjectiveFunction* of = this;
00177     std::unique_ptr< ObjectiveFunction > deleter;
00178     if( eval_type == CALCULATE )
00179     {
00180         of->clear();
00181         b = of->evaluate( ACCUMULATE, pd, OF_val, OF_FREE_EVALS_ONLY, err );
00182         if( err )
00183         {  // OF doesn't support BCD type evals, try slow method
00184             err.clear();
00185             of->clear();
00186             b = compute_patch_numerical_gradient( CALCULATE, CALCULATE, pd, OF_val, grad, err );
00187             return !MSQ_CHKERR( err ) && b;
00188         }
00189         else if( !b )
00190             return b;
00191     }
00192     else
00193     {
00194         b = this->evaluate( eval_type, pd, OF_val, OF_FREE_EVALS_ONLY, err );
00195         if( MSQ_CHKERR( err ) || !b ) return false;
00196         of      = this->clone();
00197         deleter = std::unique_ptr< ObjectiveFunction >( of );
00198     }
00199 
00200     // Determine number of layers of adjacent elements based on metric type.
00201     unsigned layers = min_patch_layers();
00202 
00203     // Create a subpatch for each free vertex and use it to evaluate the
00204     // gradient for that vertex.
00205     double flocal;
00206     PatchData subpatch;
00207     for( size_t i = 0; i < pd.num_free_vertices(); ++i )
00208     {
00209         pd.get_subpatch( i, layers, subpatch, err );
00210         MSQ_ERRZERO( err );
00211         b = of->compute_subpatch_numerical_gradient( SAVE, TEMPORARY, subpatch, flocal, grad[i], err );
00212         if( MSQ_CHKERR( err ) || !b )
00213         {
00214             of->clear();
00215             return false;
00216         }
00217     }
00218 
00219     of->clear();
00220     return true;
00221 }
00222 
00223 bool ObjectiveFunction::evaluate_with_Hessian_diagonal( EvalType type, PatchData& pd, double& value_out,
00224                                                         std::vector< Vector3D >& grad_out,
00225                                                         std::vector< SymMatrix3D >& hess_diag_out, MsqError& err )
00226 {
00227     MsqHessian hess;
00228     hess.initialize( pd, err );
00229     MSQ_ERRZERO( err );
00230     bool val = evaluate_with_Hessian( type, pd, value_out, grad_out, hess, err );
00231     MSQ_ERRZERO( err );
00232     hess_diag_out.resize( hess.size() );
00233     for( size_t i = 0; i < hess.size(); ++i )
00234         hess_diag_out[i] = hess.get_block( i, i )->upper();
00235     return val;
00236 }
00237 
00238 bool ObjectiveFunction::evaluate_with_Hessian( EvalType /*type*/, PatchData& /*pd*/, double& /*value_out*/,
00239                                                std::vector< Vector3D >& /*grad_out*/, MsqHessian& /*Hessian_out*/,
00240                                                MsqError& err )
00241 {
00242     MSQ_SETERR( err )
00243     ( "No Hessian available for this objective function.\n"
00244       "Choose either a different objective function or a "
00245       "different solver.\n",
00246       MsqError::INVALID_STATE );
00247     return false;
00248 }
00249 
00250 }  // namespace MBMesquite
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines