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