// MPS_howto.cc
// Author: Mario J. Bencomo
// last modified: 12/12/16

#include "asg_defn.hh"
#include "asg.hh"
#include "MPS_howto.hh"
#include "MPS_to_RHS.hh"

using TSOpt::add_to_pars;
using TSOpt::ExScal_MPS_Space;
using TSOpt::MPS_to_RHS;
using TSOpt::MPS_KEYS;

const char *sdoc[] = {
  " ============================================================================",
  " MPS_howto.x ",
  " ============================================================================",
  " Authors: Mario J. Bencomo",
  " ",
  " 2D acoustic variable density centered difference modeling (asg) with an example",
  " scalar (pressure) multipole source."
  "",
  " PDE:",
  "   dp/dt + beta div(v) = f_1(t) delta(x-x^*) + f_2(t) {d/dx_0 + d/dx_1} delta(x-x^*)",
  "   dv/dt + kappa grad(p) = 0",
  "",
  " Input:",
  "   kappa = bulk modulus (units GPa)",
  "   beta = buoyancy (units cm^3/g)",
  "   {f_1,f_2} = multipole coeffiecients (units GPa*m^2/s for f_1, units GPa*m^3/s for f_2)",   
  "",
  " Output:",
  "   p = pressure data (units GPa)",
  " ",
  " Typical parameter list. May be copied, edited, and used for input: either",
  " include parameters on command line (for example in Flow), or place",
  " in file <foo> and include \"par=<foo>\" on command line. Any parameter",
  " included explicitly in command line overrides parameter with same key",
  " in par file.",
  " ",
  "  Invoke single threaded execution by ",
  " \"MPS_howto.x [parameters] [standalone install]\"",
  " ",
  " or multi-threaded execution using interactive or batch MPI (for which",
  " code must be compiled with MPI enabled).",
  " ",
  " Given values are defaults; non-optional values indicated by corner brackets.",
  " ",
  " --------------------------- begin parameters ---------------------------",
  " FD info:",
  " ",
  "          order = 2           spatial half-order",
  "            cfl = 0.75        proportion of max dt/dx",
  "           cmin = 1.0         min permitted velocity (m/ms) - sanity check",
  "           cmax = 4.5         max permitted velocity (m/ms) - used in dt comp",
  "           dmin = 1.0         min permitted density (g/cm^3) - sanity check",
  "           dmax = 1.0         max permitted density (g/cm^3) - sanity check",
  " ",
  " ------------------------------------------------------------------------",
  " Source info:",
  " ",
  "      source_p = <path>    RHS source file, SU format - need filename.",
  "                           Currently driver assumes file does not exist,"
  "                           thus will generate such file as part of the"
  "                           initialization of MPS_to_RHS.",
  "      appx_ord = <2*order> Singular source approximation order. ",
  "      MPS_file = <path>    MPS file containing MPS coefficients, SU format - must exist.",
  "                           For this particular MPS space the MPS file must contain"
  "                           two traces, and source locations.",
  " ",
  " ------------------------------------------------------------------------",
  " Trace info:",
  " ",
  "          data_p = <path>    data file, SU format - must exist ",
  "                             headers establish acquisition geometry",
  "                             output for forward modeling",
  " ",
  " ------------------------------------------------------------------------",
  " Model info:",
  " ",
  "        bulkmod = <path>      input bulk modulus, ",
  "                              determines simulation spatial grid,",
  "       buoyancy = <path>      input buoyancy",
  " ",
  " ------------------------------------------------------------------------",
  " MPI info:",
  " ",
  "        mpi_np1 = 1           number of subdomains along axis 1",
  "        mpi_np2 = 1           number of subdomains along axis 2",
  "        mpi_np3 = 1           number of subdomains along axis 3",
  "        partask = 1           number of shots to execute in parallel",
  " ",
  " ------------------------------------------------------------------------",
  " Output info:",
  " ",
  " FD ouput - written to coutxxx.txt on rank xxx",
  "       printact = 0           monitor output:",
  "                            < 0 - none",
  "                              0 - announce each shot simulation",
  "                              1 - print time step index",
  "                              2 - diagnostic messages from main routines",
  "                              > 2 - much more, useful only for serious ",
  "                                 debugging",
  "        dump_pi = 0           dump parallel/dom. decomp info",
  "       dump_lda = 0           dump grid data for allocated arrays",
  "       dump_ldc = 0           dump grid data for computational arrays",
  "       dump_ldr = 0           dump grid data for receive buffers",
  "       dump_lds = 0           dump grid data for send buffers",
  "      dump_term = 0           dump trace header data",
  " ---------------------------end parameters ------------------------------",
  NULL };


IOKEY IWaveInfo::iwave_iokeys[]
= {
  {"bulkmod",  0, true,  1 },
  {"buoyancy", 1, true,  1 },
  {"source_p", 2, true,  2 },
  {"data_p",   2, false, 2 },
  {"",         0, false, 0 }
};


int xargc;
char **xargv;

int main(int argc, char ** argv) {

  std::stringstream str;

  try {

#ifdef IWAVE_USE_MPI
    int ts=0;
    MPI_Init_thread(&argc,&argv,MPI_THREAD_FUNNELED,&ts);    
#endif

    PARARRAY * pars = NULL;
    FILE * stream = NULL;
    TSOpt::IWaveEnvironment(argc, argv, 0, &pars, &stream);

#ifdef IWAVE_USE_MPI
    if( retrieveGlobalRank()==0 ){
#endif

      if(argc<2){
	pagedoc();
	exit(0);
      }

      cerr << "\n============================================================\n" 
	   << "Running MPS_howto.x\n"
	   << "============================================================\n" ;
      //ps_printall(*pars,stderr);
#ifdef IWAVE_USE_MPI
    }
#endif

    MPS_KEYS mks;
    mks.appx_ord  = "appx_ord"; //required for MPS_to_RHS ops
    mks.MPS_file  = "MPS_file"; //required for MPS_Space and MPS_to_RHS ops
    mks.grid_file = "bulkmod";  //required for MPS_Space and MPS_to_RHS ops
    mks.RHS_files.push_back("source_p"); //required for MPS_to_RHS ops

    //------------------------//
    // Initializing mps space //
    //------------------------//

    // Assuming MPS file exists
    ExScal_MPS_Space mps_sp(mks,*pars);

    str << "------------------------------\n"
	<< "Printing out info of MPS space\n"
	<< "------------------------------\n";
    mps_sp.write(str);


    //----------------------------//
    // Initializing MPS_to_RHS op //
    //----------------------------//

    // Assuming RHS files do not exist 
    // Assuming MPS file exist    
    MPS_to_RHS m2r(mks,*pars, mps_sp);

    str << "----------------------------------\n"
	<< "Printing out info of MPS_to_RHS op\n"
	<< "----------------------------------\n";
    m2r.write(str);


    //------------------------//
    // Applying MPS_to_RHS op //
    //------------------------//
    
    //setting up mps vector
    Vector<ireal> mps(mps_sp);
    RVL::AssignFilename af_mps(RVL::valparse<string>(*pars,mks.MPS_file));
    mps.eval(af_mps);

    //setting up rhs vector
    Vector<ireal> rhs(m2r.getRange());
    RVL::AssignFilename af_rhs(RVL::valparse<string>(*pars,mks.RHS_files[0]));
    rhs.eval(af_rhs);
    
    //applying op
    m2r.applyOp(mps,rhs);


    //----------------------------------//
    // Applying MPS_to_RHS + IWaveLOVOp //
    //----------------------------------//
 
    //operator composition
    TSOpt::IWaveLOVOp iwop(m2r.get_pars(),stream);
    RVL::LinCompLOVOp<float> IWOP(m2r,iwop);

    cerr << "Writing out IWOP\n";
    IWOP.write(cerr);

    //setting up domain vector
    Vector<ireal> x(IWOP.getDomain());
    Components<ireal> cx(x);
    Components<ireal> cx_m(cx[0]); //model params
    Components<ireal> cx_f(cx[1]); //source params
    RVL::AssignFilename af_bmod(RVL::valparse<string>(*pars,"bulkmod"));
    RVL::AssignFilename af_buoy(RVL::valparse<string>(*pars,"buoyancy"));
    cx_m[0].eval(af_bmod);
    cx_m[1].eval(af_buoy);
    cx_f[0].eval(af_mps);

    //setting up range vector
    Vector<ireal> y(IWOP.getRange());
    RVL::AssignFilename af_data(RVL::valparse<string>(*pars,"data_p"));
    y.eval(af_data);
    
    //applying op
    RVL::OperatorEvaluation<ireal> eval_IWOP(IWOP,x);
    y.copy(eval_IWOP.getValue());
    
    
    //-------------//
    // Writing out //
    //-------------//
    string outfile = RVL::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{
	  cerr << str.str();
	}
#ifdef IWAVE_USE_MPI
      }
      MPI_Barrier(retrieveGlobalComm());
    }
#endif

    //clean up
    ps_delete(&pars);
    iwave_fdestroy();

#ifdef IWAVE_USE_MPI    
    if( retrieveGlobalRank()==0 ){
#endif            
      cerr << "============================================================\n" 
	   << "Exiting MSP_howto.x\n"
	   << "============================================================\n\n";
#ifdef IWAVE_USE_MPI    
    }
    MPI_Finalize();
#endif

  }
  catch(RVLException &e){
    e << "Exiting with error!\n"
      << "Dumping output stream:\n"
      << str.str();

#ifdef IWAVE_USE_MPI
    if( retrieveGlobalRank()==0 ){
#endif
      e.write(cerr);
#ifdef IWAVE_USE_MPI
    }   
    MPI_Abort(MPI_COMM_WORLD,0);
#endif
    exit(1);
  }
}