//$Id: tWaveletFunctions.cc,v 0.2 2001/06/06 15:17:00 sazonov Exp $

/*
 Run program:

 tWaveletFunctions ; echo exit status = $stat
 and see if 'exit status = 0'

 or for more detailed output

 tWaveletFunctions --verbose ; echo exit status = $stat
*/

#include <stream.h>
#include <stdlib.h>

#include "CallChain.hh"
#include "UDT.hh"
#include "SequenceUDT.hh"
#include "ScalarUDT.hh"
#include "TimeSeries.hh"
#include "Statistics.hh"

#include "WaveletFunctions.hh"
#include "unittest.h"

#define NN 256          // number of data samples
#define TEST_LEVEL 4
#define TEST_ORDER 10
#define TEST_TREE  0

// the following tests are performed:
/*
   test if lifting wavelet for <float> initialized properly
   test if lifting wavelet for <double> initialized properly
   test if Daubechies wavelet for <float> initialized properly
   test if Daubechies wavelet for <double> initialized properly
   test if Meyer wavelet for <float> initialized properly
   test if Meyer wavelet for <double> initialized properly

   test if forward lifting transform for <float> computed properly
   test if forward lifting transform for <double> computed properly
   test if forward Daubechies transform for <float> computed properly
   test if forward Daubechies transform for <double> computed properly
   test if forward Meyer transform for <float> computed properly
   test if forward Meyer transform for <double> computed properly

   test if inverse lifting transform for <float> computed properly
   test if inverse lifting transform for <double> computed properly
   test if inverse Daubechies transform for <float> computed properly
   test if inverse Daubechies transform for <double> computed properly
   test if inverse Meyer transform for <float> computed properly
   test if inverse Meyer transform for <double> computed properly
*/

using namespace datacondAPI;

void ilwd_dump(CallChain &chain)
{
// dump result of CallChain in ASCII

    ILwd::LdasContainer* r=NULL;
    cout << "Results:" << endl;
    if ((r = chain.GetResults()))
    {
      cout <<" r="<<r<<endl;
      r->write(cout, ILwd::ASCII);
      cout << endl;
    }
}

UnitTest Test;

template<class T>
bool testInit(enum wavelet_t wavelet_type)
{
 bool pass = true;

 CallChain chain;

// list of parameters
 char* par[5]={"input", "order", "tree", "wavelet_type", NULL};

 int t= int(wavelet_type);

 Scalar<int> ord(10);           // wavelet order
 Scalar<int> wtr(0);            // wavelet tree type 0,1
 Scalar<int> wtp(t); // enumerator goes to CallChain as int

 Sequence<T> in(NN); // input sequence

// input data from file
 ifstream infile("d256.dat");                   // file with input data
 if (! infile )
 {
   cout <<" File d256.dat is not found!"<<endl;
   return false;
 }

 for (int i=0; i<NN && infile ; i++) infile >> in[i];
 infile.close();

 try {

    // Need to create copies of the input data because once AddSymbol
    // has been called on an object it's the CallChains responsibility
    // to delete it
    CallChain::Symbol* input = in.Clone();
    CallChain::Symbol* order = ord.Clone();
    CallChain::Symbol* wtree = wtr.Clone();
    CallChain::Symbol* wtype = wtp.Clone();

    // set names and references for input parameters
    chain.AddSymbol("input", input);
    chain.AddSymbol("order", order);
    chain.AddSymbol("tree", wtree);
    chain.AddSymbol("wavelet_type", wtype);

    // Don't need to do anything for the result to be created, it is
    // created automatically by the LHS of the assignment

    chain.AppendCallFunction("wavelet_init",  par, "output");

    chain.Execute();

// dump output in ASCII
//   ilwd_dump(chain);

// get pointer to output
    udt* out=chain.GetSymbol("output");

// another way to access output 
//    udt* out=const_cast<udt*>(chain.Peek());

    if ( udt::IsA<WaveletUDT<T> >(*out) )
    {  
      WaveletUDT<T>* wu = &(udt::Cast<WaveletUDT<T> >(*out));

// check wavelet metadata
      if ( wu->size()  != NN )   pass = false;
      if ( wu->m_order != ord.GetValue() )   pass = false;
      if ( wu->m_level != 0 )   pass = false;
      if ( wu->m_tree  != wtr.GetValue() )   pass = false;
      if ( wu->m_type  != wtp.GetValue() )   pass = false;

// compare data from output object with data from input
      for( unsigned int i=0; i<NN; i++)
        if ( (*wu)[i] != in[i] ) pass = false;

    }
    else pass = false;

  } // end of 'try'

  catch (std::invalid_argument& x)
  {
      cout << x.what() << endl;
      pass = false;
    }

  catch (std::exception& x)
  {
      cout << x.what() << endl;
      pass = false;
    }

  catch(CallChain::Exception& x)
  {
    cout << "Caught a CallChain exception: " <<x.what()<< endl;
    pass = false;
  }

  catch(...)
  {
    cout << "Caught an Unknown exception.\n";
    pass = false;
  }

  return pass;
}

