#ifndef __FILE_MOMENT_NORM_H_SEEN__
#define __FILE_MOMENT_NORM_H_SEEN__

#include "libscl.h"
#include "tree_sim.h"
#include "tree_mf.h"

class moment_norm : public scl::nleqns_base {
private:
  scl::realmat data;
  const INTEGER n_obs;
  const INTEGER mf_lags;
  const INTEGER HAC_lags;
  scl::realmat e;
  scl::realmat theta;
  INTEGER t;
  tree_moment_function tree_mf;
  scl::gmm tree_gmm;
  INTEGER d;
  scl::realmat m;
  scl::realmat W;
  REAL logdetW;
  INTEGER rankW;
  scl::realmat S;
  scl::realmat R;
  scl::realmat Z;
public:
  moment_norm(const scl::realmat& dat, INTEGER n, INTEGER mfl, INTEGER hacl)
  : data(dat), n_obs(n), mf_lags(mfl), HAC_lags(hacl)
  {
    scl::gmm tmp_gmm(&tree_mf, mf_lags, &data, n_obs, HAC_lags);
    tmp_gmm.set_correct_W_for_mean(true);
    tree_gmm.set_warning_messages(true);
    tree_gmm = tmp_gmm;
    d = tree_mf.get_d();
  }
  void set_m_theta(const scl::realmat& u, const scl::realmat& parms)
  { e = u; theta = parms; } 
  void set_t(INTEGER obs) 
  { t = obs; if (t<1 || n_obs<t) scl::error("Error, bad t"); }
  void set_x(const scl::realmat& x) 
  { 
    INTEGER xsize = x.size();
    INTEGER rows = data.nrow();
    INTEGER base = rows*(t - 1);
    if (base + xsize > n_obs*rows) scl::error("Error, pre_image norm, bad x");
    for (INTEGER i=1; i<=xsize; ++i) data[base + i] = x[i];
  }
  void set_data(const scl::realmat& dat) { data = dat; }
  scl::realmat get_data() { return data; }
  INTEGER get_d() { return d; }
  INTEGER get_n() { return n_obs; }
  bool get_f(const scl::realmat& x, scl::realmat& f)
  {
    set_x(x);
    bool rv = tree_gmm.set_data(&data);
    tree_gmm(theta,m,W,logdetW,rankW,S);
    if (tree_gmm.get_W_numerr() > 0) {
      scl::warn("Warning, moment_norm, numerr > 0");
      rv = false;
      f.resize(1,1,REAL_MAX);
      return rv;
    }
    REAL sum = 0.0;
    for (INTEGER i=1; i<=d; ++i) sum += pow(m[i]-e[i],2);
    sum = sqrt(sum);
    f.resize(1,1,sum);
    return rv;
  }
  bool get_F(const scl::realmat& x, scl::realmat& f, scl::realmat& F)
  {
    if (this->get_f(x,f)) {
      return nleqns_base::df(x,F);
    }
    else {
      return false;
    }
  }
};

#endif
