density.hpp

Go to the documentation of this file.

// Copyright (C) 2013-2015 Kasper Kristensen
// License: GPL-2

namespace density {
using namespace tmbutils;

#define TYPEDEFS(scalartype_)                   \
public:                                         \
typedef scalartype_ scalartype;                 \
typedef vector<scalartype> vectortype;          \
typedef matrix<scalartype> matrixtype;          \
typedef array<scalartype> arraytype

#define VARIANCE_NOT_YET_IMPLEMENTED            \
private:                                        \
vectortype variance(){return vectortype();}     \
public:

#define SIMULATE_NOT_YET_IMPLEMENTED            \
private:                                        \
vectortype sqrt_cov_scale(vectortype x){}       \
void simulate(vectortype &u){}                  \
vectortype simulate(){}                         \
public:

/* Add this macro to all classes with *dynamic* dimension (e.g. AR1) */
#define SIMULATE_IMPLEMENTED_UNKNOWN_SIZE       \
void simulate(vectortype &x) {                  \
  rnorm_fill(x);                                \
  x = sqrt_cov_scale(x);                        \
  x = zero_derivatives(x);                      \
}
/* Add this macro to all classes with *fixed* dimension (e.g. MVNORM) */
#define SIMULATE_IMPLEMENTED_KNOWN_SIZE(SIZE)   \
SIMULATE_IMPLEMENTED_UNKNOWN_SIZE               \
vectortype simulate() {                         \
  vectortype x(SIZE);                           \
  simulate(x);                                  \
  return x;                                     \
}

/* Utility function: The simulators should not track derivatives when
   running with AD types. A workaround is to zero out the derivatives
   after simulation (FIXME: there is a minor efficiency loss by
   tracking the derivatives in the first place...). */
template<class arraytype>
arraytype zero_derivatives(arraytype x) {
  for(int i=0; i<x.size(); i++)
    x(i) = asDouble(x(i));
  return x;
}

/* Utility function */
template<class arraytype>
void rnorm_fill(arraytype &x) {
  for(int i=0; i<x.size(); i++)
    x(i) = Rf_rnorm(0., 1.);
}

template <class scalartype_>
class MVNORM_t{
  TYPEDEFS(scalartype_);
  matrixtype Q;       /* Inverse covariance matrix */
  scalartype logdetQ; /* log-determinant of Q */
  matrixtype Sigma;   /* Keep for convenience - not used */
  matrixtype L_Sigma; /* Used by simulate() */
public:
  MVNORM_t(){}
  MVNORM_t(matrixtype Sigma_, bool use_atomic=true){
    setSigma(Sigma_, use_atomic);
  }

  matrixtype cov(){return Sigma;}