template<class T>
bool testForward(enum wavelet_t wavelet_type, char* fname, T er)
{
 bool pass = true;

 CallChain chain;

// list of parameters
 char* par[3]={"input", "level", NULL};

 Scalar<int> lev(TEST_LEVEL); 	// requested level of wavelet decomposition
 int ord = TEST_ORDER;   			// wavelet order
 int wtr = TEST_TREE;	// wavelet tree type 0,1

//input wavelet
 WaveletUDT<T>* in = new WaveletUDT<T>(NN, ord, wtr, wavelet_type);

// example of wavelet sequence from file
 Sequence<T> ex(NN);

// input data from file
 ifstream infile("d256.dat");			// file with input data
 if (! infile )
 {
   cout <<" File d256.dat is not found!"<<endl;
   return false;
 }

 for (int i=0; i<NN && infile ; i++) infile >> (*in)[i];
 infile.close();

// input data for comparision with result
 infile.open(fname); 				// file with reference result 
 if (! infile )
 {
   cout <<" File "<<fname<<" is not found!"<<endl;
   return false;
 }

 for (int i=0; i<NN && infile ; i++) infile >> ex[i];
 infile.close();

 Statistics stat;
 double rms_in = stat.rms(ex);

 if (Test.IsVerbose()) cout <<endl<<" Wavelet RMS="<<rms_in<<endl; 

 try {

    // Need to create copies of the input data because once AddSymbol
    // has been called on an object it's the CallChains responsibility
    // to delete it
    CallChain::Symbol* input = in->Clone();
    CallChain::Symbol* level = lev.Clone();

    // set names and references for input parameters
    chain.AddSymbol("input", input);
    chain.AddSymbol("level", level);

    // Don't need to do anything for the result to be created, it is
    // created automatically by the LHS of the assignment

    chain.AppendCallFunction("wavelet_decompose",  par, "output");

    chain.Execute();

// dump CallChain symbol table
//    chain.Dump();

// dump output in ASCII
//   ilwd_dump(chain);

// get pointer to output
    udt* out=chain.GetSymbol("output");

// another way to access output 
//    udt* out=const_cast<udt*>(chain.Peek());

    if ( udt::IsA<Sequence<T> >(*out) )
    {
       Sequence<T>* s = &(udt::Cast<Sequence<T> >(*out));

// compare output data samples with reference data from file

       ex -= *s;		// take difference between sequences

       double rmsd = stat.rms(ex);
       if ( rmsd > er*rms_in ) pass = false;
       double maxd = stat.max(ex);

       if (Test.IsVerbose())
       {
         cout<<" RMS difference ="<<rmsd<<endl;
         cout<<" Max difference ="<<maxd<<endl;
       }
    }
    else pass = false;

// check wavelet metadata
    if ( udt::IsA<WaveletUDT<T> >(*out) )
    {  
      WaveletUDT<T>* wu = &(udt::Cast<WaveletUDT<T> >(*out));
      if ( wu->m_order != ord )   pass = false;
      if ( wu->m_level != lev.GetValue() )   pass = false;
      if ( wu->m_tree  != wtr )   pass = false;
      if ( wu->m_type  != wavelet_type )   pass = false;
    }
    else pass = false;

  } // end of 'try'

  catch (std::invalid_argument& x)
  {
      cout << x.what() << endl;
      pass = false;
    }

  catch (std::exception& x)
  {
      cout << x.what() << endl;
      pass = false;
    }

  catch(CallChain::Exception& x)
  {
    cout << "Caught a CallChain exception: " <<x.what()<< endl;
    pass = false;
  }

  catch(...)
  {
    cout << "Caught an Unknown exception.\n";
    pass = false;
  }

  return pass;
}

