Optimizes the objective function using the Polack-Ribiere scheme. More...

#include <ConjugateGradient.hpp>

Inheritance diagram for MBMesquite::ConjugateGradient:

Collaboration diagram for MBMesquite::ConjugateGradient:

Public Member Functions
MESQUITE_EXPORT	ConjugateGradient (ObjectiveFunction *objective)
MESQUITE_EXPORT	ConjugateGradient (ObjectiveFunction *objective, MsqError &err)
virtual MESQUITE_EXPORT	~ConjugateGradient ()
virtual MESQUITE_EXPORT std::string	get_name () const
	Get string name for use in diagnostic and status output.
virtual PatchSet *	get_patch_set ()
MESQUITE_EXPORT void	set_debugging_level (int new_lev)
Protected Member Functions
virtual void	initialize (PatchData &pd, MsqError &err)
	Initialize data for smoothing process.
virtual void	optimize_vertex_positions (PatchData &pd, MsqError &err)
virtual void	initialize_mesh_iteration (PatchData &pd, MsqError &err)
virtual void	terminate_mesh_iteration (PatchData &pd, MsqError &err)
virtual void	cleanup ()
	Delete arrays initially created in initialize().
double	get_step (PatchData &pd, double f0, int &j, MsqError &err)
	Returns the step distance to take in the search direction.
Private Attributes
std::vector< Vector3D >	fGrad
	Culls the vertex list free_vertex_list.
std::vector< Vector3D >	pGrad
std::vector< Vector3D >	fNewGrad
PatchDataVerticesMemento *	pMemento
int	conjGradDebug

Detailed Description

Optimizes the objective function using the Polack-Ribiere scheme.

Definition at line 51 of file ConjugateGradient.hpp.

Constructor & Destructor Documentation

MBMesquite::ConjugateGradient::ConjugateGradient ( ObjectiveFunction * objective )

Definition at line 60 of file ConjugateGradient.cpp.

    : VertexMover( objective ), PatchSetUser( true ), pMemento( NULL ), conjGradDebug( 0 )
{
}

MBMesquite::ConjugateGradient::ConjugateGradient	(	ObjectiveFunction *	objective,
		MsqError &	err
	)

Definition at line 65 of file ConjugateGradient.cpp.

References MBMesquite::TerminationCriterion::add_iteration_limit(), MBMesquite::QualityImprover::get_inner_termination_criterion(), MBMesquite::MsqError::INVALID_STATE, MSQ_ERRRTN, MSQ_SETERR, and set_debugging_level().

    : VertexMover( objective ), PatchSetUser( true ), pMemento( NULL ), conjGradDebug( 0 )
{
    // Michael:: default to global?
    set_debugging_level( 0 );
    // set the default inner termination criterion
    TerminationCriterion* default_crit = get_inner_termination_criterion();
    if( default_crit == NULL )
    {
        MSQ_SETERR( err )
        ( "QualityImprover did not create a default inner "
          "termination criterion.",
          MsqError::INVALID_STATE );
        return;
    }
    else
    {
        default_crit->add_iteration_limit( 5 );MSQ_ERRRTN( err );
    }
}

MBMesquite::ConjugateGradient::~ConjugateGradient ( ) [virtual]

Definition at line 86 of file ConjugateGradient.cpp.

References pMemento.

{
    // Checks that cleanup() has been called.
    assert( pMemento == NULL );
}

Member Function Documentation

void MBMesquite::ConjugateGradient::cleanup ( ) [protected, virtual]

Delete arrays initially created in initialize().

Implements MBMesquite::VertexMover.

Definition at line 284 of file ConjugateGradient.cpp.

References fGrad, fNewGrad, pGrad, and pMemento.

{
    //  cout << "- Executing ConjugateGradient::iteration_end()\n";
    fGrad.clear();
    pGrad.clear();
    fNewGrad.clear();
    // pMemento->~PatchDataVerticesMemento();
    delete pMemento;
    pMemento = NULL;
}

std::string MBMesquite::ConjugateGradient::get_name ( ) const [virtual]

Get string name for use in diagnostic and status output.

Implements MBMesquite::Instruction.

Definition at line 50 of file ConjugateGradient.cpp.

{
    return "ConjugateGradient";
}

PatchSet * MBMesquite::ConjugateGradient::get_patch_set ( ) [virtual]

Reimplemented from MBMesquite::PatchSetUser.