  /* initializer via covariance matrix */
  void setSigma(matrixtype Sigma_, bool use_atomic=true){
    Sigma = Sigma_;
    scalartype logdetS;
    if(use_atomic){
      Q = atomic::matinvpd(Sigma,logdetS);
    } else {
      matrixtype I(Sigma.rows(),Sigma.cols());
      I.setIdentity();
      Eigen::LDLT<Eigen::Matrix<scalartype,Dynamic,Dynamic> > ldlt(Sigma);
      Q = ldlt.solve(I);
      vectortype D = ldlt.vectorD();
      logdetS = D.log().sum();
    }
    logdetQ = -logdetS;
  }
  scalartype Quadform(vectortype x){
    return (x*(vectortype(Q*x))).sum();
  }
  scalartype operator()(vectortype x){
    return -scalartype(.5)*logdetQ + scalartype(.5)*Quadform(x) + x.size()*scalartype(log(sqrt(2.0*M_PI)));
  }
  scalartype operator()(vectortype x, vectortype keep){
    matrix<scalartype> S = Sigma;
    vector<scalartype> not_keep = scalartype(1.0) - keep;
    for(int i = 0; i < S.rows(); i++){
      for(int j = 0; j < S.cols(); j++){
        S(i,j) = S(i,j) * keep(i) * keep(j);
      }
      S(i,i) += not_keep(i) * scalartype(1.0 / (2.0 * M_PI));
    }
    return MVNORM_t<scalartype>(S)(x * keep);
  }
  arraytype jacobian(arraytype x){
    arraytype y(x.dim);
    matrixtype m(x.size()/x.cols(),x.cols());
    for(int i=0;i<x.size();i++)m(i)=x[i];
    matrixtype mQ=m*Q;
    for(int i=0;i<x.size();i++)y[i]=mQ(i);
    return y;
  }
  int ndim(){return 1;}
  VARIANCE_NOT_YET_IMPLEMENTED
  vectortype sqrt_cov_scale(vectortype u) {
    if(L_Sigma.rows() == 0) {
      Eigen::LLT<Eigen::Matrix<scalartype,Dynamic,Dynamic> > llt(Sigma);
      L_Sigma = llt.matrixL();
    }
    vectortype ans = L_Sigma * u;
    return ans;
  }
  SIMULATE_IMPLEMENTED_KNOWN_SIZE(Sigma.rows());
};

template <class scalartype>
MVNORM_t<scalartype> MVNORM(matrix<scalartype> Sigma, bool use_atomic = true){
  return MVNORM_t<scalartype>(Sigma, use_atomic);
}

template <class scalartype_>
class UNSTRUCTURED_CORR_t : public MVNORM_t<scalartype_>{
  TYPEDEFS(scalartype_);
  UNSTRUCTURED_CORR_t(){}
  UNSTRUCTURED_CORR_t(vectortype x){
    // (n*n-n)/2=nx  ==>  n*n-n-2*nx=0 ==> n=(1+sqrt(1+8*nx))/2
    int nx=x.size();
    int n=int((1.0+sqrt(1.0+8*nx))/2.0);
    if((n*n-n)/2!=nx)std::cout << "vector does not specify an UNSTRUCTERED_CORR\n";
    matrixtype L(n,n);
    L.setIdentity();
    int i,j,k=0;
    for(i=0;i<L.rows();i++){
      for(j=0;j<L.cols();j++){
        if(i>j){L(i,j)=x[k];k++;}
      }
    }
    matrixtype llt=L*L.transpose();
    matrixtype Sigma=llt;
    for(i=0;i<Sigma.rows();i++){
      for(j=0;j<Sigma.cols();j++){
        Sigma(i,j)/=sqrt(llt(i,i)*llt(j,j));
      }
    }    
    this->setSigma(Sigma); /* Call MVNORM_t initializer */
  }
};
template <class scalartype>
UNSTRUCTURED_CORR_t<scalartype> UNSTRUCTURED_CORR(vector<scalartype> x){
  return UNSTRUCTURED_CORR_t<scalartype>(x);
}

template<class scalartype_> 
class N01{
  TYPEDEFS(scalartype_);
public:
  scalartype operator()(scalartype x){
    return x*x*.5 + log(sqrt(2.0*M_PI));
  }
  scalartype operator()(vectortype x){
    return (x*x*scalartype(.5) + scalartype(log(sqrt(2.0*M_PI))) ).sum() ;
  }
  scalartype operator()(arraytype x){
    return (x*x*scalartype(.5) + scalartype(log(sqrt(2.0*M_PI))) ).sum() ;
  }
  arraytype jacobian(arraytype x){return x;}
  int ndim(){return 1;}
  VARIANCE_NOT_YET_IMPLEMENTED
  // FIXME: 1D will not suffice for e.g. SEPARABLE(N01, OTHER);
  vectortype simulate() {
    vectortype x(1);
    x[0] = rnorm(0.0, 1.0);
    return x;
  }
  // Inplace simulate
  void simulate(vectortype &x) {
    rnorm_fill(x);
  }
  vectortype sqrt_cov_scale(vectortype u) {
    return u;
  }

};

template <class distribution>
class AR1_t{
  TYPEDEFS(typename distribution::scalartype);
private:
  scalartype phi;
  distribution MARGINAL;
public:
  AR1_t(){/*phi=phi_;MARGINAL=f_;*/}
  AR1_t(scalartype phi_, distribution f_) : phi(phi_), MARGINAL(f_) {}
  scalartype operator()(vectortype x){
    scalartype value;
    value=scalartype(0);
    int n=x.rows();
    int m=x.size()/n;
    scalartype sigma=sqrt(scalartype(1)-phi*phi); /* Steady-state standard deviation */
    value+=MARGINAL(x(0));                                       /* E.g. x0 ~ N(0,1)  */
    for(int i=1;i<n;i++)value+=MARGINAL((x(i)-x(i-1)*phi)/sigma);/* x(i)-phi*x(i-1) ~ N(0,sigma^2) */
    value+=scalartype((n-1)*m)*log(sigma);
    return value;
  }

  /* Copied vector version - apply over outermost dimension */
  scalartype operator()(arraytype x){
    scalartype value;
    value=scalartype(0);
    int n=x.cols();
    int m=x.size()/n;
    scalartype sigma=sqrt(scalartype(1)-phi*phi); /* Steady-state standard deviation */
    value+=MARGINAL(x.col(0));                                  /* E.g. x0 ~ N(0,1)  */
    for(int i=1;i<n;i++)value+=MARGINAL((x.col(i)-x.col(i-1)*phi)/sigma);/* x(i)-phi*x(i-1) ~ N(0,sigma^2) */
    value+=scalartype((n-1)*m)*log(sigma);
    return value;
  }


  arraytype jacobian(arraytype x){
    scalartype sigma=sqrt(scalartype(1)-phi*phi); /* Steady-state standard deviation */
    int n=x.cols();
    arraytype y(x.dim);
    y.setZero();
    y.col(0) = y.col(0) + MARGINAL.jacobian(x.col(0));
    for(int i=1;i<n;i++){
      //MARGINAL((x(i)-x(i-1)*phi)/sigma);
      y.col(i-1) = y.col(i-1) - (phi/sigma) * MARGINAL.jacobian((x.col(i)-x.col(i-1)*phi)/sigma);
      y.col(i) = y.col(i) +  MARGINAL.jacobian((x.col(i)-x.col(i-1)*phi)/sigma)/sigma;
    }
    return y;
  }
  int ndim(){return 1;}
  VARIANCE_NOT_YET_IMPLEMENTED

