// matrix/optimization.h
  
  // Copyright 2012  Johns Hopkins University (author: Daniel Povey)
  //
  // See ../../COPYING for clarification regarding multiple authors
  //
  // Licensed under the Apache License, Version 2.0 (the "License");
  // you may not use this file except in compliance with the License.
  // You may obtain a copy of the License at
  //
  //  http://www.apache.org/licenses/LICENSE-2.0
  //
  // THIS CODE IS PROVIDED *AS IS* BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
  // KIND, EITHER EXPRESS OR IMPLIED, INCLUDING WITHOUT LIMITATION ANY IMPLIED
  // WARRANTIES OR CONDITIONS OF TITLE, FITNESS FOR A PARTICULAR PURPOSE,
  // MERCHANTABLITY OR NON-INFRINGEMENT.
  // See the Apache 2 License for the specific language governing permissions and
  // limitations under the License.
  //
  // (*) incorporates, with permission, FFT code from his book
  // "Signal Processing with Lapped Transforms", Artech, 1992.
  
  
  
  #ifndef KALDI_MATRIX_OPTIMIZATION_H_
  #define KALDI_MATRIX_OPTIMIZATION_H_
  
  #include "matrix/kaldi-vector.h"
  #include "matrix/kaldi-matrix.h"
  
  namespace kaldi {
  
  
  /// @addtogroup matrix_optimization
  /// @{
  
  struct LinearCgdOptions {
    int32 max_iters;  //  Maximum number of iters (if >= 0).
    BaseFloat max_error;  // Maximum 2-norm of the residual A x - b (convergence
                          // test)
    // Every time the residual 2-norm decreases by this recompute_residual_factor
    // since the last time it was computed from scratch, recompute it from
    // scratch.  This helps to keep the computed residual accurate even in the
    // presence of roundoff.
    BaseFloat recompute_residual_factor;
    
    LinearCgdOptions(): max_iters(-1),
                        max_error(0.0),
                        recompute_residual_factor(0.01) { }
  };
    
  /*
    This function uses linear conjugate gradient descent to approximately solve
    the system A x = b.  The value of x at entry corresponds to the initial guess
    of x.  The algorithm continues until the number of iterations equals b.Dim(),
    or until the 2-norm of (A x - b) is <= max_error, or until the number of
    iterations equals max_iter, whichever happens sooner.  It is a requirement
    that A be positive definite.
    It returns the number of iterations that were actually executed (this is
    useful for testing purposes).
  */
  template<typename Real>
  int32 LinearCgd(const LinearCgdOptions &opts,
                  const SpMatrix<Real> &A, const VectorBase<Real> &b,
                  VectorBase<Real> *x);
  
  
  
  
  
  
  /**
     This is an implementation of L-BFGS.  It pushes responsibility for
     determining when to stop, onto the user.  There is no call-back here:
     everything is done via calls to the class itself (see the example in
     matrix-lib-test.cc).  This does not implement constrained L-BFGS, but it will
     handle constrained problems correctly as long as the function approaches
     +infinity (or -infinity for maximization problems) when it gets close to the
     bound of the constraint.  In these types of problems, you just let the
     function value be +infinity for minimization problems, or -infinity for
     maximization problems, outside these bounds).
  */
  
  struct LbfgsOptions {
    bool minimize; // if true, we're minimizing, else maximizing.
    int m; // m is the number of stored vectors L-BFGS keeps.
    float first_step_learning_rate; // The very first step of L-BFGS is
    // like gradient descent.  If you want to configure the size of that step,
    // you can do it using this variable.
    float first_step_length; // If this variable is >0.0, it overrides
    // first_step_learning_rate; on the first step we choose an approximate
    // Hessian that is the multiple of the identity that would generate this
    // step-length, or 1.0 if the gradient is zero.
    float first_step_impr; // If this variable is >0.0, it overrides
    // first_step_learning_rate; on the first step we choose an approximate
    // Hessian that is the multiple of the identity that would generate this
    // amount of objective function improvement (assuming the "real" objf
    // was linear).
    float c1; // A constant in Armijo rule = Wolfe condition i)
    float c2; // A constant in Wolfe condition ii)
    float d; // An amount > 1.0 (default 2.0) that we initially multiply or
    // divide the step length by, in the line search.
    int max_line_search_iters; // after this many iters we restart L-BFGS.
    int avg_step_length; // number of iters to avg step length over, in
    // RecentStepLength().
    
    LbfgsOptions (bool minimize = true):
        minimize(minimize),
        m(10),
        first_step_learning_rate(1.0),
        first_step_length(0.0),
        first_step_impr(0.0),
        c1(1.0e-04),
        c2(0.9),
        d(2.0),
        max_line_search_iters(50),
        avg_step_length(4) { }
  };
  
  template<typename Real>
  class OptimizeLbfgs {
   public:
    /// Initializer takes the starting value of x.
    OptimizeLbfgs(const VectorBase<Real> &x,
                  const LbfgsOptions &opts);
    
    /// This returns the value of the variable x that has the best objective
    /// function so far, and the corresponding objective function value if
    /// requested.  This would typically be called only at the end.
    const VectorBase<Real>& GetValue(Real *objf_value = NULL) const;
    
    /// This returns the value at which the function wants us
    /// to compute the objective function and gradient.
    const VectorBase<Real>& GetProposedValue() const { return new_x_; }
    
    /// Returns the average magnitude of the last n steps (but not
    /// more than the number we have stored).  Before we have taken
    /// any steps, returns +infinity.  Note: if the most recent
    /// step length was 0, it returns 0, regardless of the other
    /// step lengths.  This makes it suitable as a convergence test
    /// (else we'd generate NaN's).
    Real RecentStepLength() const;
    
    /// The user calls this function to provide the class with the
    /// function and gradient info at the point GetProposedValue().
    /// If this point is outside the constraints you can set function_value
    /// to {+infinity,-infinity} for {minimization,maximization} problems.
    /// In this case the gradient, and also the second derivative (if you call
    /// the second overloaded version of this function) will be ignored.
    void DoStep(Real function_value,
                const VectorBase<Real> &gradient);
    
    /// The user can call this version of DoStep() if it is desired to set some
    /// kind of approximate Hessian on this iteration.  Note: it is a prerequisite
    /// that diag_approx_2nd_deriv must be strictly positive (minimizing), or
    /// negative (maximizing).
    void DoStep(Real function_value,
                const VectorBase<Real> &gradient,
                const VectorBase<Real> &diag_approx_2nd_deriv);
    
   private:
    KALDI_DISALLOW_COPY_AND_ASSIGN(OptimizeLbfgs);
  
  
    // The following variable says what stage of the computation we're at.
    // Refer to Algorithm 7.5 (L-BFGS) of Nodecdal & Wright, "Numerical
    // Optimization", 2nd edition.
    // kBeforeStep means we're about to do
    /// "compute p_k <-- - H_k \delta f_k" (i.e. Algorithm 7.4).
    // kWithinStep means we're at some point within line search; note
    // that line search is iterative so we can stay in this state more
    // than one time on each iteration.
    enum ComputationState {
      kBeforeStep,
      kWithinStep, // This means we're within the step-size computation, and
      // have not yet done the 1st function evaluation.
    };
    
    inline MatrixIndexT Dim() { return x_.Dim(); }
    inline MatrixIndexT M() { return opts_.m; }
    SubVector<Real> Y(MatrixIndexT i) {
      return SubVector<Real>(data_, (i % M()) * 2); // vector y_i
    }
    SubVector<Real> S(MatrixIndexT i) {
      return SubVector<Real>(data_, (i % M()) * 2 + 1); // vector s_i
    }
    // The following are subroutines within DoStep():
    bool AcceptStep(Real function_value,
                    const VectorBase<Real> &gradient);
    void Restart(const VectorBase<Real> &x,
                 Real function_value,
                 const VectorBase<Real> &gradient);
    void ComputeNewDirection(Real function_value,
                             const VectorBase<Real> &gradient);
    void ComputeHifNeeded(const VectorBase<Real> &gradient);
    void StepSizeIteration(Real function_value,
                           const VectorBase<Real> &gradient);
    void RecordStepLength(Real s);
    
    
    LbfgsOptions opts_;
    SignedMatrixIndexT k_; // Iteration number, starts from zero.  Gets set back to zero
    // when we restart.
    
    ComputationState computation_state_;
    bool H_was_set_; // True if the user specified H_; if false,
    // we'll use a heuristic to estimate it.
  
  
    Vector<Real> x_; // current x.
    Vector<Real> new_x_; // the x proposed in the line search.
    Vector<Real> best_x_; // the x with the best objective function so far
                          // (either the same as x_ or something in the current line search.)
    Vector<Real> deriv_; // The most recently evaluated derivative-- at x_k.
    Vector<Real> temp_;
    Real f_; // The function evaluated at x_k.
    Real best_f_; // the best objective function so far.
    Real d_; // a number d > 1.0, but during an iteration we may decrease this, when
    // we switch between armijo and wolfe failures.
  
    int num_wolfe_i_failures_; // the num times we decreased step size.
    int num_wolfe_ii_failures_; // the num times we increased step size.
    enum { kWolfeI, kWolfeII, kNone } last_failure_type_; // last type of step-search
    // failure on this iter.
    
    Vector<Real> H_; // Current inverse-Hessian estimate.  May be computed by this class itself,
    // or provided by user using 2nd form of SetGradientInfo().
    Matrix<Real> data_; // dimension (m*2) x dim.  Even rows store
    // gradients y_i, odd rows store steps s_i.
    Vector<Real> rho_; // dimension m; rho_(m) = 1/(y_m^T s_m), Eq. 7.17.
  
    std::vector<Real> step_lengths_; // The step sizes we took on the last
    // (up to m) iterations; these are not stored in a rotating buffer but
    // are shifted by one each time (this is more convenient when we
    // restart, as we keep this info past restarting).
    
  
  };
    
  /// @} 
  
  
  } // end namespace kaldi
  
  
  
  #endif