CSolverStepLengthSelection.txx

#ifndef C_SOLVER_STEP_LENGTH_SELECTION_TXX
#define C_SOLVER_STEP_LENGTH_SELECTION_TXX


template < class TState >
CSolverStepLengthSelection< TState>::CSolverStepLengthSelection()
  : DefaultMinGradAllowed( 0.0001 ),
    DefaultMinDisplacementAllowed( 0.001 ),
    DefaultDecreaseConstant( 0.0001 ),
    DefaultMaxNumberOfIterations( 100 ),
    DefaultMaxNumberOfTries( 10 ),
    m_ExternallySetMinGradAllowed( false ),
    m_ExternallySetMinDisplacementAllowed( false ),
    m_ExternallySetDecreaseConstant( false ),
    m_ExternallySetMaxNumberOfIterations( false ),
    m_ExternallySetMaxNumberOfTries( false )
{
  // default setting for the parameters
  // minimum value allowed for the norm of the gradient
  m_MinGradAllowed = DefaultMinGradAllowed;

  // minimum value allowed for the displacement
  m_MinDisplacementAllowed = DefaultMinDisplacementAllowed;

  // constant for the sufficient decrease condition
  m_DecreaseConstant = DefaultDecreaseConstant;

  // maximal number of iterations
  m_MaxNumberOfIterations = DefaultMaxNumberOfIterations;

  // maximal number of tries in one line search
  m_MaxNumberOfTries = DefaultMaxNumberOfTries;

}


template < class TState >
CSolverStepLengthSelection< TState>::~CSolverStepLengthSelection()
{
}


template < class TState >
void CSolverStepLengthSelection< TState>::SetAutoConfiguration( CJSONConfiguration * combined, CJSONConfiguration * cleaned )
{
  Superclass::SetAutoConfiguration( combined, cleaned );

  Json::Value& currentConfigurationIn = this->m_CombinedJSONConfig->GetFromKey( "LineSearch", Json::nullValue );
  Json::Value& currentConfigurationOut = this->m_CleanedJSONConfig->GetFromKey( "LineSearch", Json::nullValue, CONF_NORMAL );

  SetJSONHelpForRootKey( LineSearch, "setting for the linesearch algorithm", CONF_NORMAL );

  SetJSONFromKeyDouble( currentConfigurationIn, currentConfigurationOut, MinGradAllowed, CONF_ADVANCED );
  SetJSONFromKeyDouble( currentConfigurationIn, currentConfigurationOut, MinDisplacementAllowed, CONF_ADVANCED );
  SetJSONFromKeyDouble( currentConfigurationIn, currentConfigurationOut, DecreaseConstant, CONF_EXPERT );
  SetJSONFromKeyUInt( currentConfigurationIn, currentConfigurationOut, MaxNumberOfIterations, CONF_NORMAL );
  SetJSONFromKeyUInt( currentConfigurationIn, currentConfigurationOut, MaxNumberOfTries, CONF_ADVANCED );
  SetJSONFromKeyBool( currentConfigurationIn, currentConfigurationOut, OutputStateInformation, CONF_EXPERT );
  SetJSONFromKeyUInt( currentConfigurationIn, currentConfigurationOut, OutputStateInformationFrequency, CONF_EXPERT );

  SetJSONHelpForKey( currentConfigurationIn, currentConfigurationOut, MinGradAllowed,
                     "minimum value allowed for the norm of the gradient", CONF_ADVANCED );
  SetJSONHelpForKey( currentConfigurationIn, currentConfigurationOut, MinDisplacementAllowed,
                     "minimum value allowed for the displacement", CONF_ADVANCED );
  SetJSONHelpForKey( currentConfigurationIn, currentConfigurationOut, DecreaseConstant,
                     "require sufficient decrease of energy", CONF_EXPERT );
  SetJSONHelpForKey( currentConfigurationIn, currentConfigurationOut, MaxNumberOfIterations,
                     "maximal number of iterations", CONF_NORMAL );
  SetJSONHelpForKey( currentConfigurationIn, currentConfigurationOut, MaxNumberOfTries,
                     "number of tries to find a smaller solution", CONF_ADVANCED );
  SetJSONHelpForKey( currentConfigurationIn, currentConfigurationOut, OutputStateInformation,
                     "output internal information of algorithm (for debugging)", CONF_EXPERT );
  SetJSONHelpForKey( currentConfigurationIn, currentConfigurationOut, OutputStateInformationFrequency,
                     "controls at which iterations the state information should be output", CONF_EXPERT );
}


template < class TState >
bool CSolverStepLengthSelection< TState>::StepSelectionLineSearch()
{

    /* Current values for x, f(x) and g(x) */
    ObjectiveFunctionType * f = this->GetObjectiveFunction();
    TState *x_cur = f->GetStatePointer();
    CEnergyValues f_cur = f->GetCurrentEnergy();
    f->ComputeGradient();
    TState *g_cur = f->GetGradientPointer();
    T g_norm = g_cur->SquaredNorm();

    /* Iterations */
    unsigned int it_count = 0;
    while( (g_norm > m_MinGradAllowed) && (it_count <= m_MaxNumberOfIterations) ) //add condition on dx
    {
        it_count++;

        /* compute alpha0 */
        T alpha_cur = 100/(1+g_norm);
        if ( alpha_cur > 1.0 ) alpha_cur = 1.0;

        /* get the new position x and its energy f(x) */
        ObjectiveFunctionType f_test = *f;                  //need alloc?
        TState *x_new = f_test.GetStatePointer();
        *x_new = *g_cur;
        *x_new *= -alpha_cur;
        *x_new += *x_cur;
        CEnergyValues f_new = f_test.GetCurrentEnergy();

        /* variables to store previous values */
        CEnergyValues f_prev = f_new;
        T alpha_prev = alpha_cur;
        T alpha_new = alpha_cur;

        unsigned int it_armijo = 0;
        while( (f_new > f_cur + m_DecreaseConstant * alpha_cur * g_norm) && (it_armijo <= m_MaxNumberOfTries) )
        {
            it_armijo++;

            /* compute new alpha with step length selection */
            if(it_armijo == 1) //quadratic model
            {
                alpha_new = stepLengthSelection(f_cur, g_norm, alpha_cur, f_new); // TO DO
            }
            else // cubic model
            {
                alpha_new = stepLengthSelection(f_cur, g_norm, alpha_cur, f_new, alpha_prev, f_prev); // TO DO
            }

            /* store previous values */
            f_prev = f_new;
            alpha_prev  =   alpha_cur;
            alpha_cur   =   alpha_new;

            /* calculate new value to test for Armijo */
            *x_new = *g_cur;
            *x_new *= -alpha_cur;
            *x_new += *x_cur;
            f_new = f_test.GetCurrentEnergy();

        }

        /* Update new values x, f(x) and g(x) */
        *x_cur = *x_new;        //x_cur updated
        f_cur = f_new;          //f_cur updated
        f->ComputeGradient();   //g_cur updated

    }

}

#endif