  /* Simulation */
  arraytype sqrt_cov_scale(arraytype u) {
    if(u.dim.size() == 1) {
      // FIXME: This entire class should be reconsidered with this
      // special case in mind !
      u.dim.resize(2); u.dim << 1, u.size();
    }
    int n = u.cols();
    arraytype x(u.dim);
    scalartype sigma = sqrt(1.0 - phi * phi); /* Steady-state standard deviation */
    x.col(0) = MARGINAL.sqrt_cov_scale(u.col(0));
    for(int i=1; i<n; i++)
      x.col(i) = phi * x.col(i-1) + sigma * MARGINAL.sqrt_cov_scale(u.col(i));
    return x;
  }
  vectortype sqrt_cov_scale(vectortype u) {
    vector<int> dim(2);
    dim << 1, u.size();
    arraytype u_array(u, dim);
    return sqrt_cov_scale(u_array);
  }
  void simulate(vectortype &x) {
    rnorm_fill(x);
    x = sqrt_cov_scale(x);
    x = zero_derivatives(x);
  }
  void simulate(arraytype &x) {
    rnorm_fill(x);
    x = sqrt_cov_scale(x);
    x = zero_derivatives(x);
  }

};
template <class scalartype, class distribution>
AR1_t<distribution> AR1(scalartype phi_, distribution f_){
  return AR1_t<distribution>(phi_, f_);
}
template <class scalartype>
AR1_t<N01<scalartype> > AR1(scalartype phi_){
  return AR1_t<N01<scalartype> >(phi_, N01<scalartype>());
}


template <class scalartype_>
class ARk_t{
  TYPEDEFS(scalartype_);
  //private:
  int k;
  vectortype phi;   /* [phi1,...,phik] */
  vectortype gamma; /* [gamma(1),...,gamma(k)] (note gamma(0) is 1) */
  /* Initial distribution matrices. */
  matrixtype V0;    /* kxk variance  */
  matrixtype Q0;    /* kxk precision */
  matrixtype L0;    /* kxk Cholesky Q0 = L0*L0' */
  /* gamma is found through (I-M)*gamma=phi ... */
  matrixtype M;     /* kxk   */
  matrixtype I;     /* kxk   */
  scalartype sigma;/* increment standard deviation */
  scalartype logdetQ0;
public:
  ARk_t(){/*phi=phi_;MARGINAL=f_;*/}
  ARk_t(vectortype phi_){
    phi=phi_;
    k=phi.size();
    V0.resize(k,k);Q0.resize(k,k);
    M.resize(k,k);I.resize(k,k);
    /* build M-matrix */
    M.setZero();
    int d;
    for(int i=0;i<k;i++){
      for(int j=0;j<k;j++){
        d=abs(i-j);
        if(d!=0){
          M(i,d-1)+=phi[j];
        }
      }
    }
    I.setIdentity();
    gamma=((I-M).inverse())*matrixtype(phi);
    /* Increment sd */
    sigma=sqrt(scalartype(1)-(phi*gamma).sum());
    /* build V0 matrix */
    for(int i=0;i<k;i++){
      for(int j=0;j<k;j++){
        d=abs(i-j);
        if(d==0)
          V0(i,j)=scalartype(1); 
        else 
          V0(i,j)=gamma(d-1);
      }
    }
    /* build Q0 matrix */
    Q0=V0.inverse();
    /* log determinant */
    L0=Q0.llt().matrixL(); /* L0 L0' = Q0 */
    logdetQ0=scalartype(0);
    for(int i=0;i<k;i++)logdetQ0+=scalartype(2)*log(L0(i,i));
  }
  vectortype cov(int n){
    vectortype rho(n);
    for(int i=0;i<n;i++){
      if(i==0){rho(0)=scalartype(1);}
      else if(i<=k){rho(i)=gamma(i-1);}
      else { /* youle walker */
        scalartype tmp=0;
        for(int j=0;j<k;j++)tmp+=phi[j]*rho[i-1-j];
        rho(i)=tmp;
      }
    }
    return rho;
  }

