#include "GridSpace.hh"

namespace RVL {

  std::vector<Grid *> make_asg_ctrllist(Grid const & phys, int maxoff);
  std::vector<Grid *> make_asg_statelist(Grid const & phys, int maxoff);
  std::map<std::string,int> make_asg_indices(int dim);
  float * sgcoeffs(int maxoff);

  class ASGaux: public Writeable {
  private:
    void initialize(GridDomain const & gdom,
		    int maxoff, float dt, float amp);
    ASGaux();
  public:
    GridDomain const & gdom;
    float dt;
    float amp;
    IPNT lbc;
    IPNT rbc;
    float * vdiv_alloc;
    float * vdiv;
    float ** pgrad_alloc;
    float ** pgrad;
    IPNT nls;
    IPNT nrs;
    IPNT nlsloc;
    IPNT nrsloc;
    float ** ep;
    float ** epp;
    float ** ev;
    float ** evp;
    float ** ep_alloc;
    float ** epp_alloc;
    float ** ev_alloc;
    float ** evp_alloc;
    float ** coeffs;
    int maxoff;

    ASGaux(GridDomain const & gdom,
	   int maxoff, float dt, float amp);
    ASGaux(ASGaux const & aux);
    ~ASGaux();

    ostream & write(ostream & str) const;
    
  };
    
  class ASGapplyFO: public TSFO {
  private:
    ASGaux const & aux;
    bool const & integerstep;
  public:
    ASGapplyFO(ASGaux const & _aux,
	       bool const & _integerstep)
      : aux(_aux),
	integerstep(_integerstep) {}
    ASGapplyFO(ASGapplyFO const & f)
      : aux(f.aux), integerstep(f.integerstep) {}
    ~ASGapplyFO() {}

    void operator()(TSDC & y) const;

    std::string getName() const { std::string tmp = "ASGapplyFO"; return tmp;}
  };
    
  class ASGapplyAdjFO: public TSFO {
  private:
    ASGaux const & aux;
    bool const & integerstep;
  public:
    ASGapplyAdjFO(ASGaux const & _aux,
	       bool const & _integerstep)
      : aux(_aux),
	integerstep(_integerstep) {}
    ASGapplyAdjFO(ASGapplyAdjFO const & f)
      : aux(f.aux), integerstep(f.integerstep) {}
    ~ASGapplyAdjFO() {}

    void operator()(TSDC & y) const;

    std::string getName() const { std::string tmp = "ASGapplyAdjFO"; return tmp;}
  };
  
  class ASGapplyTangentFO: public TSFO {
  private:
    ASGaux const & aux;
    bool const & integerstep;
  public:
    ASGapplyTangentFO(ASGaux const & _aux,
	       bool const & _integerstep)
      : aux(_aux),
	integerstep(_integerstep) {}
    ASGapplyTangentFO(ASGapplyTangentFO const & f)
      : aux(f.aux), integerstep(f.integerstep) {}
    ~ASGapplyTangentFO() {}

    void operator()(TSDC & y) const;

    std::string getName() const { std::string tmp = "ASGapplyTangentFO"; return tmp;}
  };

  class ASGapplyAdjTangentFO: public TSFO{
  private:
    ASGaux const & aux;
    bool const & integerstep;
  public:
    ASGapplyAdjTangentFO(ASGaux const & _aux,
	       bool const & _integerstep)
      : aux(_aux),
	integerstep(_integerstep) {}
    ASGapplyAdjTangentFO(ASGapplyAdjTangentFO const & f)
      : aux(f.aux), integerstep(f.integerstep) {}
    ~ASGapplyAdjTangentFO() {}

    void operator()(TSDC & y) const;

    std::string getName() const { std::string tmp = "ASGapplyAdjTangentFO"; return tmp;}
  };
  
  class ASGrandAccessFO: public TSRAFO<float> {
  private:
    ASGaux const & aux;
    bool const & integerstep;
    ASGrandAccessFO();
  public:
    ASGrandAccessFO(ASGaux const & _aux,
		    bool const & _integerstep)
      : aux(_aux),
	integerstep(_integerstep) {}
    ASGrandAccessFO(ASGrandAccessFO const & f)
      : aux(f.aux), integerstep(f.integerstep) {}
    ~ASGrandAccessFO() {}

    void operator()(TSDC & y) const;

    std::string getName() const { std::string tmp = "ASGrandAccessFO"; return tmp;}
  };
  
  class ASGStep: public TSStep<float> {

  private:

    GridDomain gdom;

    float dt;
    float tmin;

    ASGaux aux;
    
    mutable bool integerstep;
    
    mutable ASGapplyFO app;
    mutable ASGapplyAdjFO appAdj;
    mutable ASGapplyTangentFO appTan;
    mutable ASGapplyAdjTangentFO appAdjTan;
    mutable ASGrandAccessFO randAcc;

    TSFO & applyFO() const { return app; }
    TSFO & applyAdjFO() const { return appAdj; }
    TSFO & applyTangentFO() const { return appTan; }
    TSFO & applyAdjTangentFO() const { return appAdjTan; }
    TSRAFO<float> & randAccessFO() const { return randAcc; }
    void stepTimeFwd() const {
      this->stepTime(+0.5*dt);
      integerstep = !integerstep;
    } 
    void stepTimeBwd() const {
      this->stepTime(-0.5*dt);
      integerstep = !integerstep;
    } 
    
  protected:

    Operator<float> * clone() const { return new ASGStep(*this); }
    
  public:

    ASGStep(ASGStep const & stp)
      : gdom(stp.gdom),tmin(stp.tmin),dt(stp.dt),
	aux(stp.aux),
	app(aux,integerstep),
	appAdj(aux,integerstep),
	appTan(aux,integerstep),
	appAdjTan(aux,integerstep),
	randAcc(aux,integerstep) {

	//	randAcc(stp.randAcc) {    
      try {
	this->setTime(tmin);
      }
      catch (RVLException & e) {
	e<<"\ncalled from ASGStep constructor\n";
	throw e;
      }
    }

    ASGStep(Grid const & phys, int maxoff, float _tmin, float _dt, float amp)
      : gdom(make_asg_ctrllist(phys,maxoff),
	     make_asg_statelist(phys,maxoff),
	     make_asg_indices(phys.getDimension())),
	tmin(_tmin), dt(_dt),
	aux(gdom,maxoff,dt,amp),
	integerstep(true),
	app(aux,integerstep),
	appAdj(aux,integerstep),
	appTan(aux,integerstep),
	appAdjTan(aux,integerstep),
	randAcc(aux,integerstep){
      try {
	this->setTime(tmin);
      }
      catch (RVLException & e) {
	e<<"\ncalled from ASGStep constructor\n";
	throw e;
      }
    }

    Space<float> const & getDomain() const { return gdom; }
    Space<float> const & getRange() const { return gdom; }

    ostream & write(ostream & str) const {
      str<<"ASGStep\n";
      return str;
    }
  };

}