#include "iwop.hh"
#include "adjtest.hh"

#include "MPS_to_RHS.hh"
#include "MPS_iwop.hh"
#include "MPS_Space_Examples.hh"

#define CSQ "csq"
#define BULKMOD "bulkmod"
#define BUOYANCY "buoyancy"
#define SOURCE_P "source_p"
#define SOURCE_V0 "source_v0"
#define SOURCE_V1 "source_v1"
#define DATA_P "data_p"
#define DATA_V0 "data_v0"
#define DATA_V1 "data_v1"

#ifndef MY_TOL
#define MY_TOL 1e-3
#endif

#define VERBOSE_MJB

using RVL::valparse;
using RVL::AssignFilename;
using RVL::OperatorEvaluation;

using TSOpt::IWaveLOVOp;
using TSOpt::IWaveSpace;
using TSOpt::CanScal_MPS_Space;
using TSOpt::CanVec_MPS_Space;
using TSOpt::ExVec_MPS_Space;
using TSOpt::Scal_MPS_Space;
using TSOpt::MPS_to_RHS;
using TSOpt::MPS_IWaveLOVOp;


template< class T_MPS_Space >
void test_kernel(PARARRAY &pars, FILE &stream, string test_info){

  std::stringstream str;
  string message;

  string jobname = valparse<string>(pars,"jobname");

  try {
    str << "\n============================================================\n"
	<< "Test name: "<< jobname <<"\n"
	<< "Runnning test for MPS_iwop.\n"
	<< test_info <<"\n";
#ifdef IWAVE_USE_MPI
    str << "rank = "<< retrieveGlobalRank() <<"\n";
#endif
    str << "============================================================\n";  

    string driver_type = valparse<string>(pars,"driver");
    string src_type = valparse<string>(pars,"src_type");
    string data_type = valparse<string>(pars,"data_type");

    if( driver_type.compare("acd")!=0 &&
        driver_type.compare("asg")!=0 ){
      RVLException e;
      e << "Driver ("<< driver_type <<") is of unexpected type!\n";
      throw e;
    }
    if( src_type.compare("pressure")!=0 &&
	src_type.compare("velocity")!=0 ){
      RVLException e;
      e << "Source ("<< src_type <<") is of unexpected type!\n";
      throw e;
    }
    if( data_type.compare("pressure")!=0 &&
	data_type.compare("velocity")!=0 ){
      RVLException e;
      e << "Data ("<< data_type <<") is of unexpected type!\n";
      throw e;
    }


    TSOpt::MPS_KEYS mks;
    mks.appx_ord   = "appx_ord";
    mks.MPS_file   = "MPS_file";
    mks.delta_file = "MPS_delta";

    //setting mks.grid_file
    if( driver_type.compare("acd")==0 ){
      mks.grid_file = CSQ;
    }
    else{
      mks.grid_file = BULKMOD;
    }
    //setting mks.RHS_files
    if( src_type.compare("pressure")==0 ){
      mks.RHS_files.push_back(SOURCE_P);
    }
    else{
      mks.RHS_files.push_back(SOURCE_V0);
      mks.RHS_files.push_back(SOURCE_V1);
    }
    //setting mks.hdr_files and mks.G_files
    if( data_type.compare("pressure")==0 ){
      mks.hdr_files.push_back(DATA_P);
      mks.G_files.push_back("data_p_g");
    }
    else{
      mks.hdr_files.push_back(DATA_V0);
      mks.hdr_files.push_back(DATA_V1);
      mks.G_files.push_back("data_v0_g");
      mks.G_files.push_back("data_v1_g");
    }


    /////////////////
    // making data //
    /////////////////
#ifdef VERBOSE_MJB
    str << "\n[] Making data:\n";
#endif

    //initializing
    T_MPS_Space mps_sp(mks,pars);
    MPS_to_RHS m2r(mks,pars,mps_sp);
    IWaveLOVOp iwop(m2r.get_pars(),&stream);
    RVL::LinCompLOVOp<float> IWOP(m2r,iwop);

    //setting inputs
    Vector<ireal> x(IWOP.getDomain());
    Components<ireal> cx(x);
    Components<ireal> cx0(cx[0]);
    Components<ireal> cx1(cx[1]);

    if( driver_type.compare("acd")==0 ){
      AssignFilename af_csq(valparse<string>(pars,CSQ));
      cx0[0].eval(af_csq);
    }
    else{
      AssignFilename af_bmod(valparse<string>(pars,BULKMOD));
      AssignFilename af_buoy(valparse<string>(pars,BUOYANCY));
      cx0[0].eval(af_bmod);
      cx0[1].eval(af_buoy);
    }
    AssignFilename af_src(valparse<string>(pars,mks.MPS_file));
    cx1[0].eval(af_src);

    //setting outputs
    Vector<ireal> y(IWOP.getRange());
    Components<ireal> cy(y);

    if(data_type.compare("pressure")==0){
      AssignFilename af_data_p(valparse<string>(pars,DATA_P));
      cy[0].eval(af_data_p);
    }
    else{
      AssignFilename af_data_v0(valparse<string>(pars,DATA_V0));
      AssignFilename af_data_v1(valparse<string>(pars,DATA_V1));
      cy[0].eval(af_data_v0);
      cy[1].eval(af_data_v1);
    }

    //evaluating operator
    OperatorEvaluation<ireal> opeval(IWOP,x);
    y.copy(opeval.getValue());

#ifdef IWAVE_USE_MPI
    MPI_Barrier(retrieveGlobalComm());
#endif

    /////////////////////////////////
    // making data via convolution //
    /////////////////////////////////
#ifdef VERBOSE_MJB
    str << "\n[] Making data via convolution:\n\n";
#endif

    //initialing
    MPS_IWaveLOVOp<T_MPS_Space> mps_lovop(mks,pars,&stream);
    //str << "writing out MPS_IWaveLOVOp:\n";
    //mps_lovop.write(str);

    //output
    Vector<ireal> y_re(mps_lovop.getRange());
    Components<ireal> cy_re(y_re);

    if( data_type.compare("pressure")==0 ){
      AssignFilename af_data_p_re(valparse<string>(pars,"data_p_re"));
      cy_re[0].eval(af_data_p_re);
    }
    else{
      AssignFilename af_data_v0_re(valparse<string>(pars,"data_v0_re"));
      AssignFilename af_data_v1_re(valparse<string>(pars,"data_v1_re"));
      cy_re[0].eval(af_data_v0_re);
      cy_re[1].eval(af_data_v1_re);
    }
 
    //eval
    OperatorEvaluation<ireal> opeval_re(mps_lovop,x);
    y_re.copy(opeval_re.getValue());

#ifdef IWAVE_USE_MPI
    MPI_Barrier(retrieveGlobalComm());
#endif
    
    //residual
#ifdef VERBOSE_MJB
    str << "\n[] Residual\n\n";
#endif

    Vector<ireal> fwd_res(y_re);
    fwd_res.copy(y_re);
    fwd_res.linComb(-1,y);
    float norm_fwd_res = fwd_res.norm() / y.norm();

    bool fwd_test = (fabs(norm_fwd_res)<MY_TOL);
    if( fwd_test ) message="\n\033[1;32m[PASSED]\033[0m";
    else message="\n\033[1;31m[FAILED]\033[0m";

    message += " Fwd MPS_iwop test\n";
    str  << message;
    if(!fwd_test){
      str << "|data - G*f|/|data| = " << norm_fwd_res << "\n"
	  << "|data| = " << y.norm()    << "\n"
	  << "|G*f| = "  << y_re.norm() << "\n";
    }

#ifdef VERBOSE_MJB
      str << "|data - G*f|/|data| = " << norm_fwd_res << "\n"
	  << "|data| = " << y.norm()    << "\n"
	  << "|G*f| = "  << y_re.norm() << "\n";
#endif

#ifdef IWAVE_USE_MPI
    MPI_Barrier(retrieveGlobalComm());
#endif

    //////////////////
    // Adjoint test //
    //////////////////
#ifdef VERBOSE_MJB
    str << "\n[] Running adjoint test on MPS_IWaveLOVOp:\n\n";
#endif

    RVL::LinearRestrictOp<float> lop(mps_lovop,cx[0]);
    
    std::stringstream res;
    RVL::RVLRandomize<float> rnd(getpid(),-1.0,1.0);
    bool adj_test=RVL::AdjointTest<float>(lop,rnd,res);
    
    if(adj_test) message="\n\033[1;32m[PASSED]\033[0m";
    else message="\n\033[1;31m[FAILED]\033[0m";
    
    message += " Adjoint test for MPS_conv op\n";
    str  << message;
    if(!adj_test) str << res.str();

#ifdef VERBOSE_MJB
    str << res.str();
#endif


#ifdef IWAVE_USE_MPI
    MPI_Barrier(retrieveGlobalComm());
#endif

    // writing out
    string outfile = valparse<string>(pars,"outfile","");

#ifdef IWAVE_USE_MPI
    for( int r=0; r<retrieveGlobalSize(); r++ ){
      if( retrieveGlobalRank()==r ){
	cerr << "rank="<<r<<"\n";
#endif
	if(outfile.size()>0){
	  ofstream outf;
	  outf.open(outfile.c_str(),ios::app);
	  outf << str.str();
	  outf.close();
	}
	else{
	  cout << str.str();
	}
#ifdef IWAVE_USE_MPI
      }
      MPI_Barrier(retrieveGlobalComm());
    }
#endif
  }
  catch(RVLException &e){
    e << "ERROR from test_kernel()!\n"
      << "Dumping output stream:\n"
      << str.str();
    throw e;
  }
}