1  /*************************************************************************


2  Copyright (c) Sergey Bochkanov (ALGLIB project).


3 


4  >>> SOURCE LICENSE >>>


5  This program is free software; you can redistribute it and/or modify


6  it under the terms of the GNU General Public License as published by


7  the Free Software Foundation (www.fsf.org); either version 2 of the


8  License, or (at your option) any later version.


9 


10  This program is distributed in the hope that it will be useful,


11  but WITHOUT ANY WARRANTY; without even the implied warranty of


12  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the


13  GNU General Public License for more details.


14 


15  A copy of the GNU General Public License is available at


16  http://www.fsf.org/licensing/licenses


17  >>> END OF LICENSE >>>


18  *************************************************************************/


19  #pragma warning disable 162


20  #pragma warning disable 219


21  using System;


22 


23  public partial class alglib


24  {


25 


26 


27  /*************************************************************************


28  Portable high quality random number generator state.


29  Initialized with HQRNDRandomize() or HQRNDSeed().


30 


31  Fields:


32  S1, S2  seed values


33  V  precomputed value


34  MagicV  'magic' value used to determine whether State structure


35  was correctly initialized.


36  *************************************************************************/


37  public class hqrndstate


38  {


39  //


40  // Public declarations


41  //


42 


43  public hqrndstate()


44  {


45  _innerobj = new hqrnd.hqrndstate();


46  }


47 


48  //


49  // Although some of declarations below are public, you should not use them


50  // They are intended for internal use only


51  //


52  private hqrnd.hqrndstate _innerobj;


53  public hqrnd.hqrndstate innerobj { get { return _innerobj; } }


54  public hqrndstate(hqrnd.hqrndstate obj)


55  {


56  _innerobj = obj;


57  }


58  }


59 


60  /*************************************************************************


61  HQRNDState initialization with random values which come from standard


62  RNG.


63 


64   ALGLIB 


65  Copyright 02.12.2009 by Bochkanov Sergey


66  *************************************************************************/


67  public static void hqrndrandomize(out hqrndstate state)


68  {


69  state = new hqrndstate();


70  hqrnd.hqrndrandomize(state.innerobj);


71  return;


72  }


73 


74  /*************************************************************************


75  HQRNDState initialization with seed values


76 


77   ALGLIB 


78  Copyright 02.12.2009 by Bochkanov Sergey


79  *************************************************************************/


80  public static void hqrndseed(int s1, int s2, out hqrndstate state)


81  {


82  state = new hqrndstate();


83  hqrnd.hqrndseed(s1, s2, state.innerobj);


84  return;


85  }


86 


87  /*************************************************************************


88  This function generates random real number in (0,1),


89  not including interval boundaries


90 


91  State structure must be initialized with HQRNDRandomize() or HQRNDSeed().


92 


93   ALGLIB 


94  Copyright 02.12.2009 by Bochkanov Sergey


95  *************************************************************************/


96  public static double hqrnduniformr(hqrndstate state)


97  {


98 


99  double result = hqrnd.hqrnduniformr(state.innerobj);


100  return result;


101  }


102 


103  /*************************************************************************


104  This function generates random integer number in [0, N)


105 


106  1. N must be less than HQRNDMax1.


107  2. State structure must be initialized with HQRNDRandomize() or HQRNDSeed()


108 


109   ALGLIB 


110  Copyright 02.12.2009 by Bochkanov Sergey


111  *************************************************************************/


112  public static int hqrnduniformi(hqrndstate state, int n)


113  {


114 


115  int result = hqrnd.hqrnduniformi(state.innerobj, n);


116  return result;


117  }


118 


119  /*************************************************************************


120  Random number generator: normal numbers


121 


122  This function generates one random number from normal distribution.


123  Its performance is equal to that of HQRNDNormal2()


124 


125  State structure must be initialized with HQRNDRandomize() or HQRNDSeed().


126 


127   ALGLIB 


128  Copyright 02.12.2009 by Bochkanov Sergey


129  *************************************************************************/


130  public static double hqrndnormal(hqrndstate state)


131  {


132 


133  double result = hqrnd.hqrndnormal(state.innerobj);


134  return result;


135  }


136 


137  /*************************************************************************


138  Random number generator: random X and Y such that X^2+Y^2=1


139 


140  State structure must be initialized with HQRNDRandomize() or HQRNDSeed().


141 


142   ALGLIB 


143  Copyright 02.12.2009 by Bochkanov Sergey


144  *************************************************************************/


145  public static void hqrndunit2(hqrndstate state, out double x, out double y)


146  {


147  x = 0;


148  y = 0;


149  hqrnd.hqrndunit2(state.innerobj, ref x, ref y);


150  return;


151  }


152 


153  /*************************************************************************


154  Random number generator: normal numbers


155 


156  This function generates two independent random numbers from normal


157  distribution. Its performance is equal to that of HQRNDNormal()


158 


159  State structure must be initialized with HQRNDRandomize() or HQRNDSeed().


160 


161   ALGLIB 


162  Copyright 02.12.2009 by Bochkanov Sergey


163  *************************************************************************/


164  public static void hqrndnormal2(hqrndstate state, out double x1, out double x2)


165  {


166  x1 = 0;


167  x2 = 0;


168  hqrnd.hqrndnormal2(state.innerobj, ref x1, ref x2);


169  return;


170  }


171 


172  /*************************************************************************


173  Random number generator: exponential distribution


174 


175  State structure must be initialized with HQRNDRandomize() or HQRNDSeed().


176 


177   ALGLIB 


178  Copyright 11.08.2007 by Bochkanov Sergey


179  *************************************************************************/


180  public static double hqrndexponential(hqrndstate state, double lambdav)


181  {


182 


183  double result = hqrnd.hqrndexponential(state.innerobj, lambdav);


184  return result;


185  }


186 


187  /*************************************************************************


188  This function generates random number from discrete distribution given by


189  finite sample X.


190 


191  INPUT PARAMETERS


192  State  high quality random number generator, must be


193  initialized with HQRNDRandomize() or HQRNDSeed().


194  X  finite sample


195  N  number of elements to use, N>=1


196 


197  RESULT


198  this function returns one of the X[i] for random i=0..N1


199 


200   ALGLIB 


201  Copyright 08.11.2011 by Bochkanov Sergey


202  *************************************************************************/


203  public static double hqrnddiscrete(hqrndstate state, double[] x, int n)


204  {


205 


206  double result = hqrnd.hqrnddiscrete(state.innerobj, x, n);


207  return result;


208  }


209 


210  /*************************************************************************


211  This function generates random number from continuous distribution given


212  by finite sample X.


213 


214  INPUT PARAMETERS


215  State  high quality random number generator, must be


216  initialized with HQRNDRandomize() or HQRNDSeed().


217  X  finite sample, array[N] (can be larger, in this case only


218  leading N elements are used). THIS ARRAY MUST BE SORTED BY


219  ASCENDING.


220  N  number of elements to use, N>=1


221 


222  RESULT


223  this function returns random number from continuous distribution which


224  tries to approximate X as mush as possible. min(X)<=Result<=max(X).


225 


226   ALGLIB 


227  Copyright 08.11.2011 by Bochkanov Sergey


228  *************************************************************************/


229  public static double hqrndcontinuous(hqrndstate state, double[] x, int n)


230  {


231 


232  double result = hqrnd.hqrndcontinuous(state.innerobj, x, n);


233  return result;


234  }


235 


236  }


237  public partial class alglib


238  {


239 


240 


241  /*************************************************************************


242 


243  *************************************************************************/


244  public class kdtree


245  {


246  //


247  // Public declarations


248  //


249 


250  public kdtree()


251  {


252  _innerobj = new nearestneighbor.kdtree();


253  }


254 


255  //


256  // Although some of declarations below are public, you should not use them


257  // They are intended for internal use only


258  //


259  private nearestneighbor.kdtree _innerobj;


260  public nearestneighbor.kdtree innerobj { get { return _innerobj; } }


261  public kdtree(nearestneighbor.kdtree obj)


262  {


263  _innerobj = obj;


264  }


265  }


266 


267 


268  /*************************************************************************


269  This function serializes data structure to string.


270 


271  Important properties of s_out:


272  * it contains alphanumeric characters, dots, underscores, minus signs


273  * these symbols are grouped into words, which are separated by spaces


274  and Windowsstyle (CR+LF) newlines


275  * although serializer uses spaces and CR+LF as separators, you can


276  replace any separator character by arbitrary combination of spaces,


277  tabs, Windows or Unix newlines. It allows flexible reformatting of


278  the string in case you want to include it into text or XML file.


279  But you should not insert separators into the middle of the "words"


280  nor you should change case of letters.


281  * s_out can be freely moved between 32bit and 64bit systems, little


282  and big endian machines, and so on. You can serialize structure on


283  32bit machine and unserialize it on 64bit one (or vice versa), or


284  serialize it on SPARC and unserialize on x86. You can also


285  serialize it in C# version of ALGLIB and unserialize in C++ one,


286  and vice versa.


287  *************************************************************************/


288  public static void kdtreeserialize(kdtree obj, out string s_out)


289  {


290  alglib.serializer s = new alglib.serializer();


291  s.alloc_start();


292  nearestneighbor.kdtreealloc(s, obj.innerobj);


293  s.sstart_str();


294  nearestneighbor.kdtreeserialize(s, obj.innerobj);


295  s.stop();


296  s_out = s.get_string();


297  }


298 


299 


300  /*************************************************************************


301  This function unserializes data structure from string.


302  *************************************************************************/


303  public static void kdtreeunserialize(string s_in, out kdtree obj)


304  {


305  alglib.serializer s = new alglib.serializer();


306  obj = new kdtree();


307  s.ustart_str(s_in);


308  nearestneighbor.kdtreeunserialize(s, obj.innerobj);


309  s.stop();


310  }


311 


312  /*************************************************************************


313  KDtree creation


314 


315  This subroutine creates KDtree from set of Xvalues and optional Yvalues


316 


317  INPUT PARAMETERS


318  XY  dataset, array[0..N1,0..NX+NY1].


319  one row corresponds to one point.


320  first NX columns contain Xvalues, next NY (NY may be zero)


321  columns may contain associated Yvalues


322  N  number of points, N>=0.


323  NX  space dimension, NX>=1.


324  NY  number of optional Yvalues, NY>=0.


325  NormType norm type:


326  * 0 denotes infinitynorm


327  * 1 denotes 1norm


328  * 2 denotes 2norm (Euclidean norm)


329 


330  OUTPUT PARAMETERS


331  KDT  KDtree


332 


333 


334  NOTES


335 


336  1. KDtree creation have O(N*logN) complexity and O(N*(2*NX+NY)) memory


337  requirements.


338  2. Although KDtrees may be used with any combination of N and NX, they


339  are more efficient than bruteforce search only when N >> 4^NX. So they


340  are most useful in lowdimensional tasks (NX=2, NX=3). NX=1 is another


341  inefficient case, because simple binary search (without additional


342  structures) is much more efficient in such tasks than KDtrees.


343 


344   ALGLIB 


345  Copyright 28.02.2010 by Bochkanov Sergey


346  *************************************************************************/


347  public static void kdtreebuild(double[,] xy, int n, int nx, int ny, int normtype, out kdtree kdt)


348  {


349  kdt = new kdtree();


350  nearestneighbor.kdtreebuild(xy, n, nx, ny, normtype, kdt.innerobj);


351  return;


352  }


353  public static void kdtreebuild(double[,] xy, int nx, int ny, int normtype, out kdtree kdt)


354  {


355  int n;


356 


357  kdt = new kdtree();


358  n = ap.rows(xy);


359  nearestneighbor.kdtreebuild(xy, n, nx, ny, normtype, kdt.innerobj);


360 


361  return;


362  }


363 


364  /*************************************************************************


365  KDtree creation


366 


367  This subroutine creates KDtree from set of Xvalues, integer tags and


368  optional Yvalues


369 


370  INPUT PARAMETERS


371  XY  dataset, array[0..N1,0..NX+NY1].


372  one row corresponds to one point.


373  first NX columns contain Xvalues, next NY (NY may be zero)


374  columns may contain associated Yvalues


375  Tags  tags, array[0..N1], contains integer tags associated


376  with points.


377  N  number of points, N>=0


378  NX  space dimension, NX>=1.


379  NY  number of optional Yvalues, NY>=0.


380  NormType norm type:


381  * 0 denotes infinitynorm


382  * 1 denotes 1norm


383  * 2 denotes 2norm (Euclidean norm)


384 


385  OUTPUT PARAMETERS


386  KDT  KDtree


387 


388  NOTES


389 


390  1. KDtree creation have O(N*logN) complexity and O(N*(2*NX+NY)) memory


391  requirements.


392  2. Although KDtrees may be used with any combination of N and NX, they


393  are more efficient than bruteforce search only when N >> 4^NX. So they


394  are most useful in lowdimensional tasks (NX=2, NX=3). NX=1 is another


395  inefficient case, because simple binary search (without additional


396  structures) is much more efficient in such tasks than KDtrees.


397 


398   ALGLIB 


399  Copyright 28.02.2010 by Bochkanov Sergey


400  *************************************************************************/


401  public static void kdtreebuildtagged(double[,] xy, int[] tags, int n, int nx, int ny, int normtype, out kdtree kdt)


402  {


403  kdt = new kdtree();


404  nearestneighbor.kdtreebuildtagged(xy, tags, n, nx, ny, normtype, kdt.innerobj);


405  return;


406  }


407  public static void kdtreebuildtagged(double[,] xy, int[] tags, int nx, int ny, int normtype, out kdtree kdt)


408  {


409  int n;


410  if( (ap.rows(xy)!=ap.len(tags)))


411  throw new alglibexception("Error while calling 'kdtreebuildtagged': looks like one of arguments has wrong size");


412  kdt = new kdtree();


413  n = ap.rows(xy);


414  nearestneighbor.kdtreebuildtagged(xy, tags, n, nx, ny, normtype, kdt.innerobj);


415 


416  return;


417  }


418 


419  /*************************************************************************


420  KNN query: K nearest neighbors


421 


422  INPUT PARAMETERS


423  KDT  KDtree


424  X  point, array[0..NX1].


425  K  number of neighbors to return, K>=1


426  SelfMatch  whether selfmatches are allowed:


427  * if True, nearest neighbor may be the point itself


428  (if it exists in original dataset)


429  * if False, then only points with nonzero distance


430  are returned


431  * if not given, considered True


432 


433  RESULT


434  number of actual neighbors found (either K or N, if K>N).


435 


436  This subroutine performs query and stores its result in the internal


437  structures of the KDtree. You can use following subroutines to obtain


438  these results:


439  * KDTreeQueryResultsX() to get Xvalues


440  * KDTreeQueryResultsXY() to get X and Yvalues


441  * KDTreeQueryResultsTags() to get tag values


442  * KDTreeQueryResultsDistances() to get distances


443 


444   ALGLIB 


445  Copyright 28.02.2010 by Bochkanov Sergey


446  *************************************************************************/


447  public static int kdtreequeryknn(kdtree kdt, double[] x, int k, bool selfmatch)


448  {


449 


450  int result = nearestneighbor.kdtreequeryknn(kdt.innerobj, x, k, selfmatch);


451  return result;


452  }


453  public static int kdtreequeryknn(kdtree kdt, double[] x, int k)


454  {


455  bool selfmatch;


456 


457 


458  selfmatch = true;


459  int result = nearestneighbor.kdtreequeryknn(kdt.innerobj, x, k, selfmatch);


460 


461  return result;


462  }


463 


464  /*************************************************************************


465  RNN query: all points within Rsphere centered at X


466 


467  INPUT PARAMETERS


468  KDT  KDtree


469  X  point, array[0..NX1].


470  R  radius of sphere (in corresponding norm), R>0


471  SelfMatch  whether selfmatches are allowed:


472  * if True, nearest neighbor may be the point itself


473  (if it exists in original dataset)


474  * if False, then only points with nonzero distance


475  are returned


476  * if not given, considered True


477 


478  RESULT


479  number of neighbors found, >=0


480 


481  This subroutine performs query and stores its result in the internal


482  structures of the KDtree. You can use following subroutines to obtain


483  actual results:


484  * KDTreeQueryResultsX() to get Xvalues


485  * KDTreeQueryResultsXY() to get X and Yvalues


486  * KDTreeQueryResultsTags() to get tag values


487  * KDTreeQueryResultsDistances() to get distances


488 


489   ALGLIB 


490  Copyright 28.02.2010 by Bochkanov Sergey


491  *************************************************************************/


492  public static int kdtreequeryrnn(kdtree kdt, double[] x, double r, bool selfmatch)


493  {


494 


495  int result = nearestneighbor.kdtreequeryrnn(kdt.innerobj, x, r, selfmatch);


496  return result;


497  }


498  public static int kdtreequeryrnn(kdtree kdt, double[] x, double r)


499  {


500  bool selfmatch;


501 


502 


503  selfmatch = true;


504  int result = nearestneighbor.kdtreequeryrnn(kdt.innerobj, x, r, selfmatch);


505 


506  return result;


507  }


508 


509  /*************************************************************************


510  KNN query: approximate K nearest neighbors


511 


512  INPUT PARAMETERS


513  KDT  KDtree


514  X  point, array[0..NX1].


515  K  number of neighbors to return, K>=1


516  SelfMatch  whether selfmatches are allowed:


517  * if True, nearest neighbor may be the point itself


518  (if it exists in original dataset)


519  * if False, then only points with nonzero distance


520  are returned


521  * if not given, considered True


522  Eps  approximation factor, Eps>=0. epsapproximate nearest


523  neighbor is a neighbor whose distance from X is at


524  most (1+eps) times distance of true nearest neighbor.


525 


526  RESULT


527  number of actual neighbors found (either K or N, if K>N).


528 


529  NOTES


530  significant performance gain may be achieved only when Eps is is on


531  the order of magnitude of 1 or larger.


532 


533  This subroutine performs query and stores its result in the internal


534  structures of the KDtree. You can use following subroutines to obtain


535  these results:


536  * KDTreeQueryResultsX() to get Xvalues


537  * KDTreeQueryResultsXY() to get X and Yvalues


538  * KDTreeQueryResultsTags() to get tag values


539  * KDTreeQueryResultsDistances() to get distances


540 


541   ALGLIB 


542  Copyright 28.02.2010 by Bochkanov Sergey


543  *************************************************************************/


544  public static int kdtreequeryaknn(kdtree kdt, double[] x, int k, bool selfmatch, double eps)


545  {


546 


547  int result = nearestneighbor.kdtreequeryaknn(kdt.innerobj, x, k, selfmatch, eps);


548  return result;


549  }


550  public static int kdtreequeryaknn(kdtree kdt, double[] x, int k, double eps)


551  {


552  bool selfmatch;


553 


554 


555  selfmatch = true;


556  int result = nearestneighbor.kdtreequeryaknn(kdt.innerobj, x, k, selfmatch, eps);


557 


558  return result;


559  }


560 


561  /*************************************************************************


562  Xvalues from last query


563 


564  INPUT PARAMETERS


565  KDT  KDtree


566  X  possibly preallocated buffer. If X is too small to store


567  result, it is resized. If size(X) is enough to store


568  result, it is left unchanged.


569 


570  OUTPUT PARAMETERS


571  X  rows are filled with Xvalues


572 


573  NOTES


574  1. points are ordered by distance from the query point (first = closest)


575  2. if XY is larger than required to store result, only leading part will


576  be overwritten; trailing part will be left unchanged. So if on input


577  XY = [[A,B],[C,D]], and result is [1,2], then on exit we will get


578  XY = [[1,2],[C,D]]. This is done purposely to increase performance; if


579  you want function to resize array according to result size, use


580  function with same name and suffix 'I'.


581 


582  SEE ALSO


583  * KDTreeQueryResultsXY() X and Yvalues


584  * KDTreeQueryResultsTags() tag values


585  * KDTreeQueryResultsDistances() distances


586 


587   ALGLIB 


588  Copyright 28.02.2010 by Bochkanov Sergey


589  *************************************************************************/


590  public static void kdtreequeryresultsx(kdtree kdt, ref double[,] x)


591  {


592 


593  nearestneighbor.kdtreequeryresultsx(kdt.innerobj, ref x);


594  return;


595  }


596 


597  /*************************************************************************


598  X and Yvalues from last query


599 


600  INPUT PARAMETERS


601  KDT  KDtree


602  XY  possibly preallocated buffer. If XY is too small to store


603  result, it is resized. If size(XY) is enough to store


604  result, it is left unchanged.


605 


606  OUTPUT PARAMETERS


607  XY  rows are filled with points: first NX columns with


608  Xvalues, next NY columns  with Yvalues.


609 


610  NOTES


611  1. points are ordered by distance from the query point (first = closest)


612  2. if XY is larger than required to store result, only leading part will


613  be overwritten; trailing part will be left unchanged. So if on input


614  XY = [[A,B],[C,D]], and result is [1,2], then on exit we will get


615  XY = [[1,2],[C,D]]. This is done purposely to increase performance; if


616  you want function to resize array according to result size, use


617  function with same name and suffix 'I'.


618 


619  SEE ALSO


620  * KDTreeQueryResultsX() Xvalues


621  * KDTreeQueryResultsTags() tag values


622  * KDTreeQueryResultsDistances() distances


623 


624   ALGLIB 


625  Copyright 28.02.2010 by Bochkanov Sergey


626  *************************************************************************/


627  public static void kdtreequeryresultsxy(kdtree kdt, ref double[,] xy)


628  {


629 


630  nearestneighbor.kdtreequeryresultsxy(kdt.innerobj, ref xy);


631  return;


632  }


633 


634  /*************************************************************************


635  Tags from last query


636 


637  INPUT PARAMETERS


638  KDT  KDtree


639  Tags  possibly preallocated buffer. If X is too small to store


640  result, it is resized. If size(X) is enough to store


641  result, it is left unchanged.


642 


643  OUTPUT PARAMETERS


644  Tags  filled with tags associated with points,


645  or, when no tags were supplied, with zeros


646 


647  NOTES


648  1. points are ordered by distance from the query point (first = closest)


649  2. if XY is larger than required to store result, only leading part will


650  be overwritten; trailing part will be left unchanged. So if on input


651  XY = [[A,B],[C,D]], and result is [1,2], then on exit we will get


652  XY = [[1,2],[C,D]]. This is done purposely to increase performance; if


653  you want function to resize array according to result size, use


654  function with same name and suffix 'I'.


655 


656  SEE ALSO


657  * KDTreeQueryResultsX() Xvalues


658  * KDTreeQueryResultsXY() X and Yvalues


659  * KDTreeQueryResultsDistances() distances


660 


661   ALGLIB 


662  Copyright 28.02.2010 by Bochkanov Sergey


663  *************************************************************************/


664  public static void kdtreequeryresultstags(kdtree kdt, ref int[] tags)


665  {


666 


667  nearestneighbor.kdtreequeryresultstags(kdt.innerobj, ref tags);


668  return;


669  }


670 


671  /*************************************************************************


672  Distances from last query


673 


674  INPUT PARAMETERS


675  KDT  KDtree


676  R  possibly preallocated buffer. If X is too small to store


677  result, it is resized. If size(X) is enough to store


678  result, it is left unchanged.


679 


680  OUTPUT PARAMETERS


681  R  filled with distances (in corresponding norm)


682 


683  NOTES


684  1. points are ordered by distance from the query point (first = closest)


685  2. if XY is larger than required to store result, only leading part will


686  be overwritten; trailing part will be left unchanged. So if on input


687  XY = [[A,B],[C,D]], and result is [1,2], then on exit we will get


688  XY = [[1,2],[C,D]]. This is done purposely to increase performance; if


689  you want function to resize array according to result size, use


690  function with same name and suffix 'I'.


691 


692  SEE ALSO


693  * KDTreeQueryResultsX() Xvalues


694  * KDTreeQueryResultsXY() X and Yvalues


695  * KDTreeQueryResultsTags() tag values


696 


697   ALGLIB 


698  Copyright 28.02.2010 by Bochkanov Sergey


699  *************************************************************************/


700  public static void kdtreequeryresultsdistances(kdtree kdt, ref double[] r)


701  {


702 


703  nearestneighbor.kdtreequeryresultsdistances(kdt.innerobj, ref r);


704  return;


705  }


706 


707  /*************************************************************************


708  Xvalues from last query; 'interactive' variant for languages like Python


709  which support constructs like "X = KDTreeQueryResultsXI(KDT)" and


710  interactive mode of interpreter.


711 


712  This function allocates new array on each call, so it is significantly


713  slower than its 'noninteractive' counterpart, but it is more convenient


714  when you call it from command line.


715 


716   ALGLIB 


717  Copyright 28.02.2010 by Bochkanov Sergey


718  *************************************************************************/


719  public static void kdtreequeryresultsxi(kdtree kdt, out double[,] x)


720  {


721  x = new double[0,0];


722  nearestneighbor.kdtreequeryresultsxi(kdt.innerobj, ref x);


723  return;


724  }


725 


726  /*************************************************************************


727  XYvalues from last query; 'interactive' variant for languages like Python


728  which support constructs like "XY = KDTreeQueryResultsXYI(KDT)" and


729  interactive mode of interpreter.


730 


731  This function allocates new array on each call, so it is significantly


732  slower than its 'noninteractive' counterpart, but it is more convenient


733  when you call it from command line.


734 


735   ALGLIB 


736  Copyright 28.02.2010 by Bochkanov Sergey


737  *************************************************************************/


738  public static void kdtreequeryresultsxyi(kdtree kdt, out double[,] xy)


739  {


740  xy = new double[0,0];


741  nearestneighbor.kdtreequeryresultsxyi(kdt.innerobj, ref xy);


742  return;


743  }


744 


745  /*************************************************************************


746  Tags from last query; 'interactive' variant for languages like Python


747  which support constructs like "Tags = KDTreeQueryResultsTagsI(KDT)" and


748  interactive mode of interpreter.


749 


750  This function allocates new array on each call, so it is significantly


751  slower than its 'noninteractive' counterpart, but it is more convenient


752  when you call it from command line.


753 


754   ALGLIB 


755  Copyright 28.02.2010 by Bochkanov Sergey


756  *************************************************************************/


757  public static void kdtreequeryresultstagsi(kdtree kdt, out int[] tags)


758  {


759  tags = new int[0];


760  nearestneighbor.kdtreequeryresultstagsi(kdt.innerobj, ref tags);


761  return;


762  }


763 


764  /*************************************************************************


765  Distances from last query; 'interactive' variant for languages like Python


766  which support constructs like "R = KDTreeQueryResultsDistancesI(KDT)"


767  and interactive mode of interpreter.


768 


769  This function allocates new array on each call, so it is significantly


770  slower than its 'noninteractive' counterpart, but it is more convenient


771  when you call it from command line.


772 


773   ALGLIB 


774  Copyright 28.02.2010 by Bochkanov Sergey


775  *************************************************************************/


776  public static void kdtreequeryresultsdistancesi(kdtree kdt, out double[] r)


777  {


778  r = new double[0];


779  nearestneighbor.kdtreequeryresultsdistancesi(kdt.innerobj, ref r);


780  return;


781  }


782 


783  }


784  public partial class alglib


785  {


786  public class hqrnd


787  {


788  /*************************************************************************


789  Portable high quality random number generator state.


790  Initialized with HQRNDRandomize() or HQRNDSeed().


791 


792  Fields:


793  S1, S2  seed values


794  V  precomputed value


795  MagicV  'magic' value used to determine whether State structure


796  was correctly initialized.


797  *************************************************************************/


798  public class hqrndstate


799  {


800  public int s1;


801  public int s2;


802  public double v;


803  public int magicv;


804  };


805 


806 


807 


808 


809  public const int hqrndmax = 2147483563;


810  public const int hqrndm1 = 2147483563;


811  public const int hqrndm2 = 2147483399;


812  public const int hqrndmagic = 1634357784;


813 


814 


815  /*************************************************************************


816  HQRNDState initialization with random values which come from standard


817  RNG.


818 


819   ALGLIB 


820  Copyright 02.12.2009 by Bochkanov Sergey


821  *************************************************************************/


822  public static void hqrndrandomize(hqrndstate state)


823  {


824  hqrndseed(math.randominteger(hqrndm1), math.randominteger(hqrndm2), state);


825  }


826 


827 


828  /*************************************************************************


829  HQRNDState initialization with seed values


830 


831   ALGLIB 


832  Copyright 02.12.2009 by Bochkanov Sergey


833  *************************************************************************/


834  public static void hqrndseed(int s1,


835  int s2,


836  hqrndstate state)


837  {


838  state.s1 = s1%(hqrndm11)+1;


839  state.s2 = s2%(hqrndm21)+1;


840  state.v = (double)1/(double)hqrndmax;


841  state.magicv = hqrndmagic;


842  }


843 


844 


845  /*************************************************************************


846  This function generates random real number in (0,1),


847  not including interval boundaries


848 


849  State structure must be initialized with HQRNDRandomize() or HQRNDSeed().


850 


851   ALGLIB 


852  Copyright 02.12.2009 by Bochkanov Sergey


853  *************************************************************************/


854  public static double hqrnduniformr(hqrndstate state)


855  {


856  double result = 0;


857 


858  result = state.v*hqrndintegerbase(state);


859  return result;


860  }


861 


862 


863  /*************************************************************************


864  This function generates random integer number in [0, N)


865 


866  1. N must be less than HQRNDMax1.


867  2. State structure must be initialized with HQRNDRandomize() or HQRNDSeed()


868 


869   ALGLIB 


870  Copyright 02.12.2009 by Bochkanov Sergey


871  *************************************************************************/


872  public static int hqrnduniformi(hqrndstate state,


873  int n)


874  {


875  int result = 0;


876  int mx = 0;


877 


878 


879  //


880  // Correct handling of N's close to RNDBaseMax


881  // (avoiding skewed distributions for RNDBaseMax<>K*N)


882  //


883  alglib.ap.assert(n>0, "HQRNDUniformI: N<=0!");


884  alglib.ap.assert(n<hqrndmax1, "HQRNDUniformI: N>=RNDBaseMax1!");


885  mx = hqrndmax1(hqrndmax1)%n;


886  do


887  {


888  result = hqrndintegerbase(state)1;


889  }


890  while( result>=mx );


891  result = result%n;


892  return result;


893  }


894 


895 


896  /*************************************************************************


897  Random number generator: normal numbers


898 


899  This function generates one random number from normal distribution.


900  Its performance is equal to that of HQRNDNormal2()


901 


902  State structure must be initialized with HQRNDRandomize() or HQRNDSeed().


903 


904   ALGLIB 


905  Copyright 02.12.2009 by Bochkanov Sergey


906  *************************************************************************/


907  public static double hqrndnormal(hqrndstate state)


908  {


909  double result = 0;


910  double v1 = 0;


911  double v2 = 0;


912 


913  hqrndnormal2(state, ref v1, ref v2);


914  result = v1;


915  return result;


916  }


917 


918 


919  /*************************************************************************


920  Random number generator: random X and Y such that X^2+Y^2=1


921 


922  State structure must be initialized with HQRNDRandomize() or HQRNDSeed().


923 


924   ALGLIB 


925  Copyright 02.12.2009 by Bochkanov Sergey


926  *************************************************************************/


927  public static void hqrndunit2(hqrndstate state,


928  ref double x,


929  ref double y)


930  {


931  double v = 0;


932  double mx = 0;


933  double mn = 0;


934 


935  x = 0;


936  y = 0;


937 


938  do


939  {


940  hqrndnormal2(state, ref x, ref y);


941  }


942  while( !((double)(x)!=(double)(0)  (double)(y)!=(double)(0)) );


943  mx = Math.Max(Math.Abs(x), Math.Abs(y));


944  mn = Math.Min(Math.Abs(x), Math.Abs(y));


945  v = mx*Math.Sqrt(1+math.sqr(mn/mx));


946  x = x/v;


947  y = y/v;


948  }


949 


950 


951  /*************************************************************************


952  Random number generator: normal numbers


953 


954  This function generates two independent random numbers from normal


955  distribution. Its performance is equal to that of HQRNDNormal()


956 


957  State structure must be initialized with HQRNDRandomize() or HQRNDSeed().


958 


959   ALGLIB 


960  Copyright 02.12.2009 by Bochkanov Sergey


961  *************************************************************************/


962  public static void hqrndnormal2(hqrndstate state,


963  ref double x1,


964  ref double x2)


965  {


966  double u = 0;


967  double v = 0;


968  double s = 0;


969 


970  x1 = 0;


971  x2 = 0;


972 


973  while( true )


974  {


975  u = 2*hqrnduniformr(state)1;


976  v = 2*hqrnduniformr(state)1;


977  s = math.sqr(u)+math.sqr(v);


978  if( (double)(s)>(double)(0) && (double)(s)<(double)(1) )


979  {


980 


981  //


982  // two Sqrt's instead of one to


983  // avoid overflow when S is too small


984  //


985  s = Math.Sqrt((2*Math.Log(s)))/Math.Sqrt(s);


986  x1 = u*s;


987  x2 = v*s;


988  return;


989  }


990  }


991  }


992 


993 


994  /*************************************************************************


995  Random number generator: exponential distribution


996 


997  State structure must be initialized with HQRNDRandomize() or HQRNDSeed().


998 


999   ALGLIB 


1000  Copyright 11.08.2007 by Bochkanov Sergey


1001  *************************************************************************/


1002  public static double hqrndexponential(hqrndstate state,


1003  double lambdav)


1004  {


1005  double result = 0;


1006 


1007  alglib.ap.assert((double)(lambdav)>(double)(0), "HQRNDExponential: LambdaV<=0!");


1008  result = (Math.Log(hqrnduniformr(state))/lambdav);


1009  return result;


1010  }


1011 


1012 


1013  /*************************************************************************


1014  This function generates random number from discrete distribution given by


1015  finite sample X.


1016 


1017  INPUT PARAMETERS


1018  State  high quality random number generator, must be


1019  initialized with HQRNDRandomize() or HQRNDSeed().


1020  X  finite sample


1021  N  number of elements to use, N>=1


1022 


1023  RESULT


1024  this function returns one of the X[i] for random i=0..N1


1025 


1026   ALGLIB 


1027  Copyright 08.11.2011 by Bochkanov Sergey


1028  *************************************************************************/


1029  public static double hqrnddiscrete(hqrndstate state,


1030  double[] x,


1031  int n)


1032  {


1033  double result = 0;


1034 


1035  alglib.ap.assert(n>0, "HQRNDDiscrete: N<=0");


1036  alglib.ap.assert(n<=alglib.ap.len(x), "HQRNDDiscrete: Length(X)<N");


1037  result = x[hqrnduniformi(state, n)];


1038  return result;


1039  }


1040 


1041 


1042  /*************************************************************************


1043  This function generates random number from continuous distribution given


1044  by finite sample X.


1045 


1046  INPUT PARAMETERS


1047  State  high quality random number generator, must be


1048  initialized with HQRNDRandomize() or HQRNDSeed().


1049  X  finite sample, array[N] (can be larger, in this case only


1050  leading N elements are used). THIS ARRAY MUST BE SORTED BY


1051  ASCENDING.


1052  N  number of elements to use, N>=1


1053 


1054  RESULT


1055  this function returns random number from continuous distribution which


1056  tries to approximate X as mush as possible. min(X)<=Result<=max(X).


1057 


1058   ALGLIB 


1059  Copyright 08.11.2011 by Bochkanov Sergey


1060  *************************************************************************/


1061  public static double hqrndcontinuous(hqrndstate state,


1062  double[] x,


1063  int n)


1064  {


1065  double result = 0;


1066  double mx = 0;


1067  double mn = 0;


1068  int i = 0;


1069 


1070  alglib.ap.assert(n>0, "HQRNDContinuous: N<=0");


1071  alglib.ap.assert(n<=alglib.ap.len(x), "HQRNDContinuous: Length(X)<N");


1072  if( n==1 )


1073  {


1074  result = x[0];


1075  return result;


1076  }


1077  i = hqrnduniformi(state, n1);


1078  mn = x[i];


1079  mx = x[i+1];


1080  alglib.ap.assert((double)(mx)>=(double)(mn), "HQRNDDiscrete: X is not sorted by ascending");


1081  if( (double)(mx)!=(double)(mn) )


1082  {


1083  result = (mxmn)*hqrnduniformr(state)+mn;


1084  }


1085  else


1086  {


1087  result = mn;


1088  }


1089  return result;


1090  }


1091 


1092 


1093  /*************************************************************************


1094 


1095  L'Ecuyer, Efficient and portable combined random number generators


1096  *************************************************************************/


1097  private static int hqrndintegerbase(hqrndstate state)


1098  {


1099  int result = 0;


1100  int k = 0;


1101 


1102  alglib.ap.assert(state.magicv==hqrndmagic, "HQRNDIntegerBase: State is not correctly initialized!");


1103  k = state.s1/53668;


1104  state.s1 = 40014*(state.s1k*53668)k*12211;


1105  if( state.s1<0 )


1106  {


1107  state.s1 = state.s1+2147483563;


1108  }


1109  k = state.s2/52774;


1110  state.s2 = 40692*(state.s2k*52774)k*3791;


1111  if( state.s2<0 )


1112  {


1113  state.s2 = state.s2+2147483399;


1114  }


1115 


1116  //


1117  // Result


1118  //


1119  result = state.s1state.s2;


1120  if( result<1 )


1121  {


1122  result = result+2147483562;


1123  }


1124  return result;


1125  }


1126 


1127 


1128  }


1129  public class nearestneighbor


1130  {


1131  public class kdtree


1132  {


1133  public int n;


1134  public int nx;


1135  public int ny;


1136  public int normtype;


1137  public double[,] xy;


1138  public int[] tags;


1139  public double[] boxmin;


1140  public double[] boxmax;


1141  public int[] nodes;


1142  public double[] splits;


1143  public double[] x;


1144  public int kneeded;


1145  public double rneeded;


1146  public bool selfmatch;


1147  public double approxf;


1148  public int kcur;


1149  public int[] idx;


1150  public double[] r;


1151  public double[] buf;


1152  public double[] curboxmin;


1153  public double[] curboxmax;


1154  public double curdist;


1155  public int debugcounter;


1156  public kdtree()


1157  {


1158  xy = new double[0,0];


1159  tags = new int[0];


1160  boxmin = new double[0];


1161  boxmax = new double[0];


1162  nodes = new int[0];


1163  splits = new double[0];


1164  x = new double[0];


1165  idx = new int[0];


1166  r = new double[0];


1167  buf = new double[0];


1168  curboxmin = new double[0];


1169  curboxmax = new double[0];


1170  }


1171  };


1172 


1173 


1174 


1175 


1176  public const int splitnodesize = 6;


1177  public const int kdtreefirstversion = 0;


1178 


1179 


1180  /*************************************************************************


1181  KDtree creation


1182 


1183  This subroutine creates KDtree from set of Xvalues and optional Yvalues


1184 


1185  INPUT PARAMETERS


1186  XY  dataset, array[0..N1,0..NX+NY1].


1187  one row corresponds to one point.


1188  first NX columns contain Xvalues, next NY (NY may be zero)


1189  columns may contain associated Yvalues


1190  N  number of points, N>=0.


1191  NX  space dimension, NX>=1.


1192  NY  number of optional Yvalues, NY>=0.


1193  NormType norm type:


1194  * 0 denotes infinitynorm


1195  * 1 denotes 1norm


1196  * 2 denotes 2norm (Euclidean norm)


1197 


1198  OUTPUT PARAMETERS


1199  KDT  KDtree


1200 


1201 


1202  NOTES


1203 


1204  1. KDtree creation have O(N*logN) complexity and O(N*(2*NX+NY)) memory


1205  requirements.


1206  2. Although KDtrees may be used with any combination of N and NX, they


1207  are more efficient than bruteforce search only when N >> 4^NX. So they


1208  are most useful in lowdimensional tasks (NX=2, NX=3). NX=1 is another


1209  inefficient case, because simple binary search (without additional


1210  structures) is much more efficient in such tasks than KDtrees.


1211 


1212   ALGLIB 


1213  Copyright 28.02.2010 by Bochkanov Sergey


1214  *************************************************************************/


1215  public static void kdtreebuild(double[,] xy,


1216  int n,


1217  int nx,


1218  int ny,


1219  int normtype,


1220  kdtree kdt)


1221  {


1222  int[] tags = new int[0];


1223  int i = 0;


1224 


1225  alglib.ap.assert(n>=0, "KDTreeBuild: N<0");


1226  alglib.ap.assert(nx>=1, "KDTreeBuild: NX<1");


1227  alglib.ap.assert(ny>=0, "KDTreeBuild: NY<0");


1228  alglib.ap.assert(normtype>=0 && normtype<=2, "KDTreeBuild: incorrect NormType");


1229  alglib.ap.assert(alglib.ap.rows(xy)>=n, "KDTreeBuild: rows(X)<N");


1230  alglib.ap.assert(alglib.ap.cols(xy)>=nx+ny  n==0, "KDTreeBuild: cols(X)<NX+NY");


1231  alglib.ap.assert(apserv.apservisfinitematrix(xy, n, nx+ny), "KDTreeBuild: X contains infinite or NaN values");


1232  if( n>0 )


1233  {


1234  tags = new int[n];


1235  for(i=0; i<=n1; i++)


1236  {


1237  tags[i] = 0;


1238  }


1239  }


1240  kdtreebuildtagged(xy, tags, n, nx, ny, normtype, kdt);


1241  }


1242 


1243 


1244  /*************************************************************************


1245  KDtree creation


1246 


1247  This subroutine creates KDtree from set of Xvalues, integer tags and


1248  optional Yvalues


1249 


1250  INPUT PARAMETERS


1251  XY  dataset, array[0..N1,0..NX+NY1].


1252  one row corresponds to one point.


1253  first NX columns contain Xvalues, next NY (NY may be zero)


1254  columns may contain associated Yvalues


1255  Tags  tags, array[0..N1], contains integer tags associated


1256  with points.


1257  N  number of points, N>=0


1258  NX  space dimension, NX>=1.


1259  NY  number of optional Yvalues, NY>=0.


1260  NormType norm type:


1261  * 0 denotes infinitynorm


1262  * 1 denotes 1norm


1263  * 2 denotes 2norm (Euclidean norm)


1264 


1265  OUTPUT PARAMETERS


1266  KDT  KDtree


1267 


1268  NOTES


1269 


1270  1. KDtree creation have O(N*logN) complexity and O(N*(2*NX+NY)) memory


1271  requirements.


1272  2. Although KDtrees may be used with any combination of N and NX, they


1273  are more efficient than bruteforce search only when N >> 4^NX. So they


1274  are most useful in lowdimensional tasks (NX=2, NX=3). NX=1 is another


1275  inefficient case, because simple binary search (without additional


1276  structures) is much more efficient in such tasks than KDtrees.


1277 


1278   ALGLIB 


1279  Copyright 28.02.2010 by Bochkanov Sergey


1280  *************************************************************************/


1281  public static void kdtreebuildtagged(double[,] xy,


1282  int[] tags,


1283  int n,


1284  int nx,


1285  int ny,


1286  int normtype,


1287  kdtree kdt)


1288  {


1289  int i = 0;


1290  int j = 0;


1291  int maxnodes = 0;


1292  int nodesoffs = 0;


1293  int splitsoffs = 0;


1294  int i_ = 0;


1295  int i1_ = 0;


1296 


1297  alglib.ap.assert(n>=0, "KDTreeBuildTagged: N<0");


1298  alglib.ap.assert(nx>=1, "KDTreeBuildTagged: NX<1");


1299  alglib.ap.assert(ny>=0, "KDTreeBuildTagged: NY<0");


1300  alglib.ap.assert(normtype>=0 && normtype<=2, "KDTreeBuildTagged: incorrect NormType");


1301  alglib.ap.assert(alglib.ap.rows(xy)>=n, "KDTreeBuildTagged: rows(X)<N");


1302  alglib.ap.assert(alglib.ap.cols(xy)>=nx+ny  n==0, "KDTreeBuildTagged: cols(X)<NX+NY");


1303  alglib.ap.assert(apserv.apservisfinitematrix(xy, n, nx+ny), "KDTreeBuildTagged: X contains infinite or NaN values");


1304 


1305  //


1306  // initialize


1307  //


1308  kdt.n = n;


1309  kdt.nx = nx;


1310  kdt.ny = ny;


1311  kdt.normtype = normtype;


1312  kdt.kcur = 0;


1313 


1314  //


1315  // N=0 => quick exit


1316  //


1317  if( n==0 )


1318  {


1319  return;


1320  }


1321 


1322  //


1323  // Allocate


1324  //


1325  kdtreeallocdatasetindependent(kdt, nx, ny);


1326  kdtreeallocdatasetdependent(kdt, n, nx, ny);


1327 


1328  //


1329  // Initial fill


1330  //


1331  for(i=0; i<=n1; i++)


1332  {


1333  for(i_=0; i_<=nx1;i_++)


1334  {


1335  kdt.xy[i,i_] = xy[i,i_];


1336  }


1337  i1_ = (0)  (nx);


1338  for(i_=nx; i_<=2*nx+ny1;i_++)


1339  {


1340  kdt.xy[i,i_] = xy[i,i_+i1_];


1341  }


1342  kdt.tags[i] = tags[i];


1343  }


1344 


1345  //


1346  // Determine bounding box


1347  //


1348  for(i_=0; i_<=nx1;i_++)


1349  {


1350  kdt.boxmin[i_] = kdt.xy[0,i_];


1351  }


1352  for(i_=0; i_<=nx1;i_++)


1353  {


1354  kdt.boxmax[i_] = kdt.xy[0,i_];


1355  }


1356  for(i=1; i<=n1; i++)


1357  {


1358  for(j=0; j<=nx1; j++)


1359  {


1360  kdt.boxmin[j] = Math.Min(kdt.boxmin[j], kdt.xy[i,j]);


1361  kdt.boxmax[j] = Math.Max(kdt.boxmax[j], kdt.xy[i,j]);


1362  }


1363  }


1364 


1365  //


1366  // prepare tree structure


1367  // * MaxNodes=N because we guarantee no trivial splits, i.e.


1368  // every split will generate two nonempty boxes


1369  //


1370  maxnodes = n;


1371  kdt.nodes = new int[splitnodesize*2*maxnodes];


1372  kdt.splits = new double[2*maxnodes];


1373  nodesoffs = 0;


1374  splitsoffs = 0;


1375  for(i_=0; i_<=nx1;i_++)


1376  {


1377  kdt.curboxmin[i_] = kdt.boxmin[i_];


1378  }


1379  for(i_=0; i_<=nx1;i_++)


1380  {


1381  kdt.curboxmax[i_] = kdt.boxmax[i_];


1382  }


1383  kdtreegeneratetreerec(kdt, ref nodesoffs, ref splitsoffs, 0, n, 8);


1384  }


1385 


1386 


1387  /*************************************************************************


1388  KNN query: K nearest neighbors


1389 


1390  INPUT PARAMETERS


1391  KDT  KDtree


1392  X  point, array[0..NX1].


1393  K  number of neighbors to return, K>=1


1394  SelfMatch  whether selfmatches are allowed:


1395  * if True, nearest neighbor may be the point itself


1396  (if it exists in original dataset)


1397  * if False, then only points with nonzero distance


1398  are returned


1399  * if not given, considered True


1400 


1401  RESULT


1402  number of actual neighbors found (either K or N, if K>N).


1403 


1404  This subroutine performs query and stores its result in the internal


1405  structures of the KDtree. You can use following subroutines to obtain


1406  these results:


1407  * KDTreeQueryResultsX() to get Xvalues


1408  * KDTreeQueryResultsXY() to get X and Yvalues


1409  * KDTreeQueryResultsTags() to get tag values


1410  * KDTreeQueryResultsDistances() to get distances


1411 


1412   ALGLIB 


1413  Copyright 28.02.2010 by Bochkanov Sergey


1414  *************************************************************************/


1415  public static int kdtreequeryknn(kdtree kdt,


1416  double[] x,


1417  int k,


1418  bool selfmatch)


1419  {


1420  int result = 0;


1421 


1422  alglib.ap.assert(k>=1, "KDTreeQueryKNN: K<1!");


1423  alglib.ap.assert(alglib.ap.len(x)>=kdt.nx, "KDTreeQueryKNN: Length(X)<NX!");


1424  alglib.ap.assert(apserv.isfinitevector(x, kdt.nx), "KDTreeQueryKNN: X contains infinite or NaN values!");


1425  result = kdtreequeryaknn(kdt, x, k, selfmatch, 0.0);


1426  return result;


1427  }


1428 


1429 


1430  /*************************************************************************


1431  RNN query: all points within Rsphere centered at X


1432 


1433  INPUT PARAMETERS


1434  KDT  KDtree


1435  X  point, array[0..NX1].


1436  R  radius of sphere (in corresponding norm), R>0


1437  SelfMatch  whether selfmatches are allowed:


1438  * if True, nearest neighbor may be the point itself


1439  (if it exists in original dataset)


1440  * if False, then only points with nonzero distance


1441  are returned


1442  * if not given, considered True


1443 


1444  RESULT


1445  number of neighbors found, >=0


1446 


1447  This subroutine performs query and stores its result in the internal


1448  structures of the KDtree. You can use following subroutines to obtain


1449  actual results:


1450  * KDTreeQueryResultsX() to get Xvalues


1451  * KDTreeQueryResultsXY() to get X and Yvalues


1452  * KDTreeQueryResultsTags() to get tag values


1453  * KDTreeQueryResultsDistances() to get distances


1454 


1455   ALGLIB 


1456  Copyright 28.02.2010 by Bochkanov Sergey


1457  *************************************************************************/


1458  public static int kdtreequeryrnn(kdtree kdt,


1459  double[] x,


1460  double r,


1461  bool selfmatch)


1462  {


1463  int result = 0;


1464  int i = 0;


1465  int j = 0;


1466 


1467  alglib.ap.assert((double)(r)>(double)(0), "KDTreeQueryRNN: incorrect R!");


1468  alglib.ap.assert(alglib.ap.len(x)>=kdt.nx, "KDTreeQueryRNN: Length(X)<NX!");


1469  alglib.ap.assert(apserv.isfinitevector(x, kdt.nx), "KDTreeQueryRNN: X contains infinite or NaN values!");


1470 


1471  //


1472  // Handle special case: KDT.N=0


1473  //


1474  if( kdt.n==0 )


1475  {


1476  kdt.kcur = 0;


1477  result = 0;


1478  return result;


1479  }


1480 


1481  //


1482  // Prepare parameters


1483  //


1484  kdt.kneeded = 0;


1485  if( kdt.normtype!=2 )


1486  {


1487  kdt.rneeded = r;


1488  }


1489  else


1490  {


1491  kdt.rneeded = math.sqr(r);


1492  }


1493  kdt.selfmatch = selfmatch;


1494  kdt.approxf = 1;


1495  kdt.kcur = 0;


1496 


1497  //


1498  // calculate distance from point to current bounding box


1499  //


1500  kdtreeinitbox(kdt, x);


1501 


1502  //


1503  // call recursive search


1504  // results are returned as heap


1505  //


1506  kdtreequerynnrec(kdt, 0);


1507 


1508  //


1509  // pop from heap to generate ordered representation


1510  //


1511  // last element is not pop'ed because it is already in


1512  // its place


1513  //


1514  result = kdt.kcur;


1515  j = kdt.kcur;


1516  for(i=kdt.kcur; i>=2; i)


1517  {


1518  tsort.tagheappopi(ref kdt.r, ref kdt.idx, ref j);


1519  }


1520  return result;


1521  }


1522 


1523 


1524  /*************************************************************************


1525  KNN query: approximate K nearest neighbors


1526 


1527  INPUT PARAMETERS


1528  KDT  KDtree


1529  X  point, array[0..NX1].


1530  K  number of neighbors to return, K>=1


1531  SelfMatch  whether selfmatches are allowed:


1532  * if True, nearest neighbor may be the point itself


1533  (if it exists in original dataset)


1534  * if False, then only points with nonzero distance


1535  are returned


1536  * if not given, considered True


1537  Eps  approximation factor, Eps>=0. epsapproximate nearest


1538  neighbor is a neighbor whose distance from X is at


1539  most (1+eps) times distance of true nearest neighbor.


1540 


1541  RESULT


1542  number of actual neighbors found (either K or N, if K>N).


1543 


1544  NOTES


1545  significant performance gain may be achieved only when Eps is is on


1546  the order of magnitude of 1 or larger.


1547 


1548  This subroutine performs query and stores its result in the internal


1549  structures of the KDtree. You can use following subroutines to obtain


1550  these results:


1551  * KDTreeQueryResultsX() to get Xvalues


1552  * KDTreeQueryResultsXY() to get X and Yvalues


1553  * KDTreeQueryResultsTags() to get tag values


1554  * KDTreeQueryResultsDistances() to get distances


1555 


1556   ALGLIB 


1557  Copyright 28.02.2010 by Bochkanov Sergey


1558  *************************************************************************/


1559  public static int kdtreequeryaknn(kdtree kdt,


1560  double[] x,


1561  int k,


1562  bool selfmatch,


1563  double eps)


1564  {


1565  int result = 0;


1566  int i = 0;


1567  int j = 0;


1568 


1569  alglib.ap.assert(k>0, "KDTreeQueryAKNN: incorrect K!");


1570  alglib.ap.assert((double)(eps)>=(double)(0), "KDTreeQueryAKNN: incorrect Eps!");


1571  alglib.ap.assert(alglib.ap.len(x)>=kdt.nx, "KDTreeQueryAKNN: Length(X)<NX!");


1572  alglib.ap.assert(apserv.isfinitevector(x, kdt.nx), "KDTreeQueryAKNN: X contains infinite or NaN values!");


1573 


1574  //


1575  // Handle special case: KDT.N=0


1576  //


1577  if( kdt.n==0 )


1578  {


1579  kdt.kcur = 0;


1580  result = 0;


1581  return result;


1582  }


1583 


1584  //


1585  // Prepare parameters


1586  //


1587  k = Math.Min(k, kdt.n);


1588  kdt.kneeded = k;


1589  kdt.rneeded = 0;


1590  kdt.selfmatch = selfmatch;


1591  if( kdt.normtype==2 )


1592  {


1593  kdt.approxf = 1/math.sqr(1+eps);


1594  }


1595  else


1596  {


1597  kdt.approxf = 1/(1+eps);


1598  }


1599  kdt.kcur = 0;


1600 


1601  //


1602  // calculate distance from point to current bounding box


1603  //


1604  kdtreeinitbox(kdt, x);


1605 


1606  //


1607  // call recursive search


1608  // results are returned as heap


1609  //


1610  kdtreequerynnrec(kdt, 0);


1611 


1612  //


1613  // pop from heap to generate ordered representation


1614  //


1615  // last element is non pop'ed because it is already in


1616  // its place


1617  //


1618  result = kdt.kcur;


1619  j = kdt.kcur;


1620  for(i=kdt.kcur; i>=2; i)


1621  {


1622  tsort.tagheappopi(ref kdt.r, ref kdt.idx, ref j);


1623  }


1624  return result;


1625  }


1626 


1627 


1628  /*************************************************************************


1629  Xvalues from last query


1630 


1631  INPUT PARAMETERS


1632  KDT  KDtree


1633  X  possibly preallocated buffer. If X is too small to store


1634  result, it is resized. If size(X) is enough to store


1635  result, it is left unchanged.


1636 


1637  OUTPUT PARAMETERS


1638  X  rows are filled with Xvalues


1639 


1640  NOTES


1641  1. points are ordered by distance from the query point (first = closest)


1642  2. if XY is larger than required to store result, only leading part will


1643  be overwritten; trailing part will be left unchanged. So if on input


1644  XY = [[A,B],[C,D]], and result is [1,2], then on exit we will get


1645  XY = [[1,2],[C,D]]. This is done purposely to increase performance; if


1646  you want function to resize array according to result size, use


1647  function with same name and suffix 'I'.


1648 


1649  SEE ALSO


1650  * KDTreeQueryResultsXY() X and Yvalues


1651  * KDTreeQueryResultsTags() tag values


1652  * KDTreeQueryResultsDistances() distances


1653 


1654   ALGLIB 


1655  Copyright 28.02.2010 by Bochkanov Sergey


1656  *************************************************************************/


1657  public static void kdtreequeryresultsx(kdtree kdt,


1658  ref double[,] x)


1659  {


1660  int i = 0;


1661  int k = 0;


1662  int i_ = 0;


1663  int i1_ = 0;


1664 


1665  if( kdt.kcur==0 )


1666  {


1667  return;


1668  }


1669  if( alglib.ap.rows(x)<kdt.kcur  alglib.ap.cols(x)<kdt.nx )


1670  {


1671  x = new double[kdt.kcur, kdt.nx];


1672  }


1673  k = kdt.kcur;


1674  for(i=0; i<=k1; i++)


1675  {


1676  i1_ = (kdt.nx)  (0);


1677  for(i_=0; i_<=kdt.nx1;i_++)


1678  {


1679  x[i,i_] = kdt.xy[kdt.idx[i],i_+i1_];


1680  }


1681  }


1682  }


1683 


1684 


1685  /*************************************************************************


1686  X and Yvalues from last query


1687 


1688  INPUT PARAMETERS


1689  KDT  KDtree


1690  XY  possibly preallocated buffer. If XY is too small to store


1691  result, it is resized. If size(XY) is enough to store


1692  result, it is left unchanged.


1693 


1694  OUTPUT PARAMETERS


1695  XY  rows are filled with points: first NX columns with


1696  Xvalues, next NY columns  with Yvalues.


1697 


1698  NOTES


1699  1. points are ordered by distance from the query point (first = closest)


1700  2. if XY is larger than required to store result, only leading part will


1701  be overwritten; trailing part will be left unchanged. So if on input


1702  XY = [[A,B],[C,D]], and result is [1,2], then on exit we will get


1703  XY = [[1,2],[C,D]]. This is done purposely to increase performance; if


1704  you want function to resize array according to result size, use


1705  function with same name and suffix 'I'.


1706 


1707  SEE ALSO


1708  * KDTreeQueryResultsX() Xvalues


1709  * KDTreeQueryResultsTags() tag values


1710  * KDTreeQueryResultsDistances() distances


1711 


1712   ALGLIB 


1713  Copyright 28.02.2010 by Bochkanov Sergey


1714  *************************************************************************/


1715  public static void kdtreequeryresultsxy(kdtree kdt,


1716  ref double[,] xy)


1717  {


1718  int i = 0;


1719  int k = 0;


1720  int i_ = 0;


1721  int i1_ = 0;


1722 


1723  if( kdt.kcur==0 )


1724  {


1725  return;


1726  }


1727  if( alglib.ap.rows(xy)<kdt.kcur  alglib.ap.cols(xy)<kdt.nx+kdt.ny )


1728  {


1729  xy = new double[kdt.kcur, kdt.nx+kdt.ny];


1730  }


1731  k = kdt.kcur;


1732  for(i=0; i<=k1; i++)


1733  {


1734  i1_ = (kdt.nx)  (0);


1735  for(i_=0; i_<=kdt.nx+kdt.ny1;i_++)


1736  {


1737  xy[i,i_] = kdt.xy[kdt.idx[i],i_+i1_];


1738  }


1739  }


1740  }


1741 


1742 


1743  /*************************************************************************


1744  Tags from last query


1745 


1746  INPUT PARAMETERS


1747  KDT  KDtree


1748  Tags  possibly preallocated buffer. If X is too small to store


1749  result, it is resized. If size(X) is enough to store


1750  result, it is left unchanged.


1751 


1752  OUTPUT PARAMETERS


1753  Tags  filled with tags associated with points,


1754  or, when no tags were supplied, with zeros


1755 


1756  NOTES


1757  1. points are ordered by distance from the query point (first = closest)


1758  2. if XY is larger than required to store result, only leading part will


1759  be overwritten; trailing part will be left unchanged. So if on input


1760  XY = [[A,B],[C,D]], and result is [1,2], then on exit we will get


1761  XY = [[1,2],[C,D]]. This is done purposely to increase performance; if


1762  you want function to resize array according to result size, use


1763  function with same name and suffix 'I'.


1764 


1765  SEE ALSO


1766  * KDTreeQueryResultsX() Xvalues


1767  * KDTreeQueryResultsXY() X and Yvalues


1768  * KDTreeQueryResultsDistances() distances


1769 


1770   ALGLIB 


1771  Copyright 28.02.2010 by Bochkanov Sergey


1772  *************************************************************************/


1773  public static void kdtreequeryresultstags(kdtree kdt,


1774  ref int[] tags)


1775  {


1776  int i = 0;


1777  int k = 0;


1778 


1779  if( kdt.kcur==0 )


1780  {


1781  return;


1782  }


1783  if( alglib.ap.len(tags)<kdt.kcur )


1784  {


1785  tags = new int[kdt.kcur];


1786  }


1787  k = kdt.kcur;


1788  for(i=0; i<=k1; i++)


1789  {


1790  tags[i] = kdt.tags[kdt.idx[i]];


1791  }


1792  }


1793 


1794 


1795  /*************************************************************************


1796  Distances from last query


1797 


1798  INPUT PARAMETERS


1799  KDT  KDtree


1800  R  possibly preallocated buffer. If X is too small to store


1801  result, it is resized. If size(X) is enough to store


1802  result, it is left unchanged.


1803 


1804  OUTPUT PARAMETERS


1805  R  filled with distances (in corresponding norm)


1806 


1807  NOTES


1808  1. points are ordered by distance from the query point (first = closest)


1809  2. if XY is larger than required to store result, only leading part will


1810  be overwritten; trailing part will be left unchanged. So if on input


1811  XY = [[A,B],[C,D]], and result is [1,2], then on exit we will get


1812  XY = [[1,2],[C,D]]. This is done purposely to increase performance; if


1813  you want function to resize array according to result size, use


1814  function with same name and suffix 'I'.


1815 


1816  SEE ALSO


1817  * KDTreeQueryResultsX() Xvalues


1818  * KDTreeQueryResultsXY() X and Yvalues


1819  * KDTreeQueryResultsTags() tag values


1820 


1821   ALGLIB 


1822  Copyright 28.02.2010 by Bochkanov Sergey


1823  *************************************************************************/


1824  public static void kdtreequeryresultsdistances(kdtree kdt,


1825  ref double[] r)


1826  {


1827  int i = 0;


1828  int k = 0;


1829 


1830  if( kdt.kcur==0 )


1831  {


1832  return;


1833  }


1834  if( alglib.ap.len(r)<kdt.kcur )


1835  {


1836  r = new double[kdt.kcur];


1837  }


1838  k = kdt.kcur;


1839 


1840  //


1841  // unload norms


1842  //


1843  // Abs() call is used to handle cases with negative norms


1844  // (generated during KFN requests)


1845  //


1846  if( kdt.normtype==0 )


1847  {


1848  for(i=0; i<=k1; i++)


1849  {


1850  r[i] = Math.Abs(kdt.r[i]);


1851  }


1852  }


1853  if( kdt.normtype==1 )


1854  {


1855  for(i=0; i<=k1; i++)


1856  {


1857  r[i] = Math.Abs(kdt.r[i]);


1858  }


1859  }


1860  if( kdt.normtype==2 )


1861  {


1862  for(i=0; i<=k1; i++)


1863  {


1864  r[i] = Math.Sqrt(Math.Abs(kdt.r[i]));


1865  }


1866  }


1867  }


1868 


1869 


1870  /*************************************************************************


1871  Xvalues from last query; 'interactive' variant for languages like Python


1872  which support constructs like "X = KDTreeQueryResultsXI(KDT)" and


1873  interactive mode of interpreter.


1874 


1875  This function allocates new array on each call, so it is significantly


1876  slower than its 'noninteractive' counterpart, but it is more convenient


1877  when you call it from command line.


1878 


1879   ALGLIB 


1880  Copyright 28.02.2010 by Bochkanov Sergey


1881  *************************************************************************/


1882  public static void kdtreequeryresultsxi(kdtree kdt,


1883  ref double[,] x)


1884  {


1885  x = new double[0,0];


1886 


1887  kdtreequeryresultsx(kdt, ref x);


1888  }


1889 


1890 


1891  /*************************************************************************


1892  XYvalues from last query; 'interactive' variant for languages like Python


1893  which support constructs like "XY = KDTreeQueryResultsXYI(KDT)" and


1894  interactive mode of interpreter.


1895 


1896  This function allocates new array on each call, so it is significantly


1897  slower than its 'noninteractive' counterpart, but it is more convenient


1898  when you call it from command line.


1899 


1900   ALGLIB 


1901  Copyright 28.02.2010 by Bochkanov Sergey


1902  *************************************************************************/


1903  public static void kdtreequeryresultsxyi(kdtree kdt,


1904  ref double[,] xy)


1905  {


1906  xy = new double[0,0];


1907 


1908  kdtreequeryresultsxy(kdt, ref xy);


1909  }


1910 


1911 


1912  /*************************************************************************


1913  Tags from last query; 'interactive' variant for languages like Python


1914  which support constructs like "Tags = KDTreeQueryResultsTagsI(KDT)" and


1915  interactive mode of interpreter.


1916 


1917  This function allocates new array on each call, so it is significantly


1918  slower than its 'noninteractive' counterpart, but it is more convenient


1919  when you call it from command line.


1920 


1921   ALGLIB 


1922  Copyright 28.02.2010 by Bochkanov Sergey


1923  *************************************************************************/


1924  public static void kdtreequeryresultstagsi(kdtree kdt,


1925  ref int[] tags)


1926  {


1927  tags = new int[0];


1928 


1929  kdtreequeryresultstags(kdt, ref tags);


1930  }


1931 


1932 


1933  /*************************************************************************


1934  Distances from last query; 'interactive' variant for languages like Python


1935  which support constructs like "R = KDTreeQueryResultsDistancesI(KDT)"


1936  and interactive mode of interpreter.


1937 


1938  This function allocates new array on each call, so it is significantly


1939  slower than its 'noninteractive' counterpart, but it is more convenient


1940  when you call it from command line.


1941 


1942   ALGLIB 


1943  Copyright 28.02.2010 by Bochkanov Sergey


1944  *************************************************************************/


1945  public static void kdtreequeryresultsdistancesi(kdtree kdt,


1946  ref double[] r)


1947  {


1948  r = new double[0];


1949 


1950  kdtreequeryresultsdistances(kdt, ref r);


1951  }


1952 


1953 


1954  /*************************************************************************


1955  Serializer: allocation


1956 


1957   ALGLIB 


1958  Copyright 14.03.2011 by Bochkanov Sergey


1959  *************************************************************************/


1960  public static void kdtreealloc(alglib.serializer s,


1961  kdtree tree)


1962  {


1963 


1964  //


1965  // Header


1966  //


1967  s.alloc_entry();


1968  s.alloc_entry();


1969 


1970  //


1971  // Data


1972  //


1973  s.alloc_entry();


1974  s.alloc_entry();


1975  s.alloc_entry();


1976  s.alloc_entry();


1977  apserv.allocrealmatrix(s, tree.xy, 1, 1);


1978  apserv.allocintegerarray(s, tree.tags, 1);


1979  apserv.allocrealarray(s, tree.boxmin, 1);


1980  apserv.allocrealarray(s, tree.boxmax, 1);


1981  apserv.allocintegerarray(s, tree.nodes, 1);


1982  apserv.allocrealarray(s, tree.splits, 1);


1983  }


1984 


1985 


1986  /*************************************************************************


1987  Serializer: serialization


1988 


1989   ALGLIB 


1990  Copyright 14.03.2011 by Bochkanov Sergey


1991  *************************************************************************/


1992  public static void kdtreeserialize(alglib.serializer s,


1993  kdtree tree)


1994  {


1995 


1996  //


1997  // Header


1998  //


1999  s.serialize_int(scodes.getkdtreeserializationcode());


2000  s.serialize_int(kdtreefirstversion);


2001 


2002  //


2003  // Data


2004  //


2005  s.serialize_int(tree.n);


2006  s.serialize_int(tree.nx);


2007  s.serialize_int(tree.ny);


2008  s.serialize_int(tree.normtype);


2009  apserv.serializerealmatrix(s, tree.xy, 1, 1);


2010  apserv.serializeintegerarray(s, tree.tags, 1);


2011  apserv.serializerealarray(s, tree.boxmin, 1);


2012  apserv.serializerealarray(s, tree.boxmax, 1);


2013  apserv.serializeintegerarray(s, tree.nodes, 1);


2014  apserv.serializerealarray(s, tree.splits, 1);


2015  }


2016 


2017 


2018  /*************************************************************************


2019  Serializer: unserialization


2020 


2021   ALGLIB 


2022  Copyright 14.03.2011 by Bochkanov Sergey


2023  *************************************************************************/


2024  public static void kdtreeunserialize(alglib.serializer s,


2025  kdtree tree)


2026  {


2027  int i0 = 0;


2028  int i1 = 0;


2029 


2030 


2031  //


2032  // check correctness of header


2033  //


2034  i0 = s.unserialize_int();


2035  alglib.ap.assert(i0==scodes.getkdtreeserializationcode(), "KDTreeUnserialize: stream header corrupted");


2036  i1 = s.unserialize_int();


2037  alglib.ap.assert(i1==kdtreefirstversion, "KDTreeUnserialize: stream header corrupted");


2038 


2039  //


2040  // Unserialize data


2041  //


2042  tree.n = s.unserialize_int();


2043  tree.nx = s.unserialize_int();


2044  tree.ny = s.unserialize_int();


2045  tree.normtype = s.unserialize_int();


2046  apserv.unserializerealmatrix(s, ref tree.xy);


2047  apserv.unserializeintegerarray(s, ref tree.tags);


2048  apserv.unserializerealarray(s, ref tree.boxmin);


2049  apserv.unserializerealarray(s, ref tree.boxmax);


2050  apserv.unserializeintegerarray(s, ref tree.nodes);


2051  apserv.unserializerealarray(s, ref tree.splits);


2052  kdtreealloctemporaries(tree, tree.n, tree.nx, tree.ny);


2053  }


2054 


2055 


2056  /*************************************************************************


2057  Rearranges nodes [I1,I2) using partition in Dth dimension with S as threshold.


2058  Returns split position I3: [I1,I3) and [I3,I2) are created as result.


2059 


2060  This subroutine doesn't create tree structures, just rearranges nodes.


2061  *************************************************************************/


2062  private static void kdtreesplit(kdtree kdt,


2063  int i1,


2064  int i2,


2065  int d,


2066  double s,


2067  ref int i3)


2068  {


2069  int i = 0;


2070  int j = 0;


2071  int ileft = 0;


2072  int iright = 0;


2073  double v = 0;


2074 


2075  i3 = 0;


2076 


2077  alglib.ap.assert(kdt.n>0, "KDTreeSplit: internal error");


2078 


2079  //


2080  // split XY/Tags in two parts:


2081  // * [ILeft,IRight] is nonprocessed part of XY/Tags


2082  //


2083  // After cycle is done, we have Ileft=IRight. We deal with


2084  // this element separately.


2085  //


2086  // After this, [I1,ILeft) contains left part, and [ILeft,I2)


2087  // contains right part.


2088  //


2089  ileft = i1;


2090  iright = i21;


2091  while( ileft<iright )


2092  {


2093  if( (double)(kdt.xy[ileft,d])<=(double)(s) )


2094  {


2095 


2096  //


2097  // XY[ILeft] is on its place.


2098  // Advance ILeft.


2099  //


2100  ileft = ileft+1;


2101  }


2102  else


2103  {


2104 


2105  //


2106  // XY[ILeft,..] must be at IRight.


2107  // Swap and advance IRight.


2108  //


2109  for(i=0; i<=2*kdt.nx+kdt.ny1; i++)


2110  {


2111  v = kdt.xy[ileft,i];


2112  kdt.xy[ileft,i] = kdt.xy[iright,i];


2113  kdt.xy[iright,i] = v;


2114  }


2115  j = kdt.tags[ileft];


2116  kdt.tags[ileft] = kdt.tags[iright];


2117  kdt.tags[iright] = j;


2118  iright = iright1;


2119  }


2120  }


2121  if( (double)(kdt.xy[ileft,d])<=(double)(s) )


2122  {


2123  ileft = ileft+1;


2124  }


2125  else


2126  {


2127  iright = iright1;


2128  }


2129  i3 = ileft;


2130  }


2131 


2132 


2133  /*************************************************************************


2134  Recursive kdtree generation subroutine.


2135 


2136  PARAMETERS


2137  KDT tree


2138  NodesOffs unused part of Nodes[] which must be filled by tree


2139  SplitsOffs unused part of Splits[]


2140  I1, I2 points from [I1,I2) are processed


2141 


2142  NodesOffs[] and SplitsOffs[] must be large enough.


2143 


2144   ALGLIB 


2145  Copyright 28.02.2010 by Bochkanov Sergey


2146  *************************************************************************/


2147  private static void kdtreegeneratetreerec(kdtree kdt,


2148  ref int nodesoffs,


2149  ref int splitsoffs,


2150  int i1,


2151  int i2,


2152  int maxleafsize)


2153  {


2154  int n = 0;


2155  int nx = 0;


2156  int ny = 0;


2157  int i = 0;


2158  int j = 0;


2159  int oldoffs = 0;


2160  int i3 = 0;


2161  int cntless = 0;


2162  int cntgreater = 0;


2163  double minv = 0;


2164  double maxv = 0;


2165  int minidx = 0;


2166  int maxidx = 0;


2167  int d = 0;


2168  double ds = 0;


2169  double s = 0;


2170  double v = 0;


2171  int i_ = 0;


2172  int i1_ = 0;


2173 


2174  alglib.ap.assert(kdt.n>0, "KDTreeGenerateTreeRec: internal error");


2175  alglib.ap.assert(i2>i1, "KDTreeGenerateTreeRec: internal error");


2176 


2177  //


2178  // Generate leaf if needed


2179  //


2180  if( i2i1<=maxleafsize )


2181  {


2182  kdt.nodes[nodesoffs+0] = i2i1;


2183  kdt.nodes[nodesoffs+1] = i1;


2184  nodesoffs = nodesoffs+2;


2185  return;


2186  }


2187 


2188  //


2189  // Load values for easier access


2190  //


2191  nx = kdt.nx;


2192  ny = kdt.ny;


2193 


2194  //


2195  // select dimension to split:


2196  // * D is a dimension number


2197  //


2198  d = 0;


2199  ds = kdt.curboxmax[0]kdt.curboxmin[0];


2200  for(i=1; i<=nx1; i++)


2201  {


2202  v = kdt.curboxmax[i]kdt.curboxmin[i];


2203  if( (double)(v)>(double)(ds) )


2204  {


2205  ds = v;


2206  d = i;


2207  }


2208  }


2209 


2210  //


2211  // Select split position S using sliding midpoint rule,


2212  // rearrange points into [I1,I3) and [I3,I2)


2213  //


2214  s = kdt.curboxmin[d]+0.5*ds;


2215  i1_ = (i1)  (0);


2216  for(i_=0; i_<=i2i11;i_++)


2217  {


2218  kdt.buf[i_] = kdt.xy[i_+i1_,d];


2219  }


2220  n = i2i1;


2221  cntless = 0;


2222  cntgreater = 0;


2223  minv = kdt.buf[0];


2224  maxv = kdt.buf[0];


2225  minidx = i1;


2226  maxidx = i1;


2227  for(i=0; i<=n1; i++)


2228  {


2229  v = kdt.buf[i];


2230  if( (double)(v)<(double)(minv) )


2231  {


2232  minv = v;


2233  minidx = i1+i;


2234  }


2235  if( (double)(v)>(double)(maxv) )


2236  {


2237  maxv = v;


2238  maxidx = i1+i;


2239  }


2240  if( (double)(v)<(double)(s) )


2241  {


2242  cntless = cntless+1;


2243  }


2244  if( (double)(v)>(double)(s) )


2245  {


2246  cntgreater = cntgreater+1;


2247  }


2248  }


2249  if( cntless>0 && cntgreater>0 )


2250  {


2251 


2252  //


2253  // normal midpoint split


2254  //


2255  kdtreesplit(kdt, i1, i2, d, s, ref i3);


2256  }


2257  else


2258  {


2259 


2260  //


2261  // sliding midpoint


2262  //


2263  if( cntless==0 )


2264  {


2265 


2266  //


2267  // 1. move split to MinV,


2268  // 2. place one point to the left bin (move to I1),


2269  // others  to the right bin


2270  //


2271  s = minv;


2272  if( minidx!=i1 )


2273  {


2274  for(i=0; i<=2*kdt.nx+kdt.ny1; i++)


2275  {


2276  v = kdt.xy[minidx,i];


2277  kdt.xy[minidx,i] = kdt.xy[i1,i];


2278  kdt.xy[i1,i] = v;


2279  }


2280  j = kdt.tags[minidx];


2281  kdt.tags[minidx] = kdt.tags[i1];


2282  kdt.tags[i1] = j;


2283  }


2284  i3 = i1+1;


2285  }


2286  else


2287  {


2288 


2289  //


2290  // 1. move split to MaxV,


2291  // 2. place one point to the right bin (move to I21),


2292  // others  to the left bin


2293  //


2294  s = maxv;


2295  if( maxidx!=i21 )


2296  {


2297  for(i=0; i<=2*kdt.nx+kdt.ny1; i++)


2298  {


2299  v = kdt.xy[maxidx,i];


2300  kdt.xy[maxidx,i] = kdt.xy[i21,i];


2301  kdt.xy[i21,i] = v;


2302  }


2303  j = kdt.tags[maxidx];


2304  kdt.tags[maxidx] = kdt.tags[i21];


2305  kdt.tags[i21] = j;


2306  }


2307  i3 = i21;


2308  }


2309  }


2310 


2311  //


2312  // Generate 'split' node


2313  //


2314  kdt.nodes[nodesoffs+0] = 0;


2315  kdt.nodes[nodesoffs+1] = d;


2316  kdt.nodes[nodesoffs+2] = splitsoffs;


2317  kdt.splits[splitsoffs+0] = s;


2318  oldoffs = nodesoffs;


2319  nodesoffs = nodesoffs+splitnodesize;


2320  splitsoffs = splitsoffs+1;


2321 


2322  //


2323  // Recirsive generation:


2324  // * update CurBox


2325  // * call subroutine


2326  // * restore CurBox


2327  //


2328  kdt.nodes[oldoffs+3] = nodesoffs;


2329  v = kdt.curboxmax[d];


2330  kdt.curboxmax[d] = s;


2331  kdtreegeneratetreerec(kdt, ref nodesoffs, ref splitsoffs, i1, i3, maxleafsize);


2332  kdt.curboxmax[d] = v;


2333  kdt.nodes[oldoffs+4] = nodesoffs;


2334  v = kdt.curboxmin[d];


2335  kdt.curboxmin[d] = s;


2336  kdtreegeneratetreerec(kdt, ref nodesoffs, ref splitsoffs, i3, i2, maxleafsize);


2337  kdt.curboxmin[d] = v;


2338  }


2339 


2340 


2341  /*************************************************************************


2342  Recursive subroutine for NN queries.


2343 


2344   ALGLIB 


2345  Copyright 28.02.2010 by Bochkanov Sergey


2346  *************************************************************************/


2347  private static void kdtreequerynnrec(kdtree kdt,


2348  int offs)


2349  {


2350  double ptdist = 0;


2351  int i = 0;


2352  int j = 0;


2353  int nx = 0;


2354  int i1 = 0;


2355  int i2 = 0;


2356  int d = 0;


2357  double s = 0;


2358  double v = 0;


2359  double t1 = 0;


2360  int childbestoffs = 0;


2361  int childworstoffs = 0;


2362  int childoffs = 0;


2363  double prevdist = 0;


2364  bool todive = new bool();


2365  bool bestisleft = new bool();


2366  bool updatemin = new bool();


2367 


2368  alglib.ap.assert(kdt.n>0, "KDTreeQueryNNRec: internal error");


2369 


2370  //


2371  // Leaf node.


2372  // Process points.


2373  //


2374  if( kdt.nodes[offs]>0 )


2375  {


2376  i1 = kdt.nodes[offs+1];


2377  i2 = i1+kdt.nodes[offs];


2378  for(i=i1; i<=i21; i++)


2379  {


2380 


2381  //


2382  // Calculate distance


2383  //


2384  ptdist = 0;


2385  nx = kdt.nx;


2386  if( kdt.normtype==0 )


2387  {


2388  for(j=0; j<=nx1; j++)


2389  {


2390  ptdist = Math.Max(ptdist, Math.Abs(kdt.xy[i,j]kdt.x[j]));


2391  }


2392  }


2393  if( kdt.normtype==1 )


2394  {


2395  for(j=0; j<=nx1; j++)


2396  {


2397  ptdist = ptdist+Math.Abs(kdt.xy[i,j]kdt.x[j]);


2398  }


2399  }


2400  if( kdt.normtype==2 )


2401  {


2402  for(j=0; j<=nx1; j++)


2403  {


2404  ptdist = ptdist+math.sqr(kdt.xy[i,j]kdt.x[j]);


2405  }


2406  }


2407 


2408  //


2409  // Skip points with zero distance if selfmatches are turned off


2410  //


2411  if( (double)(ptdist)==(double)(0) && !kdt.selfmatch )


2412  {


2413  continue;


2414  }


2415 


2416  //


2417  // We CAN'T process point if Rcriterion isn't satisfied,


2418  // i.e. (RNeeded<>0) AND (PtDist>R).


2419  //


2420  if( (double)(kdt.rneeded)==(double)(0)  (double)(ptdist)<=(double)(kdt.rneeded) )


2421  {


2422 


2423  //


2424  // Rcriterion is satisfied, we must either:


2425  // * replace worst point, if (KNeeded<>0) AND (KCur=KNeeded)


2426  // (or skip, if worst point is better)


2427  // * add point without replacement otherwise


2428  //


2429  if( kdt.kcur<kdt.kneeded  kdt.kneeded==0 )


2430  {


2431 


2432  //


2433  // add current point to heap without replacement


2434  //


2435  tsort.tagheappushi(ref kdt.r, ref kdt.idx, ref kdt.kcur, ptdist, i);


2436  }


2437  else


2438  {


2439 


2440  //


2441  // New points are added or not, depending on their distance.


2442  // If added, they replace element at the top of the heap


2443  //


2444  if( (double)(ptdist)<(double)(kdt.r[0]) )


2445  {


2446  if( kdt.kneeded==1 )


2447  {


2448  kdt.idx[0] = i;


2449  kdt.r[0] = ptdist;


2450  }


2451  else


2452  {


2453  tsort.tagheapreplacetopi(ref kdt.r, ref kdt.idx, kdt.kneeded, ptdist, i);


2454  }


2455  }


2456  }


2457  }


2458  }


2459  return;


2460  }


2461 


2462  //


2463  // Simple split


2464  //


2465  if( kdt.nodes[offs]==0 )


2466  {


2467 


2468  //


2469  // Load:


2470  // * D dimension to split


2471  // * S split position


2472  //


2473  d = kdt.nodes[offs+1];


2474  s = kdt.splits[kdt.nodes[offs+2]];


2475 


2476  //


2477  // Calculate:


2478  // * ChildBestOffs child box with best chances


2479  // * ChildWorstOffs child box with worst chances


2480  //


2481  if( (double)(kdt.x[d])<=(double)(s) )


2482  {


2483  childbestoffs = kdt.nodes[offs+3];


2484  childworstoffs = kdt.nodes[offs+4];


2485  bestisleft = true;


2486  }


2487  else


2488  {


2489  childbestoffs = kdt.nodes[offs+4];


2490  childworstoffs = kdt.nodes[offs+3];


2491  bestisleft = false;


2492  }


2493 


2494  //


2495  // Navigate through childs


2496  //


2497  for(i=0; i<=1; i++)


2498  {


2499 


2500  //


2501  // Select child to process:


2502  // * ChildOffs current child offset in Nodes[]


2503  // * UpdateMin whether minimum or maximum value


2504  // of bounding box is changed on update


2505  //


2506  if( i==0 )


2507  {


2508  childoffs = childbestoffs;


2509  updatemin = !bestisleft;


2510  }


2511  else


2512  {


2513  updatemin = bestisleft;


2514  childoffs = childworstoffs;


2515  }


2516 


2517  //


2518  // Update bounding box and current distance


2519  //


2520  if( updatemin )


2521  {


2522  prevdist = kdt.curdist;


2523  t1 = kdt.x[d];


2524  v = kdt.curboxmin[d];


2525  if( (double)(t1)<=(double)(s) )


2526  {


2527  if( kdt.normtype==0 )


2528  {


2529  kdt.curdist = Math.Max(kdt.curdist, st1);


2530  }


2531  if( kdt.normtype==1 )


2532  {


2533  kdt.curdist = kdt.curdistMath.Max(vt1, 0)+st1;


2534  }


2535  if( kdt.normtype==2 )


2536  {


2537  kdt.curdist = kdt.curdistmath.sqr(Math.Max(vt1, 0))+math.sqr(st1);


2538  }


2539  }


2540  kdt.curboxmin[d] = s;


2541  }


2542  else


2543  {


2544  prevdist = kdt.curdist;


2545  t1 = kdt.x[d];


2546  v = kdt.curboxmax[d];


2547  if( (double)(t1)>=(double)(s) )


2548  {


2549  if( kdt.normtype==0 )


2550  {


2551  kdt.curdist = Math.Max(kdt.curdist, t1s);


2552  }


2553  if( kdt.normtype==1 )


2554  {


2555  kdt.curdist = kdt.curdistMath.Max(t1v, 0)+t1s;


2556  }


2557  if( kdt.normtype==2 )


2558  {


2559  kdt.curdist = kdt.curdistmath.sqr(Math.Max(t1v, 0))+math.sqr(t1s);


2560  }


2561  }


2562  kdt.curboxmax[d] = s;


2563  }


2564 


2565  //


2566  // Decide: to dive into cell or not to dive


2567  //


2568  if( (double)(kdt.rneeded)!=(double)(0) && (double)(kdt.curdist)>(double)(kdt.rneeded) )


2569  {


2570  todive = false;


2571  }


2572  else


2573  {


2574  if( kdt.kcur<kdt.kneeded  kdt.kneeded==0 )


2575  {


2576 


2577  //


2578  // KCur<KNeeded (i.e. not all points are found)


2579  //


2580  todive = true;


2581  }


2582  else


2583  {


2584 


2585  //


2586  // KCur=KNeeded, decide to dive or not to dive


2587  // using point position relative to bounding box.


2588  //


2589  todive = (double)(kdt.curdist)<=(double)(kdt.r[0]*kdt.approxf);


2590  }


2591  }


2592  if( todive )


2593  {


2594  kdtreequerynnrec(kdt, childoffs);


2595  }


2596 


2597  //


2598  // Restore bounding box and distance


2599  //


2600  if( updatemin )


2601  {


2602  kdt.curboxmin[d] = v;


2603  }


2604  else


2605  {


2606  kdt.curboxmax[d] = v;


2607  }


2608  kdt.curdist = prevdist;


2609  }


2610  return;


2611  }


2612  }


2613 


2614 


2615  /*************************************************************************


2616  Copies X[] to KDT.X[]


2617  Loads distance from X[] to bounding box.


2618  Initializes CurBox[].


2619 


2620   ALGLIB 


2621  Copyright 28.02.2010 by Bochkanov Sergey


2622  *************************************************************************/


2623  private static void kdtreeinitbox(kdtree kdt,


2624  double[] x)


2625  {


2626  int i = 0;


2627  double vx = 0;


2628  double vmin = 0;


2629  double vmax = 0;


2630 


2631  alglib.ap.assert(kdt.n>0, "KDTreeInitBox: internal error");


2632 


2633  //


2634  // calculate distance from point to current bounding box


2635  //


2636  kdt.curdist = 0;


2637  if( kdt.normtype==0 )


2638  {


2639  for(i=0; i<=kdt.nx1; i++)


2640  {


2641  vx = x[i];


2642  vmin = kdt.boxmin[i];


2643  vmax = kdt.boxmax[i];


2644  kdt.x[i] = vx;


2645  kdt.curboxmin[i] = vmin;


2646  kdt.curboxmax[i] = vmax;


2647  if( (double)(vx)<(double)(vmin) )


2648  {


2649  kdt.curdist = Math.Max(kdt.curdist, vminvx);


2650  }


2651  else


2652  {


2653  if( (double)(vx)>(double)(vmax) )


2654  {


2655  kdt.curdist = Math.Max(kdt.curdist, vxvmax);


2656  }


2657  }


2658  }


2659  }


2660  if( kdt.normtype==1 )


2661  {


2662  for(i=0; i<=kdt.nx1; i++)


2663  {


2664  vx = x[i];


2665  vmin = kdt.boxmin[i];


2666  vmax = kdt.boxmax[i];


2667  kdt.x[i] = vx;


2668  kdt.curboxmin[i] = vmin;


2669  kdt.curboxmax[i] = vmax;


2670  if( (double)(vx)<(double)(vmin) )


2671  {


2672  kdt.curdist = kdt.curdist+vminvx;


2673  }


2674  else


2675  {


2676  if( (double)(vx)>(double)(vmax) )


2677  {


2678  kdt.curdist = kdt.curdist+vxvmax;


2679  }


2680  }


2681  }


2682  }


2683  if( kdt.normtype==2 )


2684  {


2685  for(i=0; i<=kdt.nx1; i++)


2686  {


2687  vx = x[i];


2688  vmin = kdt.boxmin[i];


2689  vmax = kdt.boxmax[i];


2690  kdt.x[i] = vx;


2691  kdt.curboxmin[i] = vmin;


2692  kdt.curboxmax[i] = vmax;


2693  if( (double)(vx)<(double)(vmin) )


2694  {


2695  kdt.curdist = kdt.curdist+math.sqr(vminvx);


2696  }


2697  else


2698  {


2699  if( (double)(vx)>(double)(vmax) )


2700  {


2701  kdt.curdist = kdt.curdist+math.sqr(vxvmax);


2702  }


2703  }


2704  }


2705  }


2706  }


2707 


2708 


2709  /*************************************************************************


2710  This function allocates all datasetindependent array fields of KDTree,


2711  i.e. such array fields that their dimensions do not depend on dataset


2712  size.


2713 


2714  This function do not sets KDT.NX or KDT.NY  it just allocates arrays


2715 


2716   ALGLIB 


2717  Copyright 14.03.2011 by Bochkanov Sergey


2718  *************************************************************************/


2719  private static void kdtreeallocdatasetindependent(kdtree kdt,


2720  int nx,


2721  int ny)


2722  {


2723  alglib.ap.assert(kdt.n>0, "KDTreeAllocDatasetIndependent: internal error");


2724  kdt.x = new double[nx];


2725  kdt.boxmin = new double[nx];


2726  kdt.boxmax = new double[nx];


2727  kdt.curboxmin = new double[nx];


2728  kdt.curboxmax = new double[nx];


2729  }


2730 


2731 


2732  /*************************************************************************


2733  This function allocates all datasetdependent array fields of KDTree, i.e.


2734  such array fields that their dimensions depend on dataset size.


2735 


2736  This function do not sets KDT.N, KDT.NX or KDT.NY 


2737  it just allocates arrays.


2738 


2739   ALGLIB 


2740  Copyright 14.03.2011 by Bochkanov Sergey


2741  *************************************************************************/


2742  private static void kdtreeallocdatasetdependent(kdtree kdt,


2743  int n,


2744  int nx,


2745  int ny)


2746  {


2747  alglib.ap.assert(n>0, "KDTreeAllocDatasetDependent: internal error");


2748  kdt.xy = new double[n, 2*nx+ny];


2749  kdt.tags = new int[n];


2750  kdt.idx = new int[n];


2751  kdt.r = new double[n];


2752  kdt.x = new double[nx];


2753  kdt.buf = new double[Math.Max(n, nx)];


2754  kdt.nodes = new int[splitnodesize*2*n];


2755  kdt.splits = new double[2*n];


2756  }


2757 


2758 


2759  /*************************************************************************


2760  This function allocates temporaries.


2761 


2762  This function do not sets KDT.N, KDT.NX or KDT.NY 


2763  it just allocates arrays.


2764 


2765   ALGLIB 


2766  Copyright 14.03.2011 by Bochkanov Sergey


2767  *************************************************************************/


2768  private static void kdtreealloctemporaries(kdtree kdt,


2769  int n,


2770  int nx,


2771  int ny)


2772  {


2773  alglib.ap.assert(n>0, "KDTreeAllocTemporaries: internal error");


2774  kdt.x = new double[nx];


2775  kdt.idx = new int[n];


2776  kdt.r = new double[n];


2777  kdt.buf = new double[Math.Max(n, nx)];


2778  kdt.curboxmin = new double[nx];


2779  kdt.curboxmax = new double[nx];


2780  }


2781 


2782 


2783  }


2784  }


2785 

