GNU Octave  3.8.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
__dsearchn__.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 2007-2013 David Bateman
4 
5 This file is part of Octave.
6 
7 Octave is free software; you can redistribute it and/or modify it
8 under the terms of the GNU General Public License as published by the
9 Free Software Foundation; either version 3 of the License, or (at your
10 option) any later version.
11 
12 Octave is distributed in the hope that it will be useful, but WITHOUT
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
15 for more details.
16 
17 You should have received a copy of the GNU General Public License
18 along with Octave; see the file COPYING. If not, see
19 <http://www.gnu.org/licenses/>.
20 
21 */
22 
23 #ifdef HAVE_CONFIG_H
24 #include <config.h>
25 #endif
26 
27 #include <iostream>
28 #include <fstream>
29 #include <string>
30 
31 #include "lo-math.h"
32 
33 #include "defun-dld.h"
34 #include "error.h"
35 #include "oct-obj.h"
36 
37 DEFUN_DLD (__dsearchn__, args, ,
38  "-*- texinfo -*-\n\
39 @deftypefn {Loadable Function} {[@var{idx}, @var{d}] =} dsearch (@var{x}, @var{xi})\n\
40 Undocumented internal function.\n\
41 @end deftypefn")
42 {
43  int nargin = args.length ();
44  octave_value_list retval;
45 
46  if (nargin != 2)
47  {
48  print_usage ();
49  return retval;
50  }
51 
52  Matrix x = args(0).matrix_value ().transpose ();
53  Matrix xi = args(1).matrix_value ().transpose ();
54 
55  if (! error_state)
56  {
57  if (x.rows () != xi.rows () || x.columns () < 1)
58  error ("__dsearch__: number of rows of X and XI must match");
59  else
60  {
61  octave_idx_type n = x.rows ();
62  octave_idx_type nx = x.columns ();
63  octave_idx_type nxi = xi.columns ();
64 
65  ColumnVector idx (nxi);
66  double *pidx = idx.fortran_vec ();
67  ColumnVector dist (nxi);
68  double *pdist = dist.fortran_vec ();
69 
70 #define DIST(dd, y, yi, m) \
71  dd = 0.; \
72  for (octave_idx_type k = 0; k < m; k++) \
73  { \
74  double yd = y[k] - yi[k]; \
75  dd += yd * yd; \
76  } \
77  dd = sqrt (dd);
78 
79  const double *pxi = xi.fortran_vec ();
80  for (octave_idx_type i = 0; i < nxi; i++)
81  {
82  double d0;
83  const double *px = x.fortran_vec ();
84  DIST(d0, px, pxi, n);
85  *pidx = 1.;
86  for (octave_idx_type j = 1; j < nx; j++)
87  {
88  px += n;
89  double d;
90  DIST (d, px, pxi, n);
91  if (d < d0)
92  {
93  d0 = d;
94  *pidx = static_cast<double>(j + 1);
95  }
97  }
98 
99  *pdist++ = d0;
100  pidx++;
101  pxi += n;
102  }
103 
104  retval(1) = dist;
105  retval(0) = idx;
106  }
107  }
108 
109  return retval;
110 }
111 
112 /*
113 ## No test needed for internal helper function.
114 %!assert (1)
115 */