lcar.cpp

using Eigen::SparseMatrix;

template<class Type> 
SparseMatrix<Type> lcar_prior(size_t I, Type rho) {
    // LCAR prior is:
    //  pW + (1-p)I
    //    where W[j,k] = 
    //      n if j==k
    //      -1 if j~k
    //      0 otherwise
    // initialize sparse matrix
    SparseMatrix<Type> K(I,I); 
    // initialize diagonal to (1 - rho)
    for (size_t i = 0; i < I; i++)
        K.insert(i,i) = Type(1.) - rho;
    // set off-diagonal to -1*rho if neighbors
    // and add (count of neighbors)*rho to the diagonal
    for (size_t i1 = 0; i1 < I; i1++) {
        for (size_t i2 = 0; i2 < I; i2++) {
            if (abs(int(i1) - int(i2)) == 1) { 
                K.insert(i1,i2) = Type(-1.) * rho;
                K.coeffRef(i1,i1) += rho;
            }
        }
    }
    return K; 
}