LSST Applications  21.0.0-172-gfb10e10a+18fedfabac,22.0.0+297cba6710,22.0.0+80564b0ff1,22.0.0+8d77f4f51a,22.0.0+a28f4c53b1,22.0.0+dcf3732eb2,22.0.1-1-g7d6de66+2a20fdde0d,22.0.1-1-g8e32f31+297cba6710,22.0.1-1-geca5380+7fa3b7d9b6,22.0.1-12-g44dc1dc+2a20fdde0d,22.0.1-15-g6a90155+515f58c32b,22.0.1-16-g9282f48+790f5f2caa,22.0.1-2-g92698f7+dcf3732eb2,22.0.1-2-ga9b0f51+7fa3b7d9b6,22.0.1-2-gd1925c9+bf4f0e694f,22.0.1-24-g1ad7a390+a9625a72a8,22.0.1-25-g5bf6245+3ad8ecd50b,22.0.1-25-gb120d7b+8b5510f75f,22.0.1-27-g97737f7+2a20fdde0d,22.0.1-32-gf62ce7b1+aa4237961e,22.0.1-4-g0b3f228+2a20fdde0d,22.0.1-4-g243d05b+871c1b8305,22.0.1-4-g3a563be+32dcf1063f,22.0.1-4-g44f2e3d+9e4ab0f4fa,22.0.1-42-gca6935d93+ba5e5ca3eb,22.0.1-5-g15c806e+85460ae5f3,22.0.1-5-g58711c4+611d128589,22.0.1-5-g75bb458+99c117b92f,22.0.1-6-g1c63a23+7fa3b7d9b6,22.0.1-6-g50866e6+84ff5a128b,22.0.1-6-g8d3140d+720564cf76,22.0.1-6-gd805d02+cc5644f571,22.0.1-8-ge5750ce+85460ae5f3,master-g6e05de7fdc+babf819c66,master-g99da0e417a+8d77f4f51a,w.2021.48
LSST Data Management Base Package
Public Member Functions | List of all members
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
 
KdTreeoperator= (const KdTree &)=delete
 
 KdTree (KdTree &&)=delete
 
KdTreeoperator= (KdTree &&)=delete
 
 KdTree ()=default
 
void Initialize (ndarray::Array< T, 2, 2 > const &dt)
 Build a KD Tree to store the data for GaussianProcess. More...
 
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. More...
 
getData (int ipt, int idim) const
 Return one element of one node on the tree. More...
 
ndarray::Array< T, 1, 1 > getData (int ipt) const
 Return an entire node from the tree. More...
 
void addPoint (ndarray::Array< const T, 1, 1 > const &v)
 Add a point to the tree. More...
 
void removePoint (int dex)
 Remove a point from the tree. More...
 
int getNPoints () const
 return the number of data points stored in the tree More...
 
void getTreeNode (ndarray::Array< int, 1, 1 > const &v, int dex) const
 Return the _tree information for a given data point. More...
 

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]vthe point you are adding to the tree
Exceptions
pex::exceptions::RuntimeErrorif the branch ending in the new point is not properly constructed

Definition at line 210 of file GaussianProcess.cc.

210  {
211  if (v.getNumElements() != static_cast<ndarray::Size>(_dimensions)) {
213  "You are trying to add a point of the incorrect dimensionality to KdTree\n");
214  }
215 
216  int i, j, node, dim, dir;
217 
218  node = _findNode(v);
219  dim = _tree[node][DIMENSION] + 1;
220  if (dim == _dimensions) dim = 0;
221 
222  if (_npts == _room) {
223  ndarray::Array<T, 2, 2> dbuff = allocate(ndarray::makeVector(_npts, _dimensions));
224  ndarray::Array<int, 2, 2> tbuff = allocate(ndarray::makeVector(_npts, 4));
225 
226  dbuff.deep() = _data;
227  tbuff.deep() = _tree;
228 
229  _room += _roomStep;
230 
231  _tree = allocate(ndarray::makeVector(_room, 4));
232  _data = allocate(ndarray::makeVector(_room, _dimensions));
233 
234  for (i = 0; i < _npts; i++) {
235  for (j = 0; j < _dimensions; j++) _data[i][j] = dbuff[i][j];
236  for (j = 0; j < 4; j++) _tree[i][j] = tbuff[i][j];
237  }
238  }
239 
240  _tree[_npts][DIMENSION] = dim;
241  _tree[_npts][PARENT] = node;
242  i = _tree[node][DIMENSION];
243 
244  if (_data[node][i] > v[i]) {
245  if (_tree[node][LT] >= 0) {
247  "Trying to add to KdTree in a node that is already occupied\n");
248  }
249  _tree[node][LT] = _npts;
250  dir = LT;
251  } else {
252  if (_tree[node][GEQ] >= 0) {
254  "Trying to add to KdTree in a node that is already occupied\n");
255  }
256  _tree[node][GEQ] = _npts;
257  dir = GEQ;
258  }
259  _tree[_npts][LT] = -1;
260  _tree[_npts][GEQ] = -1;
261  for (i = 0; i < _dimensions; i++) {
262  _data[_npts][i] = v[i];
263  }
264 
265  _npts++;
266 
267  i = _walkUpTree(_tree[_npts - 1][PARENT], dir, _npts - 1);
268  if (i != _masterParent) {
269  throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError, "Adding to KdTree failed\n");
270  }
271 }
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
Reports errors that are due to events beyond the control of the program.
Definition: Runtime.h:104

