newton_base.hh

Go to the documentation of this file.

#ifndef NEWTON_BASE_HH
#define NEWTON_BASE_HH
  
#include <vector>
#include <memory>
 
#include "abstract_interface.hh"
#include "algorithm_base.hh"
 
namespace Kaskade
{
 
//------------------------------------------------------------------------------------------
 
class NewtonParameters : public IterationParameters
{
public:
  NewtonParameters(double desiredAccuracy_, int maxSteps_)
    : IterationParameters(desiredAccuracy_, maxSteps_), reuseFactorization(false)
  {
    reset();
  }
 
 
//Parameters with values that must be supplied by client
//  double desiredAccuracy;
//  int maxSteps;
 
  virtual ~NewtonParameters() {}
  
  virtual void reset() {
    doForAll(LQAction::Reset);
  }
 
  bool reuseFactorization;
 
protected:
  friend class NewtonsMethod;
  friend class SimpleNewtonMethod;
 
  void logStep() { 
    doForAll(LQAction::LogValue);
  }
  
  virtual void doForAll(LQAction::ToDo td)
  {
    normCorr.doAction(td,"|Corr|");
  }
 
  LoggedQuantity<double> normCorr;
};
 
class IdentityChart : public AbstractChart
{
public:
  void addPerturbation(AbstractFunctionSpaceElement& newIterate, 
                       AbstractFunctionSpaceElement const& perturbation, 
                       AbstractLinearization const& lin,
                       std::vector<std::shared_ptr<AbstractFunctionSpaceElement> > basis = std::vector<std::shared_ptr<AbstractFunctionSpaceElement> >()) const
  {
    newIterate = lin.getOrigin();
    newIterate += perturbation;
  }
};
 
//class TrivialAbstractChart : public AbstractChart
//{
//public:
//  void addPerturbation(AbstractFunctionSpaceElement& newIterate,
//                       AbstractFunctionSpaceElement const& perturbation,
//                       AbstractLinearization const& lin,
//                       std::vector<AbstractFunctionSpaceElement* > basis = std::vector<AbstractFunctionSpaceElement* >()) const
//  {
//    newIterate = lin.getOrigin();
//    newIterate += perturbation;
//  }
//};
 
 
class NewtonsMethod : public Algorithm
{
public:
  typedef NewtonParameters Parameters; 
 
  virtual ~NewtonsMethod() {}
 
  NewtonsMethod(AbstractNewtonDirection& l, AbstractChart& chart_, AbstractNorm& n, Parameters& p_): 
    linearSolver(l), chart(chart_), norm(n), p(p_) {}
 
  void solve(AbstractFunctional* f, AbstractFunctionSpaceElement& x);
 
  void oneStep(AbstractFunctional* f, AbstractLinearization* l, AbstractFunctionSpaceElement& x);
 
  void resolve(AbstractFunctionSpaceElement& x, AbstractLinearization const&  l) 
  { 
    if(p.reuseFactorization)
      linearSolver.simplified(x,l); 
    else
    {
      std::cout << "No Factorization available: set NewtonParameters::reuseFactorization=true" << std::endl;
    }
}
 
  virtual Parameters& getParameters() {return p;}
  void setDesiredAccuracy(double da) { assert(da>0); p.desiredAccuracy=da;}
  virtual void setDesiredRelativeAccuracy(double ra) { assert(false); }
  void resetParameters() {p.reset();}
 
  int stepsPerformed() {return step;}
 
  bool changedGrid() { return gridHasChanged; }
 
  virtual void getSearchDirection(AbstractFunctionSpaceElement& direction) 
  {
    linearSolver.ordinary(direction, *newtonLinearization);
    direction *= -1.0;    
  }
 
 
  virtual void computeTrialIterate(AbstractFunctionSpaceElement& trialIterate, AbstractFunctionSpaceElement const& direction, AbstractLinearization const& lin)
  {
          chart.addPerturbation(trialIterate, direction, lin);
  }
 
  AbstractLinearization const& getLastLinearization() { assert(newtonLinearization.get()); return *newtonLinearization; }
 
  AbstractLinearization const& getLastSimplifiedLinearization() { assert(simplifiedLinearization.get()); return *simplifiedLinearization; }
 
protected:
 
  virtual void initialize();
  virtual void finalize(int);
  virtual void logQuantities() { p.logStep();}
 
  virtual void initNewtonStep() {};
 
  virtual Convergence convergenceTest(AbstractFunctionSpaceElement const& correction, AbstractFunctionSpaceElement const& iterate) = 0;
 
  virtual void predictNextDampingFactor(AbstractFunctionSpaceElement& correction) {}
 
  virtual double& dampingFactor() { dummy=1.0; return dummy;}
  virtual double maxSteps() {return p.maxSteps;}
 
  virtual void updateIterate(AbstractFunctionSpaceElement& iterate, 
                             AbstractFunctionSpaceElement& trialIterate,
                             AbstractLinearization const& lin);
 
  virtual RegularityTest regularityTest(double scalingFactor)
  {
    return RegularityTest::Passed;
  }
 
  virtual AcceptanceTest evaluateTrialIterate(
    AbstractFunctionSpaceElement const& trialIterate, 
    AbstractFunctionSpaceElement const& correction, 
    AbstractLinearization const& lin) = 0;
 
 
  template<class T>
  class CycleDetection
  {
  public:
    CycleDetection() : hasIncreased(false) {};
 
    bool detectCycling(T& quantity)
    {
      if(quantityCycle.size()!=0)
        if(*std::min_element(quantityCycle.begin(),quantityCycle.end())<quantity) hasIncreased=true;
 
      if(hasIncreased && quantity < *std::min_element(quantityCycle.begin(),quantityCycle.end()))
      {
        quantity = *std::min_element(quantityCycle.begin(),quantityCycle.end());
        std::cout << "Detected Cycling: decreased below already accepted step size:" << quantity << std::endl;
        return true;
      }
 
      for(int i = 0; i< quantityCycle.size(); i++)
        if(std::abs(quantityCycle[i]-quantity)<std::abs(quantity)*std::numeric_limits<double>::epsilon()){
          quantity = quantityCycle[i+1];
          std::cout << "Detected Cycling: two identical step sizes:" << quantity << std::endl;
          return true;
        }
 
      if(quantityCycle.size() > 20){
        quantity=quantity/4.0;
        std::cout << "Detected Cycling: too many trials:" << quantity << std::endl;
        return true;
      }
      quantityCycle.push_back(quantity);
      return false;
    }
 
  private:
    std::vector<T> quantityCycle;
    bool hasIncreased;
  };
 
  virtual void terminationMessage(int flag);
 
  AbstractFunctional* functional;
  AbstractNewtonDirection& linearSolver;
  AbstractChart& chart;
  AbstractNorm& norm;
  Parameters& p;
 
  std::unique_ptr<AbstractFunctionSpaceElement> iterate, trialIterate, correction;
 
  std::unique_ptr<AbstractLinearization> newtonLinearization;
  std::unique_ptr<AbstractLinearization> simplifiedLinearization;
  AbstractLinearization* newtonPtr;
 
  int step;
 
private:
  int runAlgorithm();
 
  NewtonsMethod(NewtonsMethod&);
 
  double dummy;
 
 
  bool doOneStep;
  bool gridHasChanged;
};
 
 
 
}  // namespace Kaskade
#endif