Free cookie consent management tool by TermsFeed Policy Generator

source: stable/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/3.7.0/ALGLIB-3.7.0/alglibmisc.cs @ 16204

Last change on this file since 16204 was 9268, checked in by mkommend, 12 years ago

#2007: Added AlgLib 3.7.0 to the external libraries.

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