reducedalgorithm.h

Go to the documentation of this file.

#ifndef REDUCED_ALGORITHM_H_
#define REDUCED_ALGORITHM_H_

#include "optproblemcontainer.h"
#include "reducedprobleminterface.h"
#include "dopeexceptionhandler.h"
#include "outputhandler.h"
#include "controlvector.h"

#include <deal.II/lac/vector.h>

#include <iostream>
#include <assert.h>
#include <iomanip>

namespace DOpE
{

  template<typename PROBLEM, typename VECTOR>
    class ReducedAlgorithm
  {
  public:
    ReducedAlgorithm(PROBLEM* OP,
                     ReducedProblemInterface<PROBLEM,VECTOR>* S,
                     ParameterReader &param_reader,
                     DOpEExceptionHandler<VECTOR>* Except = NULL,
                     DOpEOutputHandler<VECTOR>* Output = NULL, int base_priority = 0);
    virtual ~ReducedAlgorithm();
    
    static void
      declare_params(ParameterReader &param_reader);
    
    virtual void
      ReInit()
    {
      this->GetReducedProblem()->ReInit();
      if (rem_output_)
      {
        this->GetOutputHandler()->ReInit();
      }
    }
    
    virtual int
      Solve(ControlVector<VECTOR>& q, double global_tol = -1.)=0;
    
    virtual void
      SolveForward(ControlVector<VECTOR>& q);
    
     virtual void
      CheckGrads(double c, ControlVector<VECTOR>& q,
                 ControlVector<VECTOR>& dq, unsigned int niter = 1, double eps = 1.);
     virtual void
      FirstDifferenceQuotient(double exact, double eps,
                              const ControlVector<VECTOR>& q, const ControlVector<VECTOR>& dq);
    virtual void
      CheckHessian(double c, ControlVector<VECTOR>& q,
                   ControlVector<VECTOR>& dq, unsigned int niter = 1, double eps = 1.);
    virtual void
      SecondDifferenceQuotient(double exact, double eps,
                               const ControlVector<VECTOR>& q, const ControlVector<VECTOR>& dq);
    
    DOpEExceptionHandler<VECTOR>*
      GetExceptionHandler()
    {
      return ExceptionHandler_;
    }
    DOpEOutputHandler<VECTOR>*
      GetOutputHandler()
    {
      return OutputHandler_;
        }
  protected:
    PROBLEM*
      GetProblem()
    {
      return OP_;
    }
    const PROBLEM*
      GetProblem() const
    {
      return OP_;
    }
    const  ReducedProblemInterface<PROBLEM,VECTOR>*
      GetReducedProblem() const
    {
      return Solver_;
    }
    ReducedProblemInterface<PROBLEM,VECTOR>*
      GetReducedProblem()
    {
          return Solver_;
    }
    
    int
      GetBasePriority() const
    {
      return base_priority_;
    }
    
  private:
    PROBLEM* OP_;
    ReducedProblemInterface<PROBLEM,VECTOR>* Solver_;
    DOpEExceptionHandler<VECTOR>* ExceptionHandler_;
    DOpEOutputHandler<VECTOR>* OutputHandler_;
    bool rem_exception_;
    bool rem_output_;
    int base_priority_;
  };
  
  /***************************************************************************************/
  /****************************************IMPLEMENTATION*********************************/
  /***************************************************************************************/
  using namespace dealii;

  /******************************************************/

  template<typename PROBLEM, typename VECTOR>
    void
    ReducedAlgorithm<PROBLEM, VECTOR>::declare_params(
        ParameterReader &param_reader)
    {
      DOpEOutputHandler<VECTOR>::declare_params(param_reader);
    }

  /******************************************************/

  template<typename PROBLEM, typename VECTOR>
    ReducedAlgorithm<PROBLEM, VECTOR>::ReducedAlgorithm(
        PROBLEM* OP,  ReducedProblemInterface<PROBLEM,VECTOR>* S,
        ParameterReader &param_reader, DOpEExceptionHandler<VECTOR>* Except,
        DOpEOutputHandler<VECTOR>* Output, int base_priority)
    {
      assert(OP);
      assert(S);

      OP_ = OP;
      Solver_ = S;
      if (Output == NULL)
      {
        OutputHandler_ = new DOpEOutputHandler<VECTOR>(S, param_reader);
        rem_output_ = true;
      }
      else
      {
        OutputHandler_ = Output;
        rem_output_ = false;
      }
      if (Except == NULL)
      {
        ExceptionHandler_ = new DOpEExceptionHandler<VECTOR>(OutputHandler_);
        rem_exception_ = true;
      }
      else
      {
        ExceptionHandler_ = Except;
        rem_exception_ = false;
      }
      OP_->RegisterOutputHandler(OutputHandler_);
      OP_->RegisterExceptionHandler(ExceptionHandler_);
      Solver_->RegisterOutputHandler(OutputHandler_);
      Solver_->RegisterExceptionHandler(ExceptionHandler_);

      base_priority_ = base_priority;
    }

  /******************************************************/

  template<typename PROBLEM, typename VECTOR>
    ReducedAlgorithm<PROBLEM, VECTOR>::~ReducedAlgorithm()
    {
      if (ExceptionHandler_ && rem_exception_)
      {
        delete ExceptionHandler_;
      }
      if (OutputHandler_ && rem_output_)
      {
        delete OutputHandler_;
      }
    }
  /******************************************************/

