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