template<class T>
bool testInverse(enum wavelet_t wavelet_type, char* fname, T er)
{
 bool pass = true;

 Scalar<int> lev(TEST_LEVEL);   // requested level of wavelet decomposition
 int ord = TEST_ORDER;   	// wavelet order
 int wtr = TEST_TREE;           // wavelet tree type 0,1
 int wtp = int(wavelet_type);

 CallChain chain;

// list of parameters
 char* par[6]={"input", "level", NULL};

// input sequence
 WaveletUDT<T>* in = new WaveletUDT<T>(NN, ord, wtr, wavelet_type);
 in->SetLevel(TEST_LEVEL);

// example of output sequence from file
 Sequence<T> ex(NN); // example of wavelet sequence from file

// input data from file
 ifstream infile(fname);                   // file with input data
 if (! infile )
 {
   cout <<" File "<<fname<<" is not found!"<<endl;
   return false;
 }

 for (int i=0; i<NN && infile ; i++) infile >> (*in)[i];
 infile.close();

// input data for comparision with result
 infile.open("d256.dat");                  // file with reference result
 if (! infile )
 {
   cout <<" File "<<"d256.dat"<<" is not found!"<<endl;
   return false;
 }

 for (int i=0; i<NN && infile ; i++) infile >> ex[i];
 infile.close();

 Statistics stat;
 double rms_in = stat.rms(ex);
 if (Test.IsVerbose()) cout <<endl<<" Input RMS="<<rms_in<<endl;

 try {

    // Need to create copies of the input data because once AddSymbol
    // has been called on an object it's the CallChains responsibility
    // to delete it
    CallChain::Symbol* input = (*in).Clone();
    CallChain::Symbol* level = lev.Clone();

    // set names and references for input parameters
    chain.AddSymbol("input", input);
    chain.AddSymbol("level", level);

    // Don't need to do anything for the result to be created, it is
    // created automatically by the LHS of the assignment

    chain.AppendCallFunction("wavelet_reconstruct",  par, "output");

    chain.Execute();

// dump CallChain symbol table
//    chain.Dump();

// dump output in ASCII
//   ilwd_dump(chain);

// get pointer to output
    udt* out=chain.GetSymbol("output");

// another way to access output
//    udt* out=const_cast<udt*>(chain.Peek());

    if ( udt::IsA<Sequence<T> >(*out) )
    {
       Sequence<T>* s = &(udt::Cast<Sequence<T> >(*out));

// compare output data samples with reference data from file
       ex -= *s;                // take difference between sequences

       double rmsd = stat.rms(ex);
       if ( rmsd > er*rms_in ) pass = false;
       double maxd = stat.max(ex);

       if (Test.IsVerbose())
       {
         cout<<" RMS difference ="<<rmsd<<endl;
         cout<<" Max difference ="<<maxd<<endl;
       }
    }
    else pass = false;

// check wavelet metadata
    if ( udt::IsA<WaveletUDT<T> >(*out) )
    {
      WaveletUDT<T>* wu = &(udt::Cast<WaveletUDT<T> >(*out));
      if ( wu->m_order != ord )   pass = false;
      if ( wu->m_level != 0 )   pass = false;
      if ( wu->m_tree  != wtr )   pass = false;
      if ( wu->m_type  != wtp )   pass = false;
    }
    else pass = false;

  } // end of 'try'
  catch (std::invalid_argument& x)
  {
      cout << x.what() << endl;
      pass = false;
    }

  catch (std::exception& x)
  {
      cout << x.what() << endl;
      pass = false;
    }

  catch(CallChain::Exception& x)
  {
    cout << "Caught a CallChain exception: " <<x.what()<< endl;
    pass = false;
  }

  catch(...)
  {
    cout << "Caught an Unknown exception.\n";
    pass = false;
  }

  return pass;
}