  template<typename PROBLEM, typename VECTOR>
    void
    ReducedAlgorithm<PROBLEM, VECTOR>::SolveForward(
        ControlVector<VECTOR>& q)
    {
      q.ReInit();

      //Solve j'(q) = 0
      double cost = 0.;
      std::stringstream out;

      out << "**************************************************\n";
      out << "*             Starting Forward Solver            *\n";
      out << "*   Solving : " << this->GetProblem()->GetName() << "\t*\n";
      out << "*  CDoFs : ";
      q.PrintInfos(out);
      out << "*  SDoFs : ";
      this->GetReducedProblem()->StateSizeInfo(out);
      out << "**************************************************";
      this->GetOutputHandler()->Write(out, 1 + this->GetBasePriority(), 1, 1);

      try
      {
        cost = this->GetReducedProblem()->ComputeReducedCostFunctional(q);
      }
      catch (DOpEException& e)
      {
        this->GetExceptionHandler()->HandleCriticalException(e);
      }

      out << "CostFunctional: " << cost;
      this->GetOutputHandler()->Write(out, 2 + this->GetBasePriority());

      try
      {
        this->GetReducedProblem()->ComputeReducedFunctionals(q);
      }
      catch (DOpEException& e)
      {
        this->GetExceptionHandler()->HandleCriticalException(e);
      }

    }
  /******************************************************/

  template<typename PROBLEM, typename VECTOR>
    void
    ReducedAlgorithm<PROBLEM, VECTOR>::CheckGrads(double c,
        ControlVector<VECTOR>& q, ControlVector<VECTOR>& dq, unsigned int niter,
        double eps)
    {
      q.ReInit();
      dq.ReInit();

      dq = c;

      ControlVector<VECTOR> point(q);
      point = q;
      std::stringstream out;

      ControlVector<VECTOR> gradient(q), gradient_transposed(q);

      this->GetReducedProblem()->ComputeReducedCostFunctional(point);
      this->GetReducedProblem()->ComputeReducedGradient(point, gradient,
          gradient_transposed);
      double cost_diff = gradient * dq;
      out << "Checking Gradients...." << std::endl;
      out << " Epsilon \t Exact \t Diff.Quot. \t Rel. Error ";
      this->GetOutputHandler()->Write(out, 3 + this->GetBasePriority());

      for (unsigned int i = 0; i < niter; i++)
      {
        FirstDifferenceQuotient(cost_diff, eps, q, dq);
        eps /= 10.;
      }
    }
  /******************************************************/

  template<typename PROBLEM, typename VECTOR>
    void
    ReducedAlgorithm<PROBLEM, VECTOR>::FirstDifferenceQuotient(
        double exact, double eps, const ControlVector<VECTOR>& q,
        const ControlVector<VECTOR>& dq)
    {
      ControlVector<VECTOR> point(q);
      point = q;

      std::stringstream out;

      point.add(eps, dq);

      double cost_right = 0.;
      //Differenzenquotient
      cost_right = this->GetReducedProblem()->ComputeReducedCostFunctional(
          point);

      point.add(-2. * eps, dq);

      double cost_left = 0.;
      //Differenzenquotient
      cost_left = this->GetReducedProblem()->ComputeReducedCostFunctional(
          point);

      double diffquot = (cost_right - cost_left) / (2. * eps);
      out << eps << "\t" << exact << "\t" << diffquot << "\t"
          << (exact - diffquot) / exact << std::endl;
      this->GetOutputHandler()->Write(out, 3 + this->GetBasePriority());
    }

  /******************************************************/

  template<typename PROBLEM, typename VECTOR>
    void
    ReducedAlgorithm<PROBLEM, VECTOR>::CheckHessian(double c,
        ControlVector<VECTOR>& q, ControlVector<VECTOR>& dq, unsigned int niter,
        double eps)
    {
      q.ReInit();
      dq.ReInit();

      dq = c;

      ControlVector<VECTOR> point(q);
      point = q;
      std::stringstream out;

      ControlVector<VECTOR> gradient(q), gradient_transposed(q), hessian(q),
          hessian_transposed(q);

      this->GetReducedProblem()->ComputeReducedCostFunctional(point);
      this->GetReducedProblem()->ComputeReducedGradient(point, gradient,
          gradient_transposed);

      this->GetReducedProblem()->ComputeReducedHessianVector(point, dq, hessian,
          hessian_transposed);

      double cost_diff = hessian * dq;
      out << "Checking Hessian...." << std::endl;
      out << " Epsilon \t Exact \t Diff.Quot. \t Rel. Error ";
      this->GetOutputHandler()->Write(out, 3 + this->GetBasePriority());

      for (unsigned int i = 0; i < niter; i++)
      {
        SecondDifferenceQuotient(cost_diff, eps, q, dq);
        eps /= 10.;
      }
    }
  /******************************************************/

  template<typename PROBLEM, typename VECTOR>
    void
    ReducedAlgorithm<PROBLEM, VECTOR>::SecondDifferenceQuotient(
        double exact, double eps, const ControlVector<VECTOR>& q,
        const ControlVector<VECTOR>& dq)
    {
      ControlVector<VECTOR> point(q);
      point = q;
      std::stringstream out;

      double cost_mid = 0.;
      //Differenzenquotient
      cost_mid = this->GetReducedProblem()->ComputeReducedCostFunctional(point);

      point.add(eps, dq);

      double cost_right = 0.;
      //Differenzenquotient
      cost_right = this->GetReducedProblem()->ComputeReducedCostFunctional(
          point);

      point.add(-2. * eps, dq);

      double cost_left = 0.;
      //Differenzenquotient
      cost_left = this->GetReducedProblem()->ComputeReducedCostFunctional(
          point);

      double diffquot = (cost_left - 2. * cost_mid + cost_right) / (eps * eps);

      out << eps << "\t" << exact << "\t" << diffquot << "\t"
          << (exact - diffquot) / exact << std::endl;
      this->GetOutputHandler()->Write(out, 3 + this->GetBasePriority());
    }

/******************************************************/

}
#endif