generate matrix for NNGP

[alashworth]

Issue by LuZhangstat Thursday Feb 02, 2017 at 09:11 GMT Originally opened as https://github.com/stan-dev/stan/issues/2222

---

Summary:

A function for generating sparse distance matrix and the index of its' nonzero entities to be used in the calculation of the log likelihood of Nearest-Neighbor Gaussian Process.

`#ifndef GET_NN_DATA_HPP

define GET_NN_DATA_HPP

include

include

include

include <Eigen/Dense>

include "get_dist_v.hpp"

include "getNNIndx.hpp"

include <stan/math/prim/mat.hpp>

using std::vector; using Eigen::Matrix; using Eigen::Dynamic; using Eigen::RowVectorXd; template <typename T, int neighborNum> void get_NN_data(const Matrix<T, Dynamic, Dynamic>& sorted_coords, vector& NN_indM, vector& NN_distv, vector& NN_dist_M){

/**
 *  Generate sparse distance matrix and the index of
 *  its' nonzero entities to be used in the calculation of
 *  the log likelihood of Nearest-Neighbor Gaussian Process.
 *  
 *  The workspace is allocated base on the number of Nearest 
 *  Neighbors and the number of locations
 *
 *  neighborNum:    int     number of nearest neighbors
 *  sorted_coords:  matrix  coordinates with sorted first
 *  NN_indM:        vector  index of nearest neighbors
 *  NN_distv:       vector  distances among obs and 
 *  its' neighbors
 *  NN_dist_M:      vector  distance matrix of neighbors
 *
 *
 **/

int l1 = static_cast<int>
(static_cast<double>(neighborNum) / 2 * (neighborNum - 1)+
 (sorted_coords.rows() - neighborNum) * neighborNum);

NN_indM.resize(l1);
NN_distv.resize(l1);

int l2 = static_cast<int>
(static_cast<double>(neighborNum - 2) * (neighborNum - 1) *
 (neighborNum) / 6 + (sorted_coords.rows() - neighborNum) *
 static_cast<double>(neighborNum ) / 2 * (neighborNum - 1));

NN_dist_M.resize(l2);

int iNN, iNNIndx, iNNMIndx;

for (int i = 0; i < sorted_coords.rows(); i++){

    vector<double> dist_v(i);

    dist_v = get_dist_v<double>(sorted_coords, i);

    vector<int> sortind = stan::math::sort_indices_asc(dist_v);
    getNNIndx(i, neighborNum, iNNIndx, iNN, iNNMIndx);

    for (int j = 0; j < iNN; j++){
        NN_indM[j + iNNIndx] = sortind[j] - 1;
        NN_distv[j + iNNIndx] = dist_v[sortind[j] - 1];
    }

    int k = 0;

    for (int j = 0; j < (iNN - 1); j++){
        for (int h = j + 1; h < iNN; h++){
            RowVectorXd loc1 = sorted_coords.row(NN_indM[j + iNNIndx]);
            RowVectorXd loc2 = sorted_coords.row(NN_indM[h + iNNIndx]);
            NN_dist_M[iNNMIndx + k] = stan::math::squared_distance(loc1, loc2);
            k++;
        }
    }
}
return;

}

endif

`

steps:

1. Allocate the workspace for NN_indM, NN_distv and NN_dist_M
2. In the for-loop, _get_distv saves the distance vector between ith obs and all obs with order lower than i, getNNIndx find the index of the starting point for saving data.

Here is the test code

vector<int> NN_indM(0); vector<double> NN_distv(0); vector<double> NN_dist_M(0); get_NN_data<double, 8>(sorted_coords, NN_indM, NN_distv, NN_dist_M);

Results:

Can get results in C++. No sure how to test it with stan. Haven't done the test unit.

Current Version:

v2.14.0