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

Copyright (C) 2012, 2013

A. Ronald Gallant
Post Office Box 659
Chapel Hill NC 27514-0659
USA   

This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.

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

#include "svmnt.h"
#include "svparms.h"

using namespace std;
using namespace scl;

bool svmnt::set_data(const scl::realmat* datptr)
{
  data = datptr;
  return true;
}

bool svmnt::set_particle(const realmat* partptr)
{
  particle = partptr;
  return true;
}

bool svmnt::set_sample_size(INTEGER Tmax)
{
  n = Tmax;      // validity of n checked in operator()
  if (n <= 0) error("Error, svmnt, Tmax not positive");
  return true;
}

bool svmnt::set_theta(const scl::realmat& parm)
{
  if (parm.size() == 4) {
    theta = parm;
    return true;
  }
  else {
    theta.resize(4,1);
    theta[1] = T1;  //rho
    theta[2] = T2;  //phi
    theta[3] = T3;  //sigma
    theta[4] = T4;  //beta
    error("Error, svmnt, size theta != 4");
    return false;
  }
}
  
#if defined USE_FOC_MOMENT_CONDITIONS

  realmat svmnt::operator() (INTEGER t)
  {
    if (data == 0) error("Error, svmnt, data not set");
    if (n > data.ncol()) error("Error, svmnt, n too large");
    if (particle == 0) error("Error, svmnt, particle not set");
    if (t < get_T0()) error("Error, svmnt, t too small");
    if (t > n) error("Error, svmnt, t too large");
    if (t > particle->ncol()) error("Error, svmnt, t too large");
  
    REAL rho   = theta[1];
    REAL phi   = theta[2];
    REAL sigma = theta[3];
    REAL beta  = theta[4];
  
    const realmat* part = particle;
  
    REAL xres = (*part)[t] - phi*(*part)[t-1];
    REAL yres = (*data)[t] - rho*(*data)[t-1];

    realmat mt(4,1);
  
    mt[1] = xres*(*part)[t-1];
    mt[2] = pow(sigma,2) - pow(xres,2);
    mt[3] = pow(beta*exp((*part)[t]),2) - pow(yres,2);
    mt[4] = yres*(*data)[t-1];
  
    return mt;
  }

#else

  realmat svmnt::operator() (INTEGER t)
  {
    if (t < get_T0()) error("Error, svmnt, t too small");
    if (t > n) error("Error, svmnt, t too large or n not set");
  
    REAL rho = theta[1];
    REAL phi = theta[2];
    REAL sigma = theta[3];
    REAL beta = theta[4];
  
    const realmat* part = particle;
  
    realmat mt(L+4,1);

    REAL xres = (*part)[t] - phi*(*part)[t-1];
    REAL yres = (*data)[t] - rho*(*data)[t-1];
  
    REAL vt = fabs(yres);
    REAL Vt = beta*exp( (*part)[t] );
  
    mt[1] = vt*vt - Vt*Vt;
  
    const REAL sclfac = sqrt(2.0/3.1415926535897932);
  
    for (INTEGER j=1; j<=L; ++j) {

      REAL yreslagj = (*data)[t-j] - rho*(*data)[t-1-j];
  
      REAL vs = fabs(yreslagj);
      REAL Vs = beta*exp( (*part)[t-j] );
  
      Vt *= sclfac;
      Vs *= sclfac;
  
      mt[j+1] = vt*vs - Vt*Vs;
  
    }

    mt[L+2] = yres*(*data)[t-1];
    mt[L+3] = xres*(*part)[t-1];
    mt[L+4] = pow(sigma,2) - pow(xres,2);
  
    return mt;
  }

#endif