int main(int ArgC, char** ArgV)
{
  Test.Init(ArgC, ArgV);

  if (Test.IsVerbose())
  {
    cout << "$Id: tWaveletFunctions.cc,v 0.2 2001/06/06 20:02:00 sazonov Exp $"
    << endl << endl;
  }

  Test.Check(testInit<float>(lifting))
  << "(test if lifting wavelet for <float> initialized properly)" << endl;

  Test.Check(testInit<double>(lifting))
  << "(test if lifting wavelet for <double> initialized properly)" << endl;

  Test.Check(testInit<float>(Daubechies))
  << "(test if Daubechies wavelet for <float> initialized properly)" << endl;

  Test.Check(testInit<double>(Daubechies))
  << "(test if Daubechies wavelet for <double> initialized properly)" << endl;

  Test.Check(testInit<float>(Meyer))
  << "(test if Meyer wavelet for <float> initialized properly)" << endl;

  Test.Check(testInit<double>(Meyer))
  << "(test if Meyer wavelet for <double> initialized properly)" << endl;

  Test.Check(testForward<float>(lifting, "wl256_4_10_0.dat", 1.e-9))
  << "(test if forward lifting transform for <float> computed properly)"
  << endl;

  Test.Check(testForward<double>(lifting, "wl256_4_10_0.dat",1.e-14))
  << "(test if forward lifting transform for <double> computed properly)"
  << endl;

  Test.Check(testForward<float>(Daubechies, "wd256_4_10_0.dat", 1.e-9))
  << "(test if forward Daubechies transform for <float> computed properly)" 
  << endl;

  Test.Check(testForward<double>(Daubechies, "wd256_4_10_0.dat",1.e-14))
  << "(test if forward Daubechies transform for <double> computed properly)" 
  << endl;

  Test.Check(testForward<float>(Meyer, "wm256_4_x_0.dat", 1.e-9))
  << "(test if forward Meyer transform for <float> computed properly)"
  << endl;

  Test.Check(testForward<double>(Meyer, "wm256_4_x_0.dat",1.e-14))
  << "(test if forward Meyer transform for <double> computed properly)"
  << endl;

  Test.Check(testInverse<float>(lifting, "wl256_4_10_0.dat", 1.e-7))
  << "(test if inverse lifting transform for <float> computed properly)"
  << endl;

  Test.Check(testInverse<double>(lifting, "wl256_4_10_0.dat",1.e-12))
  << "(test if inverse lifting transform for <double> computed properly)"
  << endl;

  Test.Check(testInverse<float>(Daubechies, "wd256_4_10_0.dat", 1.e-7))
  << "(test if inverse Daubechies transform for <float> computed properly)"
  << endl;

  Test.Check(testInverse<double>(Daubechies, "wd256_4_10_0.dat",1.e-11))
  << "(test if inverse Daubechies transform for <double> computed properly)"
  << endl;

  Test.Check(testInverse<float>(Meyer, "wm256_4_x_0.dat", 1.e-2))
  << "(test if inverse Meyer transform for <float> computed properly)"
  << endl;

  Test.Check(testInverse<double>(Meyer, "wm256_4_x_0.dat",1.e-2))
  << "(test if inverse Meyer transform for <double> computed properly)"
  << endl;

// all done!!
  Test.Exit();
}