Definition at line 55 of file ConjugateGradient.cpp.

{
    return PatchSetUser::get_patch_set();
}

double MBMesquite::ConjugateGradient::get_step	(	PatchData &	pd,
		double	f0,
		int &	j,
		MsqError &	err
	)		`[protected]`

Returns the step distance to take in the search direction.

Computes a distance to move vertices given an initial position and search direction (stored in data member pGrad).

Returns alp, the double which scales the search direction vector which when added to the old nodal positions yields the new nodal positions.

Todo:: Michael NOTE: ConjugateGradient::get_step's int &j is only to remain consisitent with CUBIT for an initial test. It can be removed.

Definition at line 304 of file ConjugateGradient.cpp.

References MBMesquite::arrptr(), MBMesquite::MsqError::BARRIER_VIOLATED, MBMesquite::MsqError::clear(), MBMesquite::MsqError::error_code(), MBMesquite::OFEvaluator::evaluate(), MBMesquite::VertexMover::get_objective_function_evaluator(), MBMesquite::MsqError::INVALID_MESH, MSQ_ERRZERO, MBMesquite::MSQ_MIN, MSQ_PRINT, MSQ_SETERR, MBMesquite::PatchData::num_free_vertices(), pGrad, pMemento, MBMesquite::PatchData::recreate_vertices_memento(), MBMesquite::PatchData::set_free_vertices_constrained(), and MBMesquite::PatchData::set_to_vertices_memento().

Referenced by optimize_vertex_positions().

{
    // get OF evaluator
    OFEvaluator& objFunc = get_objective_function_evaluator();

    size_t num_vertices = pd.num_free_vertices();
    // initial guess for alp
    double alp = 1.0;
    int jmax   = 100;
    double rho = 0.5;
    // feasible=false implies the mesh is not in the feasible region
    bool feasible = false;
    int found     = 0;
    // f and fnew hold the objective function value
    double f    = 0;
    double fnew = 0;
    // Counter to avoid infinitly scaling alp
    j = 0;
    // save memento
    pd.recreate_vertices_memento( pMemento, err );
    // if we must check feasiblility
    // while step takes mesh into infeasible region and ...
    while( j < jmax && !feasible && alp > MSQ_MIN )
    {
        ++j;
        pd.set_free_vertices_constrained( pMemento, arrptr( pGrad ), num_vertices, alp, err );
        feasible = objFunc.evaluate( pd, f, err );
        if( err.error_code() == err.BARRIER_VIOLATED )
            err.clear();  // barrier violation does not represent an actual error here
        MSQ_ERRZERO( err );
        // if not feasible, try a smaller alp (take smaller step)
        if( !feasible )
        {
            alp *= rho;
        }
    }  // end while ...

    // if above while ended due to j>=jmax, no valid step was found.
    if( j >= jmax )
    {
        MSQ_PRINT( 2 )( "\nFeasible Point Not Found" );
        return 0.0;
    }
    // Message::print_info("\nOriginal f %f, first new f = %f, alp = %f",f0,f,alp);
    // if new f is larger than original, our step was too large
    if( f >= f0 )
    {
        j = 0;
        while( j < jmax && found == 0 )
        {
            ++j;
            alp *= rho;
            pd.set_free_vertices_constrained( pMemento, arrptr( pGrad ), num_vertices, alp, err );
            // Get new obj value
            // if patch is now invalid, then the feasible region is  convex or
            // we have an error.  For now, we assume an error.
            if( !objFunc.evaluate( pd, f, err ) )
            {
                MSQ_SETERR( err )
                ( "Non-convex feasiblility region found.", MsqError::INVALID_MESH );
            }
            pd.set_to_vertices_memento( pMemento, err );
            MSQ_ERRZERO( err );
            // if our step has now improved the objective function value
            if( f < f0 )
            {
                found = 1;
            }
        }  //   end while j less than jmax
           // Message::print_info("\nj = %d found = %d f = %20.18f f0 = %20.18f\n",j,found,f,f0);
           // if above ended because of j>=jmax, take no step
        if( found == 0 )
        {
            // Message::print_info("alp = %10.8f, but returning zero\n",alp);
            alp = 0.0;
            return alp;
        }

        j = 0;
        // while shrinking the step improves the objFunc value further,
        // scale alp down.  Return alp, when scaling once more would
        // no longer improve the objFunc value.
        while( j < jmax )
        {
            ++j;
            alp *= rho;
            // step alp in search direction from original positions
            pd.set_free_vertices_constrained( pMemento, arrptr( pGrad ), num_vertices, alp, err );
            MSQ_ERRZERO( err );

            // get new objective function value
            if( !objFunc.evaluate( pd, fnew, err ) ) MSQ_SETERR( err )
            ( "Non-convex feasiblility region found while "
              "computing new f.",
              MsqError::INVALID_MESH );
            if( fnew < f )
            {
                f = fnew;
            }
            else
            {
                // Reset the vertices to original position
                pd.set_to_vertices_memento( pMemento, err );
                MSQ_ERRZERO( err );
                alp /= rho;
                return alp;
            }
        }
        // Reset the vertices to original position and return alp
        pd.set_to_vertices_memento( pMemento, err );
        MSQ_ERRZERO( err );
        return alp;
    }
    // else our new f was already smaller than our original
    else
    {
        j = 0;
        // check to see how large of step we can take
        while( j < jmax && found == 0 )
        {
            ++j;
            // scale alp up (rho must be less than 1)
            alp /= rho;
            // step alp in search direction from original positions
            pd.set_free_vertices_constrained( pMemento, arrptr( pGrad ), num_vertices, alp, err );
            MSQ_ERRZERO( err );

            feasible = objFunc.evaluate( pd, fnew, err );
            if( err.error_code() == err.BARRIER_VIOLATED )
                err.clear();  // evaluate() error does not represent an actual problem here
            MSQ_ERRZERO( err );
            if( !feasible )
            {
                alp *= rho;
                // Reset the vertices to original position and return alp
                pd.set_to_vertices_memento( pMemento, err );
                MSQ_ERRZERO( err );
                return alp;
            }
            if( fnew < f )
            {
                f = fnew;
            }
            else
            {
                found = 1;
                alp *= rho;
            }
        }

        // Reset the vertices to original position and return alp
        pd.set_to_vertices_memento( pMemento, err );
        MSQ_ERRZERO( err );
        return alp;
    }
}

