//+------------------------------------------------------------------+ //| alglibmisc.mqh | //| Copyright 2003-2022 Sergey Bochkanov (ALGLIB project) | //| Copyright 2012-2026, MetaQuotes Ltd. | //| www.mql5.com | //+------------------------------------------------------------------+ //| Implementation of ALGLIB library in MetaQuotes Language 5 | //| | //| The features of the library include: | //| - Linear algebra (direct algorithms, EVD, SVD) | //| - Solving systems of linear and non-linear equations | //| - Interpolation | //| - Optimization | //| - FFT (Fast Fourier Transform) | //| - Numerical integration | //| - Linear and nonlinear least-squares fitting | //| - Ordinary differential equations | //| - Computation of special functions | //| - Descriptive statistics and hypothesis testing | //| - Data analysis - classification, regression | //| - Implementing linear algebra algorithms, interpolation, etc. | //| in high-precision arithmetic (using MPFR) | //| | //| This file is free software; you can redistribute it and/or | //| modify it under the terms of the GNU General Public License as | //| published by the Free Software Foundation (www.fsf.org); either | //| version 2 of the License, or (at your option) any later version. | //| | //| This program is distributed in the hope that it will be useful, | //| but WITHOUT ANY WARRANTY; without even the implied warranty of | //| MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | //| GNU General Public License for more details. | //+------------------------------------------------------------------+ #include "ap.mqh" #include "alglibinternal.mqh" //+------------------------------------------------------------------+ //| Buffer object which is used to perform nearest neighbor requests | //| in the multithreaded mode (multiple threads working with same | //| KD-tree object). | //| This object should be created with KDTreeCreateRequestBuffer(). | //+------------------------------------------------------------------+ class CKDTreeRequestBuffer { public: int m_kneeded; int m_kcur; bool m_selfmatch; double m_rneeded; double m_approxf; double m_curdist; //--- arrays CRowInt m_idx; CRowDouble m_x; CRowDouble m_boxmin; CRowDouble m_boxmax; CRowDouble m_r; CRowDouble m_buf; CRowDouble m_curboxmin; CRowDouble m_curboxmax; //--- constructor, destructor CKDTreeRequestBuffer(void) { Init(); } ~CKDTreeRequestBuffer(void) {} void Init(void) { m_kneeded=0; m_kcur=0; m_selfmatch=false; m_rneeded=0; m_approxf=0; m_curdist=0;} //--- copy void Copy(const CKDTreeRequestBuffer &obj); }; //+------------------------------------------------------------------+ //| Copy | //+------------------------------------------------------------------+ void CKDTreeRequestBuffer::Copy(const CKDTreeRequestBuffer &obj) { m_x=obj.m_x; m_boxmin=obj.m_boxmin; m_boxmax=obj.m_boxmax; m_kneeded=obj.m_kneeded; m_rneeded=obj.m_rneeded; m_selfmatch=obj.m_selfmatch; m_approxf=obj.m_approxf; m_kcur=obj.m_kcur; m_idx=obj.m_idx; m_r=obj.m_r; m_buf=obj.m_buf; m_curboxmin=obj.m_curboxmin; m_curboxmax=obj.m_curboxmax; m_curdist=obj.m_curdist; } //+------------------------------------------------------------------+ //| This class is a shell for class CKDTree | //+------------------------------------------------------------------+ class CKDTreeRequestBufferShell { private: CKDTreeRequestBuffer m_innerobj; public: //--- constructors, destructor CKDTreeRequestBufferShell(void) {} CKDTreeRequestBufferShell(CKDTreeRequestBuffer &obj) { m_innerobj.Copy(obj); } ~CKDTreeRequestBufferShell(void) {} //--- method CKDTreeRequestBuffer *GetInnerObj(void) { return(GetPointer(m_innerobj)); } }; //+------------------------------------------------------------------+ //| KD-trees | //+------------------------------------------------------------------+ class CKDTree { public: int m_n; int m_nx; int m_ny; int m_normtype; int m_debugcounter; //--- arrays CRowInt m_tags; CRowDouble m_boxmin; CRowDouble m_boxmax; CRowInt m_nodes; CRowDouble m_splits; //--- buffer CKDTreeRequestBuffer m_innerbuf; //--- matrix CMatrixDouble m_xy; //--- constructor, destructor CKDTree(void) { m_n=0; m_nx=0; m_ny=0; m_normtype=0; m_debugcounter=0; } ~CKDTree(void) {} //--- copy void Copy(const CKDTree &obj); //--- override void operator=(const CKDTree &obj) { Copy(obj); } }; //+------------------------------------------------------------------+ //| Copy | //+------------------------------------------------------------------+ void CKDTree::Copy(const CKDTree &obj) { //--- copy variables m_n=obj.m_n; m_nx=obj.m_nx; m_ny=obj.m_ny; m_normtype=obj.m_normtype; m_innerbuf.Copy(obj.m_innerbuf); m_debugcounter=obj.m_debugcounter; //--- copy arrays m_tags=obj.m_tags; m_boxmin=obj.m_boxmin; m_boxmax=obj.m_boxmax; m_nodes=obj.m_nodes; m_splits=obj.m_splits; //--- copy matrix m_xy=obj.m_xy; } //+------------------------------------------------------------------+ //| This class is a shell for class CKDTree | //+------------------------------------------------------------------+ class CKDTreeShell { private: CKDTree m_innerobj; public: //--- constructors, destructor CKDTreeShell(void) {} CKDTreeShell(CKDTree &obj) { m_innerobj.Copy(obj); } ~CKDTreeShell(void) {} //--- method CKDTree *GetInnerObj(void) { return(GetPointer(m_innerobj)); } }; //+------------------------------------------------------------------+ //| Build KD-trees | //+------------------------------------------------------------------+ class CNearestNeighbor { public: //--- class constants static const int m_splitnodesize; static const int m_kdtreefirstversion; //--- build static void KDTreeBuild(CMatrixDouble &xy,const int n,const int nx,const int ny,const int normtype,CKDTree &kdt); static void KDTreeBuildTagged(CMatrixDouble &xy,int &tags[],const int n,const int nx,const int ny,const int normtype,CKDTree &kdt); static void KDTreeBuildTagged(CMatrixDouble &xy,CRowInt &tags,const int n,const int nx,const int ny,const int normtype,CKDTree &kdt); static void KDTreeCreateRequestBuffer(CKDTree &kdt,CKDTreeRequestBuffer &buf); static int KDTreeQueryKNN(CKDTree &kdt,double &x[],const int k,const bool selfmatch=true); static int KDTreeQueryKNN(CKDTree &kdt,CRowDouble &x,const int k,const bool selfmatch=true); static int KDTreeTsQueryKNN(CKDTree &kdt,CKDTreeRequestBuffer &buf,CRowDouble &x,int k,bool selfmatch=true); static int KDTreeQueryRNN(CKDTree &kdt,double &x[],const double r,const bool selfmatch=true); static int KDTreeQueryRNN(CKDTree &kdt,CRowDouble &x,const double r,const bool selfmatch=true); static int KDTreeQueryRNNU(CKDTree &kdt,double &x[],const double r,const bool selfmatch=true); static int KDTreeQueryRNNU(CKDTree &kdt,CRowDouble &x,const double r,const bool selfmatch=true); static int KDTreeTsQueryRNN(CKDTree &kdt,CKDTreeRequestBuffer &buf,CRowDouble &x,const double r,const bool selfmatch=true); static int KDTreeTsQueryRNNU(CKDTree &kdt,CKDTreeRequestBuffer &buf,CRowDouble &x,const double r,const bool selfmatch=true); static int KDTreeQueryAKNN(CKDTree &kdt,double &x[],int k,const bool selfmatch=true,const double eps=0); static int KDTreeQueryAKNN(CKDTree &kdt,CRowDouble &x,int k,const bool selfmatch=true,const double eps=0); static int KDTreeTsQueryAKNN(CKDTree &kdt,CKDTreeRequestBuffer &buf,CRowDouble &x,int k,const bool selfmatch=true,const double eps=0); static int KDTreeQueryBox(CKDTree &kdt,double &boxmin[],double &boxmax[]); static int KDTreeQueryBox(CKDTree &kdt,CRowDouble &boxmin,CRowDouble &boxmax); static int KDTreeTsQueryBox(CKDTree &kdt,CKDTreeRequestBuffer &buf,double &boxmin[],double &boxmax[]); static int KDTreeTsQueryBox(CKDTree &kdt,CKDTreeRequestBuffer &buf,CRowDouble &boxmin,CRowDouble &boxmax); static void KDTreeQueryResultsX(CKDTree &kdt,CMatrixDouble &x); static void KDTreeTsQueryResultsX(CKDTree &kdt,CKDTreeRequestBuffer &buf,CMatrixDouble &x); static void KDTreeQueryResultsXY(CKDTree &kdt,CMatrixDouble &xy); static void KDTreeTsQueryResultsXY(CKDTree &kdt,CKDTreeRequestBuffer &buf,CMatrixDouble &xy); static void KDTreeQueryResultsTags(CKDTree &kdt,int &tags[]); static void KDTreeQueryResultsTags(CKDTree &kdt,CRowInt &tags); static void KDTreeTsQueryResultsTags(CKDTree &kdt,CKDTreeRequestBuffer &buf,int &tags[]); static void KDTreeTsQueryResultsTags(CKDTree &kdt,CKDTreeRequestBuffer &buf,CRowInt &tags); static void KDTreeQueryResultsDistances(CKDTree &kdt,double &r[]); static void KDTreeQueryResultsDistances(CKDTree &kdt,CRowDouble &r); static void KDTreeTsQueryResultsDistances(CKDTree &kdt,CKDTreeRequestBuffer &buf,double &r[]); static void KDTreeTsQueryResultsDistances(CKDTree &kdt,CKDTreeRequestBuffer &buf,CRowDouble &r); static void KDTreeQueryResultsXI(CKDTree &kdt,CMatrixDouble &x); static void KDTreeQueryResultsXYI(CKDTree &kdt,CMatrixDouble &xy); static void KDTreeQueryResultsTagsI(CKDTree &kdt,int &tags[]); static void KDTreeQueryResultsTagsI(CKDTree &kdt,CRowInt &tags); static void KDTreeQueryResultsDistancesI(CKDTree &kdt,double &r[]); static void KDTreeQueryResultsDistancesI(CKDTree &kdt,CRowDouble &r); //--- information static void KDTreeExploreBox(CKDTree &kdt,CRowDouble &boxmin,CRowDouble &boxmax); static void KDTreeExploreNodeType(CKDTree &kdt,int node,int &nodetype); static void KDTreeExploreLeaf(CKDTree &kdt,int node,CMatrixDouble &xy,int &k); static void KDTreeExploreSplit(CKDTree &kdt,int node,int &d,double &s,int &nodele,int &nodege); //--- serialize static void KDTreeAlloc(CSerializer &s,CKDTree &tree); static void KDTreeSerialize(CSerializer &s,CKDTree &tree); static void KDTreeUnserialize(CSerializer &s,CKDTree &tree); private: static int TsQueryRNN(CKDTree &kdt,CKDTreeRequestBuffer &buf,CRowDouble &x,double r,bool selfmatch=true,bool orderedbydist=true); static bool CheckRequestBufferConsistency(CKDTree &kdt,CKDTreeRequestBuffer &buf); static void KDTreeSplit(CKDTree &kdt,const int i1,const int i2,const int d,const double s,int &i3); static void KDTreeGenerateTreeRec(CKDTree &kdt,int &nodesoffs,int &splitsoffs,const int i1,const int i2,const int maxleafsize); static void KDTreeQueryNNRec(CKDTree &kdt,const int offs); static void KDTreeQueryNNRec(CKDTree &kdt,CKDTreeRequestBuffer &buf,const int offs); static void KDTreeQueryBoxRec(CKDTree &kdt,CKDTreeRequestBuffer &buf,int offs); static void KDTreeInitBox(CKDTree &kdt,CRowDouble &x,CKDTreeRequestBuffer &buf); static void KDTreeAllocDataSetIndependent(CKDTree &kdt,const int nx,const int ny); static void KDTreeAllocDataSetDependent(CKDTree &kdt,const int n,const int nx,const int ny); }; //+------------------------------------------------------------------+ //| Initialize constants | //+------------------------------------------------------------------+ const int CNearestNeighbor::m_splitnodesize=6; const int CNearestNeighbor::m_kdtreefirstversion=0; //+------------------------------------------------------------------+ //| KD-tree creation | //| This subroutine creates KD-tree from set of X-values and optional| //| Y-values | //| INPUT PARAMETERS | //| XY - dataset, array[0..N-1,0..NX+NY-1]. | //| one row corresponds to one point. | //| first NX columns contain X-values, next NY (NY | //| may be zero) | //| columns may contain associated Y-values | //| N - number of points, N>=1 | //| NX - space dimension, NX>=1. | //| NY - number of optional Y-values, NY>=0. | //| NormType- norm type: | //| * 0 denotes infinity-norm | //| * 1 denotes 1-norm | //| * 2 denotes 2-norm (Euclidean norm) | //| OUTPUT PARAMETERS | //| KDT - KD-tree | //| NOTES | //| 1. KD-tree creation have O(N*logN) complexity and | //| O(N*(2*NX+NY)) memory requirements. | //| 2. Although KD-trees may be used with any combination of N and | //| NX, they are more efficient than brute-force search only when | //| N >> 4^NX. So they are most useful in low-dimensional tasks | //| (NX=2, NX=3). NX=1 is another inefficient case, because | //| simple binary search (without additional structures) is | //| much more efficient in such tasks than KD-trees. | //+------------------------------------------------------------------+ void CNearestNeighbor::KDTreeBuild(CMatrixDouble &xy,const int n, const int nx,const int ny, const int normtype,CKDTree &kdt) { int tags[]; //--- check if(!CAp::Assert(n>=0,__FUNCTION__+": N<0!")) return; //--- check if(!CAp::Assert(nx>=1,__FUNCTION__+": NX<1!")) return; //--- check if(!CAp::Assert(ny>=0,__FUNCTION__+": NY<0!")) return; //--- check if(!CAp::Assert(normtype>=0 && normtype<=2,__FUNCTION__+": incorrect NormType!")) return; //--- check if(!CAp::Assert((int)CAp::Rows(xy)>=n,__FUNCTION__+": rows(X)=nx+ny,__FUNCTION__+": cols(X)=1 | //| NX - space dimension, NX>=1. | //| NY - number of optional Y-values, NY>=0. | //| NormType- norm type: | //| * 0 denotes infinity-norm | //| * 1 denotes 1-norm | //| * 2 denotes 2-norm (Euclidean norm) | //| OUTPUT PARAMETERS | //| KDT - KD-tree | //| NOTES | //| 1. KD-tree creation have O(N*logN) complexity and | //| O(N*(2*NX+NY)) memory requirements. | //| 2. Although KD-trees may be used with any combination of N and | //| NX, they are more efficient than brute-force search only when | //| N >> 4^NX. So they are most useful in low-dimensional tasks | //| (NX=2, NX=3). NX=1 is another inefficient case, because simple| //| binary search (without additional structures) is much more | //| efficient in such tasks than KD-trees. | //+------------------------------------------------------------------+ void CNearestNeighbor::KDTreeBuildTagged(CMatrixDouble &xy,int &tags[], const int n,const int nx, const int ny, const int normtype,CKDTree &kdt) { CRowInt Tags=tags; KDTreeBuildTagged(xy,Tags,n,nx,ny,normtype,kdt); } //+------------------------------------------------------------------+ //| | //+------------------------------------------------------------------+ void CNearestNeighbor::KDTreeBuildTagged(CMatrixDouble &xy,CRowInt &tags, const int n,const int nx, const int ny, const int normtype,CKDTree &kdt) { int i=0; int j=0; int nodesoffs=0; int splitsoffs=0; int i_=0; int i1_=0; //--- check if(!CAp::Assert(n>=0,__FUNCTION__+": N<0!")) return; //--- check if(!CAp::Assert(nx>=1,__FUNCTION__+": NX<1!")) return; //--- check if(!CAp::Assert(ny>=0,__FUNCTION__+": NY<0!")) return; //--- check if(!CAp::Assert(normtype>=0 && normtype<=2,__FUNCTION__+": incorrect NormType!")) return; //--- check if(!CAp::Assert((int)CAp::Rows(xy)>=n,__FUNCTION__+": rows(X)=nx+ny || n==0,__FUNCTION__+": cols(X) quick exit if(n==0) return; //--- Allocate KDTreeAllocDataSetIndependent(kdt,nx,ny); KDTreeAllocDataSetDependent(kdt,n,nx,ny); KDTreeCreateRequestBuffer(kdt,kdt.m_innerbuf); //--- Initial fill for(i_=0; i_ col=xy.Col(i_); kdt.m_xy.Col(i_,col); } i1_=-nx; for(i_=nx; i_<2*nx+ny; i_++) { vector col=xy.Col(i_+i1_); kdt.m_xy.Col(i_,col); } kdt.m_tags=tags; //--- Determine bounding box kdt.m_boxmin=kdt.m_xy.Min(0)+0; kdt.m_boxmax=kdt.m_xy.Max(0)+0; kdt.m_boxmin.Resize(kdt.m_nx); kdt.m_boxmax.Resize(kdt.m_nx); //--- prepare tree structure //--- * MaxNodes=N because we guarantee no trivial splits,i.e. //--- every split will generate two non-empty boxes //--- allocation nodesoffs=0; splitsoffs=0; kdt.m_innerbuf.m_curboxmin=kdt.m_boxmin; kdt.m_innerbuf.m_curboxmax=kdt.m_boxmax; //--- function call KDTreeGenerateTreeRec(kdt,nodesoffs,splitsoffs,0,n,8); kdt.m_nodes.Resize(nodesoffs); kdt.m_splits.Resize(splitsoffs); } //+------------------------------------------------------------------+ //| This function creates buffer structure which can be used | //| to perform parallel KD-tree requests. | //| KD-tree subpackage provides two sets of request functions - ones | //| which use internal buffer of KD-tree object (these functions are | //| single-threaded because they use same buffer, which can not | //| shared between threads), and ones which use external buffer. | //| This function is used to initialize external buffer. | //| INPUT PARAMETERS | //| KDT - KD-tree which is associated with newly created buffer | //| OUTPUT PARAMETERS | //| Buf - external buffer. | //| IMPORTANT: KD-tree buffer should be used only with KD-tree object| //| which was used to initialize buffer. Any attempt | //| to use buffer with different object is dangerous - you| //| may get integrity check failure (exception) because | //| sizes of internal arrays do not fit to dimensions of | //| KD-tree structure. | //+------------------------------------------------------------------+ void CNearestNeighbor::KDTreeCreateRequestBuffer(CKDTree &kdt, CKDTreeRequestBuffer &buf) { buf.m_x.Resize(kdt.m_nx); buf.m_boxmin.Resize(kdt.m_nx); buf.m_boxmax.Resize(kdt.m_nx); buf.m_idx.Resize(kdt.m_n); buf.m_r.Resize(kdt.m_n); buf.m_buf.Resize(MathMax(kdt.m_n,kdt.m_nx)); buf.m_curboxmin.Resize(kdt.m_nx); buf.m_curboxmax.Resize(kdt.m_nx); buf.m_kcur=0; } //+------------------------------------------------------------------+ //| K-NN query: K nearest neighbors | //| INPUT PARAMETERS | //| KDT - KD-tree | //| X - point, array[0..NX-1]. | //| K - number of neighbors to return, K>=1 | //| SelfMatch - whether self-matches are allowed: | //| * if True, nearest neighbor may be the point | //| itself (if it exists in original dataset) | //| * if False, then only points with non-zero | //| distance are returned | //| * if not given, considered True | //| RESULT | //| number of actual neighbors found (either K or N, if K>N). | //| This subroutine performs query and stores its result in the | //| internal structures of the KD-tree. You can use following | //| subroutines to obtain these results: | //| * KDTreeQueryResultsX() to get X-values | //| * KDTreeQueryResultsXY() to get X- and Y-values | //| * KDTreeQueryResultsTags() to get tag values | //| * KDTreeQueryResultsDistances() to get distances | //+------------------------------------------------------------------+ int CNearestNeighbor::KDTreeQueryKNN(CKDTree &kdt,double &x[], const int k,const bool selfmatch) { //--- check if(!CAp::Assert(k>=1,__FUNCTION__+": K<1!")) return(-1); //--- check if(!CAp::Assert(CAp::Len(x)>=kdt.m_nx,__FUNCTION__+": Length(X)=1,__FUNCTION__+": K<1!")) return(-1); //--- check if(!CAp::Assert(CAp::Len(x)>=kdt.m_nx,__FUNCTION__+": Length(X)=1 | //| SelfMatch - whether self-matches are allowed: | //| * if True, nearest neighbor may be the point | //| itself (if it exists in original dataset) | //| * if False, then only points with non-zero | //| distance are returned | //| * if not given, considered True | //| RESULT | //| number of actual neighbors found (either K or N, if K>N). | //| This subroutine performs query and stores its result in the | //| internal structures of the buffer object. You can use following | //| subroutines to obtain these results (pay attention to "buf" in | //| their names): | //| * KDTreeTsQueryResultsX() to get X-values | //| * KDTreeTsQueryResultsXY() to get X- and Y-values | //| * KDTreeTsQueryResultsTags() to get tag values | //| * KDTreeTsQueryResultsDistances() to get distances | //| IMPORTANT: kd-tree buffer should be used only with KD-tree object| //| which was used to initialize buffer. Any attempt to use biffer | //| with different object is dangerous - you may get integrity check | //| failure (exception) because sizes of internal arrays do not fit | //| to dimensions of KD-tree structure. | //+------------------------------------------------------------------+ int CNearestNeighbor::KDTreeTsQueryKNN(CKDTree &kdt, CKDTreeRequestBuffer &buf, CRowDouble &x, int k, bool selfmatch=true) { //--- check if(!CAp::Assert(k>=1,__FUNCTION__+": K<1!")) return(-1); //--- check if(!CAp::Assert(CAp::Len(x)>=kdt.m_nx,__FUNCTION__+": Length(X)0| //| SelfMatch - whether self-matches are allowed: | //| * if True, nearest neighbor may be the point | //| itself (if it exists in original dataset) | //| * if False, then only points with non-zero | //| distance are returned | //| * if not given, considered True | //| RESULT | //| number of neighbors found, >=0 | //| This subroutine performs query and stores its result in the | //| internal structures of the KD-tree. You can use following | //| subroutines to obtain actual results: | //| * KDTreeQueryResultsX() to get X-values | //| * KDTreeQueryResultsXY() to get X- and Y-values | //| * KDTreeQueryResultsTags() to get tag values | //| * KDTreeQueryResultsDistances() to get distances | //+------------------------------------------------------------------+ int CNearestNeighbor::KDTreeQueryRNN(CKDTree &kdt,double &x[], const double r,const bool selfmatch=true) { CRowDouble vec=x; //--- return result return(KDTreeQueryRNN(kdt,vec,r,selfmatch)); } //+------------------------------------------------------------------+ //| | //+------------------------------------------------------------------+ int CNearestNeighbor::KDTreeQueryRNN(CKDTree &kdt,CRowDouble &x, const double r,const bool selfmatch=true) { //--- check if(!CAp::Assert((double)(r)>0.0,__FUNCTION__+": incorrect R!")) return(-1); //--- check if(!CAp::Assert((int)x.Size()>=kdt.m_nx,__FUNCTION__+": Length(X)0 | //| SelfMatch - whether self-matches are allowed: | //| * if True, nearest neighbor may be the point | //| itself (if it exists in original dataset) | //| * if False, then only points with non-zero | //| distance are returned | //| * if not given, considered True | //| RESULT | //| number of neighbors found, >=0 | //| This subroutine performs query and stores its result in the | //| internal structures of the KD-tree. You can use following | //| subroutines to obtain actual results: | //| * KDTreeQueryResultsX() to get X-values | //| * KDTreeQueryResultsXY() to get X- and Y-values | //| * KDTreeQueryResultsTags() to get tag values | //| * KDTreeQueryResultsDistances() to get distances | //| As indicated by "U" suffix, this function returns unordered | //| results. | //+------------------------------------------------------------------+ int CNearestNeighbor::KDTreeQueryRNNU(CKDTree &kdt, double &x[], double r, bool selfmatch=true) { CRowDouble vec=x; //--- return result return(KDTreeQueryRNNU(kdt,vec,r,selfmatch)); } //+------------------------------------------------------------------+ //| | //+------------------------------------------------------------------+ int CNearestNeighbor::KDTreeQueryRNNU(CKDTree &kdt, CRowDouble &x, double r, bool selfmatch=true) { //--- check if(!CAp::Assert(r>0.0,__FUNCTION__+": incorrect R!")) return(-1); //--- check if(!CAp::Assert((int)x.Size()>=kdt.m_nx,__FUNCTION__+": Length(X)0 | //| SelfMatch - whether self-matches are allowed: | //| * if True, nearest neighbor may be the point itself | //| (if it exists in original dataset) | //| * if False, then only points with non-zero distance | //| are returned | //| * if not given, considered True | //| RESULT | //| number of neighbors found, >=0 | //| This subroutine performs query and stores its result in the | //| internal structures of the buffer object. You can use following | //| subroutines to obtain these results (pay attention to "buf" in | //| their names): | //| * KDTreeTsQueryResultsX() to get X-values | //| * KDTreeTsQueryResultsXY() to get X- and Y-values | //| * KDTreeTsQueryResultsTags() to get tag values | //| * KDTreeTsQueryResultsDistances() to get distances | //| IMPORTANT: kd-tree buffer should be used only with KD-tree object| //| which was used to initialize buffer. Any attempt to | //| use biffer with different object is dangerous - you | //| may get integrity check failure (exception) because | //| sizes of internal arrays do not fit to dimensions of | //| KD-tree structure. | //+------------------------------------------------------------------+ int CNearestNeighbor::KDTreeTsQueryRNN(CKDTree &kdt, CKDTreeRequestBuffer &buf, CRowDouble &x, const double r, const bool selfmatch=true) { //--- check if(!CAp::Assert(MathIsValidNumber(r) && r>0.0,__FUNCTION__+": incorrect R!")) return(-1); //--- check if(!CAp::Assert((int)x.Size()>=kdt.m_nx,__FUNCTION__+": Length(X)0 | //| SelfMatch - whether self-matches are allowed: | //| * if True, nearest neighbor may be the point itself| //| (if it exists in original dataset) | //| * if False, then only points with non-zero distance| //| are returned | //| * if not given, considered True | //| RESULT | //| number of neighbors found, >=0 | //| This subroutine performs query and stores its result in the | //| internal structures of the buffer object. You can use following | //| subroutines to obtain these results (pay attention to "buf" in | //| their names): | //| * KDTreeTsQueryResultsX() to get X-values | //| * KDTreeTsQueryResultsXY() to get X- and Y-values | //| * KDTreeTsQueryResultsTags() to get tag values | //| * KDTreeTsQueryResultsDistances() to get distances | //| As indicated by "U" suffix, this function returns unordered | //| results. | //| IMPORTANT: kd-tree buffer should be used only with KD-tree object| //| which was used to initialize buffer. Any attempt to | //| use biffer with different object is dangerous - you | //| may get integrity check failure (exception) because | //| sizes of internal arrays do not fit to dimensions of | //| KD-tree structure. | //+------------------------------------------------------------------+ int CNearestNeighbor::KDTreeTsQueryRNNU(CKDTree &kdt, CKDTreeRequestBuffer &buf, CRowDouble &x, const double r, const bool selfmatch=true) { //--- check if(!CAp::Assert(MathIsValidNumber(r) && r>0.0,__FUNCTION__+": incorrect R!")) return(-1); //--- check if(!CAp::Assert((int)x.Size()>=kdt.m_nx,__FUNCTION__+": Length(X)=1 | //| SelfMatch - whether self-matches are allowed: | //| * if True, nearest neighbor may be the point | //| itself (if it exists in original dataset) | //| * if False, then only points with non-zero | //| distance are returned | //| * if not given, considered True | //| Eps - approximation factor, Eps>=0. eps-approximate| //| nearest neighbor is a neighbor whose distance| //| from X is at most (1+eps) times distance of | //| true nearest neighbor. | //| RESULT | //| number of actual neighbors found (either K or N, if K>N). | //| NOTES | //| significant performance gain may be achieved only when Eps is| //| on the order of magnitude of 1 or larger. | //| This subroutine performs query and stores its result in the | //| internal structures of the KD-tree. You can use following | //| these subroutines to obtain results: | //| * KDTreeQueryResultsX() to get X-values | //| * KDTreeQueryResultsXY() to get X- and Y-values | //| * KDTreeQueryResultsTags() to get tag values | //| * KDTreeQueryResultsDistances() to get distances | //+------------------------------------------------------------------+ int CNearestNeighbor::KDTreeQueryAKNN(CKDTree &kdt,double &x[], int k,const bool selfmatch=true, const double eps=0) { CRowDouble vec=x; //--- return result return(KDTreeTsQueryAKNN(kdt,kdt.m_innerbuf,vec,k,selfmatch,eps)); } //+------------------------------------------------------------------+ //| | //+------------------------------------------------------------------+ int CNearestNeighbor::KDTreeQueryAKNN(CKDTree &kdt,CRowDouble &x, int k,const bool selfmatch=true, const double eps=0) { return(KDTreeTsQueryAKNN(kdt,kdt.m_innerbuf,x,k,selfmatch,eps)); } //+------------------------------------------------------------------+ //| K-NN query: approximate K nearest neighbors, using thread-local | //| buffer. | //| You can call this function from multiple threads for same kd-tree| //| instance, assuming that different instances of buffer object are | //| passed to different threads. | //| INPUT PARAMETERS | //| KDT - KD-tree | //| Buf - request buffer object created for this particular | //| instance of kd-tree structure with | //| KDTreeCreateRequestBuffer() function. | //| X - point, array[0..NX-1]. | //| K - number of neighbors to return, K>=1 | //| SelfMatch - whether self-matches are allowed: | //| * if True, nearest neighbor may be the point itself| //| (if it exists in original dataset) | //| * if False, then only points with non-zero distance| //| are returned | //| * if not given, considered True | //| Eps - approximation factor, Eps>=0. eps-approximate nearest | //| neighbor is a neighbor whose distance from X is at | //| most (1+eps) times distance of true nearest neighbor. | //| RESULT | //| number of actual neighbors found (either K or N, if K>N). | //| NOTES | //| significant performance gain may be achieved only when Eps is | //| on the order of magnitude of 1 or larger. | //| This subroutine performs query and stores its result in the | //| internal structures of the buffer object. You can use following | //| subroutines to obtain these results (pay attention to "buf" in | //| their names): | //| * KDTreeTsQueryResultsX() to get X-values | //| * KDTreeTsQueryResultsXY() to get X- and Y-values | //| * KDTreeTsQueryResultsTags() to get tag values | //| * KDTreeTsQueryResultsDistances() to get distances | //| IMPORTANT: kd-tree buffer should be used only with KD-tree object| //| which was used to initialize buffer. Any attempt to | //| use biffer with different object is dangerous - you | //| may get integrity check failure (exception) because | //| sizes of internal arrays do not fit to dimensions of | //| KD-tree structure. | //+------------------------------------------------------------------+ int CNearestNeighbor::KDTreeTsQueryAKNN(CKDTree &kdt, CKDTreeRequestBuffer &buf, CRowDouble &x, int k, bool selfmatch=true, double eps=0) { int result=0; int i=0; int j=0; //--- check if(!CAp::Assert(k>0,__FUNCTION__+": incorrect K!")) return(-1); //--- check if(!CAp::Assert(eps>=0.0,__FUNCTION__+": incorrect Eps!")) return(-1); //--- check if(!CAp::Assert((int)x.Size()>=kdt.m_nx,__FUNCTION__+": Length(X)=2; i--) CTSort::TagHeapPopI(buf.m_r,buf.m_idx,j); return(result); } //+------------------------------------------------------------------+ //| Box query: all points within user-specified box. | //| IMPORTANT: this function can not be used in multithreaded code | //| because it uses internal temporary buffer of kd-tree | //| object, which can not be shared between multiple | //| threads. If you want to perform parallel requests, | //| use function which uses external request buffer: | //| KDTreeTsQueryBox() ("Ts" stands for "thread-safe"). | //| INPUT PARAMETERS | //| KDT - KD-tree | //| BoxMin - lower bounds, array[0..NX-1]. | //| BoxMax - upper bounds, array[0..NX-1]. | //| RESULT | //| number of actual neighbors found (in [0,N]). | //| This subroutine performs query and stores its result in the | //| internal structures of the KD-tree. You can use following | //| subroutines to obtain these results: | //| * KDTreeQueryResultsX() to get X-values | //| * KDTreeQueryResultsXY() to get X- and Y-values | //| * KDTreeQueryResultsTags() to get tag values | //| * KDTreeQueryResultsDistances() returns zeros for this request | //| NOTE: this particular query returns unordered results, because | //| there is no meaningful way of ordering points. Furthermore,| //| no 'distance' is associated with points - it is either | //| INSIDE or OUTSIDE (so request for distances will return | //| zeros). | //+------------------------------------------------------------------+ int CNearestNeighbor::KDTreeQueryBox(CKDTree &kdt, double &boxmin[], double &boxmax[]) { return(KDTreeTsQueryBox(kdt,kdt.m_innerbuf,boxmin,boxmax)); } //+------------------------------------------------------------------+ //| | //+------------------------------------------------------------------+ int CNearestNeighbor::KDTreeQueryBox(CKDTree &kdt, CRowDouble &boxmin, CRowDouble &boxmax) { return(KDTreeTsQueryBox(kdt,kdt.m_innerbuf,boxmin,boxmax)); } //+------------------------------------------------------------------+ //| Box query: all points within user-specified box, using | //| thread-local buffer. | //| You can call this function from multiple threads for same kd-tree| //| instance, assuming that different instances of buffer object are | //| passed to different threads. | //| INPUT PARAMETERS | //| KDT - KD-tree | //| Buf - request buffer object created for this particular | //| instance of kd-tree structure with | //| KDTreeCreateRequestBuffer() function. | //| BoxMin - lower bounds, array[0..NX-1]. | //| BoxMax - upper bounds, array[0..NX-1]. | //| RESULT | //| number of actual neighbors found (in [0,N]). | //| This subroutine performs query and stores its result in the | //| internal structures of the buffer object. You can use following | //| subroutines to obtain these results (pay attention to "ts" in | //| their names): | //| * KDTreeTsQueryResultsX() to get X-values | //| * KDTreeTsQueryResultsXY() to get X- and Y-values | //| * KDTreeTsQueryResultsTags() to get tag values | //| * KDTreeTsQueryResultsDistances() returns zeros for this query | //| NOTE: this particular query returns unordered results, because | //| there is no meaningful way of ordering points. Furthermore,| //| no 'distance' is associated with points - it is either | //| INSIDE or OUTSIDE (so request for distances will return | //| zeros). | //| IMPORTANT: kd-tree buffer should be used only with KD-tree object| //| which was used to initialize buffer. Any attempt to | //| use biffer with different object is dangerous - you | //| may get integrity check failure (exception) because| //| sizes of internal arrays do not fit to dimensions of | //| KD-tree structure. | //+------------------------------------------------------------------+ int CNearestNeighbor::KDTreeTsQueryBox(CKDTree &kdt, CKDTreeRequestBuffer &buf, double &boxmin[], double &boxmax[]) { int result=0; //--- check if(!CAp::Assert(CAp::Len(boxmin)>=kdt.m_nx,__FUNCTION__+": Length(BoxMin)=kdt.m_nx,__FUNCTION__+": Length(BoxMax)(double)(boxmax[j])) { buf.m_kcur=0; return(0); } //--- prepare parameters buf.m_boxmin=boxmin; buf.m_boxmax=boxmax; buf.m_curboxmin=boxmin; buf.m_curboxmax=boxmax; buf.m_kcur=0; //--- call recursive search KDTreeQueryBoxRec(kdt,buf,0); result=buf.m_kcur; return(result); } //+------------------------------------------------------------------+ //| | //+------------------------------------------------------------------+ int CNearestNeighbor::KDTreeTsQueryBox(CKDTree &kdt, CKDTreeRequestBuffer &buf, CRowDouble &boxmin, CRowDouble &boxmax) { int result=0; //--- check if(!CAp::Assert(CAp::Len(boxmin)>=kdt.m_nx,__FUNCTION__+": Length(BoxMin)=kdt.m_nx,__FUNCTION__+": Length(BoxMax)(double)(boxmax[j])) { buf.m_kcur=0; return(0); } //--- prepare parameters buf.m_boxmin=boxmin; buf.m_boxmax=boxmax; buf.m_curboxmin=boxmin; buf.m_curboxmax=boxmax; buf.m_kcur=0; //--- call recursive search KDTreeQueryBoxRec(kdt,buf,0); result=buf.m_kcur; return(result); } //+------------------------------------------------------------------+ //| X-values from last query | //| INPUT PARAMETERS | //| KDT - KD-tree | //| X - possibly pre-allocated buffer. If X is too small | //| to store result, it is resized. If size(X) is | //| enough to store result, it is left unchanged. | //| OUTPUT PARAMETERS | //| X - rows are filled with X-values | //| NOTES | //| 1. points are ordered by distance from the query point (first = | //| closest) | //| 2. if XY is larger than required to store result, only leading | //| part will be overwritten; trailing part will be left | //| unchanged. So if on input XY = [[A,B],[C,D]], and result is | //| [1,2], then on exit we will get XY = [[1,2],[C,D]]. This is | //| done purposely to increase performance; if you want function | //| to resize array according to result size, use function with | //| same name and suffix 'I'. | //| SEE ALSO | //| * KDTreeQueryResultsXY() X- and Y-values | //| * KDTreeQueryResultsTags() tag values | //| * KDTreeQueryResultsDistances() distances | //+------------------------------------------------------------------+ void CNearestNeighbor::KDTreeQueryResultsX(CKDTree &kdt,CMatrixDouble &x) { KDTreeTsQueryResultsX(kdt,kdt.m_innerbuf,x); } //+------------------------------------------------------------------+ //| X- and Y-values from last query | //| INPUT PARAMETERS | //| KDT - KD-tree | //| XY - possibly pre-allocated buffer. If XY is too small| //| to store result, it is resized. If size(XY) is | //| enough to store result, it is left unchanged. | //| OUTPUT PARAMETERS | //| XY - rows are filled with points: first NX columns | //| with X-values, next NY columns - with Y-values. | //| NOTES | //| 1. points are ordered by distance from the query point (first = | //| closest) | //| 2. if XY is larger than required to store result, only leading | //| part will be overwritten; trailing part will be left | //| unchanged. So if on input XY = [[A,B],[C,D]], and result is | //| [1,2], then on exit we will get XY = [[1,2],[C,D]]. This is | //| done purposely to increase performance; if you want function | //| to resize array according to result size, use function with | //| same name and suffix 'I'. | //| SEE ALSO | //| * KDTreeQueryResultsX() X-values | //| * KDTreeQueryResultsTags() tag values | //| * KDTreeQueryResultsDistances() distances | //+------------------------------------------------------------------+ void CNearestNeighbor::KDTreeQueryResultsXY(CKDTree &kdt,CMatrixDouble &xy) { KDTreeTsQueryResultsXY(kdt,kdt.m_innerbuf,xy); } //+------------------------------------------------------------------+ //| Tags from last query | //| INPUT PARAMETERS | //| KDT - KD-tree | //| Tags - possibly pre-allocated buffer. If X is too small | //| to store result, it is resized. If size(X) is | //| enough to store result, it is left unchanged. | //| OUTPUT PARAMETERS | //| Tags - filled with tags associated with points, | //| or, when no tags were supplied, with zeros | //| NOTES | //| 1. points are ordered by distance from the query point (first | //| = closest) | //| 2. if XY is larger than required to store result, only leading | //| part will be overwritten; trailing part will be left | //| unchanged. So if on input XY = [[A,B],[C,D]], and result is | //| [1,2], then on exit we will get XY = [[1,2],[C,D]]. This is | //| done purposely to increase performance; if you want function | //| to resize array according to result size, use function with | //| same name and suffix 'I'. | //| SEE ALSO | //| * KDTreeQueryResultsX() X-values | //| * KDTreeQueryResultsXY() X- and Y-values | //| * KDTreeQueryResultsDistances() distances | //+------------------------------------------------------------------+ void CNearestNeighbor::KDTreeQueryResultsTags(CKDTree &kdt,int &tags[]) { KDTreeTsQueryResultsTags(kdt,kdt.m_innerbuf,tags); } //+------------------------------------------------------------------+ //| | //+------------------------------------------------------------------+ void CNearestNeighbor::KDTreeQueryResultsTags(CKDTree &kdt,CRowInt &tags) { KDTreeTsQueryResultsTags(kdt,kdt.m_innerbuf,tags); } //+------------------------------------------------------------------+ //| Distances from last query | //| INPUT PARAMETERS | //| KDT - KD-tree | //| R - possibly pre-allocated buffer. If X is too small | //| to store result, it is resized. If size(X) is | //| enough to store result, it is left unchanged. | //| OUTPUT PARAMETERS | //| R - filled with distances (in corresponding norm) | //| NOTES | //| 1. points are ordered by distance from the query point (first | //| = closest) | //| 2. if XY is larger than required to store result, only leading | //| part will be overwritten; trailing part will be left | //| unchanged. So if on input XY = [[A,B],[C,D]], and result is | //| [1,2], then on exit we will get XY = [[1,2],[C,D]]. This is | //| done purposely to increase performance; if you want function | //| to resize array according to result size, use function with | //| same name and suffix 'I'. | //| SEE ALSO | //| * KDTreeQueryResultsX() X-values | //| * KDTreeQueryResultsXY() X- and Y-values | //| * KDTreeQueryResultsTags() tag values | //+------------------------------------------------------------------+ void CNearestNeighbor::KDTreeQueryResultsDistances(CKDTree &kdt, double &r[]) { KDTreeTsQueryResultsDistances(kdt,kdt.m_innerbuf,r); } //+------------------------------------------------------------------+ //| | //+------------------------------------------------------------------+ void CNearestNeighbor::KDTreeQueryResultsDistances(CKDTree &kdt, CRowDouble &r) { KDTreeTsQueryResultsDistances(kdt,kdt.m_innerbuf,r); } //+------------------------------------------------------------------+ //| X-values from last query associated with CKDTreeRequestBuffer | //| object. | //| INPUT PARAMETERS | //| KDT - KD-tree | //| Buf - request buffer object created for this particular | //| instance of kd-tree structure. | //| X - possibly pre-allocated buffer. If X is too small to| //| store result, it is resized. If size(X) is enough | //| to store result, it is left unchanged. | //| OUTPUT PARAMETERS | //| X - rows are filled with X-values | //| NOTES | //| 1. points are ordered by distance from the query point | //| (first = closest) | //| 2. if XY is larger than required to store result, only leading | //| part will be overwritten; trailing part will be left | //| unchanged. So if on input XY = [[A,B],[C,D]], and result is | //| [1,2], then on exit we will get XY = [[1,2],[C,D]]. This is | //| done purposely to increase performance; if you want function | //| to resize array according to result size, use function with | //| same name and suffix 'I'. | //| SEE ALSO | //| * KDTreeQueryResultsXY() X- and Y-values | //| * KDTreeQueryResultsTags() tag values | //| * KDTreeQueryResultsDistances() distances | //+------------------------------------------------------------------+ void CNearestNeighbor::KDTreeTsQueryResultsX(CKDTree &kdt, CKDTreeRequestBuffer &buf, CMatrixDouble &x) { int k=0; int i1_=0; //--- check if(buf.m_kcur==0) return; if((int)x.Rows()=0,__FUNCTION__+": incorrect node")) return; //--- check if(!CAp::Assert(node0) { //--- Leaf node nodetype=0; return; } if(kdt.m_nodes[node]==0) { //--- Split node nodetype=1; return; } //--- check if(!CAp::Assert(false,__FUNCTION__+": integrity check failure")) return; } //+------------------------------------------------------------------+ //| It is informational function which allows to get information | //| about leaf node. Node index is given by integer value, with 0 | //| corresponding to root node and other node indexes obtained via | //| exploration. | //| You should not expect that serialization/unserialization will | //| retain node indexes. You should keep in mind that future versions| //| of ALGLIB may introduce new node types. | //| OUTPUT VALUES: | //| XT - output buffer is reallocated (if too small) and filled by| //| XY values | //| K - number of rows in XY | //+------------------------------------------------------------------+ void CNearestNeighbor::KDTreeExploreLeaf(CKDTree &kdt, int node, CMatrixDouble &xy, int &k) { int offs=0; k=0; //--- check if(!CAp::Assert(node>=0,__FUNCTION__+": incorrect node index")) return; //--- check if(!CAp::Assert(node+10,__FUNCTION__+": incorrect node index")) return; k=kdt.m_nodes[node]; offs=kdt.m_nodes[node+1]; //--- check if(!CAp::Assert(offs>=0,__FUNCTION__+": integrity error")) return; //--- check if(!CAp::Assert((offs+k-1)<(int)CAp::Rows(kdt.m_xy),__FUNCTION__+": integrity error")) return; CApServ::RMatrixSetLengthAtLeast(xy,k,kdt.m_nx+kdt.m_ny); for(int i=0; i=0,__FUNCTION__+": incorrect node index")) return; //--- check if(!CAp::Assert(node+4=0,__FUNCTION__+": integrity failure")) return; //--- check if(!CAp::Assert(d=0,__FUNCTION__+": integrity failure")) return; //--- check if(!CAp::Assert(nodele=0,__FUNCTION__+": integrity failure")) return; //--- check if(!CAp::Assert(nodege0 | //| SelfMatch - whether self-matches are allowed: | //| * if True, nearest neighbor may be the point itself (if | //| it exists in original dataset) | //| * if False, then only points with non-zero distance are | //| returned | //| * if not given, considered True | //| RESULT | //| number of neighbors found, >=0 | //| This subroutine performs query and stores its result in the | //| internal structures of the buffer object. You can use following | //| subroutines to obtain these results (pay attention to "buf" in | //| their names): | //| * KDTreeTsQueryResultsX() to get X-values | //| * KDTreeTsQueryResultsXY() to get X- and Y-values | //| * KDTreeTsQueryResultsTags() to get tag values | //| * KDTreeTsQueryResultsDistances() to get distances | //| IMPORTANT: kd-tree buffer should be used only with KD-tree object| //| which was used to initialize buffer. Any attempt to | //| use biffer with different object is dangerous - you | //| may get integrity check failure (exception) because | //| sizes of internal arrays do not fit to dimensions of | //| KD-tree structure. | //+------------------------------------------------------------------+ int CNearestNeighbor::TsQueryRNN(CKDTree &kdt, CKDTreeRequestBuffer &buf, CRowDouble &x, double r, bool selfmatch=true, bool orderedbydist=true) { int result=0; int i=0; int j=0; //--- Handle special case: KDT.N=0 if(kdt.m_n==0) { buf.m_kcur=0; return(0); } //--- Check consistency of request buffer CheckRequestBufferConsistency(kdt,buf); //--- Prepare parameters buf.m_kneeded=0; if(kdt.m_normtype!=2) buf.m_rneeded=r; else buf.m_rneeded=CMath::Sqr(r); buf.m_selfmatch=selfmatch; buf.m_approxf=1; buf.m_kcur=0; //--- calculate distance from point to current bounding box KDTreeInitBox(kdt,x,buf); //--- call recursive search //--- results are returned as heap KDTreeQueryNNRec(kdt,buf,0); result=buf.m_kcur; //--- pop from heap to generate ordered representation //--- last element is not pop'ed because it is already in its place if(orderedbydist) { j=buf.m_kcur; for(i=buf.m_kcur; i>=2; i--) CTSort::TagHeapPopI(buf.m_r,buf.m_idx,j); } return(result); } //+------------------------------------------------------------------+ //| Rearranges nodes [I1,I2) using partition in D-th dimension with | //| S as threshold. Returns split position I3: [I1,I3) and [I3,I2) | //| are created as result. | //| This subroutine doesn't create tree structures, just rearranges | //| nodes. | //+------------------------------------------------------------------+ void CNearestNeighbor::KDTreeSplit(CKDTree &kdt,const int i1, const int i2,const int d, const double s,int &i3) { int ileft=0; int iright=0; double v=0; //--- initialization i3=0; //--- split XY/Tags in two parts: //--- * [ILeft,IRight] is non-processed part of XY/Tags //--- After cycle is done, we have Ileft=IRight. We deal with //--- this element separately. //--- After this, [I1,ILeft) contains left part, and [ILeft,I2) //--- contains right part. ileft=i1; iright=i2-1; while(ileft0,__FUNCTION__+": internal error")) return; //--- check if(!CAp::Assert(i2>i1,__FUNCTION__+": internal error")) return; //--- Generate leaf if needed if(i2-i1<=maxleafsize) { kdt.m_nodes.Set(nodesoffs+0,i2-i1); kdt.m_nodes.Set(nodesoffs+1,i1); nodesoffs=nodesoffs+2; //--- exit the function return; } //--- Load values for easier access nx=kdt.m_nx; ny=kdt.m_ny; //--- select dimension to split: //--- * D is a dimension number //--- In case bounding box has zero size, we enforce creation of the leaf node. d=0; ds=kdt.m_innerbuf.m_curboxmax[0]-kdt.m_innerbuf.m_curboxmin[0]; for(i=1; ids) { ds=v; d=i; } } if((double)(ds)==0.0) { kdt.m_nodes.Set(nodesoffs+0,i2-i1); kdt.m_nodes.Set(nodesoffs+1,i1); nodesoffs=nodesoffs+2; return; } //--- Select split position S using sliding midpoint rule, //--- rearrange points into [I1,I3) and [I3,I2) s=kdt.m_innerbuf.m_curboxmin[d]+0.5*ds; i1_=i1; for(i_=0; i_maxv) { maxv=v; maxidx=i1+i; } //--- check if(vs) cntgreater=cntgreater+1; } //--- check if(minv==maxv) { //--- In case all points has same value of D-th component //--- (MinV=MaxV) we enforce D-th dimension of bounding //--- box to become exactly zero and repeat tree construction. double v0=kdt.m_innerbuf.m_curboxmin[d]; double v1=kdt.m_innerbuf.m_curboxmax[d]; kdt.m_innerbuf.m_curboxmin.Set(d,minv); kdt.m_innerbuf.m_curboxmax.Set(d,maxv); KDTreeGenerateTreeRec(kdt,nodesoffs,splitsoffs,i1,i2,maxleafsize); kdt.m_innerbuf.m_curboxmin.Set(d,v0); kdt.m_innerbuf.m_curboxmax.Set(d,v1); return; } //--- check if(cntless>0 && cntgreater>0) { //--- normal midpoint split KDTreeSplit(kdt,i1,i2,d,s,i3); } else { //--- sliding midpoint if(cntless==0) { //--- 1. move split to MinV, //--- 2. place one point to the left bin (move to I1), //--- others - to the right bin s=minv; //--- check if(minidx!=i1) { for(i=0; i<(2*kdt.m_nx+kdt.m_ny); i++) { v=kdt.m_xy[minidx][i]; kdt.m_xy.Set(minidx,i,kdt.m_xy[i1][i]); kdt.m_xy.Set(i1,i,v); } //--- change values kdt.m_tags.Swap(minidx,i1); } i3=i1+1; } else { //--- 1. move split to MaxV, //--- 2. place one point to the right bin (move to I2-1), //--- others - to the left bin s=maxv; //--- check if(maxidx!=i2-1) { for(i=0; i<(2*kdt.m_nx+kdt.m_ny); i++) { v=kdt.m_xy[maxidx][i]; kdt.m_xy.Set(maxidx,i,kdt.m_xy[i2-1][i]); kdt.m_xy.Set(i2-1,i,v); } //--- change values kdt.m_tags.Swap(maxidx,i2-1); } i3=i2-1; } } //--- Generate 'split' node kdt.m_nodes.Set(nodesoffs,0); kdt.m_nodes.Set(nodesoffs+1,d); kdt.m_nodes.Set(nodesoffs+2,splitsoffs); kdt.m_splits.Set(splitsoffs,s); oldoffs=nodesoffs; nodesoffs+=m_splitnodesize; splitsoffs ++; //--- Recirsive generation: //--- * update CurBox //--- * call subroutine //--- * restore CurBox kdt.m_nodes.Set(oldoffs+3,nodesoffs); v=kdt.m_innerbuf.m_curboxmax[d]; kdt.m_innerbuf.m_curboxmax.Set(d,s); //--- function call KDTreeGenerateTreeRec(kdt,nodesoffs,splitsoffs,i1,i3,maxleafsize); kdt.m_innerbuf.m_curboxmax.Set(d,v); kdt.m_nodes.Set(oldoffs+4,nodesoffs); v=kdt.m_innerbuf.m_curboxmin[d]; kdt.m_innerbuf.m_curboxmin.Set(d,s); //--- function call KDTreeGenerateTreeRec(kdt,nodesoffs,splitsoffs,i3,i2,maxleafsize); kdt.m_innerbuf.m_curboxmin.Set(d,v); //--- Zero-fill unused portions of the node (avoid false warnings by Valgrind //--- about attempt to serialize uninitialized values) CAp::Assert(CNearestNeighbor::m_splitnodesize==6,__FUNCTION__+": node size has unexpectedly changed"); kdt.m_nodes.Set(oldoffs+5,0); } //+------------------------------------------------------------------+ //| Recursive subroutine for NN queries. | //+------------------------------------------------------------------+ void CNearestNeighbor::KDTreeQueryNNRec(CKDTree &kdt,const int offs) { KDTreeQueryNNRec(kdt,kdt.m_innerbuf,offs); } //+------------------------------------------------------------------+ //| Recursive subroutine for NN queries. | //+------------------------------------------------------------------+ void CNearestNeighbor::KDTreeQueryNNRec(CKDTree &kdt,CKDTreeRequestBuffer &buf,const int offs) { //--- create variables double ptdist=0; int i=0; int j=0; int nx=0; int i1=0; int i2=0; int d=0; double s=0; double v=0; double t1=0; int childbestoffs=0; int childworstoffs=0; int childoffs=0; double prevdist=0; bool todive; bool bestisleft; bool updatemin; //--- check if(!CAp::Assert(kdt.m_n>0,__FUNCTION__+": internal error")) return; //--- Leaf node. //--- Process points. if(kdt.m_nodes[offs]>0) { i1=kdt.m_nodes[offs+1]; i2=i1+kdt.m_nodes[offs]; for(i=i1; i0) AND (PtDist>R). if(buf.m_rneeded==0.0 || ptdist<=buf.m_rneeded) { //--- R-criterion is satisfied, we must either: //--- * replace worst point, if (KNeeded<>0) AND (KCur=KNeeded) //--- (or skip, if worst point is better) //--- * add point without replacement otherwise if(buf.m_kcur=s) switch(kdt.m_normtype) { case 0: buf.m_curdist=MathMax(buf.m_curdist,t1-s); break; case 1: buf.m_curdist=buf.m_curdist-MathMax(t1-v,0)+t1-s; break; case 2: buf.m_curdist=buf.m_curdist-CMath::Sqr(MathMax(t1-v,0))+CMath::Sqr(t1-s); break; } buf.m_curboxmax.Set(d,s); } //--- Decide: to dive into cell or not to dive if(buf.m_rneeded!=0.0 && buf.m_curdist>buf.m_rneeded) todive=false; else { //--- check if(buf.m_kcur0,__FUNCTION__+": internal error")) return; nx=kdt.m_nx; //--- Check that intersection of query box with bounding box is non-empty. //--- This check is performed once for Offs=0 (tree root). if(offs==0) { for(j=0; jbuf.m_curboxmax[j]) return; if(buf.m_boxmax[j]0) { i1=kdt.m_nodes[offs+1]; i2=i1+kdt.m_nodes[offs]; for(i=i1; i=buf.m_boxmin[j]); inbox=inbox && (kdt.m_xy[i][j]<=buf.m_boxmax[j]); } if(!inbox) continue; //--- Add point to unordered list buf.m_r.Set(buf.m_kcur,0.0); buf.m_idx.Set(buf.m_kcur,i); buf.m_kcur++; } return; } //--- Simple split if(kdt.m_nodes[offs]==0) { //--- Load: //--- * D dimension to split //--- * S split position d=kdt.m_nodes[offs+1]; s=kdt.m_splits[kdt.m_nodes[offs+2]]; //--- Check lower split (S is upper bound of new bounding box) if(s>=buf.m_boxmin[d]) { v=buf.m_curboxmax[d]; buf.m_curboxmax.Set(d,s); KDTreeQueryBoxRec(kdt,buf,kdt.m_nodes[offs+3]); buf.m_curboxmax.Set(d,v); } //--- Check upper split (S is lower bound of new bounding box) if(s<=buf.m_boxmax[d]) { v=buf.m_curboxmin[d]; buf.m_curboxmin.Set(d,s); KDTreeQueryBoxRec(kdt,buf,kdt.m_nodes[offs+4]); buf.m_curboxmin.Set(d,v); } return; } } //+------------------------------------------------------------------+ //| Copies X[] to KDT.X[] | //| Loads distance from X[] to bounding box. | //| Initializes CurBox[]. | //+------------------------------------------------------------------+ void CNearestNeighbor::KDTreeInitBox(CKDTree&kdt,CRowDouble &x, CKDTreeRequestBuffer&buf) { int i=0; double vx=0; double vmin=0; double vmax=0; //--- check if(!CAp::Assert(kdt.m_n>0,__FUNCTION__+": Internal error")) return; //--- calculate distance from point to current bounding box buf.m_curdist=0; if(kdt.m_normtype==0) { for(i=0; ivmax) buf.m_curdist=MathMax(buf.m_curdist,vx-vmax); } } } //--- check if(kdt.m_normtype==1) { for(i=0; ivmax) buf.m_curdist=buf.m_curdist+vx-vmax; } } } //--- check if(kdt.m_normtype==2) { for(i=0; ivmax) buf.m_curdist=buf.m_curdist+CMath::Sqr(vx-vmax); } } } } //+------------------------------------------------------------------+ //| This function allocates all dataset-independent array fields of | //| KDTree, i.e. such array fields that their dimensions do not | //| depend on dataset size. | //| This function do not sets KDT.NX or KDT.NY - it just allocates | //| arrays | //+------------------------------------------------------------------+ void CNearestNeighbor::KDTreeAllocDataSetIndependent(CKDTree&kdt, const int nx, const int ny) { //--- check if(!CAp::Assert(kdt.m_n>0,__FUNCTION__+": internal error")) return; //--- allocation kdt.m_boxmin.Resize(nx); kdt.m_boxmax.Resize(nx); } //+------------------------------------------------------------------+ //| This function allocates all dataset-dependent array fields of | //| KDTree, i.e. such array fields that their dimensions depend on | //| dataset size. | //| This function do not sets KDT.N, KDT.NX or KDT.NY - | //| it just allocates arrays. | //+------------------------------------------------------------------+ void CNearestNeighbor::KDTreeAllocDataSetDependent(CKDTree&kdt, const int n, const int nx, const int ny) { //--- check if(!CAp::Assert(n>0,__FUNCTION__+": internal error")) return; //--- allocation kdt.m_xy.Resize(n,2*nx+ny); kdt.m_tags.Resize(n); kdt.m_nodes.Resize(m_splitnodesize*2*n); kdt.m_splits.Resize(2*n); } //+------------------------------------------------------------------+ //| This function checks consistency of request buffer structure with| //| dimensions of kd-tree object. | //+------------------------------------------------------------------+ bool CNearestNeighbor::CheckRequestBufferConsistency(CKDTree&kdt,CKDTreeRequestBuffer&buf) { if(!CAp::Assert((int)buf.m_x.Size()>=kdt.m_nx,__FUNCTION__+": dimensions of CKDTreeRequestBuffer are inconsistent with CKDTree structure")) return(false); if(!CAp::Assert(CAp::Len(buf.m_idx)>=kdt.m_n,__FUNCTION__+": dimensions of CKDTreeRequestBuffer are inconsistent with CKDTree structure")) return(false); if(!CAp::Assert((int)buf.m_r.Size()>=kdt.m_n,__FUNCTION__+": dimensions of CKDTreeRequestBuffer are inconsistent with CKDTree structure")) return(false); if(!CAp::Assert((int)buf.m_buf.Size()>=MathMax(kdt.m_n,kdt.m_nx),__FUNCTION__+": dimensions of CKDTreeRequestBuffer are inconsistent with CKDTree structure")) return(false); if(!CAp::Assert((int)buf.m_curboxmin.Size()>=kdt.m_nx,__FUNCTION__+": dimensions of CKDTreeRequestBuffer are inconsistent with CKDTree structure")) return(false); if(!CAp::Assert((int)buf.m_curboxmax.Size()>=kdt.m_nx,__FUNCTION__+": dimensions of CKDTreeRequestBuffer are inconsistent with CKDTree structure")) return(false); return(true); } //+------------------------------------------------------------------+