/*! \file sigfit.cc
 * \brief fit signal by trial-signals
 * 
 * ----------------------------------------------------------------------------
 * 
 * $Id: sigfit.cc,v 1.9 2004-02-15 20:43:09 tforb Exp $
 * \author Thomas Forbriger
 * \date 28/01/2004
 * 
 * fit signal by trial-signals
 * 
 * Copyright (c) 2004 by Thomas Forbriger (BFO Schiltach) 
 * 
 * REVISIONS and CHANGES 
 *  - 28/01/2004   V1.0   Thomas Forbriger
 * 
 * ============================================================================
 */
#define SIGFIT_VERSION \
  "SIGFIT   V1.0   fit signal by trial-signals"
#define SIGFIT_CVSID \
  "$Id: sigfit.cc,v 1.9 2004-02-15 20:43:09 tforb Exp $"

#include <fstream>
#include <iostream>
#include <list>
#include <vector>
#include <tfxx/commandline.h>
#include <tfxx/error.h>
#include <sffxx.h>
#include <sffostream.h>
#include <tsxx/tsxx.h>
#include <tsxx/innerproduct.h>
#include <aff/array.h>
#include <aff/dump.h>
#include <aff/seriesoperators.h>
#include <linearxx/lapackxx.h>

typedef std::list<std::string> Tnamelist;
typedef ts::TDsfftimeseries Tbundle;
typedef Tbundle::Tseries Tseries;
typedef std::vector<Tbundle> Tbundlevec;
typedef aff::Array<Tbundle::Tseries::Tvalue> Tmatrix;

using std::cout;
using std::cerr;
using std::endl;

// structure to store commandline options
struct Options {
  bool verbose;
  double Tdate;
  bool truncate;
  std::string residualname;
  bool writeresidual;
}; // struct Options

// formatted number output
std::string formatfloat(const double &v)
{
  char seq[20];
  std::sprintf(seq, "%8.5f ", v);
  return(std::string(seq));
}