void MBMesquite::ConjugateGradient::initialize	(	PatchData &	pd,
		MsqError &	err
	)		`[protected, virtual]`

Initialize data for smoothing process.

Implements MBMesquite::VertexMover.

Definition at line 92 of file ConjugateGradient.cpp.

References MBMesquite::PatchData::create_vertices_memento(), MBMesquite::get_parallel_rank(), MBMesquite::get_parallel_size(), MSQ_DBGOUT, and pMemento.

{
    if( get_parallel_size() )
    {
        MSQ_DBGOUT( 2 ) << "\nP[" << get_parallel_rank() << "] "
                        << "o   Performing Conjugate Gradient optimization.\n";
    }
    else
    {
        MSQ_DBGOUT( 2 ) << "\no   Performing Conjugate Gradient optimization.\n";
    }
    pMemento = pd.create_vertices_memento( err );
}

void MBMesquite::ConjugateGradient::initialize_mesh_iteration	(	PatchData &	pd,
		MsqError &	err
	)		`[protected, virtual]`

Implements MBMesquite::VertexMover.

Definition at line 106 of file ConjugateGradient.cpp.

{}

void MBMesquite::ConjugateGradient::optimize_vertex_positions	(	PatchData &	pd,
		MsqError &	err
	)		`[protected, virtual]`

Performs Conjugate gradient minimization on the PatchData, pd.

Implements MBMesquite::VertexMover.

Definition at line 109 of file ConjugateGradient.cpp.

References MBMesquite::TerminationCriterion::accumulate_inner(), MBMesquite::TerminationCriterion::accumulate_patch(), MBMesquite::arrptr(), conjGradDebug, MBMesquite::divide(), fGrad, fNewGrad, MBMesquite::QualityImprover::get_inner_termination_criterion(), MBMesquite::VertexMover::get_objective_function_evaluator(), get_step(), MBMesquite::MsqError::INVALID_MESH, MBMesquite::Linf(), MBMesquite::PatchData::move_free_vertices_constrained(), MSQ_CHKERR, MSQ_DBGOUT, MSQ_ERRRTN, MSQ_FUNCTION_TIMER, MBMesquite::MSQ_MAX_CAP, MSQ_PRINT, MSQ_SETERR, MBMesquite::PatchData::num_free_vertices(), pGrad, MBMesquite::Timer::since_birth(), MBMesquite::TerminationCriterion::terminate(), and MBMesquite::OFEvaluator::update().

