// -*- mode: C++; c-indent-level: 4; c-basic-offset: 4; indent-tabs-mode: nil; -*-

#include <RcppEigen.h>

// [[Rcpp::depends(RcppEigen)]]

typedef Eigen::MappedSparseMatrix<double>         MSpMat;
typedef Eigen::SparseMatrix<double>                SpMat;
typedef Eigen::SimplicialLDLT<SpMat>              SpChol;
typedef Eigen::MatrixXd                              Mat;
typedef Eigen::VectorXd                              Vec;
typedef Eigen::Map<Mat>                             MMat;
typedef Eigen::Map<Vec>                             MVec;
typedef Eigen::VectorXi                             iVec;
typedef Eigen::Map<iVec>                           MiVec;
typedef Eigen::CholmodDecomposition<SpMat>     ChmDecomp;

using Rcpp::as;
using Rcpp::wrap;
using Rcpp::Environment;

// [[Rcpp::export]]
ChmDecomp update_predD(Environment rho) {
  // Calls to rho.get can be preceded by calls to rho.exists to make sure the name is available
  MSpMat        Zt(as<MSpMat>(rho.get("Zt")));
  MSpMat   Lambdat(as<MSpMat>(rho.get("Lambdat")));
  MiVec       Lind(as<MiVec>(rho.get("Lind")));
  MVec       theta(as<MVec>(rho.get("theta")));
  if (Lind.size() != Lambdat.nonZeros()) throw std::invalid_argument("dimension mismatch");
  int    *lipt = Lind.data();
	double *LamX = Lambdat.valuePtr(), *thpt = theta.data();
	for (int i = 0; i < Lind.size(); ++i) {
	    LamX[i] = thpt[lipt[i] - 1];
	}
  SpMat  LamtZt(Lambdat * Zt);
  SpMat  LamtZtZLam(LamtZt * LamtZt.adjoint());
  rho.assign("ZtZ",wrap(LamtZtZLam));
  return ChmDecomp(LamtZtZLam);//  ChmDecomp      L = as<ChmDecomp>(rho.get("L"));
}