MOAB: Mesh Oriented datABase  (version 5.2.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) kraftche@cae.wisc.edu
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, Vector3D* r,
00143                                         MsqError&
00144 #ifndef USE_FN_PC1
00145                                             err
00146 #endif
00147 )
00148 {
00149 #ifndef USE_FN_PC1
00150     mHessian.apply_preconditioner( z, r, err );
00151 #else
00152     for( size_t i = 0; i < preCond.size(); ++i )
00153         z[i] = preCond[i] * r[i];
00154 #endif
00155 }
00156 
00157 void TrustRegion::optimize_vertex_positions( PatchData& pd, MsqError& err )
00158 {
00159     TerminationCriterion& term = *get_inner_termination_criterion();
00160     OFEvaluator& func          = get_objective_function_evaluator();
00161 
00162     const double cg_tol        = 1e-2;
00163     const double eta_1         = 0.01;
00164     const double eta_2         = 0.90;
00165     const double tr_incr       = 10;
00166     const double tr_decr_def   = 0.25;
00167     const double tr_decr_undef = 0.25;
00168     const double tr_num_tol    = 1e-6;
00169     const int max_cg_iter      = 10000;
00170 
00171     double radius = 1000; /* delta*delta */
00172 
00173     const int nn = pd.num_free_vertices();
00174     wVect.resize( nn );
00175     Vector3D* w = arrptr( wVect );
00176     zVect.resize( nn );
00177     Vector3D* z = arrptr( zVect );
00178     dVect.resize( nn );
00179     Vector3D* d = arrptr( dVect );
00180     pVect.resize( nn );
00181     Vector3D* p = arrptr( pVect );
00182     rVect.resize( nn );
00183     Vector3D* r = arrptr( rVect );
00184 
00185     double norm_r, norm_g;
00186     double alpha, beta, kappa;
00187     double rz, rzm1;
00188     double dMp, norm_d, norm_dp1, norm_p;
00189     double obj, objn;
00190 
00191     int cg_iter;
00192     bool valid;
00193 
00194     mHess.initialize( pd, err );  // hMesh(mesh);
00195     valid = func.update( pd, obj, mGrad, mHess, err );MSQ_ERRRTN( err );
00196     if( !valid )
00197     {
00198         MSQ_SETERR( err )( "Initial objective function is not valid", MsqError::INVALID_MESH );
00199         return;
00200     }
00201     compute_preconditioner( err );MSQ_ERRRTN( err );
00202     pd.recreate_vertices_memento( mMemento, err );MSQ_ERRRTN( err );
00203 
00204     while( !term.terminate() && ( radius > 1e-20 ) )
00205     {
00206 
00207         norm_r = length_squared( arrptr( mGrad ), nn );
00208         norm_g = sqrt( norm_r );
00209 
00210         memset( d, 0, 3 * sizeof( double ) * nn );
00211         memcpy( r, arrptr( mGrad ),
00212                 nn * sizeof( Vector3D ) );  // memcpy(r, mesh->g, 3*sizeof(double)*nn);
00213         norm_g *= cg_tol;
00214 
00215         apply_preconditioner( z, r, err );MSQ_ERRRTN( err );  // prec->apply(z, r, prec, mesh);
00216         negate( p, z, nn );
00217         rz = inner( r, z, nn );
00218 
00219         dMp    = 0;
00220         norm_p = rz;
00221         norm_d = 0;
00222 
00223         cg_iter = 0;
00224         while( ( sqrt( norm_r ) > norm_g ) &&
00225 #ifdef DO_STEEP_DESC
00226                ( norm_d > tr_num_tol ) &&
00227 #endif
00228                ( cg_iter < max_cg_iter ) )
00229         {
00230             ++cg_iter;
00231 
00232             memset( w, 0, 3 * sizeof( double ) * nn );
00233             // matmul(w, mHess, p); //matmul(w, mesh, p);
00234             mHess.product( w, p );
00235 
00236             kappa = inner( p, w, nn );
00237             if( kappa <= 0.0 )
00238             {
00239                 alpha = ( sqrt( dMp * dMp + norm_p * ( radius - norm_d ) ) - dMp ) / norm_p;
00240                 plus_eq_scaled( d, alpha, p, nn );
00241                 break;
00242             }
00243 
00244             alpha = rz / kappa;
00245 
00246             norm_dp1 = norm_d + 2.0 * alpha * dMp + alpha * alpha * norm_p;
00247             if( norm_dp1 >= radius )
00248             {
00249                 alpha = ( sqrt( dMp * dMp + norm_p * ( radius - norm_d ) ) - dMp ) / norm_p;
00250                 plus_eq_scaled( d, alpha, p, nn );
00251                 break;
00252             }
00253 
00254             plus_eq_scaled( d, alpha, p, nn );
00255             plus_eq_scaled( r, alpha, w, nn );
00256             norm_r = length_squared( r, nn );
00257 
00258             apply_preconditioner( z, r, err );MSQ_ERRRTN( err );  // prec->apply(z, r, prec, mesh);
00259 
00260             rzm1 = rz;
00261             rz   = inner( r, z, nn );
00262             beta = rz / rzm1;
00263             times_eq_minus( p, beta, z, nn );
00264 
00265             dMp    = beta * ( dMp + alpha * norm_p );
00266             norm_p = rz + beta * beta * norm_p;
00267             norm_d = norm_dp1;
00268         }
00269 
00270 #ifdef DO_STEEP_DESC
00271         if( norm_d <= tr_num_tol )
00272         {
00273             norm_g    = length( arrptr( mGrad ), nn );
00274             double ll = 1.0;
00275             if( norm_g < tr_num_tol ) break;
00276             if( norm_g > radius ) ll = radius / nurm_g;
00277             for( int i = 0; i < nn; ++i )
00278                 d[i] = ll * mGrad[i];
00279         }
00280 #endif
00281 
00282         alpha = inner( arrptr( mGrad ), d, nn );  // inner(mesh->g, d, nn);
00283 
00284         memset( p, 0, 3 * sizeof( double ) * nn );
00285         // matmul(p, mHess, d); //matmul(p, mesh, d);
00286         mHess.product( p, d );
00287         beta  = 0.5 * inner( p, d, nn );
00288         kappa = alpha + beta;
00289 
00290         /* Put the new point into the locations */
00291         pd.move_free_vertices_constrained( d, nn, 1.0, err );MSQ_ERRRTN( err );
00292 
00293         valid = func.evaluate( pd, objn, err );
00294         if( err.error_code() == err.BARRIER_VIOLATED )
00295             err.clear();  // barrier violated does not represent an actual error here
00296         MSQ_ERRRTN( err );
00297 
00298         if( !valid )
00299         {
00300             /* Function not defined at trial point */
00301             radius *= tr_decr_undef;
00302             pd.set_to_vertices_memento( mMemento, err );MSQ_ERRRTN( err );
00303             continue;
00304         }
00305 
00306         if( ( fabs( kappa ) <= tr_num_tol ) && ( fabs( objn - obj ) <= tr_num_tol ) ) { kappa = 1; }
00307         else
00308         {
00309             kappa = ( objn - obj ) / kappa;
00310         }
00311 
00312         if( kappa < eta_1 )
00313         {
00314             /* Iterate is unacceptable */
00315             radius *= tr_decr_def;
00316             pd.set_to_vertices_memento( mMemento, err );MSQ_ERRRTN( err );
00317             continue;
00318         }
00319 
00320         /* Iterate is acceptable */
00321         if( kappa >= eta_2 )
00322         {
00323             /* Iterate is a very good step, increase radius */
00324             radius *= tr_incr;
00325             if( radius > 1e20 ) { radius = 1e20; }
00326         }
00327 
00328         func.update( pd, obj, mGrad, mHess, err );
00329         compute_preconditioner( err );MSQ_ERRRTN( err );
00330         pd.recreate_vertices_memento( mMemento, err );MSQ_ERRRTN( err );
00331 
00332         // checks stopping criterion
00333         term.accumulate_patch( pd, err );MSQ_ERRRTN( err );
00334         term.accumulate_inner( pd, objn, arrptr( mGrad ), err );MSQ_ERRRTN( err );
00335     }
00336 }
00337 
00338 }  // namespace MBMesquite
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines