/*************************************************************************

Copyright Rice University, 2011.
All rights reserved.

Permission is hereby granted, free of charge, to any person obtaining a
copy of this software and associated documentation files (the "Software"),
to deal in the Software without restriction, including without limitation
the rights to use, copy, modify, merge, publish, distribute, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, provided that the above copyright notice(s) and this
permission notice appear in all copies of the Software and that both the
above copyright notice(s) and this permission notice appear in supporting
documentation.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT OF THIRD PARTY
RIGHTS. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR HOLDERS INCLUDED IN THIS
NOTICE BE LIABLE FOR ANY CLAIM, OR ANY SPECIAL INDIRECT OR CONSEQUENTIAL
DAMAGES, OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR
PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS
ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
THIS SOFTWARE.

Except as contained in this notice, the name of a copyright holder shall
not be used in advertising or otherwise to promote the sale, use or other
dealings in this Software without prior written authorization of the
copyright holder.

**************************************************************************/

#include "gtest/gtest.h"
#include "rnmat.hh"
#include "VPM.hh"
#include "usempi.h"

namespace {
  
  using namespace RVL;
  using namespace RVLAlg;
  using namespace RVLUmin;

  /** Unit Test
      Purpose: testing VPM functional
  */

  /** LinOpValOp of the type f(x0,x1) = L(x0)x1, where x0 and x1 are the
      nonlinear and linear parameters resp., and L(x0) = (x0'x0)A for some
      matrix A. It will be assumed that x0 and x1 live in the same space.
  */  
 
  template<typename T>
  class Test_LinOpValOp: public LinOpValOp<T> {
    
  private:
    CartesianPowerSpace<T> dom;
    GenMat<T> A;
    
    OperatorProductDomain<T> * clonePD() const{return new Test_LinOpValOp(*this);}
    
  protected:
    //f(x0,x1) = L(x0)x1 = (x0'x0) A x1 = y
    void apply0(const Vector<T> &x0,
		const Vector<T> &x1,
		Vector<T> &y ) const {
      try{
	A.applyOp(x1,y);
	T norm2x0 = x0.normsq();
	y.scale(norm2x0);
      }
      catch(RVLException & e){
	e << "ERROR from Test_LinOpValOp::apply0\n";
	throw e;
      }
    }
    
    //L(x0)*y = (x0'x0) A* y = x1
    void applyAdj0(const Vector<T> &x0,
		   const Vector<T> &y,
		   Vector<T> & x1) const {
      try{
	A.applyAdjOp(y,x1);
	T norm2x0 = x0.normsq();
	x1.scale(norm2x0);
      }
      catch(RVLException & e){
	e << "ERROR from Test_LinOpValOp::applyAdj0\n";
	throw e;
      }
    }
    
    //df_x0(x0,x1) dx0 = 2 A x1 x0'dx0 = dy
    void applyPartialDeriv0(const Vector<T> &x0,
			    const Vector<T> &x1,
			    const Vector<T> &dx0,
			    Vector<T> dy) const {
      try{
	T tmp = dx0.inner(x0);
	tmp *= 2;
	A.applyOp(x1,dy);
	dy.scale(tmp);
      }
      catch(RVLException & e){
	e << "ERROR from Test_LinOpValOp::applyPartialDeriv0\n";
	throw e;
      }
    }
    
    //df_x0(x0,x1)* dy = 2 x0 x1' A* dy = dx0
    void applyAdjPartialDeriv0(const Vector<T> & x0,
			       const Vector<T> & x1,
			       const Vector<T> & dy,
			       Vector<T> dx0) const {
      try{
	Vector<T> tmp_x(x1);
	A.applyAdjOp(dy,tmp_x);
	T tmp_c = x1.inner(tmp_x);
	tmp_c *= 2;
	dx0.linComb(tmp_c,x0);
      }
      catch(RVLException & e){
	e << "ERROR from Test_LinOpValOp::applyAdjPartialDeriv0\n";
	throw e;
      }
    }
    
    void applyPartialDeriv20(const Vector<T> & x0,
			     const Vector<T> & x1,
			     const Vector<T> & dx00,
			     Vector<T> dy) const {
      RVLException e;
      e<<"ERROR: LinOpValOp::applyPartialDeriv20\n";
      e<<"  no implementation supplied\n";
      throw e;
    }
    
    void applyAdjPartialDeriv20(const Vector<T> & x0,
				const Vector<T> & x1,
				const Vector<T> & dy,
				Vector<T> dx00) const {
      RVLException e;
      e<<"ERROR: LinOpValOp::applyPartialDeriv20\n";
      e<<"  no implementation supplied\n";
      throw e;
    }
    
  public:
    Test_LinOpValOp(GenMat<T> const & _A) 
      : A(_A), dom(2,_A.getDomain()){}
    Test_LinOpValOp(Test_LinOpValOp<T> const & lovo)
      : A(lovo.A), dom(lovo.dom){}
    ~Test_LinOpValOp(){}
    
    ProductSpace<T> const & getProductDomain() const {
      return dom;
    }
    Space<T> const & getRange() const{
      return A.getRange();
    }
    
