Newton.h

#ifndef OPTIMIZATION_NEWTON_H
#define OPTIMIZATION_NEWTON_H

#include <KrisLibrary/math/function.h>
#include <KrisLibrary/math/root.h>
#include <KrisLibrary/math/SVDecomposition.h>
#include <vector>

namespace Math {
  template <class T> class SparseMatrixTemplate_RM;
  typedef SparseMatrixTemplate_RM<Real> SparseMatrix;
}

namespace Optimization {
using namespace Math;

class NewtonRoot
{
public:
  NewtonRoot(VectorFieldFunction* func);
  virtual ~NewtonRoot();
  bool GlobalSolve(int& iters,ConvergenceResult* res=NULL);
  ConvergenceResult Solve(int& iters);
  ConvergenceResult Solve_Sparse(int& iters);
  bool LineMinimization(const Vector& g, const Vector& p, Real *f);
  Real MaxDistance(const Vector& x);
  virtual bool SolveUnderconstrainedLS(const Matrix& A,const Vector& b,Vector& x);
  virtual bool SolveUnderconstrainedLS(const SparseMatrix& A,const Vector& b,Vector& x);
  virtual Real Merit();  //evaluates the merit function at x

  Vector x;
  VectorFieldFunction* func;
  Real tolf,tolmin,tolx;
  Real stepMax;  
  Real lambda;   
  Vector bmin,bmax; 
  bool sparse;   
  Vector bias;   
  int verbose;
  int debug;

  //temporary
  RobustSVD<Real> svd;
  Vector fx, g, p, xold;
  Matrix fJx;
};

class ConstrainedNewtonRoot : public NewtonRoot
{
public:
  ConstrainedNewtonRoot(VectorFieldFunction* func,VectorFieldFunction* c);
  bool GlobalSolve(int& iters,ConvergenceResult* res=NULL);
  ConvergenceResult SolveConstrained(int& iters);
  ConvergenceResult SolveConstrained2(int& iters);
  ConvergenceResult SolveConstrained2_Sparse(int& iters);
  ConvergenceResult SolveConstrained_SLP(int& iters);
  virtual Real Merit();

  VectorFieldFunction* c;
  Real tolc;

  //temporary
  Vector cx;
  std::vector<int> activeSetC,activeSetBound;
  Matrix A;  //solve for A*dx = rhs
  Vector rhs;
};

Real SurfaceDistance(VectorFieldFunction*C,const Vector& x);
bool SatisfiesEquality(VectorFieldFunction*C,const Vector& x,Real tol=Epsilon);
bool SatisfiesInequality(VectorFieldFunction*C,const Vector& x,Real margin=Zero);
Real InequalityMargin(VectorFieldFunction* c,const Vector& x,int* index=NULL);

} //namespace Optimization

#endif