Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/3.5.0/ALGLIB-3.5.0/alglibmisc.cs @ 7673

Last change on this file since 7673 was 7673, checked in by mkommend, 11 years ago

#1808: Added new alglib release 3.5.0 to the external libraries.

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