Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/3.6.0/ALGLIB-3.6.0/alglibmisc.cs @ 8463

Last change on this file since 8463 was 8418, checked in by abeham, 12 years ago

#1905: Added ALGLIB 3.6.0 to ExtLibs

File size: 101.9 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            int s0 = 0;
825            int s1 = 0;
826
827            s0 = math.randominteger(hqrndm1);
828            s1 = math.randominteger(hqrndm2);
829            hqrndseed(s0, s1, state);
830        }
831
832
833        /*************************************************************************
834        HQRNDState initialization with seed values
835
836          -- ALGLIB --
837             Copyright 02.12.2009 by Bochkanov Sergey
838        *************************************************************************/
839        public static void hqrndseed(int s1,
840            int s2,
841            hqrndstate state)
842        {
843            state.s1 = s1%(hqrndm1-1)+1;
844            state.s2 = s2%(hqrndm2-1)+1;
845            state.v = (double)1/(double)hqrndmax;
846            state.magicv = hqrndmagic;
847        }
848
849
850        /*************************************************************************
851        This function generates random real number in (0,1),
852        not including interval boundaries
853
854        State structure must be initialized with HQRNDRandomize() or HQRNDSeed().
855
856          -- ALGLIB --
857             Copyright 02.12.2009 by Bochkanov Sergey
858        *************************************************************************/
859        public static double hqrnduniformr(hqrndstate state)
860        {
861            double result = 0;
862
863            result = state.v*hqrndintegerbase(state);
864            return result;
865        }
866
867
868        /*************************************************************************
869        This function generates random integer number in [0, N)
870
871        1. N must be less than HQRNDMax-1.
872        2. State structure must be initialized with HQRNDRandomize() or HQRNDSeed()
873
874          -- ALGLIB --
875             Copyright 02.12.2009 by Bochkanov Sergey
876        *************************************************************************/
877        public static int hqrnduniformi(hqrndstate state,
878            int n)
879        {
880            int result = 0;
881            int mx = 0;
882
883           
884            //
885            // Correct handling of N's close to RNDBaseMax
886            // (avoiding skewed distributions for RNDBaseMax<>K*N)
887            //
888            alglib.ap.assert(n>0, "HQRNDUniformI: N<=0!");
889            alglib.ap.assert(n<hqrndmax-1, "HQRNDUniformI: N>=RNDBaseMax-1!");
890            mx = hqrndmax-1-(hqrndmax-1)%n;
891            do
892            {
893                result = hqrndintegerbase(state)-1;
894            }
895            while( result>=mx );
896            result = result%n;
897            return result;
898        }
899
900
901        /*************************************************************************
902        Random number generator: normal numbers
903
904        This function generates one random number from normal distribution.
905        Its performance is equal to that of HQRNDNormal2()
906
907        State structure must be initialized with HQRNDRandomize() or HQRNDSeed().
908
909          -- ALGLIB --
910             Copyright 02.12.2009 by Bochkanov Sergey
911        *************************************************************************/
912        public static double hqrndnormal(hqrndstate state)
913        {
914            double result = 0;
915            double v1 = 0;
916            double v2 = 0;
917
918            hqrndnormal2(state, ref v1, ref v2);
919            result = v1;
920            return result;
921        }
922
923
924        /*************************************************************************
925        Random number generator: random X and Y such that X^2+Y^2=1
926
927        State structure must be initialized with HQRNDRandomize() or HQRNDSeed().
928
929          -- ALGLIB --
930             Copyright 02.12.2009 by Bochkanov Sergey
931        *************************************************************************/
932        public static void hqrndunit2(hqrndstate state,
933            ref double x,
934            ref double y)
935        {
936            double v = 0;
937            double mx = 0;
938            double mn = 0;
939
940            x = 0;
941            y = 0;
942
943            do
944            {
945                hqrndnormal2(state, ref x, ref y);
946            }
947            while( !((double)(x)!=(double)(0) || (double)(y)!=(double)(0)) );
948            mx = Math.Max(Math.Abs(x), Math.Abs(y));
949            mn = Math.Min(Math.Abs(x), Math.Abs(y));
950            v = mx*Math.Sqrt(1+math.sqr(mn/mx));
951            x = x/v;
952            y = y/v;
953        }
954
955
956        /*************************************************************************
957        Random number generator: normal numbers
958
959        This function generates two independent random numbers from normal
960        distribution. Its performance is equal to that of HQRNDNormal()
961
962        State structure must be initialized with HQRNDRandomize() or HQRNDSeed().
963
964          -- ALGLIB --
965             Copyright 02.12.2009 by Bochkanov Sergey
966        *************************************************************************/
967        public static void hqrndnormal2(hqrndstate state,
968            ref double x1,
969            ref double x2)
970        {
971            double u = 0;
972            double v = 0;
973            double s = 0;
974
975            x1 = 0;
976            x2 = 0;
977
978            while( true )
979            {
980                u = 2*hqrnduniformr(state)-1;
981                v = 2*hqrnduniformr(state)-1;
982                s = math.sqr(u)+math.sqr(v);
983                if( (double)(s)>(double)(0) && (double)(s)<(double)(1) )
984                {
985                   
986                    //
987                    // two Sqrt's instead of one to
988                    // avoid overflow when S is too small
989                    //
990                    s = Math.Sqrt(-(2*Math.Log(s)))/Math.Sqrt(s);
991                    x1 = u*s;
992                    x2 = v*s;
993                    return;
994                }
995            }
996        }
997
998
999        /*************************************************************************
1000        Random number generator: exponential distribution
1001
1002        State structure must be initialized with HQRNDRandomize() or HQRNDSeed().
1003
1004          -- ALGLIB --
1005             Copyright 11.08.2007 by Bochkanov Sergey
1006        *************************************************************************/
1007        public static double hqrndexponential(hqrndstate state,
1008            double lambdav)
1009        {
1010            double result = 0;
1011
1012            alglib.ap.assert((double)(lambdav)>(double)(0), "HQRNDExponential: LambdaV<=0!");
1013            result = -(Math.Log(hqrnduniformr(state))/lambdav);
1014            return result;
1015        }
1016
1017
1018        /*************************************************************************
1019        This function generates  random number from discrete distribution given by
1020        finite sample X.
1021
1022        INPUT PARAMETERS
1023            State   -   high quality random number generator, must be
1024                        initialized with HQRNDRandomize() or HQRNDSeed().
1025                X   -   finite sample
1026                N   -   number of elements to use, N>=1
1027
1028        RESULT
1029            this function returns one of the X[i] for random i=0..N-1
1030
1031          -- ALGLIB --
1032             Copyright 08.11.2011 by Bochkanov Sergey
1033        *************************************************************************/
1034        public static double hqrnddiscrete(hqrndstate state,
1035            double[] x,
1036            int n)
1037        {
1038            double result = 0;
1039
1040            alglib.ap.assert(n>0, "HQRNDDiscrete: N<=0");
1041            alglib.ap.assert(n<=alglib.ap.len(x), "HQRNDDiscrete: Length(X)<N");
1042            result = x[hqrnduniformi(state, n)];
1043            return result;
1044        }
1045
1046
1047        /*************************************************************************
1048        This function generates random number from continuous  distribution  given
1049        by finite sample X.
1050
1051        INPUT PARAMETERS
1052            State   -   high quality random number generator, must be
1053                        initialized with HQRNDRandomize() or HQRNDSeed().
1054                X   -   finite sample, array[N] (can be larger, in this  case only
1055                        leading N elements are used). THIS ARRAY MUST BE SORTED BY
1056                        ASCENDING.
1057                N   -   number of elements to use, N>=1
1058
1059        RESULT
1060            this function returns random number from continuous distribution which 
1061            tries to approximate X as mush as possible. min(X)<=Result<=max(X).
1062
1063          -- ALGLIB --
1064             Copyright 08.11.2011 by Bochkanov Sergey
1065        *************************************************************************/
1066        public static double hqrndcontinuous(hqrndstate state,
1067            double[] x,
1068            int n)
1069        {
1070            double result = 0;
1071            double mx = 0;
1072            double mn = 0;
1073            int i = 0;
1074
1075            alglib.ap.assert(n>0, "HQRNDContinuous: N<=0");
1076            alglib.ap.assert(n<=alglib.ap.len(x), "HQRNDContinuous: Length(X)<N");
1077            if( n==1 )
1078            {
1079                result = x[0];
1080                return result;
1081            }
1082            i = hqrnduniformi(state, n-1);
1083            mn = x[i];
1084            mx = x[i+1];
1085            alglib.ap.assert((double)(mx)>=(double)(mn), "HQRNDDiscrete: X is not sorted by ascending");
1086            if( (double)(mx)!=(double)(mn) )
1087            {
1088                result = (mx-mn)*hqrnduniformr(state)+mn;
1089            }
1090            else
1091            {
1092                result = mn;
1093            }
1094            return result;
1095        }
1096
1097
1098        /*************************************************************************
1099
1100        L'Ecuyer, Efficient and portable combined random number generators
1101        *************************************************************************/
1102        private static int hqrndintegerbase(hqrndstate state)
1103        {
1104            int result = 0;
1105            int k = 0;
1106
1107            alglib.ap.assert(state.magicv==hqrndmagic, "HQRNDIntegerBase: State is not correctly initialized!");
1108            k = state.s1/53668;
1109            state.s1 = 40014*(state.s1-k*53668)-k*12211;
1110            if( state.s1<0 )
1111            {
1112                state.s1 = state.s1+2147483563;
1113            }
1114            k = state.s2/52774;
1115            state.s2 = 40692*(state.s2-k*52774)-k*3791;
1116            if( state.s2<0 )
1117            {
1118                state.s2 = state.s2+2147483399;
1119            }
1120           
1121            //
1122            // Result
1123            //
1124            result = state.s1-state.s2;
1125            if( result<1 )
1126            {
1127                result = result+2147483562;
1128            }
1129            return result;
1130        }
1131
1132
1133    }
1134    public class nearestneighbor
1135    {
1136        public class kdtree
1137        {
1138            public int n;
1139            public int nx;
1140            public int ny;
1141            public int normtype;
1142            public double[,] xy;
1143            public int[] tags;
1144            public double[] boxmin;
1145            public double[] boxmax;
1146            public int[] nodes;
1147            public double[] splits;
1148            public double[] x;
1149            public int kneeded;
1150            public double rneeded;
1151            public bool selfmatch;
1152            public double approxf;
1153            public int kcur;
1154            public int[] idx;
1155            public double[] r;
1156            public double[] buf;
1157            public double[] curboxmin;
1158            public double[] curboxmax;
1159            public double curdist;
1160            public int debugcounter;
1161            public kdtree()
1162            {
1163                xy = new double[0,0];
1164                tags = new int[0];
1165                boxmin = new double[0];
1166                boxmax = new double[0];
1167                nodes = new int[0];
1168                splits = new double[0];
1169                x = new double[0];
1170                idx = new int[0];
1171                r = new double[0];
1172                buf = new double[0];
1173                curboxmin = new double[0];
1174                curboxmax = new double[0];
1175            }
1176        };
1177
1178
1179
1180
1181        public const int splitnodesize = 6;
1182        public const int kdtreefirstversion = 0;
1183
1184
1185        /*************************************************************************
1186        KD-tree creation
1187
1188        This subroutine creates KD-tree from set of X-values and optional Y-values
1189
1190        INPUT PARAMETERS
1191            XY      -   dataset, array[0..N-1,0..NX+NY-1].
1192                        one row corresponds to one point.
1193                        first NX columns contain X-values, next NY (NY may be zero)
1194                        columns may contain associated Y-values
1195            N       -   number of points, N>=0.
1196            NX      -   space dimension, NX>=1.
1197            NY      -   number of optional Y-values, NY>=0.
1198            NormType-   norm type:
1199                        * 0 denotes infinity-norm
1200                        * 1 denotes 1-norm
1201                        * 2 denotes 2-norm (Euclidean norm)
1202                       
1203        OUTPUT PARAMETERS
1204            KDT     -   KD-tree
1205           
1206           
1207        NOTES
1208
1209        1. KD-tree  creation  have O(N*logN) complexity and O(N*(2*NX+NY))  memory
1210           requirements.
1211        2. Although KD-trees may be used with any combination of N  and  NX,  they
1212           are more efficient than brute-force search only when N >> 4^NX. So they
1213           are most useful in low-dimensional tasks (NX=2, NX=3). NX=1  is another
1214           inefficient case, because  simple  binary  search  (without  additional
1215           structures) is much more efficient in such tasks than KD-trees.
1216
1217          -- ALGLIB --
1218             Copyright 28.02.2010 by Bochkanov Sergey
1219        *************************************************************************/
1220        public static void kdtreebuild(double[,] xy,
1221            int n,
1222            int nx,
1223            int ny,
1224            int normtype,
1225            kdtree kdt)
1226        {
1227            int[] tags = new int[0];
1228            int i = 0;
1229
1230            alglib.ap.assert(n>=0, "KDTreeBuild: N<0");
1231            alglib.ap.assert(nx>=1, "KDTreeBuild: NX<1");
1232            alglib.ap.assert(ny>=0, "KDTreeBuild: NY<0");
1233            alglib.ap.assert(normtype>=0 && normtype<=2, "KDTreeBuild: incorrect NormType");
1234            alglib.ap.assert(alglib.ap.rows(xy)>=n, "KDTreeBuild: rows(X)<N");
1235            alglib.ap.assert(alglib.ap.cols(xy)>=nx+ny || n==0, "KDTreeBuild: cols(X)<NX+NY");
1236            alglib.ap.assert(apserv.apservisfinitematrix(xy, n, nx+ny), "KDTreeBuild: XY contains infinite or NaN values");
1237            if( n>0 )
1238            {
1239                tags = new int[n];
1240                for(i=0; i<=n-1; i++)
1241                {
1242                    tags[i] = 0;
1243                }
1244            }
1245            kdtreebuildtagged(xy, tags, n, nx, ny, normtype, kdt);
1246        }
1247
1248
1249        /*************************************************************************
1250        KD-tree creation
1251
1252        This  subroutine  creates  KD-tree  from set of X-values, integer tags and
1253        optional Y-values
1254
1255        INPUT PARAMETERS
1256            XY      -   dataset, array[0..N-1,0..NX+NY-1].
1257                        one row corresponds to one point.
1258                        first NX columns contain X-values, next NY (NY may be zero)
1259                        columns may contain associated Y-values
1260            Tags    -   tags, array[0..N-1], contains integer tags associated
1261                        with points.
1262            N       -   number of points, N>=0
1263            NX      -   space dimension, NX>=1.
1264            NY      -   number of optional Y-values, NY>=0.
1265            NormType-   norm type:
1266                        * 0 denotes infinity-norm
1267                        * 1 denotes 1-norm
1268                        * 2 denotes 2-norm (Euclidean norm)
1269
1270        OUTPUT PARAMETERS
1271            KDT     -   KD-tree
1272
1273        NOTES
1274
1275        1. KD-tree  creation  have O(N*logN) complexity and O(N*(2*NX+NY))  memory
1276           requirements.
1277        2. Although KD-trees may be used with any combination of N  and  NX,  they
1278           are more efficient than brute-force search only when N >> 4^NX. So they
1279           are most useful in low-dimensional tasks (NX=2, NX=3). NX=1  is another
1280           inefficient case, because  simple  binary  search  (without  additional
1281           structures) is much more efficient in such tasks than KD-trees.
1282
1283          -- ALGLIB --
1284             Copyright 28.02.2010 by Bochkanov Sergey
1285        *************************************************************************/
1286        public static void kdtreebuildtagged(double[,] xy,
1287            int[] tags,
1288            int n,
1289            int nx,
1290            int ny,
1291            int normtype,
1292            kdtree kdt)
1293        {
1294            int i = 0;
1295            int j = 0;
1296            int maxnodes = 0;
1297            int nodesoffs = 0;
1298            int splitsoffs = 0;
1299            int i_ = 0;
1300            int i1_ = 0;
1301
1302            alglib.ap.assert(n>=0, "KDTreeBuildTagged: N<0");
1303            alglib.ap.assert(nx>=1, "KDTreeBuildTagged: NX<1");
1304            alglib.ap.assert(ny>=0, "KDTreeBuildTagged: NY<0");
1305            alglib.ap.assert(normtype>=0 && normtype<=2, "KDTreeBuildTagged: incorrect NormType");
1306            alglib.ap.assert(alglib.ap.rows(xy)>=n, "KDTreeBuildTagged: rows(X)<N");
1307            alglib.ap.assert(alglib.ap.cols(xy)>=nx+ny || n==0, "KDTreeBuildTagged: cols(X)<NX+NY");
1308            alglib.ap.assert(apserv.apservisfinitematrix(xy, n, nx+ny), "KDTreeBuildTagged: XY contains infinite or NaN values");
1309           
1310            //
1311            // initialize
1312            //
1313            kdt.n = n;
1314            kdt.nx = nx;
1315            kdt.ny = ny;
1316            kdt.normtype = normtype;
1317            kdt.kcur = 0;
1318           
1319            //
1320            // N=0 => quick exit
1321            //
1322            if( n==0 )
1323            {
1324                return;
1325            }
1326           
1327            //
1328            // Allocate
1329            //
1330            kdtreeallocdatasetindependent(kdt, nx, ny);
1331            kdtreeallocdatasetdependent(kdt, n, nx, ny);
1332           
1333            //
1334            // Initial fill
1335            //
1336            for(i=0; i<=n-1; i++)
1337            {
1338                for(i_=0; i_<=nx-1;i_++)
1339                {
1340                    kdt.xy[i,i_] = xy[i,i_];
1341                }
1342                i1_ = (0) - (nx);
1343                for(i_=nx; i_<=2*nx+ny-1;i_++)
1344                {
1345                    kdt.xy[i,i_] = xy[i,i_+i1_];
1346                }
1347                kdt.tags[i] = tags[i];
1348            }
1349           
1350            //
1351            // Determine bounding box
1352            //
1353            for(i_=0; i_<=nx-1;i_++)
1354            {
1355                kdt.boxmin[i_] = kdt.xy[0,i_];
1356            }
1357            for(i_=0; i_<=nx-1;i_++)
1358            {
1359                kdt.boxmax[i_] = kdt.xy[0,i_];
1360            }
1361            for(i=1; i<=n-1; i++)
1362            {
1363                for(j=0; j<=nx-1; j++)
1364                {
1365                    kdt.boxmin[j] = Math.Min(kdt.boxmin[j], kdt.xy[i,j]);
1366                    kdt.boxmax[j] = Math.Max(kdt.boxmax[j], kdt.xy[i,j]);
1367                }
1368            }
1369           
1370            //
1371            // prepare tree structure
1372            // * MaxNodes=N because we guarantee no trivial splits, i.e.
1373            //   every split will generate two non-empty boxes
1374            //
1375            maxnodes = n;
1376            kdt.nodes = new int[splitnodesize*2*maxnodes];
1377            kdt.splits = new double[2*maxnodes];
1378            nodesoffs = 0;
1379            splitsoffs = 0;
1380            for(i_=0; i_<=nx-1;i_++)
1381            {
1382                kdt.curboxmin[i_] = kdt.boxmin[i_];
1383            }
1384            for(i_=0; i_<=nx-1;i_++)
1385            {
1386                kdt.curboxmax[i_] = kdt.boxmax[i_];
1387            }
1388            kdtreegeneratetreerec(kdt, ref nodesoffs, ref splitsoffs, 0, n, 8);
1389        }
1390
1391
1392        /*************************************************************************
1393        K-NN query: K nearest neighbors
1394
1395        INPUT PARAMETERS
1396            KDT         -   KD-tree
1397            X           -   point, array[0..NX-1].
1398            K           -   number of neighbors to return, K>=1
1399            SelfMatch   -   whether self-matches are allowed:
1400                            * if True, nearest neighbor may be the point itself
1401                              (if it exists in original dataset)
1402                            * if False, then only points with non-zero distance
1403                              are returned
1404                            * if not given, considered True
1405
1406        RESULT
1407            number of actual neighbors found (either K or N, if K>N).
1408
1409        This  subroutine  performs  query  and  stores  its result in the internal
1410        structures of the KD-tree. You can use  following  subroutines  to  obtain
1411        these results:
1412        * KDTreeQueryResultsX() to get X-values
1413        * KDTreeQueryResultsXY() to get X- and Y-values
1414        * KDTreeQueryResultsTags() to get tag values
1415        * KDTreeQueryResultsDistances() to get distances
1416
1417          -- ALGLIB --
1418             Copyright 28.02.2010 by Bochkanov Sergey
1419        *************************************************************************/
1420        public static int kdtreequeryknn(kdtree kdt,
1421            double[] x,
1422            int k,
1423            bool selfmatch)
1424        {
1425            int result = 0;
1426
1427            alglib.ap.assert(k>=1, "KDTreeQueryKNN: K<1!");
1428            alglib.ap.assert(alglib.ap.len(x)>=kdt.nx, "KDTreeQueryKNN: Length(X)<NX!");
1429            alglib.ap.assert(apserv.isfinitevector(x, kdt.nx), "KDTreeQueryKNN: X contains infinite or NaN values!");
1430            result = kdtreequeryaknn(kdt, x, k, selfmatch, 0.0);
1431            return result;
1432        }
1433
1434
1435        /*************************************************************************
1436        R-NN query: all points within R-sphere centered at X
1437
1438        INPUT PARAMETERS
1439            KDT         -   KD-tree
1440            X           -   point, array[0..NX-1].
1441            R           -   radius of sphere (in corresponding norm), R>0
1442            SelfMatch   -   whether self-matches are allowed:
1443                            * if True, nearest neighbor may be the point itself
1444                              (if it exists in original dataset)
1445                            * if False, then only points with non-zero distance
1446                              are returned
1447                            * if not given, considered True
1448
1449        RESULT
1450            number of neighbors found, >=0
1451
1452        This  subroutine  performs  query  and  stores  its result in the internal
1453        structures of the KD-tree. You can use  following  subroutines  to  obtain
1454        actual results:
1455        * KDTreeQueryResultsX() to get X-values
1456        * KDTreeQueryResultsXY() to get X- and Y-values
1457        * KDTreeQueryResultsTags() to get tag values
1458        * KDTreeQueryResultsDistances() to get distances
1459
1460          -- ALGLIB --
1461             Copyright 28.02.2010 by Bochkanov Sergey
1462        *************************************************************************/
1463        public static int kdtreequeryrnn(kdtree kdt,
1464            double[] x,
1465            double r,
1466            bool selfmatch)
1467        {
1468            int result = 0;
1469            int i = 0;
1470            int j = 0;
1471
1472            alglib.ap.assert((double)(r)>(double)(0), "KDTreeQueryRNN: incorrect R!");
1473            alglib.ap.assert(alglib.ap.len(x)>=kdt.nx, "KDTreeQueryRNN: Length(X)<NX!");
1474            alglib.ap.assert(apserv.isfinitevector(x, kdt.nx), "KDTreeQueryRNN: X contains infinite or NaN values!");
1475           
1476            //
1477            // Handle special case: KDT.N=0
1478            //
1479            if( kdt.n==0 )
1480            {
1481                kdt.kcur = 0;
1482                result = 0;
1483                return result;
1484            }
1485           
1486            //
1487            // Prepare parameters
1488            //
1489            kdt.kneeded = 0;
1490            if( kdt.normtype!=2 )
1491            {
1492                kdt.rneeded = r;
1493            }
1494            else
1495            {
1496                kdt.rneeded = math.sqr(r);
1497            }
1498            kdt.selfmatch = selfmatch;
1499            kdt.approxf = 1;
1500            kdt.kcur = 0;
1501           
1502            //
1503            // calculate distance from point to current bounding box
1504            //
1505            kdtreeinitbox(kdt, x);
1506           
1507            //
1508            // call recursive search
1509            // results are returned as heap
1510            //
1511            kdtreequerynnrec(kdt, 0);
1512           
1513            //
1514            // pop from heap to generate ordered representation
1515            //
1516            // last element is not pop'ed because it is already in
1517            // its place
1518            //
1519            result = kdt.kcur;
1520            j = kdt.kcur;
1521            for(i=kdt.kcur; i>=2; i--)
1522            {
1523                tsort.tagheappopi(ref kdt.r, ref kdt.idx, ref j);
1524            }
1525            return result;
1526        }
1527
1528
1529        /*************************************************************************
1530        K-NN query: approximate K nearest neighbors
1531
1532        INPUT PARAMETERS
1533            KDT         -   KD-tree
1534            X           -   point, array[0..NX-1].
1535            K           -   number of neighbors to return, K>=1
1536            SelfMatch   -   whether self-matches are allowed:
1537                            * if True, nearest neighbor may be the point itself
1538                              (if it exists in original dataset)
1539                            * if False, then only points with non-zero distance
1540                              are returned
1541                            * if not given, considered True
1542            Eps         -   approximation factor, Eps>=0. eps-approximate  nearest
1543                            neighbor  is  a  neighbor  whose distance from X is at
1544                            most (1+eps) times distance of true nearest neighbor.
1545
1546        RESULT
1547            number of actual neighbors found (either K or N, if K>N).
1548           
1549        NOTES
1550            significant performance gain may be achieved only when Eps  is  is  on
1551            the order of magnitude of 1 or larger.
1552
1553        This  subroutine  performs  query  and  stores  its result in the internal
1554        structures of the KD-tree. You can use  following  subroutines  to  obtain
1555        these results:
1556        * KDTreeQueryResultsX() to get X-values
1557        * KDTreeQueryResultsXY() to get X- and Y-values
1558        * KDTreeQueryResultsTags() to get tag values
1559        * KDTreeQueryResultsDistances() to get distances
1560
1561          -- ALGLIB --
1562             Copyright 28.02.2010 by Bochkanov Sergey
1563        *************************************************************************/
1564        public static int kdtreequeryaknn(kdtree kdt,
1565            double[] x,
1566            int k,
1567            bool selfmatch,
1568            double eps)
1569        {
1570            int result = 0;
1571            int i = 0;
1572            int j = 0;
1573
1574            alglib.ap.assert(k>0, "KDTreeQueryAKNN: incorrect K!");
1575            alglib.ap.assert((double)(eps)>=(double)(0), "KDTreeQueryAKNN: incorrect Eps!");
1576            alglib.ap.assert(alglib.ap.len(x)>=kdt.nx, "KDTreeQueryAKNN: Length(X)<NX!");
1577            alglib.ap.assert(apserv.isfinitevector(x, kdt.nx), "KDTreeQueryAKNN: X contains infinite or NaN values!");
1578           
1579            //
1580            // Handle special case: KDT.N=0
1581            //
1582            if( kdt.n==0 )
1583            {
1584                kdt.kcur = 0;
1585                result = 0;
1586                return result;
1587            }
1588           
1589            //
1590            // Prepare parameters
1591            //
1592            k = Math.Min(k, kdt.n);
1593            kdt.kneeded = k;
1594            kdt.rneeded = 0;
1595            kdt.selfmatch = selfmatch;
1596            if( kdt.normtype==2 )
1597            {
1598                kdt.approxf = 1/math.sqr(1+eps);
1599            }
1600            else
1601            {
1602                kdt.approxf = 1/(1+eps);
1603            }
1604            kdt.kcur = 0;
1605           
1606            //
1607            // calculate distance from point to current bounding box
1608            //
1609            kdtreeinitbox(kdt, x);
1610           
1611            //
1612            // call recursive search
1613            // results are returned as heap
1614            //
1615            kdtreequerynnrec(kdt, 0);
1616           
1617            //
1618            // pop from heap to generate ordered representation
1619            //
1620            // last element is non pop'ed because it is already in
1621            // its place
1622            //
1623            result = kdt.kcur;
1624            j = kdt.kcur;
1625            for(i=kdt.kcur; i>=2; i--)
1626            {
1627                tsort.tagheappopi(ref kdt.r, ref kdt.idx, ref j);
1628            }
1629            return result;
1630        }
1631
1632
1633        /*************************************************************************
1634        X-values from last query
1635
1636        INPUT PARAMETERS
1637            KDT     -   KD-tree
1638            X       -   possibly pre-allocated buffer. If X is too small to store
1639                        result, it is resized. If size(X) is enough to store
1640                        result, it is left unchanged.
1641
1642        OUTPUT PARAMETERS
1643            X       -   rows are filled with X-values
1644
1645        NOTES
1646        1. points are ordered by distance from the query point (first = closest)
1647        2. if  XY is larger than required to store result, only leading part  will
1648           be overwritten; trailing part will be left unchanged. So  if  on  input
1649           XY = [[A,B],[C,D]], and result is [1,2],  then  on  exit  we  will  get
1650           XY = [[1,2],[C,D]]. This is done purposely to increase performance;  if
1651           you want function  to  resize  array  according  to  result  size,  use
1652           function with same name and suffix 'I'.
1653
1654        SEE ALSO
1655        * KDTreeQueryResultsXY()            X- and Y-values
1656        * KDTreeQueryResultsTags()          tag values
1657        * KDTreeQueryResultsDistances()     distances
1658
1659          -- ALGLIB --
1660             Copyright 28.02.2010 by Bochkanov Sergey
1661        *************************************************************************/
1662        public static void kdtreequeryresultsx(kdtree kdt,
1663            ref double[,] x)
1664        {
1665            int i = 0;
1666            int k = 0;
1667            int i_ = 0;
1668            int i1_ = 0;
1669
1670            if( kdt.kcur==0 )
1671            {
1672                return;
1673            }
1674            if( alglib.ap.rows(x)<kdt.kcur || alglib.ap.cols(x)<kdt.nx )
1675            {
1676                x = new double[kdt.kcur, kdt.nx];
1677            }
1678            k = kdt.kcur;
1679            for(i=0; i<=k-1; i++)
1680            {
1681                i1_ = (kdt.nx) - (0);
1682                for(i_=0; i_<=kdt.nx-1;i_++)
1683                {
1684                    x[i,i_] = kdt.xy[kdt.idx[i],i_+i1_];
1685                }
1686            }
1687        }
1688
1689
1690        /*************************************************************************
1691        X- and Y-values from last query
1692
1693        INPUT PARAMETERS
1694            KDT     -   KD-tree
1695            XY      -   possibly pre-allocated buffer. If XY is too small to store
1696                        result, it is resized. If size(XY) is enough to store
1697                        result, it is left unchanged.
1698
1699        OUTPUT PARAMETERS
1700            XY      -   rows are filled with points: first NX columns with
1701                        X-values, next NY columns - with Y-values.
1702
1703        NOTES
1704        1. points are ordered by distance from the query point (first = closest)
1705        2. if  XY is larger than required to store result, only leading part  will
1706           be overwritten; trailing part will be left unchanged. So  if  on  input
1707           XY = [[A,B],[C,D]], and result is [1,2],  then  on  exit  we  will  get
1708           XY = [[1,2],[C,D]]. This is done purposely to increase performance;  if
1709           you want function  to  resize  array  according  to  result  size,  use
1710           function with same name and suffix 'I'.
1711
1712        SEE ALSO
1713        * KDTreeQueryResultsX()             X-values
1714        * KDTreeQueryResultsTags()          tag values
1715        * KDTreeQueryResultsDistances()     distances
1716
1717          -- ALGLIB --
1718             Copyright 28.02.2010 by Bochkanov Sergey
1719        *************************************************************************/
1720        public static void kdtreequeryresultsxy(kdtree kdt,
1721            ref double[,] xy)
1722        {
1723            int i = 0;
1724            int k = 0;
1725            int i_ = 0;
1726            int i1_ = 0;
1727
1728            if( kdt.kcur==0 )
1729            {
1730                return;
1731            }
1732            if( alglib.ap.rows(xy)<kdt.kcur || alglib.ap.cols(xy)<kdt.nx+kdt.ny )
1733            {
1734                xy = new double[kdt.kcur, kdt.nx+kdt.ny];
1735            }
1736            k = kdt.kcur;
1737            for(i=0; i<=k-1; i++)
1738            {
1739                i1_ = (kdt.nx) - (0);
1740                for(i_=0; i_<=kdt.nx+kdt.ny-1;i_++)
1741                {
1742                    xy[i,i_] = kdt.xy[kdt.idx[i],i_+i1_];
1743                }
1744            }
1745        }
1746
1747
1748        /*************************************************************************
1749        Tags from last query
1750
1751        INPUT PARAMETERS
1752            KDT     -   KD-tree
1753            Tags    -   possibly pre-allocated buffer. If X is too small to store
1754                        result, it is resized. If size(X) is enough to store
1755                        result, it is left unchanged.
1756
1757        OUTPUT PARAMETERS
1758            Tags    -   filled with tags associated with points,
1759                        or, when no tags were supplied, with zeros
1760
1761        NOTES
1762        1. points are ordered by distance from the query point (first = closest)
1763        2. if  XY is larger than required to store result, only leading part  will
1764           be overwritten; trailing part will be left unchanged. So  if  on  input
1765           XY = [[A,B],[C,D]], and result is [1,2],  then  on  exit  we  will  get
1766           XY = [[1,2],[C,D]]. This is done purposely to increase performance;  if
1767           you want function  to  resize  array  according  to  result  size,  use
1768           function with same name and suffix 'I'.
1769
1770        SEE ALSO
1771        * KDTreeQueryResultsX()             X-values
1772        * KDTreeQueryResultsXY()            X- and Y-values
1773        * KDTreeQueryResultsDistances()     distances
1774
1775          -- ALGLIB --
1776             Copyright 28.02.2010 by Bochkanov Sergey
1777        *************************************************************************/
1778        public static void kdtreequeryresultstags(kdtree kdt,
1779            ref int[] tags)
1780        {
1781            int i = 0;
1782            int k = 0;
1783
1784            if( kdt.kcur==0 )
1785            {
1786                return;
1787            }
1788            if( alglib.ap.len(tags)<kdt.kcur )
1789            {
1790                tags = new int[kdt.kcur];
1791            }
1792            k = kdt.kcur;
1793            for(i=0; i<=k-1; i++)
1794            {
1795                tags[i] = kdt.tags[kdt.idx[i]];
1796            }
1797        }
1798
1799
1800        /*************************************************************************
1801        Distances from last query
1802
1803        INPUT PARAMETERS
1804            KDT     -   KD-tree
1805            R       -   possibly pre-allocated buffer. If X is too small to store
1806                        result, it is resized. If size(X) is enough to store
1807                        result, it is left unchanged.
1808
1809        OUTPUT PARAMETERS
1810            R       -   filled with distances (in corresponding norm)
1811
1812        NOTES
1813        1. points are ordered by distance from the query point (first = closest)
1814        2. if  XY is larger than required to store result, only leading part  will
1815           be overwritten; trailing part will be left unchanged. So  if  on  input
1816           XY = [[A,B],[C,D]], and result is [1,2],  then  on  exit  we  will  get
1817           XY = [[1,2],[C,D]]. This is done purposely to increase performance;  if
1818           you want function  to  resize  array  according  to  result  size,  use
1819           function with same name and suffix 'I'.
1820
1821        SEE ALSO
1822        * KDTreeQueryResultsX()             X-values
1823        * KDTreeQueryResultsXY()            X- and Y-values
1824        * KDTreeQueryResultsTags()          tag values
1825
1826          -- ALGLIB --
1827             Copyright 28.02.2010 by Bochkanov Sergey
1828        *************************************************************************/
1829        public static void kdtreequeryresultsdistances(kdtree kdt,
1830            ref double[] r)
1831        {
1832            int i = 0;
1833            int k = 0;
1834
1835            if( kdt.kcur==0 )
1836            {
1837                return;
1838            }
1839            if( alglib.ap.len(r)<kdt.kcur )
1840            {
1841                r = new double[kdt.kcur];
1842            }
1843            k = kdt.kcur;
1844           
1845            //
1846            // unload norms
1847            //
1848            // Abs() call is used to handle cases with negative norms
1849            // (generated during KFN requests)
1850            //
1851            if( kdt.normtype==0 )
1852            {
1853                for(i=0; i<=k-1; i++)
1854                {
1855                    r[i] = Math.Abs(kdt.r[i]);
1856                }
1857            }
1858            if( kdt.normtype==1 )
1859            {
1860                for(i=0; i<=k-1; i++)
1861                {
1862                    r[i] = Math.Abs(kdt.r[i]);
1863                }
1864            }
1865            if( kdt.normtype==2 )
1866            {
1867                for(i=0; i<=k-1; i++)
1868                {
1869                    r[i] = Math.Sqrt(Math.Abs(kdt.r[i]));
1870                }
1871            }
1872        }
1873
1874
1875        /*************************************************************************
1876        X-values from last query; 'interactive' variant for languages like  Python
1877        which   support    constructs   like  "X = KDTreeQueryResultsXI(KDT)"  and
1878        interactive mode of interpreter.
1879
1880        This function allocates new array on each call,  so  it  is  significantly
1881        slower than its 'non-interactive' counterpart, but it is  more  convenient
1882        when you call it from command line.
1883
1884          -- ALGLIB --
1885             Copyright 28.02.2010 by Bochkanov Sergey
1886        *************************************************************************/
1887        public static void kdtreequeryresultsxi(kdtree kdt,
1888            ref double[,] x)
1889        {
1890            x = new double[0,0];
1891
1892            kdtreequeryresultsx(kdt, ref x);
1893        }
1894
1895
1896        /*************************************************************************
1897        XY-values from last query; 'interactive' variant for languages like Python
1898        which   support    constructs   like "XY = KDTreeQueryResultsXYI(KDT)" and
1899        interactive mode of interpreter.
1900
1901        This function allocates new array on each call,  so  it  is  significantly
1902        slower than its 'non-interactive' counterpart, but it is  more  convenient
1903        when you call it from command line.
1904
1905          -- ALGLIB --
1906             Copyright 28.02.2010 by Bochkanov Sergey
1907        *************************************************************************/
1908        public static void kdtreequeryresultsxyi(kdtree kdt,
1909            ref double[,] xy)
1910        {
1911            xy = new double[0,0];
1912
1913            kdtreequeryresultsxy(kdt, ref xy);
1914        }
1915
1916
1917        /*************************************************************************
1918        Tags  from  last  query;  'interactive' variant for languages like  Python
1919        which  support  constructs  like "Tags = KDTreeQueryResultsTagsI(KDT)" and
1920        interactive mode of interpreter.
1921
1922        This function allocates new array on each call,  so  it  is  significantly
1923        slower than its 'non-interactive' counterpart, but it is  more  convenient
1924        when you call it from command line.
1925
1926          -- ALGLIB --
1927             Copyright 28.02.2010 by Bochkanov Sergey
1928        *************************************************************************/
1929        public static void kdtreequeryresultstagsi(kdtree kdt,
1930            ref int[] tags)
1931        {
1932            tags = new int[0];
1933
1934            kdtreequeryresultstags(kdt, ref tags);
1935        }
1936
1937
1938        /*************************************************************************
1939        Distances from last query; 'interactive' variant for languages like Python
1940        which  support  constructs   like  "R = KDTreeQueryResultsDistancesI(KDT)"
1941        and interactive mode of interpreter.
1942
1943        This function allocates new array on each call,  so  it  is  significantly
1944        slower than its 'non-interactive' counterpart, but it is  more  convenient
1945        when you call it from command line.
1946
1947          -- ALGLIB --
1948             Copyright 28.02.2010 by Bochkanov Sergey
1949        *************************************************************************/
1950        public static void kdtreequeryresultsdistancesi(kdtree kdt,
1951            ref double[] r)
1952        {
1953            r = new double[0];
1954
1955            kdtreequeryresultsdistances(kdt, ref r);
1956        }
1957
1958
1959        /*************************************************************************
1960        Serializer: allocation
1961
1962          -- ALGLIB --
1963             Copyright 14.03.2011 by Bochkanov Sergey
1964        *************************************************************************/
1965        public static void kdtreealloc(alglib.serializer s,
1966            kdtree tree)
1967        {
1968           
1969            //
1970            // Header
1971            //
1972            s.alloc_entry();
1973            s.alloc_entry();
1974           
1975            //
1976            // Data
1977            //
1978            s.alloc_entry();
1979            s.alloc_entry();
1980            s.alloc_entry();
1981            s.alloc_entry();
1982            apserv.allocrealmatrix(s, tree.xy, -1, -1);
1983            apserv.allocintegerarray(s, tree.tags, -1);
1984            apserv.allocrealarray(s, tree.boxmin, -1);
1985            apserv.allocrealarray(s, tree.boxmax, -1);
1986            apserv.allocintegerarray(s, tree.nodes, -1);
1987            apserv.allocrealarray(s, tree.splits, -1);
1988        }
1989
1990
1991        /*************************************************************************
1992        Serializer: serialization
1993
1994          -- ALGLIB --
1995             Copyright 14.03.2011 by Bochkanov Sergey
1996        *************************************************************************/
1997        public static void kdtreeserialize(alglib.serializer s,
1998            kdtree tree)
1999        {
2000           
2001            //
2002            // Header
2003            //
2004            s.serialize_int(scodes.getkdtreeserializationcode());
2005            s.serialize_int(kdtreefirstversion);
2006           
2007            //
2008            // Data
2009            //
2010            s.serialize_int(tree.n);
2011            s.serialize_int(tree.nx);
2012            s.serialize_int(tree.ny);
2013            s.serialize_int(tree.normtype);
2014            apserv.serializerealmatrix(s, tree.xy, -1, -1);
2015            apserv.serializeintegerarray(s, tree.tags, -1);
2016            apserv.serializerealarray(s, tree.boxmin, -1);
2017            apserv.serializerealarray(s, tree.boxmax, -1);
2018            apserv.serializeintegerarray(s, tree.nodes, -1);
2019            apserv.serializerealarray(s, tree.splits, -1);
2020        }
2021
2022
2023        /*************************************************************************
2024        Serializer: unserialization
2025
2026          -- ALGLIB --
2027             Copyright 14.03.2011 by Bochkanov Sergey
2028        *************************************************************************/
2029        public static void kdtreeunserialize(alglib.serializer s,
2030            kdtree tree)
2031        {
2032            int i0 = 0;
2033            int i1 = 0;
2034
2035           
2036            //
2037            // check correctness of header
2038            //
2039            i0 = s.unserialize_int();
2040            alglib.ap.assert(i0==scodes.getkdtreeserializationcode(), "KDTreeUnserialize: stream header corrupted");
2041            i1 = s.unserialize_int();
2042            alglib.ap.assert(i1==kdtreefirstversion, "KDTreeUnserialize: stream header corrupted");
2043           
2044            //
2045            // Unserialize data
2046            //
2047            tree.n = s.unserialize_int();
2048            tree.nx = s.unserialize_int();
2049            tree.ny = s.unserialize_int();
2050            tree.normtype = s.unserialize_int();
2051            apserv.unserializerealmatrix(s, ref tree.xy);
2052            apserv.unserializeintegerarray(s, ref tree.tags);
2053            apserv.unserializerealarray(s, ref tree.boxmin);
2054            apserv.unserializerealarray(s, ref tree.boxmax);
2055            apserv.unserializeintegerarray(s, ref tree.nodes);
2056            apserv.unserializerealarray(s, ref tree.splits);
2057            kdtreealloctemporaries(tree, tree.n, tree.nx, tree.ny);
2058        }
2059
2060
2061        /*************************************************************************
2062        Rearranges nodes [I1,I2) using partition in D-th dimension with S as threshold.
2063        Returns split position I3: [I1,I3) and [I3,I2) are created as result.
2064
2065        This subroutine doesn't create tree structures, just rearranges nodes.
2066        *************************************************************************/
2067        private static void kdtreesplit(kdtree kdt,
2068            int i1,
2069            int i2,
2070            int d,
2071            double s,
2072            ref int i3)
2073        {
2074            int i = 0;
2075            int j = 0;
2076            int ileft = 0;
2077            int iright = 0;
2078            double v = 0;
2079
2080            i3 = 0;
2081
2082            alglib.ap.assert(kdt.n>0, "KDTreeSplit: internal error");
2083           
2084            //
2085            // split XY/Tags in two parts:
2086            // * [ILeft,IRight] is non-processed part of XY/Tags
2087            //
2088            // After cycle is done, we have Ileft=IRight. We deal with
2089            // this element separately.
2090            //
2091            // After this, [I1,ILeft) contains left part, and [ILeft,I2)
2092            // contains right part.
2093            //
2094            ileft = i1;
2095            iright = i2-1;
2096            while( ileft<iright )
2097            {
2098                if( (double)(kdt.xy[ileft,d])<=(double)(s) )
2099                {
2100                   
2101                    //
2102                    // XY[ILeft] is on its place.
2103                    // Advance ILeft.
2104                    //
2105                    ileft = ileft+1;
2106                }
2107                else
2108                {
2109                   
2110                    //
2111                    // XY[ILeft,..] must be at IRight.
2112                    // Swap and advance IRight.
2113                    //
2114                    for(i=0; i<=2*kdt.nx+kdt.ny-1; i++)
2115                    {
2116                        v = kdt.xy[ileft,i];
2117                        kdt.xy[ileft,i] = kdt.xy[iright,i];
2118                        kdt.xy[iright,i] = v;
2119                    }
2120                    j = kdt.tags[ileft];
2121                    kdt.tags[ileft] = kdt.tags[iright];
2122                    kdt.tags[iright] = j;
2123                    iright = iright-1;
2124                }
2125            }
2126            if( (double)(kdt.xy[ileft,d])<=(double)(s) )
2127            {
2128                ileft = ileft+1;
2129            }
2130            else
2131            {
2132                iright = iright-1;
2133            }
2134            i3 = ileft;
2135        }
2136
2137
2138        /*************************************************************************
2139        Recursive kd-tree generation subroutine.
2140
2141        PARAMETERS
2142            KDT         tree
2143            NodesOffs   unused part of Nodes[] which must be filled by tree
2144            SplitsOffs  unused part of Splits[]
2145            I1, I2      points from [I1,I2) are processed
2146           
2147        NodesOffs[] and SplitsOffs[] must be large enough.
2148
2149          -- ALGLIB --
2150             Copyright 28.02.2010 by Bochkanov Sergey
2151        *************************************************************************/
2152        private static void kdtreegeneratetreerec(kdtree kdt,
2153            ref int nodesoffs,
2154            ref int splitsoffs,
2155            int i1,
2156            int i2,
2157            int maxleafsize)
2158        {
2159            int n = 0;
2160            int nx = 0;
2161            int ny = 0;
2162            int i = 0;
2163            int j = 0;
2164            int oldoffs = 0;
2165            int i3 = 0;
2166            int cntless = 0;
2167            int cntgreater = 0;
2168            double minv = 0;
2169            double maxv = 0;
2170            int minidx = 0;
2171            int maxidx = 0;
2172            int d = 0;
2173            double ds = 0;
2174            double s = 0;
2175            double v = 0;
2176            int i_ = 0;
2177            int i1_ = 0;
2178
2179            alglib.ap.assert(kdt.n>0, "KDTreeGenerateTreeRec: internal error");
2180            alglib.ap.assert(i2>i1, "KDTreeGenerateTreeRec: internal error");
2181           
2182            //
2183            // Generate leaf if needed
2184            //
2185            if( i2-i1<=maxleafsize )
2186            {
2187                kdt.nodes[nodesoffs+0] = i2-i1;
2188                kdt.nodes[nodesoffs+1] = i1;
2189                nodesoffs = nodesoffs+2;
2190                return;
2191            }
2192           
2193            //
2194            // Load values for easier access
2195            //
2196            nx = kdt.nx;
2197            ny = kdt.ny;
2198           
2199            //
2200            // select dimension to split:
2201            // * D is a dimension number
2202            //
2203            d = 0;
2204            ds = kdt.curboxmax[0]-kdt.curboxmin[0];
2205            for(i=1; i<=nx-1; i++)
2206            {
2207                v = kdt.curboxmax[i]-kdt.curboxmin[i];
2208                if( (double)(v)>(double)(ds) )
2209                {
2210                    ds = v;
2211                    d = i;
2212                }
2213            }
2214           
2215            //
2216            // Select split position S using sliding midpoint rule,
2217            // rearrange points into [I1,I3) and [I3,I2)
2218            //
2219            s = kdt.curboxmin[d]+0.5*ds;
2220            i1_ = (i1) - (0);
2221            for(i_=0; i_<=i2-i1-1;i_++)
2222            {
2223                kdt.buf[i_] = kdt.xy[i_+i1_,d];
2224            }
2225            n = i2-i1;
2226            cntless = 0;
2227            cntgreater = 0;
2228            minv = kdt.buf[0];
2229            maxv = kdt.buf[0];
2230            minidx = i1;
2231            maxidx = i1;
2232            for(i=0; i<=n-1; i++)
2233            {
2234                v = kdt.buf[i];
2235                if( (double)(v)<(double)(minv) )
2236                {
2237                    minv = v;
2238                    minidx = i1+i;
2239                }
2240                if( (double)(v)>(double)(maxv) )
2241                {
2242                    maxv = v;
2243                    maxidx = i1+i;
2244                }
2245                if( (double)(v)<(double)(s) )
2246                {
2247                    cntless = cntless+1;
2248                }
2249                if( (double)(v)>(double)(s) )
2250                {
2251                    cntgreater = cntgreater+1;
2252                }
2253            }
2254            if( cntless>0 && cntgreater>0 )
2255            {
2256               
2257                //
2258                // normal midpoint split
2259                //
2260                kdtreesplit(kdt, i1, i2, d, s, ref i3);
2261            }
2262            else
2263            {
2264               
2265                //
2266                // sliding midpoint
2267                //
2268                if( cntless==0 )
2269                {
2270                   
2271                    //
2272                    // 1. move split to MinV,
2273                    // 2. place one point to the left bin (move to I1),
2274                    //    others - to the right bin
2275                    //
2276                    s = minv;
2277                    if( minidx!=i1 )
2278                    {
2279                        for(i=0; i<=2*kdt.nx+kdt.ny-1; i++)
2280                        {
2281                            v = kdt.xy[minidx,i];
2282                            kdt.xy[minidx,i] = kdt.xy[i1,i];
2283                            kdt.xy[i1,i] = v;
2284                        }
2285                        j = kdt.tags[minidx];
2286                        kdt.tags[minidx] = kdt.tags[i1];
2287                        kdt.tags[i1] = j;
2288                    }
2289                    i3 = i1+1;
2290                }
2291                else
2292                {
2293                   
2294                    //
2295                    // 1. move split to MaxV,
2296                    // 2. place one point to the right bin (move to I2-1),
2297                    //    others - to the left bin
2298                    //
2299                    s = maxv;
2300                    if( maxidx!=i2-1 )
2301                    {
2302                        for(i=0; i<=2*kdt.nx+kdt.ny-1; i++)
2303                        {
2304                            v = kdt.xy[maxidx,i];
2305                            kdt.xy[maxidx,i] = kdt.xy[i2-1,i];
2306                            kdt.xy[i2-1,i] = v;
2307                        }
2308                        j = kdt.tags[maxidx];
2309                        kdt.tags[maxidx] = kdt.tags[i2-1];
2310                        kdt.tags[i2-1] = j;
2311                    }
2312                    i3 = i2-1;
2313                }
2314            }
2315           
2316            //
2317            // Generate 'split' node
2318            //
2319            kdt.nodes[nodesoffs+0] = 0;
2320            kdt.nodes[nodesoffs+1] = d;
2321            kdt.nodes[nodesoffs+2] = splitsoffs;
2322            kdt.splits[splitsoffs+0] = s;
2323            oldoffs = nodesoffs;
2324            nodesoffs = nodesoffs+splitnodesize;
2325            splitsoffs = splitsoffs+1;
2326           
2327            //
2328            // Recirsive generation:
2329            // * update CurBox
2330            // * call subroutine
2331            // * restore CurBox
2332            //
2333            kdt.nodes[oldoffs+3] = nodesoffs;
2334            v = kdt.curboxmax[d];
2335            kdt.curboxmax[d] = s;
2336            kdtreegeneratetreerec(kdt, ref nodesoffs, ref splitsoffs, i1, i3, maxleafsize);
2337            kdt.curboxmax[d] = v;
2338            kdt.nodes[oldoffs+4] = nodesoffs;
2339            v = kdt.curboxmin[d];
2340            kdt.curboxmin[d] = s;
2341            kdtreegeneratetreerec(kdt, ref nodesoffs, ref splitsoffs, i3, i2, maxleafsize);
2342            kdt.curboxmin[d] = v;
2343        }
2344
2345
2346        /*************************************************************************
2347        Recursive subroutine for NN queries.
2348
2349          -- ALGLIB --
2350             Copyright 28.02.2010 by Bochkanov Sergey
2351        *************************************************************************/
2352        private static void kdtreequerynnrec(kdtree kdt,
2353            int offs)
2354        {
2355            double ptdist = 0;
2356            int i = 0;
2357            int j = 0;
2358            int nx = 0;
2359            int i1 = 0;
2360            int i2 = 0;
2361            int d = 0;
2362            double s = 0;
2363            double v = 0;
2364            double t1 = 0;
2365            int childbestoffs = 0;
2366            int childworstoffs = 0;
2367            int childoffs = 0;
2368            double prevdist = 0;
2369            bool todive = new bool();
2370            bool bestisleft = new bool();
2371            bool updatemin = new bool();
2372
2373            alglib.ap.assert(kdt.n>0, "KDTreeQueryNNRec: internal error");
2374           
2375            //
2376            // Leaf node.
2377            // Process points.
2378            //
2379            if( kdt.nodes[offs]>0 )
2380            {
2381                i1 = kdt.nodes[offs+1];
2382                i2 = i1+kdt.nodes[offs];
2383                for(i=i1; i<=i2-1; i++)
2384                {
2385                   
2386                    //
2387                    // Calculate distance
2388                    //
2389                    ptdist = 0;
2390                    nx = kdt.nx;
2391                    if( kdt.normtype==0 )
2392                    {
2393                        for(j=0; j<=nx-1; j++)
2394                        {
2395                            ptdist = Math.Max(ptdist, Math.Abs(kdt.xy[i,j]-kdt.x[j]));
2396                        }
2397                    }
2398                    if( kdt.normtype==1 )
2399                    {
2400                        for(j=0; j<=nx-1; j++)
2401                        {
2402                            ptdist = ptdist+Math.Abs(kdt.xy[i,j]-kdt.x[j]);
2403                        }
2404                    }
2405                    if( kdt.normtype==2 )
2406                    {
2407                        for(j=0; j<=nx-1; j++)
2408                        {
2409                            ptdist = ptdist+math.sqr(kdt.xy[i,j]-kdt.x[j]);
2410                        }
2411                    }
2412                   
2413                    //
2414                    // Skip points with zero distance if self-matches are turned off
2415                    //
2416                    if( (double)(ptdist)==(double)(0) && !kdt.selfmatch )
2417                    {
2418                        continue;
2419                    }
2420                   
2421                    //
2422                    // We CAN'T process point if R-criterion isn't satisfied,
2423                    // i.e. (RNeeded<>0) AND (PtDist>R).
2424                    //
2425                    if( (double)(kdt.rneeded)==(double)(0) || (double)(ptdist)<=(double)(kdt.rneeded) )
2426                    {
2427                       
2428                        //
2429                        // R-criterion is satisfied, we must either:
2430                        // * replace worst point, if (KNeeded<>0) AND (KCur=KNeeded)
2431                        //   (or skip, if worst point is better)
2432                        // * add point without replacement otherwise
2433                        //
2434                        if( kdt.kcur<kdt.kneeded || kdt.kneeded==0 )
2435                        {
2436                           
2437                            //
2438                            // add current point to heap without replacement
2439                            //
2440                            tsort.tagheappushi(ref kdt.r, ref kdt.idx, ref kdt.kcur, ptdist, i);
2441                        }
2442                        else
2443                        {
2444                           
2445                            //
2446                            // New points are added or not, depending on their distance.
2447                            // If added, they replace element at the top of the heap
2448                            //
2449                            if( (double)(ptdist)<(double)(kdt.r[0]) )
2450                            {
2451                                if( kdt.kneeded==1 )
2452                                {
2453                                    kdt.idx[0] = i;
2454                                    kdt.r[0] = ptdist;
2455                                }
2456                                else
2457                                {
2458                                    tsort.tagheapreplacetopi(ref kdt.r, ref kdt.idx, kdt.kneeded, ptdist, i);
2459                                }
2460                            }
2461                        }
2462                    }
2463                }
2464                return;
2465            }
2466           
2467            //
2468            // Simple split
2469            //
2470            if( kdt.nodes[offs]==0 )
2471            {
2472               
2473                //
2474                // Load:
2475                // * D  dimension to split
2476                // * S  split position
2477                //
2478                d = kdt.nodes[offs+1];
2479                s = kdt.splits[kdt.nodes[offs+2]];
2480               
2481                //
2482                // Calculate:
2483                // * ChildBestOffs      child box with best chances
2484                // * ChildWorstOffs     child box with worst chances
2485                //
2486                if( (double)(kdt.x[d])<=(double)(s) )
2487                {
2488                    childbestoffs = kdt.nodes[offs+3];
2489                    childworstoffs = kdt.nodes[offs+4];
2490                    bestisleft = true;
2491                }
2492                else
2493                {
2494                    childbestoffs = kdt.nodes[offs+4];
2495                    childworstoffs = kdt.nodes[offs+3];
2496                    bestisleft = false;
2497                }
2498               
2499                //
2500                // Navigate through childs
2501                //
2502                for(i=0; i<=1; i++)
2503                {
2504                   
2505                    //
2506                    // Select child to process:
2507                    // * ChildOffs      current child offset in Nodes[]
2508                    // * UpdateMin      whether minimum or maximum value
2509                    //                  of bounding box is changed on update
2510                    //
2511                    if( i==0 )
2512                    {
2513                        childoffs = childbestoffs;
2514                        updatemin = !bestisleft;
2515                    }
2516                    else
2517                    {
2518                        updatemin = bestisleft;
2519                        childoffs = childworstoffs;
2520                    }
2521                   
2522                    //
2523                    // Update bounding box and current distance
2524                    //
2525                    if( updatemin )
2526                    {
2527                        prevdist = kdt.curdist;
2528                        t1 = kdt.x[d];
2529                        v = kdt.curboxmin[d];
2530                        if( (double)(t1)<=(double)(s) )
2531                        {
2532                            if( kdt.normtype==0 )
2533                            {
2534                                kdt.curdist = Math.Max(kdt.curdist, s-t1);
2535                            }
2536                            if( kdt.normtype==1 )
2537                            {
2538                                kdt.curdist = kdt.curdist-Math.Max(v-t1, 0)+s-t1;
2539                            }
2540                            if( kdt.normtype==2 )
2541                            {
2542                                kdt.curdist = kdt.curdist-math.sqr(Math.Max(v-t1, 0))+math.sqr(s-t1);
2543                            }
2544                        }
2545                        kdt.curboxmin[d] = s;
2546                    }
2547                    else
2548                    {
2549                        prevdist = kdt.curdist;
2550                        t1 = kdt.x[d];
2551                        v = kdt.curboxmax[d];
2552                        if( (double)(t1)>=(double)(s) )
2553                        {
2554                            if( kdt.normtype==0 )
2555                            {
2556                                kdt.curdist = Math.Max(kdt.curdist, t1-s);
2557                            }
2558                            if( kdt.normtype==1 )
2559                            {
2560                                kdt.curdist = kdt.curdist-Math.Max(t1-v, 0)+t1-s;
2561                            }
2562                            if( kdt.normtype==2 )
2563                            {
2564                                kdt.curdist = kdt.curdist-math.sqr(Math.Max(t1-v, 0))+math.sqr(t1-s);
2565                            }
2566                        }
2567                        kdt.curboxmax[d] = s;
2568                    }
2569                   
2570                    //
2571                    // Decide: to dive into cell or not to dive
2572                    //
2573                    if( (double)(kdt.rneeded)!=(double)(0) && (double)(kdt.curdist)>(double)(kdt.rneeded) )
2574                    {
2575                        todive = false;
2576                    }
2577                    else
2578                    {
2579                        if( kdt.kcur<kdt.kneeded || kdt.kneeded==0 )
2580                        {
2581                           
2582                            //
2583                            // KCur<KNeeded (i.e. not all points are found)
2584                            //
2585                            todive = true;
2586                        }
2587                        else
2588                        {
2589                           
2590                            //
2591                            // KCur=KNeeded, decide to dive or not to dive
2592                            // using point position relative to bounding box.
2593                            //
2594                            todive = (double)(kdt.curdist)<=(double)(kdt.r[0]*kdt.approxf);
2595                        }
2596                    }
2597                    if( todive )
2598                    {
2599                        kdtreequerynnrec(kdt, childoffs);
2600                    }
2601                   
2602                    //
2603                    // Restore bounding box and distance
2604                    //
2605                    if( updatemin )
2606                    {
2607                        kdt.curboxmin[d] = v;
2608                    }
2609                    else
2610                    {
2611                        kdt.curboxmax[d] = v;
2612                    }
2613                    kdt.curdist = prevdist;
2614                }
2615                return;
2616            }
2617        }
2618
2619
2620        /*************************************************************************
2621        Copies X[] to KDT.X[]
2622        Loads distance from X[] to bounding box.
2623        Initializes CurBox[].
2624
2625          -- ALGLIB --
2626             Copyright 28.02.2010 by Bochkanov Sergey
2627        *************************************************************************/
2628        private static void kdtreeinitbox(kdtree kdt,
2629            double[] x)
2630        {
2631            int i = 0;
2632            double vx = 0;
2633            double vmin = 0;
2634            double vmax = 0;
2635
2636            alglib.ap.assert(kdt.n>0, "KDTreeInitBox: internal error");
2637           
2638            //
2639            // calculate distance from point to current bounding box
2640            //
2641            kdt.curdist = 0;
2642            if( kdt.normtype==0 )
2643            {
2644                for(i=0; i<=kdt.nx-1; i++)
2645                {
2646                    vx = x[i];
2647                    vmin = kdt.boxmin[i];
2648                    vmax = kdt.boxmax[i];
2649                    kdt.x[i] = vx;
2650                    kdt.curboxmin[i] = vmin;
2651                    kdt.curboxmax[i] = vmax;
2652                    if( (double)(vx)<(double)(vmin) )
2653                    {
2654                        kdt.curdist = Math.Max(kdt.curdist, vmin-vx);
2655                    }
2656                    else
2657                    {
2658                        if( (double)(vx)>(double)(vmax) )
2659                        {
2660                            kdt.curdist = Math.Max(kdt.curdist, vx-vmax);
2661                        }
2662                    }
2663                }
2664            }
2665            if( kdt.normtype==1 )
2666            {
2667                for(i=0; i<=kdt.nx-1; i++)
2668                {
2669                    vx = x[i];
2670                    vmin = kdt.boxmin[i];
2671                    vmax = kdt.boxmax[i];
2672                    kdt.x[i] = vx;
2673                    kdt.curboxmin[i] = vmin;
2674                    kdt.curboxmax[i] = vmax;
2675                    if( (double)(vx)<(double)(vmin) )
2676                    {
2677                        kdt.curdist = kdt.curdist+vmin-vx;
2678                    }
2679                    else
2680                    {
2681                        if( (double)(vx)>(double)(vmax) )
2682                        {
2683                            kdt.curdist = kdt.curdist+vx-vmax;
2684                        }
2685                    }
2686                }
2687            }
2688            if( kdt.normtype==2 )
2689            {
2690                for(i=0; i<=kdt.nx-1; i++)
2691                {
2692                    vx = x[i];
2693                    vmin = kdt.boxmin[i];
2694                    vmax = kdt.boxmax[i];
2695                    kdt.x[i] = vx;
2696                    kdt.curboxmin[i] = vmin;
2697                    kdt.curboxmax[i] = vmax;
2698                    if( (double)(vx)<(double)(vmin) )
2699                    {
2700                        kdt.curdist = kdt.curdist+math.sqr(vmin-vx);
2701                    }
2702                    else
2703                    {
2704                        if( (double)(vx)>(double)(vmax) )
2705                        {
2706                            kdt.curdist = kdt.curdist+math.sqr(vx-vmax);
2707                        }
2708                    }
2709                }
2710            }
2711        }
2712
2713
2714        /*************************************************************************
2715        This function allocates all dataset-independent array  fields  of  KDTree,
2716        i.e.  such  array  fields  that  their dimensions do not depend on dataset
2717        size.
2718
2719        This function do not sets KDT.NX or KDT.NY - it just allocates arrays
2720
2721          -- ALGLIB --
2722             Copyright 14.03.2011 by Bochkanov Sergey
2723        *************************************************************************/
2724        private static void kdtreeallocdatasetindependent(kdtree kdt,
2725            int nx,
2726            int ny)
2727        {
2728            alglib.ap.assert(kdt.n>0, "KDTreeAllocDatasetIndependent: internal error");
2729            kdt.x = new double[nx];
2730            kdt.boxmin = new double[nx];
2731            kdt.boxmax = new double[nx];
2732            kdt.curboxmin = new double[nx];
2733            kdt.curboxmax = new double[nx];
2734        }
2735
2736
2737        /*************************************************************************
2738        This function allocates all dataset-dependent array fields of KDTree, i.e.
2739        such array fields that their dimensions depend on dataset size.
2740
2741        This function do not sets KDT.N, KDT.NX or KDT.NY -
2742        it just allocates arrays.
2743
2744          -- ALGLIB --
2745             Copyright 14.03.2011 by Bochkanov Sergey
2746        *************************************************************************/
2747        private static void kdtreeallocdatasetdependent(kdtree kdt,
2748            int n,
2749            int nx,
2750            int ny)
2751        {
2752            alglib.ap.assert(n>0, "KDTreeAllocDatasetDependent: internal error");
2753            kdt.xy = new double[n, 2*nx+ny];
2754            kdt.tags = new int[n];
2755            kdt.idx = new int[n];
2756            kdt.r = new double[n];
2757            kdt.x = new double[nx];
2758            kdt.buf = new double[Math.Max(n, nx)];
2759            kdt.nodes = new int[splitnodesize*2*n];
2760            kdt.splits = new double[2*n];
2761        }
2762
2763
2764        /*************************************************************************
2765        This function allocates temporaries.
2766
2767        This function do not sets KDT.N, KDT.NX or KDT.NY -
2768        it just allocates arrays.
2769
2770          -- ALGLIB --
2771             Copyright 14.03.2011 by Bochkanov Sergey
2772        *************************************************************************/
2773        private static void kdtreealloctemporaries(kdtree kdt,
2774            int n,
2775            int nx,
2776            int ny)
2777        {
2778            alglib.ap.assert(n>0, "KDTreeAllocTemporaries: internal error");
2779            kdt.x = new double[nx];
2780            kdt.idx = new int[n];
2781            kdt.r = new double[n];
2782            kdt.buf = new double[Math.Max(n, nx)];
2783            kdt.curboxmin = new double[nx];
2784            kdt.curboxmax = new double[nx];
2785        }
2786
2787
2788    }
2789}
2790
Note: See TracBrowser for help on using the repository browser.