  scalartype operator()(vectortype x){
    scalartype value=0;
    /* Initial distribution. For i = k,...,1 the recursions for
       solving L' y = u (1-based index notation) are:

       y(i) = L(i,i)^-1 * ( u(i) - L(i+1,i) * y(i+1) - ... - L(k,i) * y(k) )

       where u(i) ~ N(0,1).
    */
    scalartype mu, sd;
    int col;
    for(int i=0; (i<k) & (i<x.size()); i++){
      mu = scalartype(0);
      col = k-1-i; /* reversed index */
      for(int j=col+1; j<k; j++) mu -= L0(j,col) * x(k-1-j);
      mu /= L0(col, col);
      sd = scalartype(1) / L0(col, col);
      value -= dnorm(x[i], mu, sd, true);
    }
    scalartype tmp;
    for(int i=k;i<x.size();i++){
      tmp=scalartype(0);
      for(int j=0;j<k;j++){
        tmp+=phi[j]*x[i-1-j];
      }
      value-=dnorm(x[i],tmp,sigma,1);
    }
    return value;
  }
  arraytype jacobian(arraytype x){
    arraytype y(x.dim);
    y.setZero();
    for(int i=0;i<k;i++)
      for(int j=0;j<k;j++)
        y.col(i)=y.col(i)+Q0(i,j)*x.col(j);
    vectortype v(k+1);
    v(0)=scalartype(1);
    for(int i=1;i<=k;i++)v[i]=-phi[i-1];
    v=v/sigma;
    for(int i=k;i<x.cols();i++){
      for(int j1=0;j1<=k;j1++){
        for(int j2=0;j2<=k;j2++){
          y.col(i-j1)=y.col(i-j1)+v[j1]*v[j2]*x.col(i-j2);
        }
      }
    }
    return y;
  }
  int ndim(){return 1;}
  VARIANCE_NOT_YET_IMPLEMENTED
  /* Simulate - copy paste from evaluation operator */
  vectortype sqrt_cov_scale(vectortype u){
    vectortype x(u.size());
    scalartype mu, sd;
    int col;
    for (int i = 0; (i < k) & (i < x.size()); i++) {
      mu = scalartype(0);
      col = k - 1 - i; /* reversed index */
      for (int j = col + 1; j < k; j++)
        mu -= L0(j, col) * x(k - 1 - j);
      mu /= L0(col, col);
      sd = scalartype(1) / L0(col, col);
      x(i) = sd * u(i) + mu;
    }
    scalartype tmp;
    for (int i = k; i < x.size(); i++) {
      tmp = scalartype(0);
      for (int j = 0; j < k; j++) {
        tmp += phi[j] * x[i - 1 - j];
      }
      x(i) = sigma * u(i) + tmp;
    }
    return x;
  }
  SIMULATE_IMPLEMENTED_UNKNOWN_SIZE
};

template <class scalartype>
ARk_t<scalartype> ARk(vector<scalartype> phi){
  return ARk_t<scalartype>(phi);
}


template <class scalartype_>
class contAR2_t{
  TYPEDEFS(scalartype_);
private:
  typedef Matrix<scalartype,2,2> matrix2x2;
  typedef Matrix<scalartype,2,1> matrix2x1;
  typedef Matrix<scalartype,4,4> matrix4x4;
  typedef Matrix<scalartype,4,1> matrix4x1;
  scalartype shape,scale,c0,c1;
  vectortype grid;
  matrix2x2 A, V0, I;
  matrix4x4 B, iB; /* B=A %x% I + I %x% A  */
  matexp<scalartype,2> expA;
  matrix4x1 vecSigma,iBvecSigma;
  vector<MVNORM_t<scalartype> > neglogdmvnorm; /* Cache the 2-dim increments */
  vector<matrix2x2 > expAdt; /* Cache matrix exponential for grid increments */
public:
  contAR2_t(){};
  contAR2_t(vectortype grid_, scalartype shape_, scalartype scale_=1){
    shape=shape_;scale=scale_;grid=grid_;
    c0=scalartype(-1);c1=scalartype(2)*shape_;
    c0=c0/(scale*scale); c1=c1/scale;
    A << scalartype(0), scalartype(1), c0, c1;
    V0 << 1,0,0,-c0;
    I.setIdentity();
    B=kronecker(I,A)+kronecker(A,I);
    iB=B.inverse();
    expA=matexp<scalartype,2>(A);
    vecSigma << 0,0,0,scalartype(-2)*c1*V0(1,1);
    iBvecSigma=iB*vecSigma;
    /* cache increment distribution N(0,V(dt)) - one for each grid point */
    neglogdmvnorm.resize(grid.size());
    neglogdmvnorm[0]=MVNORM_t<scalartype>(V0);
    for(int i=1;i<grid.size();i++)neglogdmvnorm[i]=MVNORM_t<scalartype>(V(grid(i)-grid(i-1)));
    /* cache matrix exponential */
    expAdt.resize(grid.size());
    expAdt[0]=expA(scalartype(0));
    for(int i=1;i<grid.size();i++)expAdt[i]=expA(grid(i)-grid(i-1));
  }
  /* Simple formula for matrix exponential exp(B*t) */
  matrix4x4 expB(scalartype t){
    return kronecker(expA(t),expA(t));
  }
  /* Variance as fct. of time when started deterministic */
  matrix2x2 V(scalartype t){
    matrix4x1 tmp;
    tmp=expB(t)*iBvecSigma-iBvecSigma;
    matrix2x2 ans;
    for(int i=0;i<4;i++)ans(i)=tmp(i);
    return ans;
  }
  scalartype operator()(vectortype x,vectortype dx){
    matrix2x1 y, y0;
    scalartype ans;
    y0 << x(0), dx(0);
    ans = neglogdmvnorm[0](y0);
    for(int i=1;i<grid.size();i++){
      y0 << x(i-1), dx(i-1);
      y << x(i), dx(i);
      ans += neglogdmvnorm[i](y-expAdt[i]*y0);
    }
    return ans;
  }
  /* Experiment: Implementing matrix-vector multiply - Q*x.
     To be used when creating separable extensions... 
     think best to assume that input array has dimension (2,ntime,...) 
     so that we can easily extract x(t) and multiply Q*x(t) with a 2x2 matrix */
  scalartype operator()(vectortype x){ /* x.dim=[2,n] */
    vector<int> dim(2);
    dim << 2 , x.size()/2 ;
    array<scalartype> y(x,dim);
    y=y.transpose();
    return this->operator()(y.col(0),y.col(1));
  }
  arraytype matmult(matrix2x2 Q,arraytype x){
    arraytype y(x.dim);
    y.col(0) = Q(0,0)*x.col(0)+Q(0,1)*x.col(1); /* TODO: can we subassign like this in array class? Hack: we use "y.row" for that */
    y.col(1) = Q(1,0)*x.col(0)+Q(1,1)*x.col(1);
    return y;
  }
  arraytype jacobian(arraytype x){
    arraytype y(x.dim);
    y.setZero();
    arraytype tmp(y(0).dim);
    y.col(0) = neglogdmvnorm[0].jacobian(x.col(0)); /* Time zero contrib */
    for(int i=1;i<grid.size();i++){
      /* When taking derivative of .5*(x(i)-G*x(i-1))'*Q*(x(i)-G*x(i-1)) [where G=expAdt]
         we get contributions like: 
         x(i-1): -G'*Q*(x(i)-G*x(i-1))
         x(i):   Q*(x(i)-G*x(i-1))
      */
      tmp=neglogdmvnorm[i].jacobian( x.col(i) - matmult(expAdt[i],x.col(i-1)) );
      y.col(i)=y.col(i)+tmp;
      y.col(i-1)=y.col(i-1)-matmult(expAdt[i].transpose(), tmp );
    }
    return y;
  } 
  int ndim(){return 2;} /* Number of dimensions this structure occupies in total array */
  VARIANCE_NOT_YET_IMPLEMENTED
  SIMULATE_NOT_YET_IMPLEMENTED
};
template <class scalartype, class vectortype>
contAR2_t<scalartype> contAR2(vectortype grid_, scalartype shape_, scalartype scale_=1){
  return contAR2_t<scalartype>(grid_, shape_, scale_);
}
template <class scalartype>
contAR2_t<scalartype> contAR2(scalartype shape_, scalartype scale_=1){
  return contAR2_t<scalartype>(shape_, scale_);
}

template <class scalartype_>
class GMRF_t {
  TYPEDEFS(scalartype_);
private:
  Eigen::SparseMatrix<scalartype> Q;
  scalartype logdetQ;
  int sqdist(vectortype x, vectortype x_) {
    int ans = 0;
    int tmp;
    for(int i=0; i<x.size(); i++){
      tmp = CppAD::Integer(x[i]) - CppAD::Integer(x_[i]);
      ans += tmp * tmp;
    }
    return ans;
  }
public:
  GMRF_t(){}
  GMRF_t(Eigen::SparseMatrix<scalartype> Q_, int order_=1, bool normalize=true){
    setQ(Q_, order_, normalize);
  }
  GMRF_t(arraytype x, vectortype delta, int order_=1, bool normalize=true){
    int n = x.cols();
    typedef Eigen::Triplet<scalartype> T;
    std::vector<T> tripletList;
    for(int i=0; i<n; i++){
      for(int j=0; j<n; j++){
        if (sqdist(x.col(i),x.col(j)) == 1){
          tripletList.push_back(T(i, j, scalartype(-1)));
          tripletList.push_back(T(i, i, scalartype( 1)));
        }
      }
    }
    for(int i=0; i<n; i++) {
      tripletList.push_back(T(i, i, delta[i]));
    }
    Eigen::SparseMatrix<scalartype> Q_(n, n);
    Q_.setFromTriplets(tripletList.begin(), tripletList.end());
    setQ(Q_, order_, normalize);
  }
  void setQ(Eigen::SparseMatrix<scalartype> Q_, int order=1, bool normalize=true){
    Q = Q_;
    if (normalize) {
#ifndef TMBAD_FRAMEWORK
      Eigen::SimplicialLDLT< Eigen::SparseMatrix<scalartype> > ldl(Q);
      vectortype D = ldl.vectorD();
      logdetQ = (log(D)).sum();
#else
      logdetQ = newton::log_determinant(Q);
#endif
    } else {
      logdetQ = 0;
    }
    /* Q^order */
    for(int i=1; i<order; i++){
      Q = Q * Q_;
    }
    logdetQ = scalartype(order) * logdetQ;
  }
  /* Quadratic form: x'*Q^order*x */
  scalartype Quadform(vectortype x){
    return (x * (Q * x.matrix()).array()).sum();
  }
  scalartype operator()(vectortype x){
    return -scalartype(.5) * logdetQ + scalartype(.5) * Quadform(x) + x.size() * scalartype(log(sqrt(2.0 * M_PI)));
  }
  /* jacobian */
  arraytype jacobian(arraytype x){
    arraytype y(x.dim);
    matrixtype m(x.size()/x.cols(),x.cols());
    for(int i=0; i<x.size(); i++) m(i) = x[i];
    matrixtype mQ = m * Q;
    for(int i=0; i<x.size(); i++) y[i] = mQ(i);
    return y;
  }
  int ndim() { return 1; }
  vectortype variance() {
    int n = Q.rows();
    vectortype ans(n);
    matrixtype C = invertSparseMatrix(Q);
    for(int i=0; i<n; i++) ans[i] = C(i,i);
    return ans;
  }
  /* Simulation */
  Eigen::SparseMatrix<scalartype> L;
  Eigen::PermutationMatrix<Dynamic,Dynamic> Pinv;
  vectortype sqrt_cov_scale(vectortype u) {
    if(L.rows() == 0) {
      Eigen::SimplicialLLT<Eigen::SparseMatrix<scalartype> > solver(Q);
      L = solver.matrixL();
      Pinv = solver.permutationPinv();
    }
    // L*L^T = P*Q*Pinv => Q^-1 = A*A^T where A:=P^-1*L^T^-1
    matrixtype x = L.transpose().template triangularView<Eigen::Upper>().solve(u.matrix());
    x = Pinv * x;
    return x.vec();
  }
  SIMULATE_IMPLEMENTED_KNOWN_SIZE(Q.rows())
};
template <class scalartype>
GMRF_t<scalartype> GMRF(Eigen::SparseMatrix<scalartype> Q, int order, bool normalize=true) {
  return GMRF_t<scalartype>(Q, order, normalize);
}
template <class scalartype, class arraytype >
GMRF_t<scalartype> GMRF(arraytype x, vector<scalartype> delta, int order=1, bool normalize=true) {
  return GMRF_t<scalartype>(x, delta, order, normalize);
}
template <class scalartype, class arraytype >
GMRF_t<scalartype> GMRF(arraytype x, scalartype delta, int order=1, bool normalize=true) {
  vector<scalartype> d(x.cols());
  for(int i=0; i<d.size(); i++) d[i] = delta;
  return GMRF_t<scalartype>(x, d, order, normalize);
}
template <class scalartype>
GMRF_t<scalartype> GMRF(Eigen::SparseMatrix<scalartype> Q, bool normalize = true) {
  return GMRF_t<scalartype>(Q, 1, normalize);
}

template <class distribution>
class SCALE_t{
  TYPEDEFS(typename distribution::scalartype);
private:
  distribution f;
  scalartype scale;
public:
  SCALE_t(){}
  SCALE_t(distribution f_, scalartype scale_){scale=scale_;f=f_;}
  scalartype operator()(arraytype x){
    scalartype ans=f(x/scale);
    ans+=x.size()*log(scale);
    return ans;
  }
  scalartype operator()(vectortype x){
    scalartype ans=f(x/scale);
    ans+=x.size()*log(scale);
    return ans;
  }
  arraytype jacobian(arraytype x){
    return f.jacobian(x/scale)/scale;    
  }
  int ndim(){return f.ndim();}
  vectortype variance(){
    return (scale*scale)*f.variance();
  }
  vectortype sqrt_cov_scale(vectortype u){
    return scale * f.sqrt_cov_scale(u);
  }
  SIMULATE_IMPLEMENTED_UNKNOWN_SIZE
};
template <class scalartype, class distribution>
SCALE_t<distribution> SCALE(distribution f_, scalartype scale_){
  return SCALE_t<distribution>(f_,scale_);
}

template <class distribution>
class VECSCALE_t{
  TYPEDEFS(typename distribution::scalartype);
private:
  distribution f;
  vectortype scale;
public:
  VECSCALE_t(){}
  VECSCALE_t(distribution f_, vectortype scale_){scale=scale_;f=f_;}
  scalartype operator()(arraytype x){
    // assert that x.size()==scale.size()
    scalartype ans=f(x/scale);
    ans+=(log(scale)).sum();
    return ans;
  }
  scalartype operator()(vectortype x){
    // assert that x.size()==scale.size()
    scalartype ans=f(x/scale);
    ans+=(log(scale)).sum();
    return ans;
  }
  arraytype jacobian(arraytype x){
    // assert that x.rows()==scale.size()
    arraytype y(x);
    for(int i=0;i<y.cols();i++)y.col(i)=y.col(i)/scale[i];
    y=f.jacobian(y);
    for(int i=0;i<y.cols();i++)y.col(i)=y.col(i)/scale[i];
    return y;
  }
  int ndim(){return f.ndim();}
  VARIANCE_NOT_YET_IMPLEMENTED
  vectortype sqrt_cov_scale(vectortype u){
    return scale * f.sqrt_cov_scale(u);
  }
  SIMULATE_IMPLEMENTED_UNKNOWN_SIZE
};
template <class vectortype, class distribution>
VECSCALE_t<distribution> VECSCALE(distribution f_, vectortype scale_){
  return VECSCALE_t<distribution>(f_,scale_);
}


//template <class scalartype, class vectortype, class arraytype, class distribution1, class distribution2>
template <class distribution1, class distribution2>
class SEPARABLE_t{
  TYPEDEFS(typename distribution1::scalartype);
private:
  distribution1 f;
  distribution2 g;
public:
  SEPARABLE_t(){}
  SEPARABLE_t(distribution1 f_, distribution2 g_){f=f_;g=g_;}
  /*
    Example: x.dim=[n1,n2,n3].
    Apply f on outer dimension (n3) and rotate:
    [n3,n1,n2]
    Apply g on new outer dimension (n2) and rotate back:
    [n1,n2,n3]
   */
  arraytype jacobian(arraytype x){
    int n=f.ndim();
    x=f.jacobian(x);
    x=x.rotate(n);
    x=g.jacobian(x);
    x=x.rotate(-n);
    return x;
  }
  /* Create zero vector corresponding to the last n dimensions of dimension-vector d */
  arraytype zeroVector(vector<int> d, int n){
    int m=1;
    vector<int> revd=d.reverse();
    vector<int> revnewdim(n);
    for(int i=0;i<n;i++){m=m*revd[i];revnewdim[i]=revd[i];}
    vectortype x(m);
    x.setZero();
    return arraytype(x,revnewdim.reverse());
  }
  scalartype operator()(arraytype x){
    if(this->ndim() != x.dim.size())std::cout << "Wrong dimension in SEPARABLE_t\n";
    /* Calculate quadform */
    arraytype y(x.dim);
    y=jacobian(x);
    y=x*y; /* pointwise */
    scalartype q=scalartype(.5)*(y.sum()); 
    /* Add normalizing constant */
    int n=f.ndim();
    arraytype zf=zeroVector(x.dim,n);
    q+=f(zf)*(scalartype(x.size())/scalartype(zf.size()));
    x=x.rotate(n);
    int m=g.ndim();
    arraytype zg=zeroVector(x.dim,m);
    q+=g(zg)*(scalartype(x.size())/scalartype(zg.size()));
    q-=log(sqrt(2.0*M_PI))*(zf.size()*zg.size());
    /* done */
    return q;
  }
  int ndim(){return f.ndim()+g.ndim();}
  VARIANCE_NOT_YET_IMPLEMENTED