◆ 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]neighdexthis is where the indices of the nearest neighbor points will be stored
[out]ddthis is where the distances to the nearest neighbors will be stored
[in]vthe point whose neighbors you want to find
[in]n_nnthe 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.

136  {
137  if (n_nn > _npts) {
139  "Asked for more neighbors than kd tree contains\n");
140  }
141 
142  if (n_nn <= 0) {
144  "Asked for zero or a negative number of neighbors\n");
145  }
146 
147  if (neighdex.getNumElements() != static_cast<ndarray::Size>(n_nn)) {
149  "Size of neighdex does not equal n_nn in KdTree.findNeighbors\n");
150  }
151 
152  if (dd.getNumElements() != static_cast<ndarray::Size>(n_nn)) {
154  "Size of dd does not equal n_nn in KdTree.findNeighbors\n");
155  }
156 
157  int i, start;
158 
159  _neighborCandidates = allocate(ndarray::makeVector(n_nn));
160  _neighborDistances = allocate(ndarray::makeVector(n_nn));
161  _neighborsFound = 0;
162  _neighborsWanted = n_nn;
163 
164  for (i = 0; i < n_nn; i++) _neighborDistances[i] = -1.0;
165 
166  start = _findNode(v);
167 
168  _neighborDistances[0] = _distance(v, _data[start]);
169  _neighborCandidates[0] = start;
170  _neighborsFound = 1;
171 
172  for (i = 1; i < 4; i++) {
173  if (_tree[start][i] >= 0) {
174  _lookForNeighbors(v, _tree[start][i], start);
175  }
176  }
177 
178  for (i = 0; i < n_nn; i++) {
179  neighdex[i] = _neighborCandidates[i];
180  dd[i] = _neighborDistances[i];
181  }
182 }

◆ 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]iptthe 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.

200  {
201  if (ipt >= _npts || ipt < 0) {
203  "Asked for point that does not exist in KdTree\n");
204  }
205 
206  return _data[ipt];
207 }

◆ 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]iptthe index of the node to return
[in]idimthe index of the dimension to return

Definition at line 185 of file GaussianProcess.cc.

185  {
186  if (ipt >= _npts || ipt < 0) {
188  "Asked for point that does not exist in KdTree\n");
189  }
190 
191  if (idim >= _dimensions || idim < 0) {
193  "KdTree does not have points of that many dimensions\n");
194  }
195 
196  return _data[ipt][idim];
197 }

◆ getNPoints()

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

return the number of data points stored in the tree

Definition at line 274 of file GaussianProcess.cc.

274  {
275  return _npts;
276 }

◆ 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]vthe array in which to store the entry from _tree
[in]dexthe index of the node whose information you are requesting

Definition at line 279 of file GaussianProcess.cc.

279  {
280  if (dex < 0 || dex >= _npts) {
282  "Asked for tree information on point that does not exist\n");
283  }
284 
285  if (v.getNumElements() != 4) {
287  "Need to pass a 4-element ndarray into KdTree.getTreeNode()\n");
288  }
289 
290  v[0] = _tree[dex][DIMENSION];
291  v[1] = _tree[dex][LT];
292  v[2] = _tree[dex][GEQ];
293  v[3] = _tree[dex][PARENT];
294 }

◆ 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]dtan array, the rows of which are the data points (dt[i][j] is the jth component of the ith data point)
Exceptions
pex::exceptions::RuntimeErrorif the tree is not properly constructed

Definition at line 104 of file GaussianProcess.cc.

104  {
105  int i;
106 
107  _npts = dt.template getSize<0>();
108  _dimensions = dt.template getSize<1>();
109 
110  // a buffer to use when first building the tree
111  _inn = allocate(ndarray::makeVector(_npts));
112 
113  _roomStep = 5000;
114  _room = _npts;
115 
116  _data = allocate(ndarray::makeVector(_room, _dimensions));
117 
118  _data.deep() = dt;
119 
120  _tree = allocate(ndarray::makeVector(_room, 4));
121 
122  for (i = 0; i < _npts; i++) {
123  _inn[i] = i;
124  }
125 
126  _organize(_inn, _npts, -1, -1);
127 
128  i = _testTree();
129  if (i == 0) {
130  throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError, "Failed to properly initialize KdTree\n");
131  }
132 }

◆ 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]dexthe index of the point you want to remove from the tree
Exceptions
pex::exceptions::RuntimeErrorif the entire tree is not poperly constructed after the point has been removed

Definition at line 582 of file GaussianProcess.cc.

582  {
583  if (target < 0 || target >= _npts) {
585  "You are trying to remove a point that doesn't exist from KdTree\n");
586  }
587 
588  if (_npts == 1) {
590  "There is only one point left in this KdTree. You cannot call removePoint\n");
591  }
592 
593  int nl, nr, i, j, k, side;
594  int root;
595 
596  nl = 0;
597  nr = 0;
598  if (_tree[target][LT] >= 0) {
599  nl++;
600  _count(_tree[target][LT], &nl);
601  }
602 
603  if (_tree[target][GEQ] >= 0) {
604  nr++;
605  _count(_tree[target][GEQ], &nr);
606  }
607 
608  if (nl == 0 && nr == 0) {
609  k = _tree[target][PARENT];
610 
611  if (_tree[k][LT] == target)
612  _tree[k][LT] = -1;
613  else if (_tree[k][GEQ] == target)
614  _tree[k][GEQ] = -1;
615 
616  } // if target is terminal
617  else if ((nl == 0 && nr > 0) || (nr == 0 && nl > 0)) {
618  if (nl == 0)
619  side = GEQ;
620  else
621  side = LT;
622 
623  k = _tree[target][PARENT];
624  if (k >= 0) {
625  if (_tree[k][LT] == target) {
626  _tree[k][LT] = _tree[target][side];
627  _tree[_tree[k][LT]][PARENT] = k;
628 
629  } else {
630  _tree[k][GEQ] = _tree[target][side];
631  _tree[_tree[k][GEQ]][PARENT] = k;
632  }
633  } else {
634  _masterParent = _tree[target][side];
635  _tree[_tree[target][side]][PARENT] = -1;
636  }
637  } // if only one side is populated
638  else {
639  if (nl > nr)
640  side = LT;
641  else
642  side = GEQ;
643 
644  k = _tree[target][PARENT];
645  if (k < 0) {
646  _masterParent = _tree[target][side];
647  _tree[_masterParent][PARENT] = -1;
648  } else {
649  if (_tree[k][LT] == target) {
650  _tree[k][LT] = _tree[target][side];
651  _tree[_tree[k][LT]][PARENT] = k;
652  } else {
653  _tree[k][GEQ] = _tree[target][side];
654  _tree[_tree[k][GEQ]][PARENT] = k;
655  }
656  }
657 
658  root = _tree[target][3 - side];
659 
660  _descend(root);
661  } // if both sides are populated
662 
663  if (target < _npts - 1) {
664  for (i = target + 1; i < _npts; i++) {
665  for (j = 0; j < 4; j++) _tree[i - 1][j] = _tree[i][j];
666 
667  for (j = 0; j < _dimensions; j++) _data[i - 1][j] = _data[i][j];
668  }
669 
670  for (i = 0; i < _npts; i++) {
671  for (j = 1; j < 4; j++) {
672  if (_tree[i][j] > target) _tree[i][j]--;
673  }
674  }
675 
676  if (_masterParent > target) _masterParent--;
677  }
678  _npts--;
679 
680  i = _testTree();
681  if (i == 0) {
682  throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError, "Subtracting from KdTree failed\n");
683  }
684 }
Key< Flag > const & target

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