int main(int iargc, char* argv[])
{

  // define usage information
  char usage_text[]=
  {
    SIGFIT_VERSION "\n"
    "usage: sigfit [-v] [-Tdate v] [-truncate]" "\n"
    "              signal trial [trial ...]" "\n"
    "   or: sigfit --help|-h" "\n"
  };

  // define full help text
  char help_text[]=
  {
    SIGFIT_CVSID "\n"
    "\n"
    "-v           be verbose" "\n"
    "-Sramp       add a ramp to the set of trial signals" "\n"
    "-Sconst      add a constant to the set of trial signals" "\n"
    "-Tdate v     tolerance for comparison of date of first sample" "\n"
    "             v give the tolerance in units of the sampling interval" "\n"
    "-truncate    truncate all series to the same number of samples" "\n"
    "-residual f  write waveforms containing residual to file \"f\"" "\n"   
    "\n"
    "signal       signal to by explained by a linear combination" "\n"
    "             of trial signals" "\n"
    "trial        any set of trial signals"
  };

  // define commandline options
  using namespace tfxx::cmdline;
  static Declare options[]= 
  {
    // 0: print help
    {"help",arg_no,"-"},
    // 1: verbose mode
    {"v",arg_no,"-"},
    // 2: date tolerance
    {"Tdate",arg_yes,"0."},
    // 3: truncate to a common number of samples
    {"truncate",arg_no,"-"},
    // 4: output SFF file
    {"residual",arg_yes,"-"},
    {NULL}
  };

  // no arguments? print usage...
  if (iargc<2) 
  {
    cerr << usage_text << endl;
    exit(0);
  }

  // collect options from commandline
  Commandline cmdline(iargc, argv, options);

  // help requested? print full help text...
  if (cmdline.optset(0))
  {
    cerr << usage_text << endl;
    cerr << help_text << endl;
    exit(0);
  }

/*
  // dummy operation: print option settings
  for (int iopt=0; iopt<2; iopt++)
  {
    cout << "option: '" << options[iopt].opt_string << "'" << endl;
    if (cmdline.optset(iopt)) {  cout << "  option was set"; }
    else { cout << "option was not set"; }
    cout << endl;
    cout << "  argument (string): '" << cmdline.string_arg(iopt) << "'" << endl;
    cout << "     argument (int): '" << cmdline.int_arg(iopt) << "'" << endl;
    cout << "    argument (long): '" << cmdline.long_arg(iopt) << "'" << endl;
    cout << "   argument (float): '" << cmdline.float_arg(iopt) << "'" << endl;
    cout << "  argument (double): '" << cmdline.double_arg(iopt) << "'" << endl;
    cout << "    argument (bool): '";
    if (cmdline.bool_arg(iopt))
    { cout << "true"; } else { cout << "false"; }
    cout << "'" << endl;
  }
  while (cmdline.extra()) { cout << cmdline.next() << endl; }

  // dummy operation: print rest of command line
  while (cmdline.extra()) { cout << cmdline.next() << endl; }
*/

  // read options
  Options opt;
  opt.verbose=cmdline.optset(1);
  opt.Tdate=cmdline.double_arg(2);
  opt.truncate=cmdline.optset(3);
  opt.writeresidual=cmdline.optset(4);
  opt.residualname=cmdline.string_arg(4);

  /*----------------------------------------------------------------------*/

  // read file names
  TFXX_assert(cmdline.extra(),"ERROR: missing signal file name!");
  std::string signalname=cmdline.next();

  TFXX_assert(cmdline.extra(),"ERROR: missing trial signal file name!");
  Tnamelist namelist;
  while (cmdline.extra()) { namelist.push_back(cmdline.next()); }

  // read signal file
  Tbundle signal;
  {
    if (opt.verbose)
    { cout << "read signal file \"" << signalname << "\"" << endl; }
    std::ifstream is(signalname.c_str());
    sff::FileHeader fileheader(is);
    if (opt.verbose)
    { cout << "  read trace" << endl; }
    sff::InputWaveform<Tseries> inputwaveform(is);
    sff::TraceHeader traceheader=inputwaveform.header();
    signal.header=traceheader.wid2();
    signal=inputwaveform.series();
    bool last=traceheader.last();
    std::string wid2line=signal.header.line();
    if (opt.verbose)
    { 
      cout << "  " << wid2line.substr(0,69) << endl; 
      if (!last)
      {
        cout << "there are more traces in that files" << endl 
          << "ignoring them..." << endl;
      }
    }
  }
  
  // read trial files
  if (opt.verbose)
  { cout << namelist.size() << " trial data files to read" << endl; }
  Tbundlevec bundlevec;
  for (Tnamelist::const_iterator i=namelist.begin(); i!=namelist.end(); i++)
  {
    if (opt.verbose)
    { cout << "read trial data file \"" << *i << "\"" << endl; }
    std::ifstream is(i->c_str());
    sff::FileHeader fileheader(is);
    bool last=false;
    while(!last)
    {
      if (opt.verbose)
      { cout << "  read trial signal #" << bundlevec.size()+1 << endl; }
      sff::InputWaveform<Tseries> inputwaveform(is);
      sff::TraceHeader traceheader=inputwaveform.header();
      Tbundle bundle;
      bundle.header=traceheader.wid2();
      bundle=inputwaveform.series();
      bundlevec.push_back(bundle);
      last=traceheader.last();
      std::string wid2line=bundle.header.line();
      if (opt.verbose)
      { cout << "  " << wid2line.substr(0,69) << endl; }
    }
  }

  /*----------------------------------------------------------------------*/

  // truncate length if requested
  if (opt.truncate)
  {
    if (opt.verbose) { cout << "truncate signals if necessary..." << endl; }
    long int n=signal.last();
    for (Tbundlevec::const_iterator i=bundlevec.begin(); 
         i!=bundlevec.end(); i++)
    { n=n<i->last() ? n : i->last(); }
    if (signal.last()>n)
    {
      signal.setlastindex(n);
      signal.header.nsamples=signal.size();
      if (opt.verbose) { cout << "truncate signal to "
        << signal.header.nsamples << " samples" << endl; }
    }
    for (Tbundlevec::iterator i=bundlevec.begin(); 
         i!=bundlevec.end(); i++)
    { 
      if (i->last()>n)
      {
        i->setlastindex(n);
        i->header.nsamples=signal.size();
        if (opt.verbose) { cout << "truncate trial signal" << endl
          << i->header.line() << endl; }
      }
    }
  }

  /*----------------------------------------------------------------------*/

  // check header consistency
  sff::WID2compare compare(sff::Fnsamples | sff::Fdt | sff::Fdate);
  compare.setdatetolerance(opt.Tdate);
  if (opt.verbose) { cout << "checking consistency..." << endl; }
  for (Tbundlevec::const_iterator i=bundlevec.begin(); i!=bundlevec.end(); i++)
  {
    if (!compare (i->header,signal.header))
    {
      cerr << "ERROR: header signature mismatch:" << endl;
      cerr << "signal:" << endl;
      cerr << signal.header.line();
      cerr << "trial signal:" << endl;
      cerr << i->header.line();
      TFXX_abort("baling out...");
    }
  }

  /*----------------------------------------------------------------------*/
  
  // set up system of linear equations
  if (opt.verbose) { cout << "set up system of linear equations..." << endl; }
  int N=bundlevec.size();
  if (opt.verbose) { cout << "system is of size " << N << "x" << N << endl; }
  Tmatrix Matrix(N,N), rhs(N), coeff(N);
  for (int i=1; i<=N; ++i)
  {
    for (int k=i; k<=N; ++k)
    {
      Matrix(i,k)=ts::innerproduct(bundlevec[i-1], bundlevec[k-1]);
      Matrix(k,i)=Matrix(i,k);
    }
    rhs(i)=ts::innerproduct(bundlevec[i-1], signal);
  }

  // solve
  coeff=linear::lapack::dposv(Matrix,rhs);

  // set up synthetics
  Tbundle synthetics;
  synthetics=Tseries(signal.shape());
  synthetics.series()=0.;
  for (int i=1; i<=N; ++i)
  { synthetics += coeff(i) * bundlevec[i-1]; }
  synthetics.header=signal.header;
  synthetics.header.channel="synt";

  // set up residual
  Tbundle residual;
  residual=Tseries(signal.shape());
  residual=signal-synthetics;
  residual.header=signal.header;
  residual.header.channel="diff";

  // write waveforms
  if (opt.writeresidual)
  {
    if (opt.verbose) 
    {
      cout << "write residual waveform to " 
      << opt.residualname << endl;
    }
    std::ofstream ofs(opt.residualname.c_str());
    sff::SFFostream<Tseries> os(ofs);
    os << signal;
    os << synthetics;
    os << residual;
    for (Tbundlevec::const_iterator i=bundlevec.begin(); 
         i!=bundlevec.end(); i++)
    { os << *i; }
    for (int i=1; i<=N; ++i)
    {
      Tbundle psyn=bundlevec[i-1];
      psyn=bundlevec[i-1].series()*coeff(i);
      char chan[6];
      std::sprintf(chan,"s%1.1d", i);
      psyn.header.channel=std::string(chan);
      os << psyn;
    }
  }
  
  // calculate rms values
  double signalrms=ts::rms(signal);
  double residualrms=ts::rms(residual);
  double explained=1.-(residualrms/signalrms);
  
  std::string timestring(signal.header.date.timestring());
  cout << "         time: " << timestring.substr(4,21) << endl;
  cout << "      channel: " << signal.header.channel << endl;
  cout << "      station: " << signal.header.station << endl;
  cout << "   instrument: " << signal.header.instype << endl;
  cout << "    signalrms: " << signalrms << endl;
  cout << "  residualrms: " << residualrms << endl;
  cout << "explained rms: " << explained << endl;
  cout << " coefficients: ";
  for (int i=1; i<=N; ++i)
  { cout << coeff(i) << "  "; }
  cout << endl;
  cout << " coefficients x rms: ";
  for (int i=1; i<=N; ++i)
  { 
    cout << coeff(i)*ts::rms(bundlevec[i-1]) << "  "; 
  }
  cout << endl;

  // reportline
  cout << signal.header.station << " "
    << signal.header.channel << " "
    << signal.header.instype << " "
    << timestring << " ";
  cout << formatfloat(explained);
  cout << formatfloat(signalrms);
  cout << formatfloat(residualrms);
  for (int i=1; i<=N; ++i)
  { 
    cout << formatfloat(coeff(i));
  }
  for (int i=1; i<=N; ++i)
  { 
    cout << formatfloat(coeff(i)*ts::rms(bundlevec[i-1])); 
  }
  cout << endl;
}

/* ----- END OF sigfit.cc ----- */