  /* For parallel accumulation:
     ==========================
     Copied operator() above and added extra argument "i" to divide the accumulation
     in chunks. The evaluation of
         operator()(x)
     is equivalent to summing up
         operator()(x,i)
     with i running through the _outer_dimension_ of x.
  */
  scalartype operator()(arraytype x, int i){
    if(this->ndim() != x.dim.size())std::cout << "Wrong dimension in SEPARABLE_t\n";
    /* Calculate quadform */
    arraytype y(x.dim);
    y=jacobian(x);
    y=x*y; /* pointwise */
    scalartype q=scalartype(.5)*(y.col(i).sum()); 
    /* Add normalizing constant */
    if(i==0){
      int n=f.ndim();
      arraytype zf=zeroVector(x.dim,n);
      q+=f(zf)*(scalartype(x.size())/scalartype(zf.size()));
      x=x.rotate(n);
      int m=g.ndim();
      arraytype zg=zeroVector(x.dim,m);
      q+=g(zg)*(scalartype(x.size())/scalartype(zg.size()));
      q-=log(sqrt(2.0*M_PI))*(zf.size()*zg.size());
    }
    /* done */
    return q;
  }

  arraytype sqrt_cov_scale(arraytype u) {
    vector<int> u_dim = u.dim;
    vector<int> f_dim = u_dim.tail(f.ndim());
    vector<int> g_dim = u_dim.head(g.ndim());
    int f_size = f_dim.prod();
    int g_size = g_dim.prod();
    // Collapse f dimension to a single dimension:
    vector<int> new_dim(g.ndim() + 1);
    new_dim << g_dim, f_size;
    u.setdim(new_dim);
    for(int i=0; i<f_size; i++) {
      u.col(i) = g.sqrt_cov_scale( u.col(i) );
    }
    u = u.rotate(1); // u.dim = (f_size, g_dim)
    // Collapse g dimension to a single dimension:
    new_dim.resize(f.ndim() + 1);
    new_dim << f_dim, g_size;
    u.setdim(new_dim);
    for(int i=0; i<g_size; i++) {
      u.col(i) = f.sqrt_cov_scale( u.col(i) );
    }
    u = u.rotate(1);
    u.setdim(u_dim);
    return u;
  }

