GNU Octave  4.2.1
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
__lin_interpn__.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 2007-2017 Alexander Barth
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 #if defined (HAVE_CONFIG_H)
24 # include "config.h"
25 #endif
26 
27 #include "lo-ieee.h"
28 #include "dNDArray.h"
29 #include "oct-locbuf.h"
30 
31 #include "defun.h"
32 #include "error.h"
33 #include "ovl.h"
34 
35 // equivalent to isvector.m
36 
37 template <typename T>
38 bool
39 isvector (const T& array)
40 {
41  const dim_vector dv = array.dims ();
42  return dv.ndims () == 2 && (dv(0) == 1 || dv(1) == 1);
43 }
44 
45 // lookup a value in a sorted table (lookup.m)
46 template <typename T>
48 lookup (const T *x, octave_idx_type n, T y)
49 {
51 
52  if (x[0] < x[n-1])
53  {
54  // increasing x
55 
56  if (y > x[n-1] || y < x[0])
57  return -1;
58 
59 #if defined (EXHAUSTIF)
60  for (j = 0; j < n - 1; j++)
61  {
62  if (x[j] <= y && y <= x[j+1])
63  return j;
64  }
65 #else
66  octave_idx_type j0 = 0;
67  octave_idx_type j1 = n - 1;
68 
69  while (true)
70  {
71  j = (j0+j1)/2;
72 
73  if (y <= x[j+1])
74  {
75  if (x[j] <= y)
76  return j;
77 
78  j1 = j;
79  }
80 
81  if (x[j] <= y)
82  j0 = j;
83  }
84 #endif
85  }
86  else
87  {
88  // decreasing x
89  // previous code with x -> -x and y -> -y
90 
91  if (y > x[0] || y < x[n-1])
92  return -1;
93 
94 #if defined (EXHAUSTIF)
95  for (j = 0; j < n - 1; j++)
96  {
97  if (x[j+1] <= y && y <= x[j])
98  return j;
99  }
100 #else
101  octave_idx_type j0 = 0;
102  octave_idx_type j1 = n - 1;
103 
104  while (true)
105  {
106  j = (j0+j1)/2;
107 
108  if (y >= x[j+1])
109  {
110  if (x[j] >= y)
111  return j;
112 
113  j1 = j;
114  }
115 
116  if (x[j] >= y)
117  j0 = j;
118  }
119 #endif
120  }
121 }
122 
123 // n-dimensional linear interpolation
124 
125 template <typename T>
126 void
128  octave_idx_type Ni, T extrapval, const T **x,
129  const T *v, const T **y, T *vi)
130 {
131  bool out = false;
132  int bit;
133 
134  OCTAVE_LOCAL_BUFFER (T, coef, 2*n);
136 
137  // loop over all points
138  for (octave_idx_type m = 0; m < Ni; m++)
139  {
140  // loop over all dimensions
141  for (int i = 0; i < n; i++)
142  {
143  index[i] = lookup (x[i], size[i], y[i][m]);
144  out = index[i] == -1;
145 
146  if (out)
147  break;
148  else
149  {
150  octave_idx_type j = index[i];
151  coef[2*i+1] = (y[i][m] - x[i][j])/(x[i][j+1] - x[i][j]);
152  coef[2*i] = 1 - coef[2*i+1];
153  }
154  }
155 
156  if (out)
157  vi[m] = extrapval;
158  else
159  {
160  vi[m] = 0;
161 
162  // loop over all corners of hypercube (1<<n = 2^n)
163  for (int i = 0; i < (1 << n); i++)
164  {
165  T c = 1;
166  octave_idx_type l = 0;
167 
168  // loop over all dimensions
169  for (int j = 0; j < n; j++)
170  {
171  // test if the jth bit in i is set
172  bit = i >> j & 1;
173  l += scale[j] * (index[j] + bit);
174  c *= coef[2*j+bit];
175  }
176 
177  vi[m] += c * v[l];
178  }
179  }
180  }
181 }
182 
183 template <typename T, typename M>
185 lin_interpn (int n, M *X, const M V, M *Y)
186 {
188 
189  M Vi = M (Y[0].dims ());
190 
191  OCTAVE_LOCAL_BUFFER (const T *, y, n);
193 
194  for (int i = 0; i < n; i++)
195  {
196  y[i] = Y[i].data ();
197  size[i] = V.dims ()(i);
198  }
199 
200  OCTAVE_LOCAL_BUFFER (const T *, x, n);
202 
203  const T *v = V.data ();
204  T *vi = Vi.fortran_vec ();
205  octave_idx_type Ni = Vi.numel ();
206 
207  T extrapval = octave_NA;
208 
209  // offset in memory of each dimension
210 
211  scale[0] = 1;
212 
213  for (int i = 1; i < n; i++)
214  scale[i] = scale[i-1] * size[i-1];
215 
216  // tests if X[0] is a vector, if yes, assume that all elements of X are
217  // in the ndgrid format.
218 
219  if (! isvector (X[0]))
220  {
221  for (int i = 0; i < n; i++)
222  {
223  if (X[i].dims () != V.dims ())
224  error ("interpn: incompatible size of argument number %d", i+1);
225 
226  M tmp = M (dim_vector (size[i], 1));
227 
228  for (octave_idx_type j = 0; j < size[i]; j++)
229  tmp(j) = X[i](scale[i]*j);
230 
231  X[i] = tmp;
232  }
233  }
234 
235  for (int i = 0; i < n; i++)
236  {
237  if (! isvector (X[i]) && X[i].numel () != size[i])
238  error ("interpn: incompatible size of argument number %d", i+1);
239 
240  x[i] = X[i].data ();
241  }
242 
243  lin_interpn (n, size, scale, Ni, extrapval, x, v, y, vi);
244 
245  retval = Vi;
246 
247  return retval;
248 }
249 
250 // Perform @var{n}-dimensional interpolation. Each element of then
251 // @var{n}-dimensional array @var{v} represents a value at a location
252 // given by the parameters @var{x1}, @var{x2},...,@var{xn}. The parameters
253 // @var{x1}, @var{x2}, @dots{}, @var{xn} are either @var{n}-dimensional
254 // arrays of the same size as the array @var{v} in the \"ndgrid\" format
255 // or vectors. The parameters @var{y1}, @var{y2}, @dots{}, @var{yn} are
256 // all @var{n}-dimensional arrays of the same size and represent the
257 // points at which the array @var{vi} is interpolated.
258 //
259 //This function only performs linear interpolation.
260 
261 DEFUN (__lin_interpn__, args, ,
262  doc: /* -*- texinfo -*-
263 @deftypefn {} {@var{vi} =} __lin_interpn__ (@var{x1}, @var{x2}, @dots{}, @var{xn}, @var{v}, @var{y1}, @var{y2}, @dots{}, @var{yn})
264 Undocumented internal function.
265 @end deftypefn */)
266 {
267  int nargin = args.length ();
268 
269  if (nargin < 2 || nargin % 2 == 0)
270  print_usage ();
271 
273 
274  // dimension of the problem
275  int n = (nargin-1)/2;
276 
277  if (args(n).is_single_type ())
278  {
281 
282  const FloatNDArray V = args(n).float_array_value ();
283 
284  for (int i = 0; i < n; i++)
285  {
286  X[i] = args(i).float_array_value ();
287  Y[i] = args(n+i+1).float_array_value ();
288 
289  if (Y[0].dims () != Y[i].dims ())
290  error ("interpn: incompatible size of argument number %d", n+i+2);
291  }
292 
293  retval = lin_interpn<float, FloatNDArray> (n, X, V, Y);
294  }
295  else
296  {
299 
300  const NDArray V = args(n).array_value ();
301 
302  for (int i = 0; i < n; i++)
303  {
304  X[i] = args(i).array_value ();
305  Y[i] = args(n+i+1).array_value ();
306 
307  if (Y[0].dims () != Y[i].dims ())
308  error ("interpn: incompatible size of argument number %d", n+i+2);
309  }
310 
311  retval = lin_interpn<double, NDArray> (n, X, V, Y);
312  }
313 
314  return retval;
315 }
316 
317 /*
318 ## No test needed for internal helper function.
319 %!assert (1)
320 */
F77_RET_T const F77_INT const F77_INT const F77_INT const F77_DBLE const F77_DBLE F77_INT F77_DBLE * V
OCTINTERP_API void print_usage(void)
Definition: defun.cc:52
#define DEFUN(name, args_name, nargout_name, doc)
Definition: defun.h:46
void error(const char *fmt,...)
Definition: error.cc:570
octave_idx_type lookup(const T *x, octave_idx_type n, T y)
F77_RET_T const F77_INT const F77_INT const F77_INT const F77_DBLE const F77_DBLE F77_INT & M
OCTAVE_EXPORT octave_value_list return the value of the option it must match the dimension of the state and the relative tolerance must also be a vector of the same length tem it must match the dimension of the state and the absolute tolerance must also be a vector of the same length The local error test applied at each integration step is xample roup calculate Y_a and Y _d item Given Y
Definition: DASPK-opts.cc:739
JNIEnv void * args
Definition: ov-java.cc:67
bool isvector(const T &array)
nd deftypefn *octave_map m
Definition: ov-struct.cc:2058
then the function must return scalars which will be concatenated into the return array(s).If code
Definition: cellfun.cc:398
int nargin
Definition: graphics.cc:10115
double tmp
Definition: data.cc:6300
octave_value retval
Definition: data.cc:6294
#define octave_NA
Definition: lo-ieee.h:36
the sparsity preserving column transformation such that that defines the pivoting threshold can be given in which case it defines the c
Definition: lu.cc:138
the exceeded dimensions are set to if fewer subscripts than dimensions are the exceeding dimensions are merged into the final requested dimension For consider the following dims
Definition: sub2ind.cc:255
=val(i)}if ode{val(i)}occurs in table i
Definition: lookup.cc:239
octave_idx_type ndims(void) const
Number of dimensions.
Definition: dim-vector.h:301
the element is set to zero In other the statement xample y
Definition: data.cc:5342
void scale(Matrix &m, double x, double y, double z)
Definition: graphics.cc:5150
#define OCTAVE_LOCAL_BUFFER(T, buf, size)
Definition: oct-locbuf.h:200
OCTAVE_EXPORT octave_value_list any number nd example oindent prints the prompt xample Pick a any number!nd example oindent and waits for the user to enter a value The string entered by the user is evaluated as an so it may be a literal a variable or any other valid Octave code The number of return their size
Definition: input.cc:871
Vector representing the dimensions (size) of an Array.
Definition: dim-vector.h:87
dim_vector dv
Definition: sub2ind.cc:263
void lin_interpn(int n, const octave_idx_type *size, const octave_idx_type *scale, octave_idx_type Ni, T extrapval, const T **x, const T *v, const T **y, T *vi)
F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_REAL const F77_REAL F77_REAL &F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_DBLE F77_DBLE &F77_RET_T const F77_REAL F77_REAL &F77_RET_T F77_REAL F77_REAL &F77_RET_T F77_DBLE F77_DBLE &F77_RET_T const F77_DBLE * x