{
    // pd.reorder();

    MSQ_FUNCTION_TIMER( "ConjugateGradient::optimize_vertex_positions" );

    Timer c_timer;
    size_t num_vert = pd.num_free_vertices();
    if( num_vert < 1 )
    {
        MSQ_DBGOUT( 1 ) << "\nEmpty free vertex list in ConjugateGradient\n";
        return;
    }
    /*
        //zero out arrays
      int zero_loop=0;
      while(zero_loop<arraySize){
        fGrad[zero_loop].set(0,0,0);
        pGrad[zero_loop].set(0,0,0);
        fNewGrad[zero_loop].set(0,0,0);
        ++zero_loop;
      }
    */

    // get OF evaluator
    OFEvaluator& objFunc = get_objective_function_evaluator();

    size_t ind;
    // Michael cull list:  possibly set soft_fixed flags here

    // MsqFreeVertexIndexIterator free_iter(pd, err);  MSQ_ERRRTN(err);

    double f = 0;
    // Michael, this isn't equivalent to CUBIT because we only want to check
    // the objective function value of the 'bad' elements
    // if invalid initial patch set an error.
    bool temp_bool = objFunc.update( pd, f, fGrad, err );
    assert( fGrad.size() == num_vert );
    if( MSQ_CHKERR( err ) ) return;
    if( !temp_bool )
    {
        MSQ_SETERR( err )
        ( "Conjugate Gradient not able to get valid gradient "
          "and function values on intial patch.",
          MsqError::INVALID_MESH );
        return;
    }
    double grad_norm = MSQ_MAX_CAP;

    if( conjGradDebug > 0 )
    {
        MSQ_PRINT( 2 )( "\nCG's DEGUB LEVEL = %i \n", conjGradDebug );
        grad_norm = Linf( arrptr( fGrad ), fGrad.size() );
        MSQ_PRINT( 2 )( "\nCG's FIRST VALUE = %f,grad_norm = %f", f, grad_norm );
        MSQ_PRINT( 2 )( "\n   TIME %f", c_timer.since_birth() );
        grad_norm = MSQ_MAX_CAP;
    }

    // Initializing pGrad (search direction).
    pGrad.resize( fGrad.size() );
    for( ind = 0; ind < num_vert; ++ind )
        pGrad[ind] = ( -fGrad[ind] );

    int j      = 0;            // total nb of step size changes ... not used much
    int i      = 0;            // iteration counter
    unsigned m = 0;            //
    double alp = MSQ_MAX_CAP;  // alp: scale factor of search direction
                               // we know inner_criterion is false because it was checked in
                               // loop_over_mesh before being sent here.
    TerminationCriterion* term_crit = get_inner_termination_criterion();

    // while ((i<maxIteration && alp>stepBound && grad_norm>normGradientBound)
    //     && !inner_criterion){
    while( !term_crit->terminate() )
    {
        ++i;
        // std::cout<<"\Michael delete i = "<<i;
        int k = 0;
        alp   = get_step( pd, f, k, err );
        j += k;
        if( conjGradDebug > 2 )
        {
            MSQ_PRINT( 2 )( "\n  Alp initial, alp = %20.18f", alp );
        }

        // if alp == 0, revert to steepest descent search direction
        if( alp == 0 )
        {
            for( m = 0; m < num_vert; ++m )
            {
                pGrad[m] = ( -fGrad[m] );
            }
            alp = get_step( pd, f, k, err );
            j += k;
            if( conjGradDebug > 1 )
            {
                MSQ_PRINT( 2 )( "\n CG's search direction reset." );
                if( conjGradDebug > 2 ) MSQ_PRINT( 2 )( "\n  Alp was zero, alp = %20.18f", alp );
            }
        }
        if( alp != 0 )
        {
            pd.move_free_vertices_constrained( arrptr( pGrad ), num_vert, alp, err );MSQ_ERRRTN( err );

            if( !objFunc.update( pd, f, fNewGrad, err ) )
            {
                MSQ_SETERR( err )
                ( "Error inside Conjugate Gradient, vertices moved "
                  "making function value invalid.",
                  MsqError::INVALID_MESH );
                return;
            }
            assert( fNewGrad.size() == (unsigned)num_vert );

            if( conjGradDebug > 0 )
            {
                grad_norm = Linf( arrptr( fNewGrad ), num_vert );
                MSQ_PRINT( 2 )
                ( "\nCG's VALUE = %f,  iter. = %i,  grad_norm = %f,  alp = %f", f, i, grad_norm, alp );
                MSQ_PRINT( 2 )( "\n   TIME %f", c_timer.since_birth() );
            }
            double s11 = 0;
            double s12 = 0;
            double s22 = 0;
            // free_iter.reset();
            // while (free_iter.next()) {
            //  m=free_iter.value();
            for( m = 0; m < num_vert; ++m )
            {
                s11 += fGrad[m] % fGrad[m];
                s12 += fGrad[m] % fNewGrad[m];
                s22 += fNewGrad[m] % fNewGrad[m];
            }

            // Steepest Descent (takes 2-3 times as long as P-R)
            // double bet=0;

            // Fletcher-Reeves (takes twice as long as P-R)
            // double bet = s22/s11;

            // Polack-Ribiere
            double bet;
            if( !divide( s22 - s12, s11, bet ) ) return;  // gradient is zero
            // free_iter.reset();
            // while (free_iter.next()) {
            //  m=free_iter.value();
            for( m = 0; m < num_vert; ++m )
            {
                pGrad[m] = ( -fNewGrad[m] + ( bet * pGrad[m] ) );
                fGrad[m] = fNewGrad[m];
            }
            if( conjGradDebug > 2 )
            {
                MSQ_PRINT( 2 )
                ( " \nSEARCH DIRECTION INFINITY NORM = %e", Linf( arrptr( fNewGrad ), num_vert ) );
            }

        }  // end if on alp == 0

        term_crit->accumulate_patch( pd, err );MSQ_ERRRTN( err );
        term_crit->accumulate_inner( pd, f, arrptr( fGrad ), err );MSQ_ERRRTN( err );
    }  // end while
    if( conjGradDebug > 0 )
    {
        MSQ_PRINT( 2 )( "\nConjugate Gradient complete i=%i ", i );
        MSQ_PRINT( 2 )( "\n-  FINAL value = %f, alp=%4.2e grad_norm=%4.2e", f, alp, grad_norm );
        MSQ_PRINT( 2 )( "\n   FINAL TIME %f", c_timer.since_birth() );
    }
}

