MOAB: Mesh Oriented datABase  (version 5.2.1)
ConjugateGradient.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   ConjugateGradient.cpp
00029   \brief
00030 
00031   The Conjugate Gradient class is a concrete vertex mover
00032   which performs conjugate gradient minimizaiton.
00033 
00034   \author Michael Brewer
00035   \date   2002-06-19
00036 */
00037 
00038 #include "ConjugateGradient.hpp"
00039 #include <cmath>
00040 #include "MsqDebug.hpp"
00041 #include "MsqTimer.hpp"
00042 //#include "MsqFreeVertexIndexIterator.hpp"
00043 
00044 namespace MBMesquite
00045 {
00046 
00047 extern int get_parallel_rank();
00048 extern int get_parallel_size();
00049 
00050 std::string ConjugateGradient::get_name() const
00051 {
00052     return "ConjugateGradient";
00053 }
00054 
00055 PatchSet* ConjugateGradient::get_patch_set()
00056 {
00057     return PatchSetUser::get_patch_set();
00058 }
00059 
00060 ConjugateGradient::ConjugateGradient( ObjectiveFunction* objective )
00061     : VertexMover( objective ), PatchSetUser( true ), pMemento( NULL ), conjGradDebug( 0 )
00062 {
00063 }
00064 
00065 ConjugateGradient::ConjugateGradient( ObjectiveFunction* objective, MsqError& err )
00066     : VertexMover( objective ), PatchSetUser( true ), pMemento( NULL ), conjGradDebug( 0 )
00067 {
00068     // Michael:: default to global?
00069     set_debugging_level( 0 );
00070     // set the default inner termination criterion
00071     TerminationCriterion* default_crit = get_inner_termination_criterion();
00072     if( default_crit == NULL )
00073     {
00074         MSQ_SETERR( err )
00075         ( "QualityImprover did not create a default inner "
00076           "termination criterion.",
00077           MsqError::INVALID_STATE );
00078         return;
00079     }
00080     else
00081     {
00082         default_crit->add_iteration_limit( 5 );MSQ_ERRRTN( err );
00083     }
00084 }
00085 
00086 ConjugateGradient::~ConjugateGradient()
00087 {
00088     // Checks that cleanup() has been called.
00089     assert( pMemento == NULL );
00090 }
00091 
00092 void ConjugateGradient::initialize( PatchData& pd, MsqError& err )
00093 {
00094     if( get_parallel_size() )
00095     {
00096         MSQ_DBGOUT( 2 ) << "\nP[" << get_parallel_rank() << "] "
00097                         << "o   Performing Conjugate Gradient optimization.\n";
00098     }
00099     else
00100     {
00101         MSQ_DBGOUT( 2 ) << "\no   Performing Conjugate Gradient optimization.\n";
00102     }
00103     pMemento = pd.create_vertices_memento( err );
00104 }
00105 
00106 void ConjugateGradient::initialize_mesh_iteration( PatchData& /*pd*/, MsqError& /*err*/ ) {}
00107 
00108 /*!Performs Conjugate gradient minimization on the PatchData, pd.*/
00109 void ConjugateGradient::optimize_vertex_positions( PatchData& pd, MsqError& err )
00110 {
00111     // pd.reorder();
00112 
00113     MSQ_FUNCTION_TIMER( "ConjugateGradient::optimize_vertex_positions" );
00114 
00115     Timer c_timer;
00116     size_t num_vert = pd.num_free_vertices();
00117     if( num_vert < 1 )
00118     {
00119         MSQ_DBGOUT( 1 ) << "\nEmpty free vertex list in ConjugateGradient\n";
00120         return;
00121     }
00122     /*
00123         //zero out arrays
00124       int zero_loop=0;
00125       while(zero_loop<arraySize){
00126         fGrad[zero_loop].set(0,0,0);
00127         pGrad[zero_loop].set(0,0,0);
00128         fNewGrad[zero_loop].set(0,0,0);
00129         ++zero_loop;
00130       }
00131     */
00132 
00133     // get OF evaluator
00134     OFEvaluator& objFunc = get_objective_function_evaluator();
00135 
00136     size_t ind;
00137     // Michael cull list:  possibly set soft_fixed flags here
00138 
00139     // MsqFreeVertexIndexIterator free_iter(pd, err);  MSQ_ERRRTN(err);
00140 
00141     double f = 0;
00142     // Michael, this isn't equivalent to CUBIT because we only want to check
00143     // the objective function value of the 'bad' elements
00144     // if invalid initial patch set an error.
00145     bool temp_bool = objFunc.update( pd, f, fGrad, err );
00146     assert( fGrad.size() == num_vert );
00147     if( MSQ_CHKERR( err ) ) return;
00148     if( !temp_bool )
00149     {
00150         MSQ_SETERR( err )
00151         ( "Conjugate Gradient not able to get valid gradient "
00152           "and function values on intial patch.",
00153           MsqError::INVALID_MESH );
00154         return;
00155     }
00156     double grad_norm = MSQ_MAX_CAP;
00157 
00158     if( conjGradDebug > 0 )
00159     {
00160         MSQ_PRINT( 2 )( "\nCG's DEGUB LEVEL = %i \n", conjGradDebug );
00161         grad_norm = Linf( arrptr( fGrad ), fGrad.size() );
00162         MSQ_PRINT( 2 )( "\nCG's FIRST VALUE = %f,grad_norm = %f", f, grad_norm );
00163         MSQ_PRINT( 2 )( "\n   TIME %f", c_timer.since_birth() );
00164         grad_norm = MSQ_MAX_CAP;
00165     }
00166 
00167     // Initializing pGrad (search direction).
00168     pGrad.resize( fGrad.size() );
00169     for( ind = 0; ind < num_vert; ++ind )
00170         pGrad[ind] = ( -fGrad[ind] );
00171 
00172     int j      = 0;            // total nb of step size changes ... not used much
00173     int i      = 0;            // iteration counter
00174     unsigned m = 0;            //
00175     double alp = MSQ_MAX_CAP;  // alp: scale factor of search direction
00176                                // we know inner_criterion is false because it was checked in
00177                                // loop_over_mesh before being sent here.
00178     TerminationCriterion* term_crit = get_inner_termination_criterion();
00179 
00180     // while ((i<maxIteration && alp>stepBound && grad_norm>normGradientBound)
00181     //     && !inner_criterion){
00182     while( !term_crit->terminate() )
00183     {
00184         ++i;
00185         // std::cout<<"\Michael delete i = "<<i;
00186         int k = 0;
00187         alp   = get_step( pd, f, k, err );
00188         j += k;
00189         if( conjGradDebug > 2 ) { MSQ_PRINT( 2 )( "\n  Alp initial, alp = %20.18f", alp ); }
00190 
00191         // if alp == 0, revert to steepest descent search direction
00192         if( alp == 0 )
00193         {
00194             for( m = 0; m < num_vert; ++m )
00195             {
00196                 pGrad[m] = ( -fGrad[m] );
00197             }
00198             alp = get_step( pd, f, k, err );
00199             j += k;
00200             if( conjGradDebug > 1 )
00201             {
00202                 MSQ_PRINT( 2 )( "\n CG's search direction reset." );
00203                 if( conjGradDebug > 2 ) MSQ_PRINT( 2 )( "\n  Alp was zero, alp = %20.18f", alp );
00204             }
00205         }
00206         if( alp != 0 )
00207         {
00208             pd.move_free_vertices_constrained( arrptr( pGrad ), num_vert, alp, err );MSQ_ERRRTN( err );
00209 
00210             if( !objFunc.update( pd, f, fNewGrad, err ) )
00211             {
00212                 MSQ_SETERR( err )
00213                 ( "Error inside Conjugate Gradient, vertices moved "
00214                   "making function value invalid.",
00215                   MsqError::INVALID_MESH );
00216                 return;
00217             }
00218             assert( fNewGrad.size() == (unsigned)num_vert );
00219 
00220             if( conjGradDebug > 0 )
00221             {
00222                 grad_norm = Linf( arrptr( fNewGrad ), num_vert );
00223                 MSQ_PRINT( 2 )
00224                 ( "\nCG's VALUE = %f,  iter. = %i,  grad_norm = %f,  alp = %f", f, i, grad_norm, alp );
00225                 MSQ_PRINT( 2 )( "\n   TIME %f", c_timer.since_birth() );
00226             }
00227             double s11 = 0;
00228             double s12 = 0;
00229             double s22 = 0;
00230             // free_iter.reset();
00231             // while (free_iter.next()) {
00232             //  m=free_iter.value();
00233             for( m = 0; m < num_vert; ++m )
00234             {
00235                 s11 += fGrad[m] % fGrad[m];
00236                 s12 += fGrad[m] % fNewGrad[m];
00237                 s22 += fNewGrad[m] % fNewGrad[m];
00238             }
00239 
00240             // Steepest Descent (takes 2-3 times as long as P-R)
00241             // double bet=0;
00242 
00243             // Fletcher-Reeves (takes twice as long as P-R)
00244             // double bet = s22/s11;
00245 
00246             // Polack-Ribiere
00247             double bet;
00248             if( !divide( s22 - s12, s11, bet ) ) return;  // gradient is zero
00249             // free_iter.reset();
00250             // while (free_iter.next()) {
00251             //  m=free_iter.value();
00252             for( m = 0; m < num_vert; ++m )
00253             {
00254                 pGrad[m] = ( -fNewGrad[m] + ( bet * pGrad[m] ) );
00255                 fGrad[m] = fNewGrad[m];
00256             }
00257             if( conjGradDebug > 2 )
00258             {
00259                 MSQ_PRINT( 2 )
00260                 ( " \nSEARCH DIRECTION INFINITY NORM = %e", Linf( arrptr( fNewGrad ), num_vert ) );
00261             }
00262 
00263         }  // end if on alp == 0
00264 
00265         term_crit->accumulate_patch( pd, err );MSQ_ERRRTN( err );
00266         term_crit->accumulate_inner( pd, f, arrptr( fGrad ), err );MSQ_ERRRTN( err );
00267     }  // end while
00268     if( conjGradDebug > 0 )
00269     {
00270         MSQ_PRINT( 2 )( "\nConjugate Gradient complete i=%i ", i );
00271         MSQ_PRINT( 2 )( "\n-  FINAL value = %f, alp=%4.2e grad_norm=%4.2e", f, alp, grad_norm );
00272         MSQ_PRINT( 2 )( "\n   FINAL TIME %f", c_timer.since_birth() );
00273     }
00274 }
00275 
00276 void ConjugateGradient::terminate_mesh_iteration( PatchData& /*pd*/, MsqError& /*err*/ )
00277 {
00278     //  cout << "- Executing ConjugateGradient::iteration_complete()\n";
00279 }
00280 
00281 void ConjugateGradient::cleanup()
00282 {
00283     //  cout << "- Executing ConjugateGradient::iteration_end()\n";
00284     fGrad.clear();
00285     pGrad.clear();
00286     fNewGrad.clear();
00287     // pMemento->~PatchDataVerticesMemento();
00288     delete pMemento;
00289     pMemento = NULL;
00290 }
00291 
00292 //! Computes a distance to move vertices given an initial position and search direction (stored in
00293 //! data member pGrad).
00294 /*!Returns alp, the double which scales the search direction vector
00295   which when added to the old nodal positions yields the new nodal
00296   positions.*/
00297 /*!\todo Michael NOTE:  ConjugateGradient::get_step's int &j is only
00298   to remain consisitent with CUBIT for an initial test.  It can be
00299   removed.*/
00300 
00301 double ConjugateGradient::get_step( PatchData& pd, double f0, int& j, MsqError& err )
00302 {
00303     // get OF evaluator
00304     OFEvaluator& objFunc = get_objective_function_evaluator();
00305 
00306     size_t num_vertices = pd.num_free_vertices();
00307     // initial guess for alp
00308     double alp = 1.0;
00309     int jmax   = 100;
00310     double rho = 0.5;
00311     // feasible=false implies the mesh is not in the feasible region
00312     bool feasible = false;
00313     int found     = 0;
00314     // f and fnew hold the objective function value
00315     double f    = 0;
00316     double fnew = 0;
00317     // Counter to avoid infinitly scaling alp
00318     j = 0;
00319     // save memento
00320     pd.recreate_vertices_memento( pMemento, err );
00321     // if we must check feasiblility
00322     // while step takes mesh into infeasible region and ...
00323     while( j < jmax && !feasible && alp > MSQ_MIN )
00324     {
00325         ++j;
00326         pd.set_free_vertices_constrained( pMemento, arrptr( pGrad ), num_vertices, alp, err );
00327         feasible = objFunc.evaluate( pd, f, err );
00328         if( err.error_code() == err.BARRIER_VIOLATED )
00329             err.clear();  // barrier violation does not represent an actual error here
00330         MSQ_ERRZERO( err );
00331         // if not feasible, try a smaller alp (take smaller step)
00332         if( !feasible ) { alp *= rho; }
00333     }  // end while ...
00334 
00335     // if above while ended due to j>=jmax, no valid step was found.
00336     if( j >= jmax )
00337     {
00338         MSQ_PRINT( 2 )( "\nFeasible Point Not Found" );
00339         return 0.0;
00340     }
00341     // Message::print_info("\nOriginal f %f, first new f = %f, alp = %f",f0,f,alp);
00342     // if new f is larger than original, our step was too large
00343     if( f >= f0 )
00344     {
00345         j = 0;
00346         while( j < jmax && found == 0 )
00347         {
00348             ++j;
00349             alp *= rho;
00350             pd.set_free_vertices_constrained( pMemento, arrptr( pGrad ), num_vertices, alp, err );
00351             // Get new obj value
00352             // if patch is now invalid, then the feasible region is  convex or
00353             // we have an error.  For now, we assume an error.
00354             if( !objFunc.evaluate( pd, f, err ) )
00355             {
00356                 MSQ_SETERR( err )
00357                 ( "Non-convex feasiblility region found.", MsqError::INVALID_MESH );
00358             }
00359             pd.set_to_vertices_memento( pMemento, err );
00360             MSQ_ERRZERO( err );
00361             // if our step has now improved the objective function value
00362             if( f < f0 ) { found = 1; }
00363         }  //   end while j less than jmax
00364            // Message::print_info("\nj = %d found = %d f = %20.18f f0 = %20.18f\n",j,found,f,f0);
00365            // if above ended because of j>=jmax, take no step
00366         if( found == 0 )
00367         {
00368             // Message::print_info("alp = %10.8f, but returning zero\n",alp);
00369             alp = 0.0;
00370             return alp;
00371         }
00372 
00373         j = 0;
00374         // while shrinking the step improves the objFunc value further,
00375         // scale alp down.  Return alp, when scaling once more would
00376         // no longer improve the objFunc value.
00377         while( j < jmax )
00378         {
00379             ++j;
00380             alp *= rho;
00381             // step alp in search direction from original positions
00382             pd.set_free_vertices_constrained( pMemento, arrptr( pGrad ), num_vertices, alp, err );
00383             MSQ_ERRZERO( err );
00384 
00385             // get new objective function value
00386             if( !objFunc.evaluate( pd, fnew, err ) ) MSQ_SETERR( err )
00387             ( "Non-convex feasiblility region found while "
00388               "computing new f.",
00389               MsqError::INVALID_MESH );
00390             if( fnew < f ) { f = fnew; }
00391             else
00392             {
00393                 // Reset the vertices to original position
00394                 pd.set_to_vertices_memento( pMemento, err );
00395                 MSQ_ERRZERO( err );
00396                 alp /= rho;
00397                 return alp;
00398             }
00399         }
00400         // Reset the vertices to original position and return alp
00401         pd.set_to_vertices_memento( pMemento, err );
00402         MSQ_ERRZERO( err );
00403         return alp;
00404     }
00405     // else our new f was already smaller than our original
00406     else
00407     {
00408         j = 0;
00409         // check to see how large of step we can take
00410         while( j < jmax && found == 0 )
00411         {
00412             ++j;
00413             // scale alp up (rho must be less than 1)
00414             alp /= rho;
00415             // step alp in search direction from original positions
00416             pd.set_free_vertices_constrained( pMemento, arrptr( pGrad ), num_vertices, alp, err );
00417             MSQ_ERRZERO( err );
00418 
00419             feasible = objFunc.evaluate( pd, fnew, err );
00420             if( err.error_code() == err.BARRIER_VIOLATED )
00421                 err.clear();  // evaluate() error does not represent an actual problem here
00422             MSQ_ERRZERO( err );
00423             if( !feasible )
00424             {
00425                 alp *= rho;
00426                 // Reset the vertices to original position and return alp
00427                 pd.set_to_vertices_memento( pMemento, err );
00428                 MSQ_ERRZERO( err );
00429                 return alp;
00430             }
00431             if( fnew < f ) { f = fnew; }
00432             else
00433             {
00434                 found = 1;
00435                 alp *= rho;
00436             }
00437         }
00438 
00439         // Reset the vertices to original position and return alp
00440         pd.set_to_vertices_memento( pMemento, err );
00441         MSQ_ERRZERO( err );
00442         return alp;
00443     }
00444 }
00445 
00446 /*!Quadratic one-dimensional line search.*/
00447 /*
00448 double ConjugateGradient::get_step(PatchData &pd,double f0,int &j,
00449                                    MsqError &err){
00450   const double CGOLD = 0.3819660;
00451   const double ZEPS = 1.0e-10;
00452   int n=pd.num_free_vertices();
00453   MsqVertex* vertices=pd.get_vertex_array(err);
00454   double a,b,d,etemp,fb,fu,fv,fw,fx,p,q,r,tol,tol1,tol2,u,v,w,x,xm;
00455   double e=0.0;
00456   d=0.0;
00457   tol=.001;
00458   int iter, maxiter;
00459   maxiter=100;
00460   a=0;
00461   b=.125;
00462   int m=0;
00463   fb=f0-1.0;
00464   iter=0;
00465   //find b such that a b 'should' bracket the min
00466   while (fb<=f0 && iter<maxiter){
00467     ++iter;
00468     b*=2.0;
00469     for(m=0;m<n;++m){
00470       mCoord[m]=mCoord[m] + (b*pGrad[m]);
00471       vertices[m]=(mCoord[m]);
00472     }
00473     fb=objFunc->evaluate(pd,err);
00474   }
00475   iter=0;
00476   x=w=v=(b/2.0);
00477   for(m=0;m<n;++m){
00478     mCoord[m]=mCoord[m] + (x*pGrad[m]);
00479     vertices[m]=(mCoord[m]);
00480   }
00481   fw=fv=fx=objFunc->evaluate(pd,err);
00482   for(iter=0;iter<maxiter;++iter){
00483       //Message::print_info("a=%f,b=%f,x=%f,iter=%i\n",a,b,x,iter);
00484     xm=(a+b)*.5;
00485     tol2=2.0*(tol1=tol*fabs(x)+ZEPS);
00486     if(fabs(x-xm)<= (tol2-0.5*(b-a))){
00487       return x;
00488     }
00489     if(fabs(e)>tol1){
00490       r=(x-w)*(fx-fv);
00491       q=(x-v)*(fx-fw);
00492       p=(x-v)*q-(x-w)*r;
00493       q=2.0*(q-r);
00494       if(q>0.0)
00495         p=-p;
00496       q=fabs(q);
00497       etemp=e;
00498       e=d;
00499       if(fabs(p)>=fabs(0.5*q*etemp)||(p<=q*(a-x))||(p>=q*(b-x))){
00500         d=CGOLD*(e=(x>=xm?a-x:b-x));
00501       }
00502       else{
00503         d=p/q;
00504         u=x+d;
00505         if(u-a<tol2||b-u<tol2)
00506         {
00507           if(tol1<0.0)
00508             d=x-xm;
00509           else
00510             d=xm-x;
00511         }
00512       }
00513     }
00514 
00515     else{
00516       d=CGOLD*(e=(x>=xm?a-x:b-x));
00517     }
00518     if(tol<0.0)
00519       u=(fabs(d)>=tol1?x+d:x-d);
00520     else
00521       u=(fabs(d)>=tol1?x+d:x+d);
00522     for(m=0;m<n;++m){
00523       mCoord[m]=mCoord[m] + (u*pGrad[m]);
00524       vertices[m]=(mCoord[m]);
00525     }
00526     fu=objFunc->evaluate(pd,err);
00527     if(fu<fx){
00528       if(u>=x)
00529         a=x;
00530       else
00531         b=x;
00532       v=w;
00533       w=x;
00534       x=u;
00535       fv=fw;
00536       fw=fx;
00537       fx=fu;
00538     }
00539     else{
00540       if(u<x)
00541         a=u;
00542       else
00543         b=u;
00544       if(fu<=fw||w==x){
00545         v=w;
00546         w=u;
00547         fv=fw;
00548         fw=fu;
00549       }
00550       else if (fu<=fv||v==x||v==w){
00551         v=u;
00552         fv=fu;
00553       }
00554     }
00555   }
00556   for(m=0;m<n;++m){
00557     vertices[m]=(mCoord[m]);
00558   }
00559     //PRINT_WARNING("TOO MANY ITERATIONS IN QUADRATIC LINE SEARCH");
00560   return x;
00561 }
00562 */
00563 
00564 }  // namespace MBMesquite
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines