#include "gtest/gtest.h"
#include "acd_defn.hh"
#include "grid.h"
#include "iwsim.hh"
#include "adjsteptest.hh"

//#define GTEST_VERBOSE

IOKEY IWaveInfo::iwave_iokeys[]
= {
  {"csq",    0, true,  true },
  {"data",   1, false, true },
  {"source", 1, true,  false},
  {"movie",  1, false, false},
  {"init",   1, true,  false},
  {"",       0, false, false}
};

namespace {

  using RVL::parse;
  using RVL::RVLException;
  using TSOpt::IWaveEnvironment;
  using TSOpt::IWaveTree;
  using TSOpt::IWaveSim;
  using TSOpt::TASK_RELN;
  using TSOpt::IOTask;

#include "hfile3d.h"
#include "rdot.h"

  class ACDStepTest : public ::testing::Test {
  public:

    string hfile;
    string dfile;
    IWaveInfo ic;
    ACDStepTest(): ic() {

      hfile = "csq_4layer.rsf";
      dfile = "csq_4layer.rsf@";

      // create simple par file
      ofstream ps("parfile");
      ps<<"INPUT DATA FOR IWAVE\n";
      ps<<"------------------------------------------------------------------------\n";
      ps<<"FD:\n";
      ps<<"\n";
      ps<<"         order = 4           scheme half-order\n";
      ps<<"           cfl = 0.5        cfl number - frac of max stable\n";
      ps<<"          cmin = 1.0         min velocity - checked\n";
      ps<<"          cmax = 5.0         max velocity - checked\n";
      ps<<"      max_step = 0           1 = set adaptively, 0 = use standard cfl from cmax\n";
      ps<<"         fpeak = 0.010       nominal central frequency \n";
      ps<<"\n";
      ps<<"------------------------------------------------------------------------\n";
      ps<<"Model info:\n";
      ps<<"\n";
      ps<<"           csq = csq_4layer.rsf\n";
      ps<<"        csq_d1 = csq_4layer.rsf\n";
      ps<<"        csq_d2 = csq_4layer.rsf\n";
      ps<<"        csq_b1 = csq_4layer.rsf\n";
      ps<<"        csq_b2 = csq_4layer.rsf\n";
      ps<<"\n";
      ps<<"------------------------------------------------------------------------\n";
      ps<<"MPI info:\n";
      ps<<"\n";
      ps<<"       mpi_np1 = 1      n_doms along axis 1\n";
      ps<<"       mpi_np2 = 1      n_doms along axis 2\n";
      ps<<"       mpi_np3 = 1      n_doms along axis 3\n";
      ps<<"       partask = 1      task parallelization\n";
      ps<<"\n";
      ps<<"------------------------------------------------------------------------\n";
      ps<<"Source info:\n";
      ps<<"\n";
      ps<<"        source = wavelet_fake.su\n";
      ps<<"       sampord = 1             sampling order\n";
      ps<<"\n";
      ps<<"------------------------------------------------------------------------\n";
      ps<<"Trace info:\n";
      ps<<"\n";
      ps<<"            data = data_fake.su    output data file\n";
      ps<<"\"        movie1 = p\"\n";
      ps<<"\"     moviestep = 100\"\n";
      ps<<"\n";
      ps<<"------------------------------------------------------------------------\n";
      ps<<"Output info:\n";
      ps<<"\n";
      ps<<"     printact = 1           per-time-step verbosity level\n";
      ps<<"                            0 - none\n";
      ps<<"                            1 - time step index\n";
      ps<<"                            2 - internal time step info\n";
      ps<<"                            > 5: dump everything\n";
      ps<<"      dump_pi = 1           dump parallel/dom. decomp info\n";
      ps<<"     dump_lda = 1           dump grid data for allocated arrays\n";
      ps<<"     dump_ldc = 1           dump grid data for computational arrays\n";
      ps<<"     dump_lds = 1           dump grid data for send arrays\n";
      ps<<"     dump_ldr = 1           dump grid data for receive arrays\n";
      ps<<"    dump_term = 0           dump terminator data\n";
      ps<<"    dump_pars = 0           print parameter table in IWaveOp\n";
      ps<<"   dump_steps = 0           print major steps in IWaveOp\n";
      ps.flush();
      ps.close();

      grid g;
      g.dim=2;
      g.gdim=2;
      g.axes[0].n=416;
      g.axes[1].n=800;
      g.axes[2].n=4;
      g.axes[0].d=25.0;
      g.axes[1].d=25.0;
      g.axes[2].d=1.0;
      g.axes[0].o=0.0;
      g.axes[1].o=0.0;
      g.axes[2].o=0.0;
      g.axes[0].id = 0;
      g.axes[1].id = 1;
      g.axes[2].id = 2;

      float val = 1.5;
      create_hfile(hfile, dfile, g, val);
    }
  };
    
  TEST_F(ACDStepTest, acd_tsf_order0test) {
    try {
      // fake command line environment
      int argc = 2;
      char * argvv = new char[128];
      char ** argv = new char*[2];
      argv[0]=&(argvv[0]); argv[1]=&(argvv[65]);
      strcpy(argv[1],"par=parfile");
      PARARRAY * par = NULL;
      FILE * stream = NULL;
      IWaveEnvironment(argc, argv, 0, &par, &stream);
      delete [] argv;
      delete [] argvv;
            
      // build order zero IWaveTree, check
      int order=0;
      IWaveTree * fwd = new IWaveTree(*par,stream,ic,order);
      IWaveTree * adj = new IWaveTree(*par,stream,ic,order);
      std::vector<RDOM *> rdfwd, rdadj;
            
      IWAVE * w = (fwd->getStateArray()[0]);
            
      //      FD_MODEL *specs = (FD_MODEL *)(w->model.specs);
            
      rdfwd = fwd->getRDOMArray();
      rdadj = adj->getRDOMArray();
            
      for (int i=0; i<rdfwd.size(); i++) {
	rdom_rand(rdfwd[i]);
	rdom_rand(rdadj[i]);
      }
      //      acd_timestep(rdfwd, true,  0, specs->fdpars);
      acd_timestep(rdfwd, true,  0, w->model.specs);
      try {
	acd_timestep(rdadj, false, 0, w->model.specs);
      }
      catch (RVLException & e) {
	std::string ref = "Error: acd_timestep().  iw.size() = 1 has no adjoint!\n";
	std::stringstream err;
	e.write(err);
	EXPECT_EQ(ref,err.str());
      }

      delete fwd;
      delete adj;
      ps_delete(&par);
      fclose(stream);
    }
    catch (RVLException & e) {
	  
      e.write(cerr);
      exit(1);
    }
  }    

  TEST_F(ACDStepTest, acd_tsf_deriv1_adjtest_iwave) {
    try {

      // fake command line environment
      int argc = 2;
      char * argvv = new char[128];
      char ** argv = new char*[2];
      argv[0]=&(argvv[0]); argv[1]=&(argvv[65]);
      strcpy(argv[1],"par=parfile");
      PARARRAY * par = NULL;
      FILE * stream = NULL;
      IWaveEnvironment(argc, argv, 0, &par, &stream);
      delete [] argv;
      delete [] argvv;
      srand(time(NULL));
      
      // build order one IWaveTree, check
      int order=1;
      IWaveTree * fwd = new IWaveTree(*par,stream,ic,order);
      IWaveTree * adj = new IWaveTree(*par,stream,ic,order);
      std::vector<RDOM *> rdfwd, rdadj;

      IWAVE * w = (fwd->getStateArray()[0]);

      // FD_MODEL *specs = (FD_MODEL *)(w->model.specs);
        
      rdfwd = fwd->getRDOMArray();
      rdadj = adj->getRDOMArray();
      for (int i=0; i<rdfwd.size(); i++) {
        //rdom_rand(rdfwd[i]);
        //rdom_rand(rdadj[i]);
          iwave_rdom_rand(fwd->getStateArray()[i]);
          iwave_rdom_rand(adj->getStateArray()[i]);
      }
      // initialize rdadj[1], D_CSQ to zero for adjoint operator
      //rdom_copy(rdfwd[0],rdadj[0]);
      iwave_rdom_copy(fwd->getStateArray()[0],adj->getStateArray()[0]);
      ra_a_zero(&(rdadj[1]->_s[D_CSQ]));
        
        IWaveTree * fwdcp = new IWaveTree(*par,stream,ic,order);
        IWaveTree * adjcp = new IWaveTree(*par,stream,ic,order);

      for (int i=0; i<rdfwd.size();i++) {
        iwave_rdom_copy(fwd->getStateArray()[i],fwdcp->getStateArray()[i]);
        iwave_rdom_copy(adj->getStateArray()[i],adjcp->getStateArray()[i]);
      }

        std::vector<RDOM *> rdfwdcp, rdadjcp;

        
        rdfwdcp = fwdcp->getRDOMArray();
        rdadjcp = adjcp->getRDOMArray();

      ra_a_swap(&(rdadj[1]->_s[D_UC]),&(rdadj[1]->_s[D_UP]));

      acd_timestep(rdfwd, true,  0, w->model.specs);
      acd_timestep(rdadj, false, 0, w->model.specs);
        
      ireal yn = acd_rdom_norm(rdadjcp);
      ireal axn = acd_rdom_norm(rdfwd);
      ireal axy = acd_rdom_inner(rdfwd, rdadjcp, true);
      ireal xaty = acd_rdom_inner(rdadj, rdfwdcp, false);
      ireal adjrl = abs(axy - xaty)/axn/yn;
        
#ifdef GTEST_VERBOSE
      cout << "        ||y|| = " << yn << endl;
      cout << "       ||Ax|| = " << axn << endl;
      cout <<  "< Ax,    y > = " << axy<< endl;
      cout << "<  x, A^Ty > = " << xaty<< endl;
      cout << "< Ax, y> - < x, A^Ty> \n";
      cout << "--------------------- = " << abs(axy - xaty)/axn/yn << endl;
      cout << "       |Ax||y|        \n";
#endif
        
      ireal eps = 200 * numeric_limits<float>::epsilon();
      EXPECT_LE(adjrl, eps);
        
        delete fwdcp;
        delete adjcp;
      delete fwd;
      delete adj;
      ps_delete(&par);
      fclose(stream);
    }
    catch (RVLException & e) {
      e.write(cerr);
      exit(1);
    }
  }
    
    TEST_F(ACDStepTest, acd_tsf_deriv1_adjtest) {
        try {
            
            // fake command line environment
            int argc = 2;
            char * argvv = new char[128];
            char ** argv = new char*[2];
            argv[0]=&(argvv[0]); argv[1]=&(argvv[65]);
            strcpy(argv[1],"par=parfile");
            PARARRAY * par = NULL;
            FILE * stream = NULL;
            IWaveEnvironment(argc, argv, 0, &par, &stream);
            delete [] argv;
            delete [] argvv;
            srand(time(NULL));
            
            // build order one IWaveTree, check
            int order=1;
            int err;
            IWaveTree * fwd = new IWaveTree(*par,stream,ic,order);
            IWaveTree * adj = new IWaveTree(*par,stream,ic,order);
            std::vector<RDOM *> rdfwd, rdadj;

      cerr<<"1\n";
        
            
            IWAVE * w = (fwd->getStateArray()[0]);
            
            // FD_MODEL *specs = (FD_MODEL *)(w->model.specs);
            
            rdfwd = fwd->getRDOMArray();
            rdadj = adj->getRDOMArray();
	    cerr<<"2\n";
        

            for (int i=0; i<rdfwd.size(); i++) {
                //rdom_rand(rdfwd[i]);
                //rdom_rand(rdadj[i]);
                iwave_rdom_rand(fwd->getStateArray()[i]);
                iwave_rdom_rand(adj->getStateArray()[i]);
            }
      cerr<<"3\n";
        

            // initialize rdadj[1], D_CSQ to zero for adjoint operator
            //rdom_copy(rdfwd[0],rdadj[0]);
            iwave_rdom_copy(fwd->getStateArray()[0],adj->getStateArray()[0]);
            ra_a_zero(&(rdadj[1]->_s[D_CSQ]));
      cerr<<"4\n";
        

      //            int ndim = (w->model).g.dim;
            
            IPNT dgs[RDOM_MAX_NARR], dge[RDOM_MAX_NARR];    /*< computational domain */
            IPNT dgsa[RDOM_MAX_NARR], dgea[RDOM_MAX_NARR];  /*< allocated domain */
            for (int i=0; i<RDOM_MAX_NARR; i++) {
	      IASN(dgs[i],IPNT_1); IASN(dgsa[i],IPNT_1);
	      IASN(dge[i],IPNT_0); IASN(dgea[i],IPNT_0);
	      rd_a_gse(&((w->model).ld_c),i,dgs[i],dge[i]);
              rd_gse(&((w->model).ld_c),i,dgsa[i],dgea[i]);
            }

            std::vector<RDOM > fwdcp(rdfwd.size());
            std::vector<RDOM > adjcp(rdadj.size());
            for (int i=0; i<rdfwd.size();i++) {
                /*-declare computational domain---------------------------------------------*/

                err = rd_a_create(&(fwdcp[i]), RDOM_MAX_NARR, dgs, dge);
                err = err || rd_a_create(&(adjcp[i]), RDOM_MAX_NARR, dgs, dge);
                
                if ( err ){
                    RVLException e;
                    e<<"Error: ACDStepTest - acdfwd\n";
                    e<<"  error from rd_a_create = "<<err<<"\n";
                    throw e;
                }
                
                rdom_copy(rdfwd[i],&fwdcp[i]);
                rdom_copy(rdadj[i],&adjcp[i]);
            }
            ra_a_swap(&(rdadj[1]->_s[D_UC]),&(rdadj[1]->_s[D_UP]));
            
            acd_timestep(rdfwd, true,  0, w->model.specs);
            acd_timestep(rdadj, false, 0, w->model.specs);
            
            
            ireal yn = acd_rdom_norm(adjcp);
            ireal axn = acd_rdom_norm(rdfwd);
            ireal axy = acd_rdom_inner(rdfwd, adjcp, true);
            ireal xaty = acd_rdom_inner(rdadj, fwdcp, false);
            ireal adjrl = abs(axy - xaty)/axn/yn;
            
#ifdef GTEST_VERBOSE
            cout << "        ||y|| = " << yn << endl;
            cout << "       ||Ax|| = " << axn << endl;
            cout <<  "< Ax,    y > = " << axy<< endl;
            cout << "<  x, A^Ty > = " << xaty<< endl;
            cout << "< Ax, y> - < x, A^Ty> \n";
            cout << "--------------------- = " << abs(axy - xaty)/axn/yn << endl;
            cout << "       |Ax||y|        \n";
#endif
            
            ireal eps = 200 * numeric_limits<float>::epsilon();
            EXPECT_LE(adjrl, eps);
            for (int i=0; i<rdfwd.size();i++) {
                rd_a_destroy(&(fwdcp[i]));
                rd_a_destroy(&(adjcp[i]));
            }
            
            delete fwd;
            delete adj;
            ps_delete(&par);
            fclose(stream);
        }
        catch (RVLException & e) {
            e.write(cerr);
            exit(1);
        }
    }    

  /*    unknown error - disabled 15.04.15
  TEST_F(ACDStepTest, acd_tsf_deriv2_adjtest) {
    try {
            
      // fake command line environment
      int argc = 2;
      char * argvv = new char[128];
      char ** argv = new char*[2];
      argv[0]=&(argvv[0]); argv[1]=&(argvv[65]);
      strcpy(argv[1],"par=parfile");
      PARARRAY * par = NULL;
      FILE * stream = NULL;
      IWaveEnvironment(argc, argv, 0, &par, &stream);
      delete [] argv;
      delete [] argvv;
        
      srand(time(NULL));
      // build order one IWaveTree, check
      int order=2;
      int err;
      IWaveTree * fwd = new IWaveTree(*par,stream,ic,order);
      IWaveTree * adj = new IWaveTree(*par,stream,ic,order);
      std::vector<RDOM *> rdfwd, rdadj;
            
      IWAVE * w = (fwd->getStateArray()[0]);
            
      // FD_MODEL *specs = (FD_MODEL *)(w->model.specs);
            
      rdfwd = fwd->getRDOMArray();
      rdadj = adj->getRDOMArray();
            
      for (int i=0; i<rdfwd.size(); i++) {
	    //rdom_rand(rdfwd[i]);
	    //rdom_rand(rdadj[i]);
          iwave_rdom_rand(fwd->getStateArray()[i]);
          iwave_rdom_rand(adj->getStateArray()[i]);
      }

        iwave_rdom_copy(fwd->getStateArray()[0],adj->getStateArray()[0]);
        iwave_rdom_copy(fwd->getStateArray()[1],adj->getStateArray()[1]);

      // initialize rdadj[1], D_CSQ to zero for adjoint operator
      ra_a_zero(&(rdadj[2]->_s[D_CSQ]));
      ra_a_zero(&(rdadj[2]->_s[D_UC]));
      ra_a_zero(&(rdadj[3]->_s[D_CSQ]));
            
      //      int ndim = (w->model).g.dim;
            
      IPNT dgs[RDOM_MAX_NARR], dge[RDOM_MAX_NARR];    // computational domain 
      IPNT dgsa[RDOM_MAX_NARR], dgea[RDOM_MAX_NARR];  // allocated domain 
      for (int i=0; i<RDOM_MAX_NARR; i++) {
	rd_a_gse(&((w->model).ld_a),i,dgs[i],dge[i]);
	rd_gse(&((w->model).ld_a),i,dgsa[i],dgea[i]);
      }
      std::vector<RDOM > fwdcp(rdfwd.size());
      std::vector<RDOM > adjcp(rdadj.size());
      for (int i=0; i<rdfwd.size();i++) {
      err = rd_a_create(&(fwdcp[i]), RDOM_MAX_NARR, dgs, dge);
      err = err || rd_a_create(&(adjcp[i]), RDOM_MAX_NARR, dgs, dge);
                
	if ( err ){
	  RVLException e;
	  e<<"Error: ACDStepTest - acdfwd\n";
	  e<<"  error from rd_a_create = "<<err<<"\n";
	  throw e;
	}
	rdom_copy(rdfwd[i],&fwdcp[i]);
	rdom_copy(rdadj[i],&adjcp[i]);
      }
            
      ra_a_swap(&(rdadj[2]->_s[D_UC]),&(rdadj[2]->_s[D_UP]));
      ra_a_swap(&(rdadj[3]->_s[D_UC]),&(rdadj[3]->_s[D_UP]));

      acd_timestep(rdfwd, true,  0, w->model.specs);
      acd_timestep(rdadj, false, 0, w->model.specs);
            
       // ireal yn = acd_rdom_norm(adjcp);

      ireal yn = sqrt(acd_rdom_inner_csq(adjcp[2],0)+acd_rdom_inner_csq(adjcp[3],0));
      ireal axn = sqrt(acd_rdom_inner_csq(rdfwd[2],0)+acd_rdom_inner_csq(rdfwd[3],0));
      ireal axy = acd_rdom_inner(rdfwd, adjcp, true);
      ireal xaty = acd_rdom_inner(rdadj, fwdcp, false);
      ireal adjrl = abs(axy - xaty)/axn/yn;

#ifdef GTEST_VERBOSE
      cout << "        ||y|| = " << yn << endl;
      cout << "       ||Ax|| = " << axn << endl;
      cout <<  "< Ax,    y > = " << axy<< endl;
      cout << "<  x, A^Ty > = " << xaty<< endl;
      cout << "< Ax, y> - < x, A^Ty> \n";
      cout << "--------------------- = " << abs(axy - xaty)/axn/yn << endl;
      cout << "       |Ax||y|        \n";
#endif
            
      ireal eps = 200 * numeric_limits<float>::epsilon();
      EXPECT_LE(adjrl, eps);
      for (int i=0; i<rdfwd.size();i++) {
          rd_a_destroy(&(fwdcp[i]));
          rd_a_destroy(&(adjcp[i]));
      }
            
      delete fwd;
      delete adj;
      ps_delete(&par);
      fclose(stream);
    }
    catch (RVLException & e) {
      e.write(cerr);
      exit(1);
    }
  }
	*/
}

#include "main.h"