generate matrix for NNGP

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_dist_v 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

Start by locating the existing get_dist_v and getNNIndx implementations and the Stan-facing entry point for this proposed function. Reproduce the provided C++ example, then add a unit test covering NN_indM, NN_distv, and NN_dist_M and confirm the function can be tested through Stan.