  void simulate(arraytype &x) {
    rnorm_fill(x);
    x = sqrt_cov_scale(x);
    x = zero_derivatives(x);
  }

};

template <class distribution1, class distribution2>
SEPARABLE_t<distribution1,distribution2> SEPARABLE(distribution1 f_, distribution2 g_){
  return SEPARABLE_t<distribution1,distribution2>(f_,g_);
}

template <class distribution>
class PROJ_t{
  TYPEDEFS(typename distribution::scalartype);
private:
  distribution f;
  bool initialized;
public:
  vector<int> proj;
  vector<int> cproj; /* complementary proj _sorted_ */
  int n,nA,nB;
  matrixtype Q;  /* Full precision */
  MVNORM_t<scalartype> dmvnorm; /* mean zero gaussian with covariance Q_BB */
  PROJ_t(){}
  PROJ_t(distribution f_, vector<int> proj_){
    f=f_;
    proj=proj_;
    initialized=false;
  }
  void initialize(int n_){
    if(!initialized){
      n=n_;
      nA=proj.size();
      nB=n-nA;
      cproj.resize(nB);
      vector<int> mark(n);
      mark.setZero();
      for(int i=0;i<nA;i++)mark[proj[i]]=1;
      int k=0;
      for(int i=0;i<n;i++)if(!mark[i])cproj[k++]=i;
      // Full precision
      //matrixtype I(n,n);
      //I.setIdentity();
      //vectortype v(I);

      k=0;
      vectortype v(n*n);
      for(int i=0;i<n;i++)
        for(int j=0;j<n;j++)
          v(k++)=scalartype(i==j);

      vector<int> dim(2);
      dim << n,n;
      arraytype a(v,dim);
      a=f.jacobian(a);
      Q.resize(n,n);
      for(int i=0;i<n*n;i++)Q(i)=a[i];
      // Get Q_BB
      matrixtype QBB(nB,nB);
      for(int i=0;i<nB;i++)
        for(int j=0;j<nB;j++)
          QBB(i,j)=Q(cproj[i],cproj[j]);
      dmvnorm=MVNORM_t<scalartype>(QBB);
    }
    initialized=true;
  }
  vectortype projB(vectortype x){
    vectortype y(nB);
    for(int i=0;i<nB;i++)y[i]=x[cproj[i]];
    return y;
  }
  vectortype setZeroB(vectortype x){
    for(int i=0;i<nB;i++)x[cproj[i]]=scalartype(0);
    return x;    
  }
  scalartype operator()(vectortype x){
    initialize(x.size());
    x=setZeroB(x);
    vector<int> dim(1);
    dim << x.size();
    arraytype xa(x,dim);
    vectortype y=projB(f.jacobian(xa));
    // f(x_A,0) - dmvnorm(y_B) + 2*dmvnorm(0)
    return f(xa) - dmvnorm(y) + 2*dmvnorm(y*scalartype(0));
  }
  /* array versions  */
  arraytype projB(arraytype x){
    vectortype z((x.size()/n)*nB);
    vector<int> dim(x.dim);
    dim[dim.size()-1]=nB;
    arraytype y(z,dim);
    for(int i=0;i<nB;i++)y.col(i)=x.col(cproj[i]);
    return y;
  }
  arraytype setZeroB(arraytype x){
    for(int i=0;i<nB;i++)x.col(cproj[i])=x.col(0)*scalartype(0);
    return x;    
  }
  arraytype jacobian(arraytype x){
    initialize(x.dim[x.dim.size()-1]);
    arraytype xa=setZeroB(x);
    arraytype y=projB(f.jacobian(xa));
    // WRONG: ----> return f.jacobian(xa) - dmvnorm.jacobian(y);
    // y=P*Q*Z*x  so should be  (P*Q*Z)' * dmvnorm.jacobian(y).
    // Note: only P is not symmetric.
    arraytype tmp=f.jacobian(xa);
    arraytype tmp0=tmp*scalartype(0);
    arraytype tmp2=dmvnorm.jacobian(y);
    // apply P'
    for(int i=0;i<nB;i++){ 
      tmp0.col(cproj[i])=tmp2.col(i);
    }
    // apply Q'(=Q)
    tmp0=f.jacobian(tmp0);
    // apply Z'(=Z)
    tmp0=setZeroB(tmp0);
    // Done:
    return tmp-tmp0;
  }
  int ndim(){return f.ndim();}
  VARIANCE_NOT_YET_IMPLEMENTED
  SIMULATE_NOT_YET_IMPLEMENTED
};

template <class distribution>
PROJ_t<distribution> PROJ(distribution f_, vector<int> i){
  return PROJ_t<distribution>(f_,i);
}


#undef TYPEDEFS

} // End namespace