MOAB: Mesh Oriented datABase  (version 5.4.1)
TrustRegion.cpp
Go to the documentation of this file.
00001 /* *****************************************************************
00002     MESQUITE -- The Mesh Quality Improvement Toolkit
00003 
00004     Copyright 2007 Sandia National Laboratories.  Developed at the
00005     University of Wisconsin--Madison under SNL contract number
00006     624796.  The U.S. Government and the University of Wisconsin
00007     retain certain rights to 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     (2008) [email protected]
00024 
00025   ***************************************************************** */
00026 
00027 /** \file TrustRegion.cpp
00028  *  \brief Port Todd Munson's trust region solver to Mesquite
00029  *  \author Jason Kraftcheck (Mesquite Port)
00030  */
00031 
00032 #include "Mesquite.hpp"
00033 #include "TrustRegion.hpp"
00034 #include "MsqDebug.hpp"
00035 #include "MsqError.hpp"
00036 #include "PatchData.hpp"
00037 
00038 #define USE_FN_PC1    // Use 1st preconditioner from Todd's code \
00039                       // (alternate is whatever is in MsqHessian already)
00040 #undef DO_STEEP_DESC  // Jason's apparently broken hack to fall back to \
00041                       // steepest descent search direction
00042 
00043 namespace MBMesquite
00044 {
00045 
00046 // Force std::vector to release allocated memory
00047 template < typename T >
00048 static inline void free_vector( std::vector< T >& v )
00049 {
00050     std::vector< T > temp;
00051     temp.swap( v );
00052 }
00053 
00054 std::string TrustRegion::get_name() const
00055 {
00056     return "TrustRegion";
00057 }
00058 
00059 PatchSet* TrustRegion::get_patch_set()
00060 {
00061     return PatchSetUser::get_patch_set();
00062 }
00063 
00064 TrustRegion::TrustRegion( ObjectiveFunction* of ) : VertexMover( of ), PatchSetUser( true ), mMemento( 0 ) {}
00065 
00066 TrustRegion::~TrustRegion()
00067 {
00068     delete mMemento;
00069     mMemento = 0;
00070 }
00071 
00072 void TrustRegion::initialize( PatchData& pd, MsqError& err )
00073 {
00074     mMemento = pd.create_vertices_memento( err );MSQ_CHKERR( err );
00075 }
00076 
00077 void TrustRegion::initialize_mesh_iteration( PatchData& /*pd*/, MsqError& /*err*/ ) {}
00078 
00079 void TrustRegion::terminate_mesh_iteration( PatchData& /*pd*/, MsqError& /*err*/ ) {}
00080 
00081 void TrustRegion::cleanup()
00082 {
00083     // release Memento
00084     delete mMemento;
00085     mMemento = 0;
00086     // release temporary array memory
00087     mHess.clear();
00088     free_vector( mGrad );
00089     free_vector( wVect );
00090     free_vector( zVect );
00091     free_vector( dVect );
00092     free_vector( pVect );
00093     free_vector( rVect );
00094     free_vector( preCond );
00095 }
00096 
00097 static inline void negate( Vector3D* out, const Vector3D* in, size_t nn )
00098 {
00099     for( size_t i = 0; i < nn; ++i )
00100         out[i] = -in[i];
00101 }
00102 
00103 // Do v += s * x, where v and x are arrays of length n
00104 static inline void plus_eq_scaled( Vector3D* v, double s, const Vector3D* x, size_t n )
00105 {
00106     Vector3D* end = v + n;
00107     for( ; v != end; ++v, ++x )
00108         *v += s * *x;
00109 }
00110 
00111 // Do v = s*v - x, where v and x are arrays of length n
00112 static inline void times_eq_minus( Vector3D* v, double s, const Vector3D* x, size_t n )
00113 {
00114     Vector3D* end = v + n;
00115     for( ; v != end; ++v, ++x )
00116     {
00117         *v *= s;
00118         *v -= *x;
00119     }
00120 }
00121 
00122 void TrustRegion::compute_preconditioner( MsqError&
00123 #ifndef USE_FN_PC1
00124                                               err
00125 #endif
00126 )
00127 {
00128 #ifndef USE_FN_PC1
00129     mHessian.calculate_preconditioner( err );
00130 #else
00131     double dia;
00132     preCond.resize( mHess.size() );
00133     for( size_t i = 0; i < mHess.size(); ++i )
00134     {
00135         const Matrix3D& m = *mHess.get_block( i, i );
00136         dia               = m[0][0] + m[1][1] + m[2][2];
00137         preCond[i]        = dia < DBL_EPSILON ? 1.0 : 1.0 / dia;
00138     }
00139 #endif
00140 }
00141 
00142 void TrustRegion::apply_preconditioner( Vector3D* z,
00143                                         Vector3D* r,
00144                                         MsqError&
00145 #ifndef USE_FN_PC1
00146                                             err
00147 #endif
00148 )
00149 {
00150 #ifndef USE_FN_PC1
00151     mHessian.apply_preconditioner( z, r, err );
00152 #else
00153     for( size_t i = 0; i < preCond.size(); ++i )
00154         z[i] = preCond[i] * r[i];
00155 #endif
00156 }
00157 
00158 void TrustRegion::optimize_vertex_positions( PatchData& pd, MsqError& err )
00159 {
00160     TerminationCriterion& term = *get_inner_termination_criterion();
00161     OFEvaluator& func          = get_objective_function_evaluator();
00162 
00163     const double cg_tol        = 1e-2;
00164     const double eta_1         = 0.01;
00165     const double eta_2         = 0.90;
00166     const double tr_incr       = 10;
00167     const double tr_decr_def   = 0.25;
00168     const double tr_decr_undef = 0.25;
00169     const double tr_num_tol    = 1e-6;
00170     const int max_cg_iter      = 10000;
00171 
00172     double radius = 1000; /* delta*delta */
00173 
00174     const int nn = pd.num_free_vertices();
00175     wVect.resize( nn );
00176     Vector3D* w = arrptr( wVect );
00177     zVect.resize( nn );
00178     Vector3D* z = arrptr( zVect );
00179     dVect.resize( nn );
00180     Vector3D* d = arrptr( dVect );
00181     pVect.resize( nn );
00182     Vector3D* p = arrptr( pVect );
00183     rVect.resize( nn );
00184     Vector3D* r = arrptr( rVect );
00185 
00186     double norm_r, norm_g;
00187     double alpha, beta, kappa;
00188     double rz, rzm1;
00189     double dMp, norm_d, norm_dp1, norm_p;
00190     double obj, objn;
00191 
00192     int cg_iter;
00193     bool valid;
00194 
00195     mHess.initialize( pd, err );  // hMesh(mesh);
00196     valid = func.update( pd, obj, mGrad, mHess, err );MSQ_ERRRTN( err );
00197     if( !valid )
00198     {
00199         MSQ_SETERR( err )( "Initial objective function is not valid", MsqError::INVALID_MESH );
00200         return;
00201     }
00202     compute_preconditioner( err );MSQ_ERRRTN( err );
00203     pd.recreate_vertices_memento( mMemento, err );MSQ_ERRRTN( err );
00204 
00205     while( !term.terminate() && ( radius > 1e-20 ) )
00206     {
00207 
00208         norm_r = length_squared( arrptr( mGrad ), nn );
00209         norm_g = sqrt( norm_r );
00210 
00211         memset( d, 0, 3 * sizeof( double ) * nn );
00212         memcpy( r, arrptr( mGrad ),
00213                 nn * sizeof( Vector3D ) );  // memcpy(r, mesh->g, 3*sizeof(double)*nn);
00214         norm_g *= cg_tol;
00215 
00216         apply_preconditioner( z, r, err );MSQ_ERRRTN( err );  // prec->apply(z, r, prec, mesh);
00217         negate( p, z, nn );
00218         rz = inner( r, z, nn );
00219 
00220         dMp    = 0;
00221         norm_p = rz;
00222         norm_d = 0;
00223 
00224         cg_iter = 0;
00225         while( ( sqrt( norm_r ) > norm_g ) &&
00226 #ifdef DO_STEEP_DESC
00227                ( norm_d > tr_num_tol ) &&
00228 #endif
00229                ( cg_iter < max_cg_iter ) )
00230         {
00231             ++cg_iter;
00232 
00233             memset( w, 0, 3 * sizeof( double ) * nn );
00234             // matmul(w, mHess, p); //matmul(w, mesh, p);
00235             mHess.product( w, p );
00236 
00237             kappa = inner( p, w, nn );
00238             if( kappa <= 0.0 )
00239             {
00240                 alpha = ( sqrt( dMp * dMp + norm_p * ( radius - norm_d ) ) - dMp ) / norm_p;
00241                 plus_eq_scaled( d, alpha, p, nn );
00242                 break;
00243             }
00244 
00245             alpha = rz / kappa;
00246 
00247             norm_dp1 = norm_d + 2.0 * alpha * dMp + alpha * alpha * norm_p;
00248             if( norm_dp1 >= radius )
00249             {
00250                 alpha = ( sqrt( dMp * dMp + norm_p * ( radius - norm_d ) ) - dMp ) / norm_p;
00251                 plus_eq_scaled( d, alpha, p, nn );
00252                 break;
00253             }
00254 
00255             plus_eq_scaled( d, alpha, p, nn );
00256             plus_eq_scaled( r, alpha, w, nn );
00257             norm_r = length_squared( r, nn );
00258 
00259             apply_preconditioner( z, r, err );MSQ_ERRRTN( err );  // prec->apply(z, r, prec, mesh);
00260 
00261             rzm1 = rz;
00262             rz   = inner( r, z, nn );
00263             beta = rz / rzm1;
00264             times_eq_minus( p, beta, z, nn );
00265 
00266             dMp    = beta * ( dMp + alpha * norm_p );
00267             norm_p = rz + beta * beta * norm_p;
00268             norm_d = norm_dp1;
00269         }
00270 
00271 #ifdef DO_STEEP_DESC
00272         if( norm_d <= tr_num_tol )
00273         {
00274             norm_g    = length( arrptr( mGrad ), nn );
00275             double ll = 1.0;
00276             if( norm_g < tr_num_tol ) break;
00277             if( norm_g > radius ) ll = radius / nurm_g;
00278             for( int i = 0; i < nn; ++i )
00279                 d[i] = ll * mGrad[i];
00280         }
00281 #endif
00282 
00283         alpha = inner( arrptr( mGrad ), d, nn );  // inner(mesh->g, d, nn);
00284 
00285         memset( p, 0, 3 * sizeof( double ) * nn );
00286         // matmul(p, mHess, d); //matmul(p, mesh, d);
00287         mHess.product( p, d );
00288         beta  = 0.5 * inner( p, d, nn );
00289         kappa = alpha + beta;
00290 
00291         /* Put the new point into the locations */
00292         pd.move_free_vertices_constrained( d, nn, 1.0, err );MSQ_ERRRTN( err );
00293 
00294         valid = func.evaluate( pd, objn, err );
00295         if( err.error_code() == err.BARRIER_VIOLATED )
00296             err.clear();  // barrier violated does not represent an actual error here
00297         MSQ_ERRRTN( err );
00298 
00299         if( !valid )
00300         {
00301             /* Function not defined at trial point */
00302             radius *= tr_decr_undef;
00303             pd.set_to_vertices_memento( mMemento, err );MSQ_ERRRTN( err );
00304             continue;
00305         }
00306 
00307         if( ( fabs( kappa ) <= tr_num_tol ) && ( fabs( objn - obj ) <= tr_num_tol ) )
00308         {
00309             kappa = 1;
00310         }
00311         else
00312         {
00313             kappa = ( objn - obj ) / kappa;
00314         }
00315 
00316         if( kappa < eta_1 )
00317         {
00318             /* Iterate is unacceptable */
00319             radius *= tr_decr_def;
00320             pd.set_to_vertices_memento( mMemento, err );MSQ_ERRRTN( err );
00321             continue;
00322         }
00323 
00324         /* Iterate is acceptable */
00325         if( kappa >= eta_2 )
00326         {
00327             /* Iterate is a very good step, increase radius */
00328             radius *= tr_incr;
00329             if( radius > 1e20 )
00330             {
00331                 radius = 1e20;
00332             }
00333         }
00334 
00335         func.update( pd, obj, mGrad, mHess, err );
00336         compute_preconditioner( err );MSQ_ERRRTN( err );
00337         pd.recreate_vertices_memento( mMemento, err );MSQ_ERRRTN( err );
00338 
00339         // checks stopping criterion
00340         term.accumulate_patch( pd, err );MSQ_ERRRTN( err );
00341         term.accumulate_inner( pd, objn, arrptr( mGrad ), err );MSQ_ERRRTN( err );
00342     }
00343 }
00344 
00345 }  // namespace MBMesquite
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines