Free cookie consent management tool by TermsFeed Policy Generator

source: branches/OaaS/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/3.4.0/ALGLIB-3.4.0/alglibmisc.cs @ 8614

Last change on this file since 8614 was 7294, checked in by gkronber, 13 years ago

#1733 updated alglib sources to version 3.4.0 of alglib and updated project and plugin references

File size: 95.5 KB
Line 
1/*************************************************************************
2Copyright (c) Sergey Bochkanov (ALGLIB project).
3
4>>> SOURCE LICENSE >>>
5This program is free software; you can redistribute it and/or modify
6it under the terms of the GNU General Public License as published by
7the Free Software Foundation (www.fsf.org); either version 2 of the
8License, or (at your option) any later version.
9
10This program is distributed in the hope that it will be useful,
11but WITHOUT ANY WARRANTY; without even the implied warranty of
12MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13GNU General Public License for more details.
14
15A copy of the GNU General Public License is available at
16http://www.fsf.org/licensing/licenses
17>>> END OF LICENSE >>>
18*************************************************************************/
19#pragma warning disable 162
20#pragma warning disable 219
21using System;
22
23public 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 HQRNDMax-1.
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}
188public partial class alglib
189{
190
191
192    /*************************************************************************
193
194    *************************************************************************/
195    public class kdtree
196    {
197        //
198        // Public declarations
199        //
200
201        public kdtree()
202        {
203            _innerobj = new nearestneighbor.kdtree();
204        }
205
206        //
207        // Although some of declarations below are public, you should not use them
208        // They are intended for internal use only
209        //
210        private nearestneighbor.kdtree _innerobj;
211        public nearestneighbor.kdtree innerobj { get { return _innerobj; } }
212        public kdtree(nearestneighbor.kdtree obj)
213        {
214            _innerobj = obj;
215        }
216    }
217
218
219    /*************************************************************************
220    This function serializes data structure to string.
221
222    Important properties of s_out:
223    * it contains alphanumeric characters, dots, underscores, minus signs
224    * these symbols are grouped into words, which are separated by spaces
225      and Windows-style (CR+LF) newlines
226    * although  serializer  uses  spaces and CR+LF as separators, you can
227      replace any separator character by arbitrary combination of spaces,
228      tabs, Windows or Unix newlines. It allows flexible reformatting  of
229      the  string  in  case you want to include it into text or XML file.
230      But you should not insert separators into the middle of the "words"
231      nor you should change case of letters.
232    * s_out can be freely moved between 32-bit and 64-bit systems, little
233      and big endian machines, and so on. You can serialize structure  on
234      32-bit machine and unserialize it on 64-bit one (or vice versa), or
235      serialize  it  on  SPARC  and  unserialize  on  x86.  You  can also
236      serialize  it  in  C# version of ALGLIB and unserialize in C++ one,
237      and vice versa.
238    *************************************************************************/
239    public static void kdtreeserialize(kdtree obj, out string s_out)
240    {
241        alglib.serializer s = new alglib.serializer();
242        s.alloc_start();
243        nearestneighbor.kdtreealloc(s, obj.innerobj);
244        s.sstart_str();
245        nearestneighbor.kdtreeserialize(s, obj.innerobj);
246        s.stop();
247        s_out = s.get_string();
248    }
249
250
251    /*************************************************************************
252    This function unserializes data structure from string.
253    *************************************************************************/
254    public static void kdtreeunserialize(string s_in, out kdtree obj)
255    {
256        alglib.serializer s = new alglib.serializer();
257        obj = new kdtree();
258        s.ustart_str(s_in);
259        nearestneighbor.kdtreeunserialize(s, obj.innerobj);
260        s.stop();
261    }
262
263    /*************************************************************************
264    KD-tree creation
265
266    This subroutine creates KD-tree from set of X-values and optional Y-values
267
268    INPUT PARAMETERS
269        XY      -   dataset, array[0..N-1,0..NX+NY-1].
270                    one row corresponds to one point.
271                    first NX columns contain X-values, next NY (NY may be zero)
272                    columns may contain associated Y-values
273        N       -   number of points, N>=1
274        NX      -   space dimension, NX>=1.
275        NY      -   number of optional Y-values, NY>=0.
276        NormType-   norm type:
277                    * 0 denotes infinity-norm
278                    * 1 denotes 1-norm
279                    * 2 denotes 2-norm (Euclidean norm)
280
281    OUTPUT PARAMETERS
282        KDT     -   KD-tree
283
284
285    NOTES
286
287    1. KD-tree  creation  have O(N*logN) complexity and O(N*(2*NX+NY))  memory
288       requirements.
289    2. Although KD-trees may be used with any combination of N  and  NX,  they
290       are more efficient than brute-force search only when N >> 4^NX. So they
291       are most useful in low-dimensional tasks (NX=2, NX=3). NX=1  is another
292       inefficient case, because  simple  binary  search  (without  additional
293       structures) is much more efficient in such tasks than KD-trees.
294
295      -- ALGLIB --
296         Copyright 28.02.2010 by Bochkanov Sergey
297    *************************************************************************/
298    public static void kdtreebuild(double[,] xy, int n, int nx, int ny, int normtype, out kdtree kdt)
299    {
300        kdt = new kdtree();
301        nearestneighbor.kdtreebuild(xy, n, nx, ny, normtype, kdt.innerobj);
302        return;
303    }
304    public static void kdtreebuild(double[,] xy, int nx, int ny, int normtype, out kdtree kdt)
305    {
306        int n;
307
308        kdt = new kdtree();
309        n = ap.rows(xy);
310        nearestneighbor.kdtreebuild(xy, n, nx, ny, normtype, kdt.innerobj);
311
312        return;
313    }
314
315    /*************************************************************************
316    KD-tree creation
317
318    This  subroutine  creates  KD-tree  from set of X-values, integer tags and
319    optional Y-values
320
321    INPUT PARAMETERS
322        XY      -   dataset, array[0..N-1,0..NX+NY-1].
323                    one row corresponds to one point.
324                    first NX columns contain X-values, next NY (NY may be zero)
325                    columns may contain associated Y-values
326        Tags    -   tags, array[0..N-1], contains integer tags associated
327                    with points.
328        N       -   number of points, N>=1
329        NX      -   space dimension, NX>=1.
330        NY      -   number of optional Y-values, NY>=0.
331        NormType-   norm type:
332                    * 0 denotes infinity-norm
333                    * 1 denotes 1-norm
334                    * 2 denotes 2-norm (Euclidean norm)
335
336    OUTPUT PARAMETERS
337        KDT     -   KD-tree
338
339    NOTES
340
341    1. KD-tree  creation  have O(N*logN) complexity and O(N*(2*NX+NY))  memory
342       requirements.
343    2. Although KD-trees may be used with any combination of N  and  NX,  they
344       are more efficient than brute-force search only when N >> 4^NX. So they
345       are most useful in low-dimensional tasks (NX=2, NX=3). NX=1  is another
346       inefficient case, because  simple  binary  search  (without  additional
347       structures) is much more efficient in such tasks than KD-trees.
348
349      -- ALGLIB --
350         Copyright 28.02.2010 by Bochkanov Sergey
351    *************************************************************************/
352    public static void kdtreebuildtagged(double[,] xy, int[] tags, int n, int nx, int ny, int normtype, out kdtree kdt)
353    {
354        kdt = new kdtree();
355        nearestneighbor.kdtreebuildtagged(xy, tags, n, nx, ny, normtype, kdt.innerobj);
356        return;
357    }
358    public static void kdtreebuildtagged(double[,] xy, int[] tags, int nx, int ny, int normtype, out kdtree kdt)
359    {
360        int n;
361        if( (ap.rows(xy)!=ap.len(tags)))
362            throw new alglibexception("Error while calling 'kdtreebuildtagged': looks like one of arguments has wrong size");
363        kdt = new kdtree();
364        n = ap.rows(xy);
365        nearestneighbor.kdtreebuildtagged(xy, tags, n, nx, ny, normtype, kdt.innerobj);
366
367        return;
368    }
369
370    /*************************************************************************
371    K-NN query: K nearest neighbors
372
373    INPUT PARAMETERS
374        KDT         -   KD-tree
375        X           -   point, array[0..NX-1].
376        K           -   number of neighbors to return, K>=1
377        SelfMatch   -   whether self-matches are allowed:
378                        * if True, nearest neighbor may be the point itself
379                          (if it exists in original dataset)
380                        * if False, then only points with non-zero distance
381                          are returned
382                        * if not given, considered True
383
384    RESULT
385        number of actual neighbors found (either K or N, if K>N).
386
387    This  subroutine  performs  query  and  stores  its result in the internal
388    structures of the KD-tree. You can use  following  subroutines  to  obtain
389    these results:
390    * KDTreeQueryResultsX() to get X-values
391    * KDTreeQueryResultsXY() to get X- and Y-values
392    * KDTreeQueryResultsTags() to get tag values
393    * KDTreeQueryResultsDistances() to get distances
394
395      -- ALGLIB --
396         Copyright 28.02.2010 by Bochkanov Sergey
397    *************************************************************************/
398    public static int kdtreequeryknn(kdtree kdt, double[] x, int k, bool selfmatch)
399    {
400
401        int result = nearestneighbor.kdtreequeryknn(kdt.innerobj, x, k, selfmatch);
402        return result;
403    }
404    public static int kdtreequeryknn(kdtree kdt, double[] x, int k)
405    {
406        bool selfmatch;
407
408
409        selfmatch = true;
410        int result = nearestneighbor.kdtreequeryknn(kdt.innerobj, x, k, selfmatch);
411
412        return result;
413    }
414
415    /*************************************************************************
416    R-NN query: all points within R-sphere centered at X
417
418    INPUT PARAMETERS
419        KDT         -   KD-tree
420        X           -   point, array[0..NX-1].
421        R           -   radius of sphere (in corresponding norm), R>0
422        SelfMatch   -   whether self-matches are allowed:
423                        * if True, nearest neighbor may be the point itself
424                          (if it exists in original dataset)
425                        * if False, then only points with non-zero distance
426                          are returned
427                        * if not given, considered True
428
429    RESULT
430        number of neighbors found, >=0
431
432    This  subroutine  performs  query  and  stores  its result in the internal
433    structures of the KD-tree. You can use  following  subroutines  to  obtain
434    actual results:
435    * KDTreeQueryResultsX() to get X-values
436    * KDTreeQueryResultsXY() to get X- and Y-values
437    * KDTreeQueryResultsTags() to get tag values
438    * KDTreeQueryResultsDistances() to get distances
439
440      -- ALGLIB --
441         Copyright 28.02.2010 by Bochkanov Sergey
442    *************************************************************************/
443    public static int kdtreequeryrnn(kdtree kdt, double[] x, double r, bool selfmatch)
444    {
445
446        int result = nearestneighbor.kdtreequeryrnn(kdt.innerobj, x, r, selfmatch);
447        return result;
448    }
449    public static int kdtreequeryrnn(kdtree kdt, double[] x, double r)
450    {
451        bool selfmatch;
452
453
454        selfmatch = true;
455        int result = nearestneighbor.kdtreequeryrnn(kdt.innerobj, x, r, selfmatch);
456
457        return result;
458    }
459
460    /*************************************************************************
461    K-NN query: approximate K nearest neighbors
462
463    INPUT PARAMETERS
464        KDT         -   KD-tree
465        X           -   point, array[0..NX-1].
466        K           -   number of neighbors to return, K>=1
467        SelfMatch   -   whether self-matches are allowed:
468                        * if True, nearest neighbor may be the point itself
469                          (if it exists in original dataset)
470                        * if False, then only points with non-zero distance
471                          are returned
472                        * if not given, considered True
473        Eps         -   approximation factor, Eps>=0. eps-approximate  nearest
474                        neighbor  is  a  neighbor  whose distance from X is at
475                        most (1+eps) times distance of true nearest neighbor.
476
477    RESULT
478        number of actual neighbors found (either K or N, if K>N).
479
480    NOTES
481        significant performance gain may be achieved only when Eps  is  is  on
482        the order of magnitude of 1 or larger.
483
484    This  subroutine  performs  query  and  stores  its result in the internal
485    structures of the KD-tree. You can use  following  subroutines  to  obtain
486    these results:
487    * KDTreeQueryResultsX() to get X-values
488    * KDTreeQueryResultsXY() to get X- and Y-values
489    * KDTreeQueryResultsTags() to get tag values
490    * KDTreeQueryResultsDistances() to get distances
491
492      -- ALGLIB --
493         Copyright 28.02.2010 by Bochkanov Sergey
494    *************************************************************************/
495    public static int kdtreequeryaknn(kdtree kdt, double[] x, int k, bool selfmatch, double eps)
496    {
497
498        int result = nearestneighbor.kdtreequeryaknn(kdt.innerobj, x, k, selfmatch, eps);
499        return result;
500    }
501    public static int kdtreequeryaknn(kdtree kdt, double[] x, int k, double eps)
502    {
503        bool selfmatch;
504
505
506        selfmatch = true;
507        int result = nearestneighbor.kdtreequeryaknn(kdt.innerobj, x, k, selfmatch, eps);
508
509        return result;
510    }
511
512    /*************************************************************************
513    X-values from last query
514
515    INPUT PARAMETERS
516        KDT     -   KD-tree
517        X       -   possibly pre-allocated buffer. If X is too small to store
518                    result, it is resized. If size(X) is enough to store
519                    result, it is left unchanged.
520
521    OUTPUT PARAMETERS
522        X       -   rows are filled with X-values
523
524    NOTES
525    1. points are ordered by distance from the query point (first = closest)
526    2. if  XY is larger than required to store result, only leading part  will
527       be overwritten; trailing part will be left unchanged. So  if  on  input
528       XY = [[A,B],[C,D]], and result is [1,2],  then  on  exit  we  will  get
529       XY = [[1,2],[C,D]]. This is done purposely to increase performance;  if
530       you want function  to  resize  array  according  to  result  size,  use
531       function with same name and suffix 'I'.
532
533    SEE ALSO
534    * KDTreeQueryResultsXY()            X- and Y-values
535    * KDTreeQueryResultsTags()          tag values
536    * KDTreeQueryResultsDistances()     distances
537
538      -- ALGLIB --
539         Copyright 28.02.2010 by Bochkanov Sergey
540    *************************************************************************/
541    public static void kdtreequeryresultsx(kdtree kdt, ref double[,] x)
542    {
543
544        nearestneighbor.kdtreequeryresultsx(kdt.innerobj, ref x);
545        return;
546    }
547
548    /*************************************************************************
549    X- and Y-values from last query
550
551    INPUT PARAMETERS
552        KDT     -   KD-tree
553        XY      -   possibly pre-allocated buffer. If XY is too small to store
554                    result, it is resized. If size(XY) is enough to store
555                    result, it is left unchanged.
556
557    OUTPUT PARAMETERS
558        XY      -   rows are filled with points: first NX columns with
559                    X-values, next NY columns - with Y-values.
560
561    NOTES
562    1. points are ordered by distance from the query point (first = closest)
563    2. if  XY is larger than required to store result, only leading part  will
564       be overwritten; trailing part will be left unchanged. So  if  on  input
565       XY = [[A,B],[C,D]], and result is [1,2],  then  on  exit  we  will  get
566       XY = [[1,2],[C,D]]. This is done purposely to increase performance;  if
567       you want function  to  resize  array  according  to  result  size,  use
568       function with same name and suffix 'I'.
569
570    SEE ALSO
571    * KDTreeQueryResultsX()             X-values
572    * KDTreeQueryResultsTags()          tag values
573    * KDTreeQueryResultsDistances()     distances
574
575      -- ALGLIB --
576         Copyright 28.02.2010 by Bochkanov Sergey
577    *************************************************************************/
578    public static void kdtreequeryresultsxy(kdtree kdt, ref double[,] xy)
579    {
580
581        nearestneighbor.kdtreequeryresultsxy(kdt.innerobj, ref xy);
582        return;
583    }
584
585    /*************************************************************************
586    Tags from last query
587
588    INPUT PARAMETERS
589        KDT     -   KD-tree
590        Tags    -   possibly pre-allocated buffer. If X is too small to store
591                    result, it is resized. If size(X) is enough to store
592                    result, it is left unchanged.
593
594    OUTPUT PARAMETERS
595        Tags    -   filled with tags associated with points,
596                    or, when no tags were supplied, with zeros
597
598    NOTES
599    1. points are ordered by distance from the query point (first = closest)
600    2. if  XY is larger than required to store result, only leading part  will
601       be overwritten; trailing part will be left unchanged. So  if  on  input
602       XY = [[A,B],[C,D]], and result is [1,2],  then  on  exit  we  will  get
603       XY = [[1,2],[C,D]]. This is done purposely to increase performance;  if
604       you want function  to  resize  array  according  to  result  size,  use
605       function with same name and suffix 'I'.
606
607    SEE ALSO
608    * KDTreeQueryResultsX()             X-values
609    * KDTreeQueryResultsXY()            X- and Y-values
610    * KDTreeQueryResultsDistances()     distances
611
612      -- ALGLIB --
613         Copyright 28.02.2010 by Bochkanov Sergey
614    *************************************************************************/
615    public static void kdtreequeryresultstags(kdtree kdt, ref int[] tags)
616    {
617
618        nearestneighbor.kdtreequeryresultstags(kdt.innerobj, ref tags);
619        return;
620    }
621
622    /*************************************************************************
623    Distances from last query
624
625    INPUT PARAMETERS
626        KDT     -   KD-tree
627        R       -   possibly pre-allocated buffer. If X is too small to store
628                    result, it is resized. If size(X) is enough to store
629                    result, it is left unchanged.
630
631    OUTPUT PARAMETERS
632        R       -   filled with distances (in corresponding norm)
633
634    NOTES
635    1. points are ordered by distance from the query point (first = closest)
636    2. if  XY is larger than required to store result, only leading part  will
637       be overwritten; trailing part will be left unchanged. So  if  on  input
638       XY = [[A,B],[C,D]], and result is [1,2],  then  on  exit  we  will  get
639       XY = [[1,2],[C,D]]. This is done purposely to increase performance;  if
640       you want function  to  resize  array  according  to  result  size,  use
641       function with same name and suffix 'I'.
642
643    SEE ALSO
644    * KDTreeQueryResultsX()             X-values
645    * KDTreeQueryResultsXY()            X- and Y-values
646    * KDTreeQueryResultsTags()          tag values
647
648      -- ALGLIB --
649         Copyright 28.02.2010 by Bochkanov Sergey
650    *************************************************************************/
651    public static void kdtreequeryresultsdistances(kdtree kdt, ref double[] r)
652    {
653
654        nearestneighbor.kdtreequeryresultsdistances(kdt.innerobj, ref r);
655        return;
656    }
657
658    /*************************************************************************
659    X-values from last query; 'interactive' variant for languages like  Python
660    which   support    constructs   like  "X = KDTreeQueryResultsXI(KDT)"  and
661    interactive mode of interpreter.
662
663    This function allocates new array on each call,  so  it  is  significantly
664    slower than its 'non-interactive' counterpart, but it is  more  convenient
665    when you call it from command line.
666
667      -- ALGLIB --
668         Copyright 28.02.2010 by Bochkanov Sergey
669    *************************************************************************/
670    public static void kdtreequeryresultsxi(kdtree kdt, out double[,] x)
671    {
672        x = new double[0,0];
673        nearestneighbor.kdtreequeryresultsxi(kdt.innerobj, ref x);
674        return;
675    }
676
677    /*************************************************************************
678    XY-values from last query; 'interactive' variant for languages like Python
679    which   support    constructs   like "XY = KDTreeQueryResultsXYI(KDT)" and
680    interactive mode of interpreter.
681
682    This function allocates new array on each call,  so  it  is  significantly
683    slower than its 'non-interactive' counterpart, but it is  more  convenient
684    when you call it from command line.
685
686      -- ALGLIB --
687         Copyright 28.02.2010 by Bochkanov Sergey
688    *************************************************************************/
689    public static void kdtreequeryresultsxyi(kdtree kdt, out double[,] xy)
690    {
691        xy = new double[0,0];
692        nearestneighbor.kdtreequeryresultsxyi(kdt.innerobj, ref xy);
693        return;
694    }
695
696    /*************************************************************************
697    Tags  from  last  query;  'interactive' variant for languages like  Python
698    which  support  constructs  like "Tags = KDTreeQueryResultsTagsI(KDT)" and
699    interactive mode of interpreter.
700
701    This function allocates new array on each call,  so  it  is  significantly
702    slower than its 'non-interactive' counterpart, but it is  more  convenient
703    when you call it from command line.
704
705      -- ALGLIB --
706         Copyright 28.02.2010 by Bochkanov Sergey
707    *************************************************************************/
708    public static void kdtreequeryresultstagsi(kdtree kdt, out int[] tags)
709    {
710        tags = new int[0];
711        nearestneighbor.kdtreequeryresultstagsi(kdt.innerobj, ref tags);
712        return;
713    }
714
715    /*************************************************************************
716    Distances from last query; 'interactive' variant for languages like Python
717    which  support  constructs   like  "R = KDTreeQueryResultsDistancesI(KDT)"
718    and interactive mode of interpreter.
719
720    This function allocates new array on each call,  so  it  is  significantly
721    slower than its 'non-interactive' counterpart, but it is  more  convenient
722    when you call it from command line.
723
724      -- ALGLIB --
725         Copyright 28.02.2010 by Bochkanov Sergey
726    *************************************************************************/
727    public static void kdtreequeryresultsdistancesi(kdtree kdt, out double[] r)
728    {
729        r = new double[0];
730        nearestneighbor.kdtreequeryresultsdistancesi(kdt.innerobj, ref r);
731        return;
732    }
733
734}
735public partial class alglib
736{
737    public class hqrnd
738    {
739        /*************************************************************************
740        Portable high quality random number generator state.
741        Initialized with HQRNDRandomize() or HQRNDSeed().
742
743        Fields:
744            S1, S2      -   seed values
745            V           -   precomputed value
746            MagicV      -   'magic' value used to determine whether State structure
747                            was correctly initialized.
748        *************************************************************************/
749        public class hqrndstate
750        {
751            public int s1;
752            public int s2;
753            public double v;
754            public int magicv;
755        };
756
757
758
759
760        public const int hqrndmax = 2147483563;
761        public const int hqrndm1 = 2147483563;
762        public const int hqrndm2 = 2147483399;
763        public const int hqrndmagic = 1634357784;
764
765
766        /*************************************************************************
767        HQRNDState  initialization  with  random  values  which come from standard
768        RNG.
769
770          -- ALGLIB --
771             Copyright 02.12.2009 by Bochkanov Sergey
772        *************************************************************************/
773        public static void hqrndrandomize(hqrndstate state)
774        {
775            hqrndseed(math.randominteger(hqrndm1), math.randominteger(hqrndm2), state);
776        }
777
778
779        /*************************************************************************
780        HQRNDState initialization with seed values
781
782          -- ALGLIB --
783             Copyright 02.12.2009 by Bochkanov Sergey
784        *************************************************************************/
785        public static void hqrndseed(int s1,
786            int s2,
787            hqrndstate state)
788        {
789            state.s1 = s1%(hqrndm1-1)+1;
790            state.s2 = s2%(hqrndm2-1)+1;
791            state.v = (double)1/(double)hqrndmax;
792            state.magicv = hqrndmagic;
793        }
794
795
796        /*************************************************************************
797        This function generates random real number in (0,1),
798        not including interval boundaries
799
800        State structure must be initialized with HQRNDRandomize() or HQRNDSeed().
801
802          -- ALGLIB --
803             Copyright 02.12.2009 by Bochkanov Sergey
804        *************************************************************************/
805        public static double hqrnduniformr(hqrndstate state)
806        {
807            double result = 0;
808
809            result = state.v*hqrndintegerbase(state);
810            return result;
811        }
812
813
814        /*************************************************************************
815        This function generates random integer number in [0, N)
816
817        1. N must be less than HQRNDMax-1.
818        2. State structure must be initialized with HQRNDRandomize() or HQRNDSeed()
819
820          -- ALGLIB --
821             Copyright 02.12.2009 by Bochkanov Sergey
822        *************************************************************************/
823        public static int hqrnduniformi(hqrndstate state,
824            int n)
825        {
826            int result = 0;
827            int mx = 0;
828
829           
830            //
831            // Correct handling of N's close to RNDBaseMax
832            // (avoiding skewed distributions for RNDBaseMax<>K*N)
833            //
834            ap.assert(n>0, "HQRNDUniformI: N<=0!");
835            ap.assert(n<hqrndmax-1, "HQRNDUniformI: N>=RNDBaseMax-1!");
836            mx = hqrndmax-1-(hqrndmax-1)%n;
837            do
838            {
839                result = hqrndintegerbase(state)-1;
840            }
841            while( result>=mx );
842            result = result%n;
843            return result;
844        }
845
846
847        /*************************************************************************
848        Random number generator: normal numbers
849
850        This function generates one random number from normal distribution.
851        Its performance is equal to that of HQRNDNormal2()
852
853        State structure must be initialized with HQRNDRandomize() or HQRNDSeed().
854
855          -- ALGLIB --
856             Copyright 02.12.2009 by Bochkanov Sergey
857        *************************************************************************/
858        public static double hqrndnormal(hqrndstate state)
859        {
860            double result = 0;
861            double v1 = 0;
862            double v2 = 0;
863
864            hqrndnormal2(state, ref v1, ref v2);
865            result = v1;
866            return result;
867        }
868
869
870        /*************************************************************************
871        Random number generator: random X and Y such that X^2+Y^2=1
872
873        State structure must be initialized with HQRNDRandomize() or HQRNDSeed().
874
875          -- ALGLIB --
876             Copyright 02.12.2009 by Bochkanov Sergey
877        *************************************************************************/
878        public static void hqrndunit2(hqrndstate state,
879            ref double x,
880            ref double y)
881        {
882            double v = 0;
883            double mx = 0;
884            double mn = 0;
885
886            x = 0;
887            y = 0;
888
889            do
890            {
891                hqrndnormal2(state, ref x, ref y);
892            }
893            while( !((double)(x)!=(double)(0) | (double)(y)!=(double)(0)) );
894            mx = Math.Max(Math.Abs(x), Math.Abs(y));
895            mn = Math.Min(Math.Abs(x), Math.Abs(y));
896            v = mx*Math.Sqrt(1+math.sqr(mn/mx));
897            x = x/v;
898            y = y/v;
899        }
900
901
902        /*************************************************************************
903        Random number generator: normal numbers
904
905        This function generates two independent random numbers from normal
906        distribution. Its performance is equal to that of HQRNDNormal()
907
908        State structure must be initialized with HQRNDRandomize() or HQRNDSeed().
909
910          -- ALGLIB --
911             Copyright 02.12.2009 by Bochkanov Sergey
912        *************************************************************************/
913        public static void hqrndnormal2(hqrndstate state,
914            ref double x1,
915            ref double x2)
916        {
917            double u = 0;
918            double v = 0;
919            double s = 0;
920
921            x1 = 0;
922            x2 = 0;
923
924            while( true )
925            {
926                u = 2*hqrnduniformr(state)-1;
927                v = 2*hqrnduniformr(state)-1;
928                s = math.sqr(u)+math.sqr(v);
929                if( (double)(s)>(double)(0) & (double)(s)<(double)(1) )
930                {
931                   
932                    //
933                    // two Sqrt's instead of one to
934                    // avoid overflow when S is too small
935                    //
936                    s = Math.Sqrt(-(2*Math.Log(s)))/Math.Sqrt(s);
937                    x1 = u*s;
938                    x2 = v*s;
939                    return;
940                }
941            }
942        }
943
944
945        /*************************************************************************
946        Random number generator: exponential distribution
947
948        State structure must be initialized with HQRNDRandomize() or HQRNDSeed().
949
950          -- ALGLIB --
951             Copyright 11.08.2007 by Bochkanov Sergey
952        *************************************************************************/
953        public static double hqrndexponential(hqrndstate state,
954            double lambdav)
955        {
956            double result = 0;
957
958            ap.assert((double)(lambdav)>(double)(0), "HQRNDExponential: LambdaV<=0!");
959            result = -(Math.Log(hqrnduniformr(state))/lambdav);
960            return result;
961        }
962
963
964        /*************************************************************************
965
966        L'Ecuyer, Efficient and portable combined random number generators
967        *************************************************************************/
968        private static int hqrndintegerbase(hqrndstate state)
969        {
970            int result = 0;
971            int k = 0;
972
973            ap.assert(state.magicv==hqrndmagic, "HQRNDIntegerBase: State is not correctly initialized!");
974            k = state.s1/53668;
975            state.s1 = 40014*(state.s1-k*53668)-k*12211;
976            if( state.s1<0 )
977            {
978                state.s1 = state.s1+2147483563;
979            }
980            k = state.s2/52774;
981            state.s2 = 40692*(state.s2-k*52774)-k*3791;
982            if( state.s2<0 )
983            {
984                state.s2 = state.s2+2147483399;
985            }
986           
987            //
988            // Result
989            //
990            result = state.s1-state.s2;
991            if( result<1 )
992            {
993                result = result+2147483562;
994            }
995            return result;
996        }
997
998
999    }
1000    public class nearestneighbor
1001    {
1002        public class kdtree
1003        {
1004            public int n;
1005            public int nx;
1006            public int ny;
1007            public int normtype;
1008            public double[,] xy;
1009            public int[] tags;
1010            public double[] boxmin;
1011            public double[] boxmax;
1012            public int[] nodes;
1013            public double[] splits;
1014            public double[] x;
1015            public int kneeded;
1016            public double rneeded;
1017            public bool selfmatch;
1018            public double approxf;
1019            public int kcur;
1020            public int[] idx;
1021            public double[] r;
1022            public double[] buf;
1023            public double[] curboxmin;
1024            public double[] curboxmax;
1025            public double curdist;
1026            public int debugcounter;
1027            public kdtree()
1028            {
1029                xy = new double[0,0];
1030                tags = new int[0];
1031                boxmin = new double[0];
1032                boxmax = new double[0];
1033                nodes = new int[0];
1034                splits = new double[0];
1035                x = new double[0];
1036                idx = new int[0];
1037                r = new double[0];
1038                buf = new double[0];
1039                curboxmin = new double[0];
1040                curboxmax = new double[0];
1041            }
1042        };
1043
1044
1045
1046
1047        public const int splitnodesize = 6;
1048        public const int kdtreefirstversion = 0;
1049
1050
1051        /*************************************************************************
1052        KD-tree creation
1053
1054        This subroutine creates KD-tree from set of X-values and optional Y-values
1055
1056        INPUT PARAMETERS
1057            XY      -   dataset, array[0..N-1,0..NX+NY-1].
1058                        one row corresponds to one point.
1059                        first NX columns contain X-values, next NY (NY may be zero)
1060                        columns may contain associated Y-values
1061            N       -   number of points, N>=1
1062            NX      -   space dimension, NX>=1.
1063            NY      -   number of optional Y-values, NY>=0.
1064            NormType-   norm type:
1065                        * 0 denotes infinity-norm
1066                        * 1 denotes 1-norm
1067                        * 2 denotes 2-norm (Euclidean norm)
1068                       
1069        OUTPUT PARAMETERS
1070            KDT     -   KD-tree
1071           
1072           
1073        NOTES
1074
1075        1. KD-tree  creation  have O(N*logN) complexity and O(N*(2*NX+NY))  memory
1076           requirements.
1077        2. Although KD-trees may be used with any combination of N  and  NX,  they
1078           are more efficient than brute-force search only when N >> 4^NX. So they
1079           are most useful in low-dimensional tasks (NX=2, NX=3). NX=1  is another
1080           inefficient case, because  simple  binary  search  (without  additional
1081           structures) is much more efficient in such tasks than KD-trees.
1082
1083          -- ALGLIB --
1084             Copyright 28.02.2010 by Bochkanov Sergey
1085        *************************************************************************/
1086        public static void kdtreebuild(double[,] xy,
1087            int n,
1088            int nx,
1089            int ny,
1090            int normtype,
1091            kdtree kdt)
1092        {
1093            int[] tags = new int[0];
1094            int i = 0;
1095
1096            ap.assert(n>=1, "KDTreeBuild: N<1!");
1097            ap.assert(nx>=1, "KDTreeBuild: NX<1!");
1098            ap.assert(ny>=0, "KDTreeBuild: NY<0!");
1099            ap.assert(normtype>=0 & normtype<=2, "KDTreeBuild: incorrect NormType!");
1100            ap.assert(ap.rows(xy)>=n, "KDTreeBuild: rows(X)<N!");
1101            ap.assert(ap.cols(xy)>=nx+ny, "KDTreeBuild: cols(X)<NX+NY!");
1102            ap.assert(apserv.apservisfinitematrix(xy, n, nx+ny), "KDTreeBuild: X contains infinite or NaN values!");
1103            tags = new int[n];
1104            for(i=0; i<=n-1; i++)
1105            {
1106                tags[i] = 0;
1107            }
1108            kdtreebuildtagged(xy, tags, n, nx, ny, normtype, kdt);
1109        }
1110
1111
1112        /*************************************************************************
1113        KD-tree creation
1114
1115        This  subroutine  creates  KD-tree  from set of X-values, integer tags and
1116        optional Y-values
1117
1118        INPUT PARAMETERS
1119            XY      -   dataset, array[0..N-1,0..NX+NY-1].
1120                        one row corresponds to one point.
1121                        first NX columns contain X-values, next NY (NY may be zero)
1122                        columns may contain associated Y-values
1123            Tags    -   tags, array[0..N-1], contains integer tags associated
1124                        with points.
1125            N       -   number of points, N>=1
1126            NX      -   space dimension, NX>=1.
1127            NY      -   number of optional Y-values, NY>=0.
1128            NormType-   norm type:
1129                        * 0 denotes infinity-norm
1130                        * 1 denotes 1-norm
1131                        * 2 denotes 2-norm (Euclidean norm)
1132
1133        OUTPUT PARAMETERS
1134            KDT     -   KD-tree
1135
1136        NOTES
1137
1138        1. KD-tree  creation  have O(N*logN) complexity and O(N*(2*NX+NY))  memory
1139           requirements.
1140        2. Although KD-trees may be used with any combination of N  and  NX,  they
1141           are more efficient than brute-force search only when N >> 4^NX. So they
1142           are most useful in low-dimensional tasks (NX=2, NX=3). NX=1  is another
1143           inefficient case, because  simple  binary  search  (without  additional
1144           structures) is much more efficient in such tasks than KD-trees.
1145
1146          -- ALGLIB --
1147             Copyright 28.02.2010 by Bochkanov Sergey
1148        *************************************************************************/
1149        public static void kdtreebuildtagged(double[,] xy,
1150            int[] tags,
1151            int n,
1152            int nx,
1153            int ny,
1154            int normtype,
1155            kdtree kdt)
1156        {
1157            int i = 0;
1158            int j = 0;
1159            int maxnodes = 0;
1160            int nodesoffs = 0;
1161            int splitsoffs = 0;
1162            int i_ = 0;
1163            int i1_ = 0;
1164
1165            ap.assert(n>=1, "KDTreeBuildTagged: N<1!");
1166            ap.assert(nx>=1, "KDTreeBuildTagged: NX<1!");
1167            ap.assert(ny>=0, "KDTreeBuildTagged: NY<0!");
1168            ap.assert(normtype>=0 & normtype<=2, "KDTreeBuildTagged: incorrect NormType!");
1169            ap.assert(ap.rows(xy)>=n, "KDTreeBuildTagged: rows(X)<N!");
1170            ap.assert(ap.cols(xy)>=nx+ny, "KDTreeBuildTagged: cols(X)<NX+NY!");
1171            ap.assert(apserv.apservisfinitematrix(xy, n, nx+ny), "KDTreeBuildTagged: X contains infinite or NaN values!");
1172           
1173            //
1174            // initialize
1175            //
1176            kdt.n = n;
1177            kdt.nx = nx;
1178            kdt.ny = ny;
1179            kdt.normtype = normtype;
1180           
1181            //
1182            // Allocate
1183            //
1184            kdtreeallocdatasetindependent(kdt, nx, ny);
1185            kdtreeallocdatasetdependent(kdt, n, nx, ny);
1186           
1187            //
1188            // Initial fill
1189            //
1190            for(i=0; i<=n-1; i++)
1191            {
1192                for(i_=0; i_<=nx-1;i_++)
1193                {
1194                    kdt.xy[i,i_] = xy[i,i_];
1195                }
1196                i1_ = (0) - (nx);
1197                for(i_=nx; i_<=2*nx+ny-1;i_++)
1198                {
1199                    kdt.xy[i,i_] = xy[i,i_+i1_];
1200                }
1201                kdt.tags[i] = tags[i];
1202            }
1203           
1204            //
1205            // Determine bounding box
1206            //
1207            for(i_=0; i_<=nx-1;i_++)
1208            {
1209                kdt.boxmin[i_] = kdt.xy[0,i_];
1210            }
1211            for(i_=0; i_<=nx-1;i_++)
1212            {
1213                kdt.boxmax[i_] = kdt.xy[0,i_];
1214            }
1215            for(i=1; i<=n-1; i++)
1216            {
1217                for(j=0; j<=nx-1; j++)
1218                {
1219                    kdt.boxmin[j] = Math.Min(kdt.boxmin[j], kdt.xy[i,j]);
1220                    kdt.boxmax[j] = Math.Max(kdt.boxmax[j], kdt.xy[i,j]);
1221                }
1222            }
1223           
1224            //
1225            // prepare tree structure
1226            // * MaxNodes=N because we guarantee no trivial splits, i.e.
1227            //   every split will generate two non-empty boxes
1228            //
1229            maxnodes = n;
1230            kdt.nodes = new int[splitnodesize*2*maxnodes];
1231            kdt.splits = new double[2*maxnodes];
1232            nodesoffs = 0;
1233            splitsoffs = 0;
1234            for(i_=0; i_<=nx-1;i_++)
1235            {
1236                kdt.curboxmin[i_] = kdt.boxmin[i_];
1237            }
1238            for(i_=0; i_<=nx-1;i_++)
1239            {
1240                kdt.curboxmax[i_] = kdt.boxmax[i_];
1241            }
1242            kdtreegeneratetreerec(kdt, ref nodesoffs, ref splitsoffs, 0, n, 8);
1243           
1244            //
1245            // Set current query size to 0
1246            //
1247            kdt.kcur = 0;
1248        }
1249
1250
1251        /*************************************************************************
1252        K-NN query: K nearest neighbors
1253
1254        INPUT PARAMETERS
1255            KDT         -   KD-tree
1256            X           -   point, array[0..NX-1].
1257            K           -   number of neighbors to return, K>=1
1258            SelfMatch   -   whether self-matches are allowed:
1259                            * if True, nearest neighbor may be the point itself
1260                              (if it exists in original dataset)
1261                            * if False, then only points with non-zero distance
1262                              are returned
1263                            * if not given, considered True
1264
1265        RESULT
1266            number of actual neighbors found (either K or N, if K>N).
1267
1268        This  subroutine  performs  query  and  stores  its result in the internal
1269        structures of the KD-tree. You can use  following  subroutines  to  obtain
1270        these results:
1271        * KDTreeQueryResultsX() to get X-values
1272        * KDTreeQueryResultsXY() to get X- and Y-values
1273        * KDTreeQueryResultsTags() to get tag values
1274        * KDTreeQueryResultsDistances() to get distances
1275
1276          -- ALGLIB --
1277             Copyright 28.02.2010 by Bochkanov Sergey
1278        *************************************************************************/
1279        public static int kdtreequeryknn(kdtree kdt,
1280            double[] x,
1281            int k,
1282            bool selfmatch)
1283        {
1284            int result = 0;
1285
1286            ap.assert(k>=1, "KDTreeQueryKNN: K<1!");
1287            ap.assert(ap.len(x)>=kdt.nx, "KDTreeQueryKNN: Length(X)<NX!");
1288            ap.assert(apserv.isfinitevector(x, kdt.nx), "KDTreeQueryKNN: X contains infinite or NaN values!");
1289            result = kdtreequeryaknn(kdt, x, k, selfmatch, 0.0);
1290            return result;
1291        }
1292
1293
1294        /*************************************************************************
1295        R-NN query: all points within R-sphere centered at X
1296
1297        INPUT PARAMETERS
1298            KDT         -   KD-tree
1299            X           -   point, array[0..NX-1].
1300            R           -   radius of sphere (in corresponding norm), R>0
1301            SelfMatch   -   whether self-matches are allowed:
1302                            * if True, nearest neighbor may be the point itself
1303                              (if it exists in original dataset)
1304                            * if False, then only points with non-zero distance
1305                              are returned
1306                            * if not given, considered True
1307
1308        RESULT
1309            number of neighbors found, >=0
1310
1311        This  subroutine  performs  query  and  stores  its result in the internal
1312        structures of the KD-tree. You can use  following  subroutines  to  obtain
1313        actual results:
1314        * KDTreeQueryResultsX() to get X-values
1315        * KDTreeQueryResultsXY() to get X- and Y-values
1316        * KDTreeQueryResultsTags() to get tag values
1317        * KDTreeQueryResultsDistances() to get distances
1318
1319          -- ALGLIB --
1320             Copyright 28.02.2010 by Bochkanov Sergey
1321        *************************************************************************/
1322        public static int kdtreequeryrnn(kdtree kdt,
1323            double[] x,
1324            double r,
1325            bool selfmatch)
1326        {
1327            int result = 0;
1328            int i = 0;
1329            int j = 0;
1330
1331            ap.assert((double)(r)>(double)(0), "KDTreeQueryRNN: incorrect R!");
1332            ap.assert(ap.len(x)>=kdt.nx, "KDTreeQueryRNN: Length(X)<NX!");
1333            ap.assert(apserv.isfinitevector(x, kdt.nx), "KDTreeQueryRNN: X contains infinite or NaN values!");
1334           
1335            //
1336            // Prepare parameters
1337            //
1338            kdt.kneeded = 0;
1339            if( kdt.normtype!=2 )
1340            {
1341                kdt.rneeded = r;
1342            }
1343            else
1344            {
1345                kdt.rneeded = math.sqr(r);
1346            }
1347            kdt.selfmatch = selfmatch;
1348            kdt.approxf = 1;
1349            kdt.kcur = 0;
1350           
1351            //
1352            // calculate distance from point to current bounding box
1353            //
1354            kdtreeinitbox(kdt, x);
1355           
1356            //
1357            // call recursive search
1358            // results are returned as heap
1359            //
1360            kdtreequerynnrec(kdt, 0);
1361           
1362            //
1363            // pop from heap to generate ordered representation
1364            //
1365            // last element is not pop'ed because it is already in
1366            // its place
1367            //
1368            result = kdt.kcur;
1369            j = kdt.kcur;
1370            for(i=kdt.kcur; i>=2; i--)
1371            {
1372                tsort.tagheappopi(ref kdt.r, ref kdt.idx, ref j);
1373            }
1374            return result;
1375        }
1376
1377
1378        /*************************************************************************
1379        K-NN query: approximate K nearest neighbors
1380
1381        INPUT PARAMETERS
1382            KDT         -   KD-tree
1383            X           -   point, array[0..NX-1].
1384            K           -   number of neighbors to return, K>=1
1385            SelfMatch   -   whether self-matches are allowed:
1386                            * if True, nearest neighbor may be the point itself
1387                              (if it exists in original dataset)
1388                            * if False, then only points with non-zero distance
1389                              are returned
1390                            * if not given, considered True
1391            Eps         -   approximation factor, Eps>=0. eps-approximate  nearest
1392                            neighbor  is  a  neighbor  whose distance from X is at
1393                            most (1+eps) times distance of true nearest neighbor.
1394
1395        RESULT
1396            number of actual neighbors found (either K or N, if K>N).
1397           
1398        NOTES
1399            significant performance gain may be achieved only when Eps  is  is  on
1400            the order of magnitude of 1 or larger.
1401
1402        This  subroutine  performs  query  and  stores  its result in the internal
1403        structures of the KD-tree. You can use  following  subroutines  to  obtain
1404        these results:
1405        * KDTreeQueryResultsX() to get X-values
1406        * KDTreeQueryResultsXY() to get X- and Y-values
1407        * KDTreeQueryResultsTags() to get tag values
1408        * KDTreeQueryResultsDistances() to get distances
1409
1410          -- ALGLIB --
1411             Copyright 28.02.2010 by Bochkanov Sergey
1412        *************************************************************************/
1413        public static int kdtreequeryaknn(kdtree kdt,
1414            double[] x,
1415            int k,
1416            bool selfmatch,
1417            double eps)
1418        {
1419            int result = 0;
1420            int i = 0;
1421            int j = 0;
1422
1423            ap.assert(k>0, "KDTreeQueryAKNN: incorrect K!");
1424            ap.assert((double)(eps)>=(double)(0), "KDTreeQueryAKNN: incorrect Eps!");
1425            ap.assert(ap.len(x)>=kdt.nx, "KDTreeQueryAKNN: Length(X)<NX!");
1426            ap.assert(apserv.isfinitevector(x, kdt.nx), "KDTreeQueryAKNN: X contains infinite or NaN values!");
1427           
1428            //
1429            // Prepare parameters
1430            //
1431            k = Math.Min(k, kdt.n);
1432            kdt.kneeded = k;
1433            kdt.rneeded = 0;
1434            kdt.selfmatch = selfmatch;
1435            if( kdt.normtype==2 )
1436            {
1437                kdt.approxf = 1/math.sqr(1+eps);
1438            }
1439            else
1440            {
1441                kdt.approxf = 1/(1+eps);
1442            }
1443            kdt.kcur = 0;
1444           
1445            //
1446            // calculate distance from point to current bounding box
1447            //
1448            kdtreeinitbox(kdt, x);
1449           
1450            //
1451            // call recursive search
1452            // results are returned as heap
1453            //
1454            kdtreequerynnrec(kdt, 0);
1455           
1456            //
1457            // pop from heap to generate ordered representation
1458            //
1459            // last element is non pop'ed because it is already in
1460            // its place
1461            //
1462            result = kdt.kcur;
1463            j = kdt.kcur;
1464            for(i=kdt.kcur; i>=2; i--)
1465            {
1466                tsort.tagheappopi(ref kdt.r, ref kdt.idx, ref j);
1467            }
1468            return result;
1469        }
1470
1471
1472        /*************************************************************************
1473        X-values from last query
1474
1475        INPUT PARAMETERS
1476            KDT     -   KD-tree
1477            X       -   possibly pre-allocated buffer. If X is too small to store
1478                        result, it is resized. If size(X) is enough to store
1479                        result, it is left unchanged.
1480
1481        OUTPUT PARAMETERS
1482            X       -   rows are filled with X-values
1483
1484        NOTES
1485        1. points are ordered by distance from the query point (first = closest)
1486        2. if  XY is larger than required to store result, only leading part  will
1487           be overwritten; trailing part will be left unchanged. So  if  on  input
1488           XY = [[A,B],[C,D]], and result is [1,2],  then  on  exit  we  will  get
1489           XY = [[1,2],[C,D]]. This is done purposely to increase performance;  if
1490           you want function  to  resize  array  according  to  result  size,  use
1491           function with same name and suffix 'I'.
1492
1493        SEE ALSO
1494        * KDTreeQueryResultsXY()            X- and Y-values
1495        * KDTreeQueryResultsTags()          tag values
1496        * KDTreeQueryResultsDistances()     distances
1497
1498          -- ALGLIB --
1499             Copyright 28.02.2010 by Bochkanov Sergey
1500        *************************************************************************/
1501        public static void kdtreequeryresultsx(kdtree kdt,
1502            ref double[,] x)
1503        {
1504            int i = 0;
1505            int k = 0;
1506            int i_ = 0;
1507            int i1_ = 0;
1508
1509            if( kdt.kcur==0 )
1510            {
1511                return;
1512            }
1513            if( ap.rows(x)<kdt.kcur | ap.cols(x)<kdt.nx )
1514            {
1515                x = new double[kdt.kcur, kdt.nx];
1516            }
1517            k = kdt.kcur;
1518            for(i=0; i<=k-1; i++)
1519            {
1520                i1_ = (kdt.nx) - (0);
1521                for(i_=0; i_<=kdt.nx-1;i_++)
1522                {
1523                    x[i,i_] = kdt.xy[kdt.idx[i],i_+i1_];
1524                }
1525            }
1526        }
1527
1528
1529        /*************************************************************************
1530        X- and Y-values from last query
1531
1532        INPUT PARAMETERS
1533            KDT     -   KD-tree
1534            XY      -   possibly pre-allocated buffer. If XY is too small to store
1535                        result, it is resized. If size(XY) is enough to store
1536                        result, it is left unchanged.
1537
1538        OUTPUT PARAMETERS
1539            XY      -   rows are filled with points: first NX columns with
1540                        X-values, next NY columns - with Y-values.
1541
1542        NOTES
1543        1. points are ordered by distance from the query point (first = closest)
1544        2. if  XY is larger than required to store result, only leading part  will
1545           be overwritten; trailing part will be left unchanged. So  if  on  input
1546           XY = [[A,B],[C,D]], and result is [1,2],  then  on  exit  we  will  get
1547           XY = [[1,2],[C,D]]. This is done purposely to increase performance;  if
1548           you want function  to  resize  array  according  to  result  size,  use
1549           function with same name and suffix 'I'.
1550
1551        SEE ALSO
1552        * KDTreeQueryResultsX()             X-values
1553        * KDTreeQueryResultsTags()          tag values
1554        * KDTreeQueryResultsDistances()     distances
1555
1556          -- ALGLIB --
1557             Copyright 28.02.2010 by Bochkanov Sergey
1558        *************************************************************************/
1559        public static void kdtreequeryresultsxy(kdtree kdt,
1560            ref double[,] xy)
1561        {
1562            int i = 0;
1563            int k = 0;
1564            int i_ = 0;
1565            int i1_ = 0;
1566
1567            if( kdt.kcur==0 )
1568            {
1569                return;
1570            }
1571            if( ap.rows(xy)<kdt.kcur | ap.cols(xy)<kdt.nx+kdt.ny )
1572            {
1573                xy = new double[kdt.kcur, kdt.nx+kdt.ny];
1574            }
1575            k = kdt.kcur;
1576            for(i=0; i<=k-1; i++)
1577            {
1578                i1_ = (kdt.nx) - (0);
1579                for(i_=0; i_<=kdt.nx+kdt.ny-1;i_++)
1580                {
1581                    xy[i,i_] = kdt.xy[kdt.idx[i],i_+i1_];
1582                }
1583            }
1584        }
1585
1586
1587        /*************************************************************************
1588        Tags from last query
1589
1590        INPUT PARAMETERS
1591            KDT     -   KD-tree
1592            Tags    -   possibly pre-allocated buffer. If X is too small to store
1593                        result, it is resized. If size(X) is enough to store
1594                        result, it is left unchanged.
1595
1596        OUTPUT PARAMETERS
1597            Tags    -   filled with tags associated with points,
1598                        or, when no tags were supplied, with zeros
1599
1600        NOTES
1601        1. points are ordered by distance from the query point (first = closest)
1602        2. if  XY is larger than required to store result, only leading part  will
1603           be overwritten; trailing part will be left unchanged. So  if  on  input
1604           XY = [[A,B],[C,D]], and result is [1,2],  then  on  exit  we  will  get
1605           XY = [[1,2],[C,D]]. This is done purposely to increase performance;  if
1606           you want function  to  resize  array  according  to  result  size,  use
1607           function with same name and suffix 'I'.
1608
1609        SEE ALSO
1610        * KDTreeQueryResultsX()             X-values
1611        * KDTreeQueryResultsXY()            X- and Y-values
1612        * KDTreeQueryResultsDistances()     distances
1613
1614          -- ALGLIB --
1615             Copyright 28.02.2010 by Bochkanov Sergey
1616        *************************************************************************/
1617        public static void kdtreequeryresultstags(kdtree kdt,
1618            ref int[] tags)
1619        {
1620            int i = 0;
1621            int k = 0;
1622
1623            if( kdt.kcur==0 )
1624            {
1625                return;
1626            }
1627            if( ap.len(tags)<kdt.kcur )
1628            {
1629                tags = new int[kdt.kcur];
1630            }
1631            k = kdt.kcur;
1632            for(i=0; i<=k-1; i++)
1633            {
1634                tags[i] = kdt.tags[kdt.idx[i]];
1635            }
1636        }
1637
1638
1639        /*************************************************************************
1640        Distances from last query
1641
1642        INPUT PARAMETERS
1643            KDT     -   KD-tree
1644            R       -   possibly pre-allocated buffer. If X is too small to store
1645                        result, it is resized. If size(X) is enough to store
1646                        result, it is left unchanged.
1647
1648        OUTPUT PARAMETERS
1649            R       -   filled with distances (in corresponding norm)
1650
1651        NOTES
1652        1. points are ordered by distance from the query point (first = closest)
1653        2. if  XY is larger than required to store result, only leading part  will
1654           be overwritten; trailing part will be left unchanged. So  if  on  input
1655           XY = [[A,B],[C,D]], and result is [1,2],  then  on  exit  we  will  get
1656           XY = [[1,2],[C,D]]. This is done purposely to increase performance;  if
1657           you want function  to  resize  array  according  to  result  size,  use
1658           function with same name and suffix 'I'.
1659
1660        SEE ALSO
1661        * KDTreeQueryResultsX()             X-values
1662        * KDTreeQueryResultsXY()            X- and Y-values
1663        * KDTreeQueryResultsTags()          tag values
1664
1665          -- ALGLIB --
1666             Copyright 28.02.2010 by Bochkanov Sergey
1667        *************************************************************************/
1668        public static void kdtreequeryresultsdistances(kdtree kdt,
1669            ref double[] r)
1670        {
1671            int i = 0;
1672            int k = 0;
1673
1674            if( kdt.kcur==0 )
1675            {
1676                return;
1677            }
1678            if( ap.len(r)<kdt.kcur )
1679            {
1680                r = new double[kdt.kcur];
1681            }
1682            k = kdt.kcur;
1683           
1684            //
1685            // unload norms
1686            //
1687            // Abs() call is used to handle cases with negative norms
1688            // (generated during KFN requests)
1689            //
1690            if( kdt.normtype==0 )
1691            {
1692                for(i=0; i<=k-1; i++)
1693                {
1694                    r[i] = Math.Abs(kdt.r[i]);
1695                }
1696            }
1697            if( kdt.normtype==1 )
1698            {
1699                for(i=0; i<=k-1; i++)
1700                {
1701                    r[i] = Math.Abs(kdt.r[i]);
1702                }
1703            }
1704            if( kdt.normtype==2 )
1705            {
1706                for(i=0; i<=k-1; i++)
1707                {
1708                    r[i] = Math.Sqrt(Math.Abs(kdt.r[i]));
1709                }
1710            }
1711        }
1712
1713
1714        /*************************************************************************
1715        X-values from last query; 'interactive' variant for languages like  Python
1716        which   support    constructs   like  "X = KDTreeQueryResultsXI(KDT)"  and
1717        interactive mode of interpreter.
1718
1719        This function allocates new array on each call,  so  it  is  significantly
1720        slower than its 'non-interactive' counterpart, but it is  more  convenient
1721        when you call it from command line.
1722
1723          -- ALGLIB --
1724             Copyright 28.02.2010 by Bochkanov Sergey
1725        *************************************************************************/
1726        public static void kdtreequeryresultsxi(kdtree kdt,
1727            ref double[,] x)
1728        {
1729            x = new double[0,0];
1730
1731            kdtreequeryresultsx(kdt, ref x);
1732        }
1733
1734
1735        /*************************************************************************
1736        XY-values from last query; 'interactive' variant for languages like Python
1737        which   support    constructs   like "XY = KDTreeQueryResultsXYI(KDT)" and
1738        interactive mode of interpreter.
1739
1740        This function allocates new array on each call,  so  it  is  significantly
1741        slower than its 'non-interactive' counterpart, but it is  more  convenient
1742        when you call it from command line.
1743
1744          -- ALGLIB --
1745             Copyright 28.02.2010 by Bochkanov Sergey
1746        *************************************************************************/
1747        public static void kdtreequeryresultsxyi(kdtree kdt,
1748            ref double[,] xy)
1749        {
1750            xy = new double[0,0];
1751
1752            kdtreequeryresultsxy(kdt, ref xy);
1753        }
1754
1755
1756        /*************************************************************************
1757        Tags  from  last  query;  'interactive' variant for languages like  Python
1758        which  support  constructs  like "Tags = KDTreeQueryResultsTagsI(KDT)" and
1759        interactive mode of interpreter.
1760
1761        This function allocates new array on each call,  so  it  is  significantly
1762        slower than its 'non-interactive' counterpart, but it is  more  convenient
1763        when you call it from command line.
1764
1765          -- ALGLIB --
1766             Copyright 28.02.2010 by Bochkanov Sergey
1767        *************************************************************************/
1768        public static void kdtreequeryresultstagsi(kdtree kdt,
1769            ref int[] tags)
1770        {
1771            tags = new int[0];
1772
1773            kdtreequeryresultstags(kdt, ref tags);
1774        }
1775
1776
1777        /*************************************************************************
1778        Distances from last query; 'interactive' variant for languages like Python
1779        which  support  constructs   like  "R = KDTreeQueryResultsDistancesI(KDT)"
1780        and interactive mode of interpreter.
1781
1782        This function allocates new array on each call,  so  it  is  significantly
1783        slower than its 'non-interactive' counterpart, but it is  more  convenient
1784        when you call it from command line.
1785
1786          -- ALGLIB --
1787             Copyright 28.02.2010 by Bochkanov Sergey
1788        *************************************************************************/
1789        public static void kdtreequeryresultsdistancesi(kdtree kdt,
1790            ref double[] r)
1791        {
1792            r = new double[0];
1793
1794            kdtreequeryresultsdistances(kdt, ref r);
1795        }
1796
1797
1798        /*************************************************************************
1799        Serializer: allocation
1800
1801          -- ALGLIB --
1802             Copyright 14.03.2011 by Bochkanov Sergey
1803        *************************************************************************/
1804        public static void kdtreealloc(alglib.serializer s,
1805            kdtree tree)
1806        {
1807           
1808            //
1809            // Header
1810            //
1811            s.alloc_entry();
1812            s.alloc_entry();
1813           
1814            //
1815            // Data
1816            //
1817            s.alloc_entry();
1818            s.alloc_entry();
1819            s.alloc_entry();
1820            s.alloc_entry();
1821            apserv.allocrealmatrix(s, tree.xy, -1, -1);
1822            apserv.allocintegerarray(s, tree.tags, -1);
1823            apserv.allocrealarray(s, tree.boxmin, -1);
1824            apserv.allocrealarray(s, tree.boxmax, -1);
1825            apserv.allocintegerarray(s, tree.nodes, -1);
1826            apserv.allocrealarray(s, tree.splits, -1);
1827        }
1828
1829
1830        /*************************************************************************
1831        Serializer: serialization
1832
1833          -- ALGLIB --
1834             Copyright 14.03.2011 by Bochkanov Sergey
1835        *************************************************************************/
1836        public static void kdtreeserialize(alglib.serializer s,
1837            kdtree tree)
1838        {
1839           
1840            //
1841            // Header
1842            //
1843            s.serialize_int(scodes.getkdtreeserializationcode());
1844            s.serialize_int(kdtreefirstversion);
1845           
1846            //
1847            // Data
1848            //
1849            s.serialize_int(tree.n);
1850            s.serialize_int(tree.nx);
1851            s.serialize_int(tree.ny);
1852            s.serialize_int(tree.normtype);
1853            apserv.serializerealmatrix(s, tree.xy, -1, -1);
1854            apserv.serializeintegerarray(s, tree.tags, -1);
1855            apserv.serializerealarray(s, tree.boxmin, -1);
1856            apserv.serializerealarray(s, tree.boxmax, -1);
1857            apserv.serializeintegerarray(s, tree.nodes, -1);
1858            apserv.serializerealarray(s, tree.splits, -1);
1859        }
1860
1861
1862        /*************************************************************************
1863        Serializer: unserialization
1864
1865          -- ALGLIB --
1866             Copyright 14.03.2011 by Bochkanov Sergey
1867        *************************************************************************/
1868        public static void kdtreeunserialize(alglib.serializer s,
1869            kdtree tree)
1870        {
1871            int i0 = 0;
1872            int i1 = 0;
1873
1874           
1875            //
1876            // check correctness of header
1877            //
1878            i0 = s.unserialize_int();
1879            ap.assert(i0==scodes.getkdtreeserializationcode(), "KDTreeUnserialize: stream header corrupted");
1880            i1 = s.unserialize_int();
1881            ap.assert(i1==kdtreefirstversion, "KDTreeUnserialize: stream header corrupted");
1882           
1883            //
1884            // Unserialize data
1885            //
1886            tree.n = s.unserialize_int();
1887            tree.nx = s.unserialize_int();
1888            tree.ny = s.unserialize_int();
1889            tree.normtype = s.unserialize_int();
1890            apserv.unserializerealmatrix(s, ref tree.xy);
1891            apserv.unserializeintegerarray(s, ref tree.tags);
1892            apserv.unserializerealarray(s, ref tree.boxmin);
1893            apserv.unserializerealarray(s, ref tree.boxmax);
1894            apserv.unserializeintegerarray(s, ref tree.nodes);
1895            apserv.unserializerealarray(s, ref tree.splits);
1896            kdtreealloctemporaries(tree, tree.n, tree.nx, tree.ny);
1897        }
1898
1899
1900        /*************************************************************************
1901        Rearranges nodes [I1,I2) using partition in D-th dimension with S as threshold.
1902        Returns split position I3: [I1,I3) and [I3,I2) are created as result.
1903
1904        This subroutine doesn't create tree structures, just rearranges nodes.
1905        *************************************************************************/
1906        private static void kdtreesplit(kdtree kdt,
1907            int i1,
1908            int i2,
1909            int d,
1910            double s,
1911            ref int i3)
1912        {
1913            int i = 0;
1914            int j = 0;
1915            int ileft = 0;
1916            int iright = 0;
1917            double v = 0;
1918
1919            i3 = 0;
1920
1921           
1922            //
1923            // split XY/Tags in two parts:
1924            // * [ILeft,IRight] is non-processed part of XY/Tags
1925            //
1926            // After cycle is done, we have Ileft=IRight. We deal with
1927            // this element separately.
1928            //
1929            // After this, [I1,ILeft) contains left part, and [ILeft,I2)
1930            // contains right part.
1931            //
1932            ileft = i1;
1933            iright = i2-1;
1934            while( ileft<iright )
1935            {
1936                if( (double)(kdt.xy[ileft,d])<=(double)(s) )
1937                {
1938                   
1939                    //
1940                    // XY[ILeft] is on its place.
1941                    // Advance ILeft.
1942                    //
1943                    ileft = ileft+1;
1944                }
1945                else
1946                {
1947                   
1948                    //
1949                    // XY[ILeft,..] must be at IRight.
1950                    // Swap and advance IRight.
1951                    //
1952                    for(i=0; i<=2*kdt.nx+kdt.ny-1; i++)
1953                    {
1954                        v = kdt.xy[ileft,i];
1955                        kdt.xy[ileft,i] = kdt.xy[iright,i];
1956                        kdt.xy[iright,i] = v;
1957                    }
1958                    j = kdt.tags[ileft];
1959                    kdt.tags[ileft] = kdt.tags[iright];
1960                    kdt.tags[iright] = j;
1961                    iright = iright-1;
1962                }
1963            }
1964            if( (double)(kdt.xy[ileft,d])<=(double)(s) )
1965            {
1966                ileft = ileft+1;
1967            }
1968            else
1969            {
1970                iright = iright-1;
1971            }
1972            i3 = ileft;
1973        }
1974
1975
1976        /*************************************************************************
1977        Recursive kd-tree generation subroutine.
1978
1979        PARAMETERS
1980            KDT         tree
1981            NodesOffs   unused part of Nodes[] which must be filled by tree
1982            SplitsOffs  unused part of Splits[]
1983            I1, I2      points from [I1,I2) are processed
1984           
1985        NodesOffs[] and SplitsOffs[] must be large enough.
1986
1987          -- ALGLIB --
1988             Copyright 28.02.2010 by Bochkanov Sergey
1989        *************************************************************************/
1990        private static void kdtreegeneratetreerec(kdtree kdt,
1991            ref int nodesoffs,
1992            ref int splitsoffs,
1993            int i1,
1994            int i2,
1995            int maxleafsize)
1996        {
1997            int n = 0;
1998            int nx = 0;
1999            int ny = 0;
2000            int i = 0;
2001            int j = 0;
2002            int oldoffs = 0;
2003            int i3 = 0;
2004            int cntless = 0;
2005            int cntgreater = 0;
2006            double minv = 0;
2007            double maxv = 0;
2008            int minidx = 0;
2009            int maxidx = 0;
2010            int d = 0;
2011            double ds = 0;
2012            double s = 0;
2013            double v = 0;
2014            int i_ = 0;
2015            int i1_ = 0;
2016
2017            ap.assert(i2>i1, "KDTreeGenerateTreeRec: internal error");
2018           
2019            //
2020            // Generate leaf if needed
2021            //
2022            if( i2-i1<=maxleafsize )
2023            {
2024                kdt.nodes[nodesoffs+0] = i2-i1;
2025                kdt.nodes[nodesoffs+1] = i1;
2026                nodesoffs = nodesoffs+2;
2027                return;
2028            }
2029           
2030            //
2031            // Load values for easier access
2032            //
2033            nx = kdt.nx;
2034            ny = kdt.ny;
2035           
2036            //
2037            // select dimension to split:
2038            // * D is a dimension number
2039            //
2040            d = 0;
2041            ds = kdt.curboxmax[0]-kdt.curboxmin[0];
2042            for(i=1; i<=nx-1; i++)
2043            {
2044                v = kdt.curboxmax[i]-kdt.curboxmin[i];
2045                if( (double)(v)>(double)(ds) )
2046                {
2047                    ds = v;
2048                    d = i;
2049                }
2050            }
2051           
2052            //
2053            // Select split position S using sliding midpoint rule,
2054            // rearrange points into [I1,I3) and [I3,I2)
2055            //
2056            s = kdt.curboxmin[d]+0.5*ds;
2057            i1_ = (i1) - (0);
2058            for(i_=0; i_<=i2-i1-1;i_++)
2059            {
2060                kdt.buf[i_] = kdt.xy[i_+i1_,d];
2061            }
2062            n = i2-i1;
2063            cntless = 0;
2064            cntgreater = 0;
2065            minv = kdt.buf[0];
2066            maxv = kdt.buf[0];
2067            minidx = i1;
2068            maxidx = i1;
2069            for(i=0; i<=n-1; i++)
2070            {
2071                v = kdt.buf[i];
2072                if( (double)(v)<(double)(minv) )
2073                {
2074                    minv = v;
2075                    minidx = i1+i;
2076                }
2077                if( (double)(v)>(double)(maxv) )
2078                {
2079                    maxv = v;
2080                    maxidx = i1+i;
2081                }
2082                if( (double)(v)<(double)(s) )
2083                {
2084                    cntless = cntless+1;
2085                }
2086                if( (double)(v)>(double)(s) )
2087                {
2088                    cntgreater = cntgreater+1;
2089                }
2090            }
2091            if( cntless>0 & cntgreater>0 )
2092            {
2093               
2094                //
2095                // normal midpoint split
2096                //
2097                kdtreesplit(kdt, i1, i2, d, s, ref i3);
2098            }
2099            else
2100            {
2101               
2102                //
2103                // sliding midpoint
2104                //
2105                if( cntless==0 )
2106                {
2107                   
2108                    //
2109                    // 1. move split to MinV,
2110                    // 2. place one point to the left bin (move to I1),
2111                    //    others - to the right bin
2112                    //
2113                    s = minv;
2114                    if( minidx!=i1 )
2115                    {
2116                        for(i=0; i<=2*kdt.nx+kdt.ny-1; i++)
2117                        {
2118                            v = kdt.xy[minidx,i];
2119                            kdt.xy[minidx,i] = kdt.xy[i1,i];
2120                            kdt.xy[i1,i] = v;
2121                        }
2122                        j = kdt.tags[minidx];
2123                        kdt.tags[minidx] = kdt.tags[i1];
2124                        kdt.tags[i1] = j;
2125                    }
2126                    i3 = i1+1;
2127                }
2128                else
2129                {
2130                   
2131                    //
2132                    // 1. move split to MaxV,
2133                    // 2. place one point to the right bin (move to I2-1),
2134                    //    others - to the left bin
2135                    //
2136                    s = maxv;
2137                    if( maxidx!=i2-1 )
2138                    {
2139                        for(i=0; i<=2*kdt.nx+kdt.ny-1; i++)
2140                        {
2141                            v = kdt.xy[maxidx,i];
2142                            kdt.xy[maxidx,i] = kdt.xy[i2-1,i];
2143                            kdt.xy[i2-1,i] = v;
2144                        }
2145                        j = kdt.tags[maxidx];
2146                        kdt.tags[maxidx] = kdt.tags[i2-1];
2147                        kdt.tags[i2-1] = j;
2148                    }
2149                    i3 = i2-1;
2150                }
2151            }
2152           
2153            //
2154            // Generate 'split' node
2155            //
2156            kdt.nodes[nodesoffs+0] = 0;
2157            kdt.nodes[nodesoffs+1] = d;
2158            kdt.nodes[nodesoffs+2] = splitsoffs;
2159            kdt.splits[splitsoffs+0] = s;
2160            oldoffs = nodesoffs;
2161            nodesoffs = nodesoffs+splitnodesize;
2162            splitsoffs = splitsoffs+1;
2163           
2164            //
2165            // Recirsive generation:
2166            // * update CurBox
2167            // * call subroutine
2168            // * restore CurBox
2169            //
2170            kdt.nodes[oldoffs+3] = nodesoffs;
2171            v = kdt.curboxmax[d];
2172            kdt.curboxmax[d] = s;
2173            kdtreegeneratetreerec(kdt, ref nodesoffs, ref splitsoffs, i1, i3, maxleafsize);
2174            kdt.curboxmax[d] = v;
2175            kdt.nodes[oldoffs+4] = nodesoffs;
2176            v = kdt.curboxmin[d];
2177            kdt.curboxmin[d] = s;
2178            kdtreegeneratetreerec(kdt, ref nodesoffs, ref splitsoffs, i3, i2, maxleafsize);
2179            kdt.curboxmin[d] = v;
2180        }
2181
2182
2183        /*************************************************************************
2184        Recursive subroutine for NN queries.
2185
2186          -- ALGLIB --
2187             Copyright 28.02.2010 by Bochkanov Sergey
2188        *************************************************************************/
2189        private static void kdtreequerynnrec(kdtree kdt,
2190            int offs)
2191        {
2192            double ptdist = 0;
2193            int i = 0;
2194            int j = 0;
2195            int nx = 0;
2196            int i1 = 0;
2197            int i2 = 0;
2198            int d = 0;
2199            double s = 0;
2200            double v = 0;
2201            double t1 = 0;
2202            int childbestoffs = 0;
2203            int childworstoffs = 0;
2204            int childoffs = 0;
2205            double prevdist = 0;
2206            bool todive = new bool();
2207            bool bestisleft = new bool();
2208            bool updatemin = new bool();
2209
2210           
2211            //
2212            // Leaf node.
2213            // Process points.
2214            //
2215            if( kdt.nodes[offs]>0 )
2216            {
2217                i1 = kdt.nodes[offs+1];
2218                i2 = i1+kdt.nodes[offs];
2219                for(i=i1; i<=i2-1; i++)
2220                {
2221                   
2222                    //
2223                    // Calculate distance
2224                    //
2225                    ptdist = 0;
2226                    nx = kdt.nx;
2227                    if( kdt.normtype==0 )
2228                    {
2229                        for(j=0; j<=nx-1; j++)
2230                        {
2231                            ptdist = Math.Max(ptdist, Math.Abs(kdt.xy[i,j]-kdt.x[j]));
2232                        }
2233                    }
2234                    if( kdt.normtype==1 )
2235                    {
2236                        for(j=0; j<=nx-1; j++)
2237                        {
2238                            ptdist = ptdist+Math.Abs(kdt.xy[i,j]-kdt.x[j]);
2239                        }
2240                    }
2241                    if( kdt.normtype==2 )
2242                    {
2243                        for(j=0; j<=nx-1; j++)
2244                        {
2245                            ptdist = ptdist+math.sqr(kdt.xy[i,j]-kdt.x[j]);
2246                        }
2247                    }
2248                   
2249                    //
2250                    // Skip points with zero distance if self-matches are turned off
2251                    //
2252                    if( (double)(ptdist)==(double)(0) & !kdt.selfmatch )
2253                    {
2254                        continue;
2255                    }
2256                   
2257                    //
2258                    // We CAN'T process point if R-criterion isn't satisfied,
2259                    // i.e. (RNeeded<>0) AND (PtDist>R).
2260                    //
2261                    if( (double)(kdt.rneeded)==(double)(0) | (double)(ptdist)<=(double)(kdt.rneeded) )
2262                    {
2263                       
2264                        //
2265                        // R-criterion is satisfied, we must either:
2266                        // * replace worst point, if (KNeeded<>0) AND (KCur=KNeeded)
2267                        //   (or skip, if worst point is better)
2268                        // * add point without replacement otherwise
2269                        //
2270                        if( kdt.kcur<kdt.kneeded | kdt.kneeded==0 )
2271                        {
2272                           
2273                            //
2274                            // add current point to heap without replacement
2275                            //
2276                            tsort.tagheappushi(ref kdt.r, ref kdt.idx, ref kdt.kcur, ptdist, i);
2277                        }
2278                        else
2279                        {
2280                           
2281                            //
2282                            // New points are added or not, depending on their distance.
2283                            // If added, they replace element at the top of the heap
2284                            //
2285                            if( (double)(ptdist)<(double)(kdt.r[0]) )
2286                            {
2287                                if( kdt.kneeded==1 )
2288                                {
2289                                    kdt.idx[0] = i;
2290                                    kdt.r[0] = ptdist;
2291                                }
2292                                else
2293                                {
2294                                    tsort.tagheapreplacetopi(ref kdt.r, ref kdt.idx, kdt.kneeded, ptdist, i);
2295                                }
2296                            }
2297                        }
2298                    }
2299                }
2300                return;
2301            }
2302           
2303            //
2304            // Simple split
2305            //
2306            if( kdt.nodes[offs]==0 )
2307            {
2308               
2309                //
2310                // Load:
2311                // * D  dimension to split
2312                // * S  split position
2313                //
2314                d = kdt.nodes[offs+1];
2315                s = kdt.splits[kdt.nodes[offs+2]];
2316               
2317                //
2318                // Calculate:
2319                // * ChildBestOffs      child box with best chances
2320                // * ChildWorstOffs     child box with worst chances
2321                //
2322                if( (double)(kdt.x[d])<=(double)(s) )
2323                {
2324                    childbestoffs = kdt.nodes[offs+3];
2325                    childworstoffs = kdt.nodes[offs+4];
2326                    bestisleft = true;
2327                }
2328                else
2329                {
2330                    childbestoffs = kdt.nodes[offs+4];
2331                    childworstoffs = kdt.nodes[offs+3];
2332                    bestisleft = false;
2333                }
2334               
2335                //
2336                // Navigate through childs
2337                //
2338                for(i=0; i<=1; i++)
2339                {
2340                   
2341                    //
2342                    // Select child to process:
2343                    // * ChildOffs      current child offset in Nodes[]
2344                    // * UpdateMin      whether minimum or maximum value
2345                    //                  of bounding box is changed on update
2346                    //
2347                    if( i==0 )
2348                    {
2349                        childoffs = childbestoffs;
2350                        updatemin = !bestisleft;
2351                    }
2352                    else
2353                    {
2354                        updatemin = bestisleft;
2355                        childoffs = childworstoffs;
2356                    }
2357                   
2358                    //
2359                    // Update bounding box and current distance
2360                    //
2361                    if( updatemin )
2362                    {
2363                        prevdist = kdt.curdist;
2364                        t1 = kdt.x[d];
2365                        v = kdt.curboxmin[d];
2366                        if( (double)(t1)<=(double)(s) )
2367                        {
2368                            if( kdt.normtype==0 )
2369                            {
2370                                kdt.curdist = Math.Max(kdt.curdist, s-t1);
2371                            }
2372                            if( kdt.normtype==1 )
2373                            {
2374                                kdt.curdist = kdt.curdist-Math.Max(v-t1, 0)+s-t1;
2375                            }
2376                            if( kdt.normtype==2 )
2377                            {
2378                                kdt.curdist = kdt.curdist-math.sqr(Math.Max(v-t1, 0))+math.sqr(s-t1);
2379                            }
2380                        }
2381                        kdt.curboxmin[d] = s;
2382                    }
2383                    else
2384                    {
2385                        prevdist = kdt.curdist;
2386                        t1 = kdt.x[d];
2387                        v = kdt.curboxmax[d];
2388                        if( (double)(t1)>=(double)(s) )
2389                        {
2390                            if( kdt.normtype==0 )
2391                            {
2392                                kdt.curdist = Math.Max(kdt.curdist, t1-s);
2393                            }
2394                            if( kdt.normtype==1 )
2395                            {
2396                                kdt.curdist = kdt.curdist-Math.Max(t1-v, 0)+t1-s;
2397                            }
2398                            if( kdt.normtype==2 )
2399                            {
2400                                kdt.curdist = kdt.curdist-math.sqr(Math.Max(t1-v, 0))+math.sqr(t1-s);
2401                            }
2402                        }
2403                        kdt.curboxmax[d] = s;
2404                    }
2405                   
2406                    //
2407                    // Decide: to dive into cell or not to dive
2408                    //
2409                    if( (double)(kdt.rneeded)!=(double)(0) & (double)(kdt.curdist)>(double)(kdt.rneeded) )
2410                    {
2411                        todive = false;
2412                    }
2413                    else
2414                    {
2415                        if( kdt.kcur<kdt.kneeded | kdt.kneeded==0 )
2416                        {
2417                           
2418                            //
2419                            // KCur<KNeeded (i.e. not all points are found)
2420                            //
2421                            todive = true;
2422                        }
2423                        else
2424                        {
2425                           
2426                            //
2427                            // KCur=KNeeded, decide to dive or not to dive
2428                            // using point position relative to bounding box.
2429                            //
2430                            todive = (double)(kdt.curdist)<=(double)(kdt.r[0]*kdt.approxf);
2431                        }
2432                    }
2433                    if( todive )
2434                    {
2435                        kdtreequerynnrec(kdt, childoffs);
2436                    }
2437                   
2438                    //
2439                    // Restore bounding box and distance
2440                    //
2441                    if( updatemin )
2442                    {
2443                        kdt.curboxmin[d] = v;
2444                    }
2445                    else
2446                    {
2447                        kdt.curboxmax[d] = v;
2448                    }
2449                    kdt.curdist = prevdist;
2450                }
2451                return;
2452            }
2453        }
2454
2455
2456        /*************************************************************************
2457        Copies X[] to KDT.X[]
2458        Loads distance from X[] to bounding box.
2459        Initializes CurBox[].
2460
2461          -- ALGLIB --
2462             Copyright 28.02.2010 by Bochkanov Sergey
2463        *************************************************************************/
2464        private static void kdtreeinitbox(kdtree kdt,
2465            double[] x)
2466        {
2467            int i = 0;
2468            double vx = 0;
2469            double vmin = 0;
2470            double vmax = 0;
2471
2472           
2473            //
2474            // calculate distance from point to current bounding box
2475            //
2476            kdt.curdist = 0;
2477            if( kdt.normtype==0 )
2478            {
2479                for(i=0; i<=kdt.nx-1; i++)
2480                {
2481                    vx = x[i];
2482                    vmin = kdt.boxmin[i];
2483                    vmax = kdt.boxmax[i];
2484                    kdt.x[i] = vx;
2485                    kdt.curboxmin[i] = vmin;
2486                    kdt.curboxmax[i] = vmax;
2487                    if( (double)(vx)<(double)(vmin) )
2488                    {
2489                        kdt.curdist = Math.Max(kdt.curdist, vmin-vx);
2490                    }
2491                    else
2492                    {
2493                        if( (double)(vx)>(double)(vmax) )
2494                        {
2495                            kdt.curdist = Math.Max(kdt.curdist, vx-vmax);
2496                        }
2497                    }
2498                }
2499            }
2500            if( kdt.normtype==1 )
2501            {
2502                for(i=0; i<=kdt.nx-1; i++)
2503                {
2504                    vx = x[i];
2505                    vmin = kdt.boxmin[i];
2506                    vmax = kdt.boxmax[i];
2507                    kdt.x[i] = vx;
2508                    kdt.curboxmin[i] = vmin;
2509                    kdt.curboxmax[i] = vmax;
2510                    if( (double)(vx)<(double)(vmin) )
2511                    {
2512                        kdt.curdist = kdt.curdist+vmin-vx;
2513                    }
2514                    else
2515                    {
2516                        if( (double)(vx)>(double)(vmax) )
2517                        {
2518                            kdt.curdist = kdt.curdist+vx-vmax;
2519                        }
2520                    }
2521                }
2522            }
2523            if( kdt.normtype==2 )
2524            {
2525                for(i=0; i<=kdt.nx-1; i++)
2526                {
2527                    vx = x[i];
2528                    vmin = kdt.boxmin[i];
2529                    vmax = kdt.boxmax[i];
2530                    kdt.x[i] = vx;
2531                    kdt.curboxmin[i] = vmin;
2532                    kdt.curboxmax[i] = vmax;
2533                    if( (double)(vx)<(double)(vmin) )
2534                    {
2535                        kdt.curdist = kdt.curdist+math.sqr(vmin-vx);
2536                    }
2537                    else
2538                    {
2539                        if( (double)(vx)>(double)(vmax) )
2540                        {
2541                            kdt.curdist = kdt.curdist+math.sqr(vx-vmax);
2542                        }
2543                    }
2544                }
2545            }
2546        }
2547
2548
2549        /*************************************************************************
2550        This function allocates all dataset-independent array  fields  of  KDTree,
2551        i.e.  such  array  fields  that  their dimensions do not depend on dataset
2552        size.
2553
2554        This function do not sets KDT.NX or KDT.NY - it just allocates arrays
2555
2556          -- ALGLIB --
2557             Copyright 14.03.2011 by Bochkanov Sergey
2558        *************************************************************************/
2559        private static void kdtreeallocdatasetindependent(kdtree kdt,
2560            int nx,
2561            int ny)
2562        {
2563            kdt.x = new double[nx];
2564            kdt.boxmin = new double[nx];
2565            kdt.boxmax = new double[nx];
2566            kdt.curboxmin = new double[nx];
2567            kdt.curboxmax = new double[nx];
2568        }
2569
2570
2571        /*************************************************************************
2572        This function allocates all dataset-dependent array fields of KDTree, i.e.
2573        such array fields that their dimensions depend on dataset size.
2574
2575        This function do not sets KDT.N, KDT.NX or KDT.NY -
2576        it just allocates arrays.
2577
2578          -- ALGLIB --
2579             Copyright 14.03.2011 by Bochkanov Sergey
2580        *************************************************************************/
2581        private static void kdtreeallocdatasetdependent(kdtree kdt,
2582            int n,
2583            int nx,
2584            int ny)
2585        {
2586            kdt.xy = new double[n, 2*nx+ny];
2587            kdt.tags = new int[n];
2588            kdt.idx = new int[n];
2589            kdt.r = new double[n];
2590            kdt.x = new double[nx];
2591            kdt.buf = new double[Math.Max(n, nx)];
2592            kdt.nodes = new int[splitnodesize*2*n];
2593            kdt.splits = new double[2*n];
2594        }
2595
2596
2597        /*************************************************************************
2598        This function allocates temporaries.
2599
2600        This function do not sets KDT.N, KDT.NX or KDT.NY -
2601        it just allocates arrays.
2602
2603          -- ALGLIB --
2604             Copyright 14.03.2011 by Bochkanov Sergey
2605        *************************************************************************/
2606        private static void kdtreealloctemporaries(kdtree kdt,
2607            int n,
2608            int nx,
2609            int ny)
2610        {
2611            kdt.x = new double[nx];
2612            kdt.idx = new int[n];
2613            kdt.r = new double[n];
2614            kdt.buf = new double[Math.Max(n, nx)];
2615            kdt.curboxmin = new double[nx];
2616            kdt.curboxmax = new double[nx];
2617        }
2618
2619
2620    }
2621}
2622
Note: See TracBrowser for help on using the repository browser.