#ifndef __FILE_PREIMAGE_EQNS_H_SEEN__
#define __FILE_PREIMAGE_EQNS_H_SEEN__

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

class preimage_eqns : public scl::nleqns_base {
private:
  scl::realmat data;
  const INTEGER n_obs;
  const INTEGER mf_lags;
  const INTEGER HAC_lags;
  scl::realmat z;
  scl::realmat rho;
  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:
  preimage_eqns (const scl::realmat& dat, INTEGER n, INTEGER ml, INTEGER hl)
  : data(dat), n_obs(n), mf_lags(ml), HAC_lags(hl)
  {
    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 update_z_rho(const scl::realmat& u, const scl::realmat& parms)
  { z = u; rho = parms; } 
  void update_data(const scl::realmat& x) 
  {
    INTEGER r = data.nrow();
    INTEGER last = r*n_obs;
    for (INTEGER i=0; i<d; ++i) {
      data[last - i] = x[d - i]; 
    }
  }
  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)
  {
    update_data(x);
    bool rv = tree_gmm.set_data(&data);
    tree_gmm(rho,m,W,logdetW,rankW,S);
    if (cholesky(W,R) != d) rv = false;
    Z = sqrt(n_obs)*R*m;
    f = Z - z;
    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
