lsst::afw::math::KdTree< T > Class Template Reference

The data for GaussianProcess is stored in a KD tree to facilitate nearest-neighbor searches. More...

#include <GaussianProcess.h>

Public Member Functions
	KdTree (const KdTree &)=delete
KdTree &	operator= (const KdTree &)=delete
	KdTree (KdTree &&)=delete
KdTree &	operator= (KdTree &&)=delete
	KdTree ()=default
void	Initialize (ndarray::Array< T, 2, 2 > const &dt)
	Build a KD Tree to store the data for GaussianProcess.
void	findNeighbors (ndarray::Array< int, 1, 1 > neighdex, ndarray::Array< double, 1, 1 > dd, ndarray::Array< const T, 1, 1 > const &v, int n_nn) const
	Find the nearest neighbors of a point.
T	getData (int ipt, int idim) const
	Return one element of one node on the tree.
ndarray::Array< T, 1, 1 >	getData (int ipt) const
	Return an entire node from the tree.
void	addPoint (ndarray::Array< const T, 1, 1 > const &v)
	Add a point to the tree.
void	removePoint (int dex)
	Remove a point from the tree.
int	getNPoints () const
	return the number of data points stored in the tree
void	getTreeNode (ndarray::Array< int, 1, 1 > const &v, int dex) const
	Return the _tree information for a given data point.

Detailed Description

template<typename T>
class lsst::afw::math::KdTree< T >

The data for GaussianProcess is stored in a KD tree to facilitate nearest-neighbor searches.

Note: I have removed the ability to arbitrarily specify a distance function. The KD Tree nearest neighbor search algorithm only makes sense in the case of Euclidean distances, so I have forced KdTree to use Euclidean distances.

Definition at line 224 of file GaussianProcess.h.

Constructor & Destructor Documentation

KdTree() [1/3]

template<typename T >

lsst::afw::math::KdTree< T >::KdTree ( const KdTree< T > & )

delete

KdTree() [2/3]

template<typename T >

lsst::afw::math::KdTree< T >::KdTree ( KdTree< T > && )

delete

KdTree() [3/3]

template<typename T >

lsst::afw::math::KdTree< T >::KdTree ( )

default

Member Function Documentation

addPoint()

template<typename T >

void lsst::afw::math::KdTree< T >::addPoint

(

ndarray::Array< const T, 1, 1 > const &

v

)

Add a point to the tree.

Allot more space in _tree and data if needed.

Parameters

[in]

v

the point you are adding to the tree

Exceptions

pex::exceptions::RuntimeError

if the branch ending in the new point is not properly constructed

Definition at line 210 of file GaussianProcess.cc.

                                                             {
    if (v.getNumElements() != static_cast<ndarray::Size>(_dimensions)) {
        throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
                          "You are trying to add a point of the incorrect dimensionality to KdTree\n");
    }
 
    int i, j, node, dim, dir;
 
    node = _findNode(v);
    dim = _tree[node][DIMENSION] + 1;
    if (dim == _dimensions) dim = 0;
 
    if (_npts == _room) {
        ndarray::Array<T, 2, 2> dbuff = allocate(ndarray::makeVector(_npts, _dimensions));
        ndarray::Array<int, 2, 2> tbuff = allocate(ndarray::makeVector(_npts, 4));
 
        dbuff.deep() = _data;
        tbuff.deep() = _tree;
 
        _room += _roomStep;
 
        _tree = allocate(ndarray::makeVector(_room, 4));
        _data = allocate(ndarray::makeVector(_room, _dimensions));
 
        for (i = 0; i < _npts; i++) {
            for (j = 0; j < _dimensions; j++) _data[i][j] = dbuff[i][j];
            for (j = 0; j < 4; j++) _tree[i][j] = tbuff[i][j];
        }
    }
 
    _tree[_npts][DIMENSION] = dim;
    _tree[_npts][PARENT] = node;
    i = _tree[node][DIMENSION];
 
    if (_data[node][i] > v[i]) {
        if (_tree[node][LT] >= 0) {
            throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
                              "Trying to add to KdTree in a node that is already occupied\n");
        }
        _tree[node][LT] = _npts;
        dir = LT;
    } else {
        if (_tree[node][GEQ] >= 0) {
            throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
                              "Trying to add to KdTree in a node that is already occupied\n");
        }
        _tree[node][GEQ] = _npts;
        dir = GEQ;
    }
    _tree[_npts][LT] = -1;
    _tree[_npts][GEQ] = -1;
    for (i = 0; i < _dimensions; i++) {
        _data[_npts][i] = v[i];
    }
 
    _npts++;
 
    i = _walkUpTree(_tree[_npts - 1][PARENT], dir, _npts - 1);
    if (i != _masterParent) {
        throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError, "Adding to KdTree failed\n");
    }
}

findNeighbors()

template<typename T >

void lsst::afw::math::KdTree< T >::findNeighbors	(	ndarray::Array< int, 1, 1 >	neighdex,
		ndarray::Array< double, 1, 1 >	dd,
		ndarray::Array< const T, 1, 1 > const &	v,
		int	n_nn ) const

Find the nearest neighbors of a point.

Parameters

[out]	neighdex	this is where the indices of the nearest neighbor points will be stored
[out]	dd	this is where the distances to the nearest neighbors will be stored
[in]	v	the point whose neighbors you want to find
[in]	n_nn	the number of nearest neighbors you want to find

neighbors will be returned in ascending order of distance

note that distance is forced to be the Euclidean distance

Definition at line 135 of file GaussianProcess.cc.

                                                                                    {
    if (n_nn > _npts) {
        throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
                          "Asked for more neighbors than kd tree contains\n");
    }
 
    if (n_nn <= 0) {
        throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
                          "Asked for zero or a negative number of neighbors\n");
    }
 
    if (neighdex.getNumElements() != static_cast<ndarray::Size>(n_nn)) {
        throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
                          "Size of neighdex does not equal n_nn in KdTree.findNeighbors\n");
    }
 
    if (dd.getNumElements() != static_cast<ndarray::Size>(n_nn)) {
        throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
                          "Size of dd does not equal n_nn in KdTree.findNeighbors\n");
    }
 
    int i, start;
 
    _neighborCandidates = allocate(ndarray::makeVector(n_nn));
    _neighborDistances = allocate(ndarray::makeVector(n_nn));
    _neighborsFound = 0;
    _neighborsWanted = n_nn;
 
    for (i = 0; i < n_nn; i++) _neighborDistances[i] = -1.0;
 
    start = _findNode(v);
 
    _neighborDistances[0] = _distance(v, _data[start]);
    _neighborCandidates[0] = start;
    _neighborsFound = 1;
 
    for (i = 1; i < 4; i++) {
        if (_tree[start][i] >= 0) {
            _lookForNeighbors(v, _tree[start][i], start);
        }
    }
 
    for (i = 0; i < n_nn; i++) {
        neighdex[i] = _neighborCandidates[i];
        dd[i] = _neighborDistances[i];
    }
}

getData() [1/2]

template<typename T >

ndarray::Array< T, 1, 1 > lsst::afw::math::KdTree< T >::getData

(

int

ipt

)

const

Return an entire node from the tree.

Parameters

[in]

ipt

the index of the node to return

I currently have this as a return-by-value method. When I tried it as a return-by-reference, the compiler gave me

warning: returning reference to local temporary object

Based on my reading of Stack Overflow, this is because ndarray was implicitly creating a new ndarray::Array<T,1,1> object and passing a reference thereto. It is unclear to me whether or not this object would be destroyed once the call to getData was complete.

The code still compiled, ran, and passed the unit tests, but the above behavior seemed to me like it could be dangerous (and, because ndarray was still creating a new object, it did not seem like we were saving any time), so I reverted to return-by-value.

Definition at line 200 of file GaussianProcess.cc.

                                                      {
    if (ipt >= _npts || ipt < 0) {
        throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
                          "Asked for point that does not exist in KdTree\n");
    }
 
    return _data[ipt];
}

getData() [2/2]

template<typename T >

T lsst::afw::math::KdTree< T >::getData	(	int	ipt,
		int	idim ) const

Return one element of one node on the tree.

Parameters

[in]	ipt	the index of the node to return
[in]	idim	the index of the dimension to return

Definition at line 185 of file GaussianProcess.cc.

                                            {
    if (ipt >= _npts || ipt < 0) {
        throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
                          "Asked for point that does not exist in KdTree\n");
    }
 
    if (idim >= _dimensions || idim < 0) {
        throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
                          "KdTree does not have points of that many dimensions\n");
    }
 
    return _data[ipt][idim];
}

getNPoints()

template<typename T >

int lsst::afw::math::KdTree< T >::getNPoints

(

)

const

return the number of data points stored in the tree

Definition at line 274 of file GaussianProcess.cc.

                                {
    return _npts;
}

getTreeNode()

template<typename T >

void lsst::afw::math::KdTree< T >::getTreeNode	(	ndarray::Array< int, 1, 1 > const &	v,
		int	dex ) const

Return the _tree information for a given data point.

Parameters

[out]	v	the array in which to store the entry from _tree
[in]	dex	the index of the node whose information you are requesting

Definition at line 279 of file GaussianProcess.cc.

                                                                           {
    if (dex < 0 || dex >= _npts) {
        throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
                          "Asked for tree information on point that does not exist\n");
    }
 
    if (v.getNumElements() != 4) {
        throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
                          "Need to pass a 4-element ndarray into KdTree.getTreeNode()\n");
    }
 
    v[0] = _tree[dex][DIMENSION];
    v[1] = _tree[dex][LT];
    v[2] = _tree[dex][GEQ];
    v[3] = _tree[dex][PARENT];
}

Initialize()

template<typename T >

void lsst::afw::math::KdTree< T >::Initialize

(

ndarray::Array< T, 2, 2 > const &

dt

)

Build a KD Tree to store the data for GaussianProcess.

Parameters

[in]

dt

an array, the rows of which are the data points (dt[i][j] is the jth component of the ith data point)

Exceptions

pex::exceptions::RuntimeError

if the tree is not properly constructed

Definition at line 104 of file GaussianProcess.cc.

                                                          {
    int i;
 
    _npts = dt.template getSize<0>();
    _dimensions = dt.template getSize<1>();
 
    // a buffer to use when first building the tree
    _inn = allocate(ndarray::makeVector(_npts));
 
    _roomStep = 5000;
    _room = _npts;
 
    _data = allocate(ndarray::makeVector(_room, _dimensions));
 
    _data.deep() = dt;
 
    _tree = allocate(ndarray::makeVector(_room, 4));
 
    for (i = 0; i < _npts; i++) {
        _inn[i] = i;
    }
 
    _organize(_inn, _npts, -1, -1);
 
    i = _testTree();
    if (i == 0) {
        throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError, "Failed to properly initialize KdTree\n");
    }
}

operator=() [1/2]

template<typename T >

KdTree & lsst::afw::math::KdTree< T >::operator= ( const KdTree< T > & )

delete

operator=() [2/2]

template<typename T >

KdTree & lsst::afw::math::KdTree< T >::operator= ( KdTree< T > && )

delete

removePoint()

template<typename T >

void lsst::afw::math::KdTree< T >::removePoint

(

int

dex

)

Remove a point from the tree.

Reorganize what remains so that the tree remains self-consistent

Parameters

[in]

dex

the index of the point you want to remove from the tree

Exceptions

pex::exceptions::RuntimeError

if the entire tree is not poperly constructed after the point has been removed

Definition at line 582 of file GaussianProcess.cc.

                                      {
    if (target < 0 || target >= _npts) {
        throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
                          "You are trying to remove a point that doesn't exist from KdTree\n");
    }
 
    if (_npts == 1) {
        throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
                          "There is only one point left in this KdTree.  You cannot call removePoint\n");
    }
 
    int nl, nr, i, j, k, side;
    int root;
 
    nl = 0;
    nr = 0;
    if (_tree[target][LT] >= 0) {
        nl++;
        _count(_tree[target][LT], &nl);
    }
 
    if (_tree[target][GEQ] >= 0) {
        nr++;
        _count(_tree[target][GEQ], &nr);
    }
 
    if (nl == 0 && nr == 0) {
        k = _tree[target][PARENT];
 
        if (_tree[k][LT] == target)
            _tree[k][LT] = -1;
        else if (_tree[k][GEQ] == target)
            _tree[k][GEQ] = -1;
 
    }  // if target is terminal
    else if ((nl == 0 && nr > 0) || (nr == 0 && nl > 0)) {
        if (nl == 0)
            side = GEQ;
        else
            side = LT;
 
        k = _tree[target][PARENT];
        if (k >= 0) {
            if (_tree[k][LT] == target) {
                _tree[k][LT] = _tree[target][side];
                _tree[_tree[k][LT]][PARENT] = k;
 
            } else {
                _tree[k][GEQ] = _tree[target][side];
                _tree[_tree[k][GEQ]][PARENT] = k;
            }
        } else {
            _masterParent = _tree[target][side];
            _tree[_tree[target][side]][PARENT] = -1;
        }
    }  // if only one side is populated
    else {
        if (nl > nr)
            side = LT;
        else
            side = GEQ;
 
        k = _tree[target][PARENT];
        if (k < 0) {
            _masterParent = _tree[target][side];
            _tree[_masterParent][PARENT] = -1;
        } else {
            if (_tree[k][LT] == target) {
                _tree[k][LT] = _tree[target][side];
                _tree[_tree[k][LT]][PARENT] = k;
            } else {
                _tree[k][GEQ] = _tree[target][side];
                _tree[_tree[k][GEQ]][PARENT] = k;
            }
        }
 
        root = _tree[target][3 - side];
 
        _descend(root);
    }  // if both sides are populated
 
    if (target < _npts - 1) {
        for (i = target + 1; i < _npts; i++) {
            for (j = 0; j < 4; j++) _tree[i - 1][j] = _tree[i][j];
 
            for (j = 0; j < _dimensions; j++) _data[i - 1][j] = _data[i][j];
        }
 
        for (i = 0; i < _npts; i++) {
            for (j = 1; j < 4; j++) {
                if (_tree[i][j] > target) _tree[i][j]--;
            }
        }
 
        if (_masterParent > target) _masterParent--;
    }
    _npts--;
 
    i = _testTree();
    if (i == 0) {
        throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError, "Subtracting from KdTree failed\n");
    }
}

---

The documentation for this class was generated from the following files:

• /j/snowflake/release/lsstsw/stack/lsst-scipipe-9.0.0/Linux64/afw/g9726552aa6+10f999ec6a/include/lsst/afw/math/GaussianProcess.h
• /j/snowflake/release/lsstsw/stack/lsst-scipipe-9.0.0/Linux64/afw/g9726552aa6+10f999ec6a/src/math/GaussianProcess.cc