MESQUITE_EXPORT void MBMesquite::ConjugateGradient::set_debugging_level ( int new_lev ) [inline]

Just for debugging purposes or for obtaining more data during the optimization process.

Definition at line 65 of file ConjugateGradient.hpp.

References conjGradDebug.

Referenced by ConjugateGradient(), run_solution_mesh_optimizer(), SphericalGeometryTest::test_cg_mesh_cond_sphere(), and PlanarGeometryTest::test_plane_tri_xz().

    {
        conjGradDebug = new_lev;
    }

void MBMesquite::ConjugateGradient::terminate_mesh_iteration	(	PatchData &	pd,
		MsqError &	err
	)		`[protected, virtual]`

Implements MBMesquite::VertexMover.

Definition at line 279 of file ConjugateGradient.cpp.

{
    //  cout << "- Executing ConjugateGradient::iteration_complete()\n";
}

Member Data Documentation

int MBMesquite::ConjugateGradient::conjGradDebug [private]

Definition at line 93 of file ConjugateGradient.hpp.

Referenced by optimize_vertex_positions(), and set_debugging_level().

std::vector< Vector3D > MBMesquite::ConjugateGradient::fGrad [private]

Culls the vertex list free_vertex_list.

Definition at line 90 of file ConjugateGradient.hpp.

Referenced by cleanup(), and optimize_vertex_positions().

std::vector< Vector3D > MBMesquite::ConjugateGradient::fNewGrad [private]

Definition at line 90 of file ConjugateGradient.hpp.

Referenced by cleanup(), and optimize_vertex_positions().

std::vector< Vector3D > MBMesquite::ConjugateGradient::pGrad [private]

Definition at line 90 of file ConjugateGradient.hpp.

Referenced by cleanup(), get_step(), and optimize_vertex_positions().

PatchDataVerticesMemento* MBMesquite::ConjugateGradient::pMemento [private]

Definition at line 91 of file ConjugateGradient.hpp.

Referenced by cleanup(), get_step(), initialize(), and ~ConjugateGradient().

List of all members.

The documentation for this class was generated from the following files:

Public Member Functions

Protected Member Functions

Private Attributes

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation