#include "GridSpace.hh"
#include "asg_statespace.hh"
#include "gtest/gtest.h"

namespace {

  class GridTest : public ::testing::Test {
  protected:
    size_t n1;
    size_t n2;
    size_t n3;
    float o1;
    float o2;
    float o3;
    float d1;
    float d2;
    float d3;
    int id1;
    int id2;
    int id3;
    
    GridTest()
      : n1(10), n2(20), n3(5),
	o1(0.0), o2(10.0), o3(5.0),
	d1(20.0), d2(10.0), d3(10.0),
	id1(0), id2(1), id3(2)
    {}
  };

  TEST_F(GridTest, axis_constructor) {
    RVL::Axis a1(n1,d1,o1,id1);
    std::stringstream t;
    a1.write(t);
    EXPECT_EQ(t.str(),"n=10 d=20 o=0 id=0 atype=0 gs0=0 ge0=9 gs=0 ge=9\n");
  }

  TEST_F(GridTest, axis_compare) {
    RVL::Axis a1(n1,d1,o1,id1);
    RVL::Axis a2(n2,d2,o2,id2);
    bool res = (a1==a2);
    EXPECT_EQ(0,res);
  }

  TEST_F(GridTest,grid_constructor_2D) {
    size_t ndim=2;
    RVL::Grid g(ndim);
    {
      RVL::Axis a(n1,d1,o1,id1);
      g.addAxis(a);
    }
    {
      RVL::Axis a(n2,d2,o2,id2);
      g.addAxis(a);
    }
    IPNT gs0;
    IPNT ge0;
    IPNT gs;
    IPNT ge;
    g.getAxisLimits(gs0,ge0);
    g.getSubAxisLimits(gs,ge);
    EXPECT_EQ(gs0[0],0);
    EXPECT_EQ(ge0[0],9);
    EXPECT_EQ(gs0[1],1);
    EXPECT_EQ(ge0[1],20);
    EXPECT_EQ(gs0[0],gs[0]);
    EXPECT_EQ(ge0[0],ge[0]);
    EXPECT_EQ(gs0[1],gs[1]);
    EXPECT_EQ(ge0[1],ge[1]);
    gs[0]=2; ge[0]=7;
    gs[1]=3; ge[1]=18;
    g.setSubAxisLimits(gs,ge);
    g.getSubAxisLimits(gs,ge);
    EXPECT_EQ(2,gs[0]);
    EXPECT_EQ(7,ge[0]);
    EXPECT_EQ(3,gs[1]);
    EXPECT_EQ(18,ge[1]);
  }

  TEST_F(GridTest,griddc_constructor_2D) {
    size_t ndim=2;
    RVL::Grid g(ndim);
    {
      RVL::Axis a(n1,d1,o1,id1);
      g.addAxis(a);
    }
    {
      RVL::Axis a(n2,d2,o2,id2);
      g.addAxis(a);
    }
    RVL::GridDC d(g);
    RVL::RVLRandomize<float> rnd(getpid(),-1.0f,1.0f);
    rnd(d);
    IPNT gs0;
    IPNT ge0;
    size_t stride = g.getAxis(0).getLength();
    g.getAxisLimits(gs0,ge0);
    float ** d2d = d.getData2D();
    d2d -= gs0[1];
    for (size_t i1=gs0[1]; i1<=ge0[1]; i1++) 
      d2d[i1] -= gs0[0];
    for (size_t i1=gs0[1]; i1<=ge0[1]; i1++) {
      for (size_t i0=gs0[0]; i0<=ge0[0]; i0++) {
	EXPECT_EQ(d2d[i1][i0],d.getData()[i0-gs0[0] + (i1-gs0[1])*stride]);
      }
    }    
  }

  TEST_F(GridTest,griddc_constructor_2D_global) {
    size_t ndim=2;
    RVL::Grid g(ndim);
    {
      RVL::Axis a(n1,d1,o1,id1);
      g.addAxis(a);
    }
    {
      RVL::Axis a(n2,d2,o2,id2);
      g.addAxis(a);
    }
    RVL::GridDC d(g);
    RVL::RVLRandomize<float> rnd(getpid(),-1.0f,1.0f);
    rnd(d);
    IPNT gs0;
    IPNT ge0;
    g.getAxisLimits(gs0,ge0);
    int stride = ge0[0]-gs0[0]+1;
    float ** d2d = d.getData2DGlobal();
    for (size_t i1=gs0[1]; i1<=ge0[1]; i1++) {
      for (size_t i0=gs0[0]; i0<=ge0[0]; i0++) {
	EXPECT_EQ(d2d[i1][i0],d.getData()[i0-gs0[0] + (i1-gs0[1])*stride]);
      }
    }    
  }

  TEST_F(GridTest,griddc_constructor_3D) {
    size_t ndim=3;
    RVL::Grid g(ndim);
    {
      RVL::Axis a(n1,d1,o1,id1);
      g.addAxis(a);
    }
    {
      RVL::Axis a(n2,d2,o2,id2);
      g.addAxis(a);
    }
    {
      RVL::Axis a(n3,d3,o3,id3);
      g.addAxis(a);
    }
    //    g.write(cerr);
    //    cerr<<"getdatasize="<<getDataSize(g)<<"\n";
    RVL::GridDC d(g);
    RVL::RVLRandomize<float> rnd(getpid(),-1.0f,1.0f);
    rnd(d);
    IPNT gs0;
    IPNT ge0;
    g.getAxisLimits(gs0,ge0);
    int stride0 = ge0[0]-gs0[0]+1;
    int stride1 = ge0[1]-gs0[1]+1;
    float *** d3d = d.getData3D();
    for (size_t i2=gs0[2]; i2<=ge0[2]; i2++) {
      //      cerr<<"i2="<<i2<<"\n";
      //      cerr<<"d3d[i2-gs0[2]]="<<d3d[i2-gs0[2]]<<"\n";
      for (size_t i1=gs0[1]; i1<=ge0[1]; i1++) {
	for (size_t i0=gs0[0]; i0<=ge0[0]; i0++) {
	  size_t iii = i0-gs0[0] + (i1-gs0[1])*stride0
	    + (i2-gs0[2])*stride1*stride0;
	  EXPECT_EQ(d3d[i2-gs0[2]][i1-gs0[1]][i0-gs0[0]],d.getData()[iii]);
	}
      }    
    }
    float *** d3ds = d.getData3DGlobal();
    for (size_t i2=gs0[2]; i2<=ge0[2]; i2++) {
      //      cerr<<"i2="<<i2<<"\n";
      //      cerr<<"d3ds[i2]="<<d3ds[i2]<<"\n";
      for (size_t i1=gs0[1]; i1<=ge0[1]; i1++) {
	for (size_t i0=gs0[0]; i0<=ge0[0]; i0++) {
	  size_t iii = i0-gs0[0] + (i1-gs0[1])*stride0
	    + (i2-gs0[2])*stride1*stride0;
	  EXPECT_EQ(d3ds[i2][i1][i0],d.getData()[iii]);
	}
      }    
    }
  }

  TEST_F(GridTest,gridspace_build_const_vector) {
    size_t ndim=2;
    RVL::Grid g(ndim);
    {
      RVL::Axis a(n1,d1,o1,id1);
      g.addAxis(a);
    }
    {
      RVL::Axis a(n2,d2,o2,id2);
      g.addAxis(a);
    }
    RVL::GridSpace sp(g);
    RVL::Vector<float> x(sp);
    RVL::RVLAssignConst<float> ac(1.0f);
    x.eval(ac);
    EXPECT_EQ(x.inner(x),40000.0f);
  }

  TEST_F(GridTest,build_griddomain_2D) {
    size_t ndim=2;
    RVL::Grid phys(ndim);
    {
      RVL::Axis a(201,20.0,0.0,0);
      phys.addAxis(a);
    }
    {
      RVL::Axis a(401,20.0,0.0,0);
      phys.addAxis(a);
    }
    IPNT gs0; IPNT ge0; IPNT gs; IPNT ge;
    int maxoff=2;
    std::vector<RVL::Grid *> ctrllist
      =RVL::make_asg_ctrllist(phys,maxoff);
    std::vector<RVL::Grid *> statelist
      =RVL::make_asg_statelist(phys,maxoff);
    std::map<std::string,int> ind
      =RVL::make_asg_indices(phys.getDimension());
    RVL::GridDomain dom(ctrllist,statelist,ind);
    RVL::ProductSpace<float> const & ssp
      = dynamic_cast<RVL::ProductSpace<float> const &>(dom[1]);
    // get p0
    RVL::GridSpace const & arr2
      = dynamic_cast<RVL::GridSpace const &>(ssp[dom.get_ind()["p0"]]);
    arr2.getGrid().getAxisLimits(gs0,ge0);
    arr2.getGrid().getSubAxisLimits(gs,ge);
    // cerr<<"dim 0:  gs = "<<gs[0]<<"  ge = "<<ge[0] <<" gs0 = "<<gs0[0]
    //	<<"  ge0 = "<<ge0[0]<<"\n";
    // cerr<<"dim 1:  gs = "<<gs[1]<<"  ge = "<<ge[1] <<" gs0 = "<<gs0[1]
    //	<<"  ge0 = "<<ge0[1]<<"\n";
    EXPECT_EQ(gs0[0],-1);
    EXPECT_EQ(ge0[0],201);
    EXPECT_EQ(gs0[1],1);
    EXPECT_EQ(ge0[1],399);
    EXPECT_EQ(gs[0],1);
    EXPECT_EQ(ge[0],199);
    EXPECT_EQ(gs[1],1);
    EXPECT_EQ(ge[1],399);
    
  }
}