/*======================================================================

  This file is part of the elastix software.

  Copyright (c) University Medical Center Utrecht. All rights reserved.
  See src/CopyrightElastix.txt or http://elastix.isi.uu.nl/legal.php for
  details.

     This software is distributed WITHOUT ANY WARRANTY; without even
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
     PURPOSE. See the above copyright notices for more information.

======================================================================*/

#ifndef __itkAdaptiveStochasticGradientDescentOptimizer_cxx
#define __itkAdaptiveStochasticGradientDescentOptimizer_cxx

#include "itkAdaptiveStochasticGradientDescentOptimizer.h"
#include "vnl/vnl_math.h"
#include "itkSigmoidImageFilter.h"

namespace itk
{

  /**
   * ************************* Constructor ************************
   */

  AdaptiveStochasticGradientDescentOptimizer
    ::AdaptiveStochasticGradientDescentOptimizer()
  {
    this->m_UseAdaptiveStepSizes = true;
    this->m_SigmoidMax = 1.0;
    this->m_SigmoidMin = -0.8;
    this->m_SigmoidScale = 1e-8;

  } // end Constructor


  /**
   * ************************** UpdateCurrentTime ********************
   *
   * This function computes the input for the Compute_a function.
   */

  void AdaptiveStochasticGradientDescentOptimizer
    ::UpdateCurrentTime( void )
  {
    typedef itk::Function::Sigmoid<double, double> SigmoidType;

    if ( this->m_UseAdaptiveStepSizes )
    {
      if ( this->GetCurrentIteration() > 0 )
      {
        /** Make sigmoid function
         * Compute beta such that sigmoid(0)=0
         * We assume Max>0, min<0 */
        SigmoidType sigmoid;
        sigmoid.SetOutputMaximum( this->GetSigmoidMax() );
        sigmoid.SetOutputMinimum( this->GetSigmoidMin() );
        sigmoid.SetAlpha( this->GetSigmoidScale() );
        const double beta = this->GetSigmoidScale() *
           vcl_log( - this->GetSigmoidMax() / this->GetSigmoidMin() );
        sigmoid.SetBeta( beta );

        /** Formula (2) in Cruz */
        const double inprod = inner_product(
          this->m_PreviousGradient, this->GetGradient() );
        this->m_CurrentTime += sigmoid(-inprod);
        this->m_CurrentTime = vnl_math_max( 0.0, this->m_CurrentTime );
      }

      /** Save for next iteration */
      this->m_PreviousGradient = this->GetGradient();
    }
    else
    {
      /** Almost Robbins-Monro: time = time + E_0.
       * If you want the parameter estimation but no adaptive stuff,
       * this may be use useful:  */
      this->m_CurrentTime += ( this->GetSigmoidMax() + this->GetSigmoidMin() ) / 2.0;
    }

  } // end UpdateCurrentTime


} // end namespace itk

#endif // end #ifndef __itkAdaptiveStochasticGradientDescentOptimizer_cxx