    ostream & write(ostream & str) const{
      str<<"Test_LinOpValOp operator\n";
      return str;
    }
    
  };

  class VPMTest: public ::testing::Test {

  protected:
    
    int dim;
    float rtol;
    float nrtol;
    float Delta;
    int maxcount;
    string output; 
    
    //typedef typename ScalarFieldTraits<T>::AbsType atype;

    VPMTest() {
      
      PlantSeeds(19490615);
      dim = 3;
      rtol = 1.0e-5;
      nrtol = 1.0e-5;
      Delta = numeric_limits<float>::max();
      maxcount = 100;
      
    }

    virtual ~VPMTest() {
      // You can do clean-up work that doesn't throw exceptions here.
      //      unlink(fname.c_str());
    }

    // If the constructor and destructor are not enough for setting up
    // and cleaning up each test, you can define the following methods:

    virtual void SetUp() {
      // Code here will be called immediately after the constructor (right
      // before each test).
    }

    virtual void TearDown() {
      // Code here will be called immediately after each test (right
      // before the destructor).
    }

    // Objects declared here can be used by all tests in the test case for Grid.
  };

  TEST_F(VPMTest,SimpleLOVO) {

    try {
      
#ifdef VERBOSE
      cout<<" ****************************************************\n";
      cout<<" *            BEGIN UMin UNIT TEST "<<endl;
      cout<<" ****************************************************\n\n"; 
      cout<<"\n1. Construct random GenMat A\n";
#endif
      RnSpace<float> Rn(dim);
      GenMat<float> A(Rn,Rn);
      for (int i=0;i<dim;i++) 
	A.getElement(i,i)+=ScalarFieldTraits<float>::One();    

#ifdef VERBOSE
      cout<<"printing out matrix A:\n";
      for (int i=0; i<dim; i++){
	for (int j=0; j<dim; j++){
	  cout<<"A["<<i<<","<<j<<"]="<<A.getElement(i,j)<<"\n";
	}
      } 
      cout<<"\n";
#endif
    
      RVLRandomize<float> rnd;
      LocalVector<float> x0(Rn);  x0.eval(rnd);
      LocalVector<float> sol(Rn); sol.eval(rnd);
      LocalVector<float> d(Rn);   d.zero();
      float norm2x0 = x0.normsq();
      d.linComb(norm2x0,sol);

#ifdef VERBOSE    
      cout<<"printing out x0:\n";  x0.write(cout);  cout<<"\n";
      cout<<"printing out sol:\n"; sol.write(cout); cout<<"\n";
      cout<<"printing out d:\n";   d.write(cout);   cout<<"\n";
    
      cout<<"\n2. Build Test_LinOpValOp and VPM functional\n";
#endif

      Test_LinOpValOp<float> LOVO(A);
      CGNEPolicyData<float> cgnepd(rtol,
				   nrtol,
				   Delta,
				   maxcount);
    
      VPM<float,CGNEPolicy<float>,CGNEPolicyData<float> > vpmfun(LOVO,d,cgnepd,cout);
    
#ifdef VERBOSE    
      cout<<"\n3. Run VPM functional\n";
#endif
      FunctionalEvaluation<float> feval(vpmfun,x0);
      float val = feval.getValue();
      VPM<float,CGNEPolicy<float>,CGNEPolicyData<float> > const & fcp = 
	dynamic_cast< VPM<float,CGNEPolicy<float>,CGNEPolicyData<float> > const &>(feval.getFunctional());
      LocalVector<float> const & est = fcp.getLSSoln();
    
#ifdef VERBOSE
      cout<<"printing out est vector:\n";
      est.write(cout);
      cout<<"\n";
      cout<<"VPM functional value: "<< val <<"\n"; 
#endif

      EXPECT_NEAR(0.0,val,1.0e-9);
     
      //computing difference between sol and est
      float sol_norm = sol.norm();
      sol.linComb(-1.0,est);

      float res_norm = sol.norm();
      EXPECT_NEAR(0.0,res_norm,1.0e-5);
      //cout<<"residual norm, sol-est = "<< res_norm << "\n";

      res_norm /= sol_norm;
      EXPECT_NEAR(0.0,res_norm,1.0e-5);
      //cout<<"relative norm, sol-est / norm(sol) = " <<res_norm<<"\n";
    }
    catch (RVLException & e) {
      e.write(cerr);
      exit(1);
    }
  };
} // namespace

int main(int argc, char **argv) {
#ifdef IWAVE_USE_MPI
  int ts=0;
  MPI_Init_thread(&argc,&argv,MPI_THREAD_FUNNELED,&ts);    
#endif
  try {
    ::testing::InitGoogleTest(&argc, argv);
#ifdef IWAVE_USE_MPI
    //    storeGlobalComm(MPI_COMM_WORLD);
    //    storeComm(MPI_COMM_WORLD);
#endif
    int res = RUN_ALL_TESTS();
#ifdef IWAVE_USE_MPI
    MPI_Finalize();
#endif
    return res;
  }
  catch (RVLException & e) {
    e.write(cerr);
#ifdef IWAVE_USE_MPI
    MPI_Abort(MPI_COMM_WORLD,0);
    MPI_Finalize();
#endif
    exit(1);
  }
}