Mercurial > hg > octave-nkf
annotate libinterp/corefcn/__lin_interpn__.cc @ 15880:df1aceb8f0bc ss-3-7-1
snapshot version 3.7.1
* configure.ac (AC_INIT): Set version to 3.7.1.
(OCTAVE_RELEASE_DATE): Set to 2013-01-02.
author | John W. Eaton <jwe@octave.org> |
---|---|
date | Wed, 02 Jan 2013 15:16:08 -0500 |
parents | 2fc554ffbc28 |
children | d63878346099 |
rev | line source |
---|---|
6702 | 1 /* |
2 | |
14138
72c96de7a403
maint: update copyright notices for 2012
John W. Eaton <jwe@octave.org>
parents:
12805
diff
changeset
|
3 Copyright (C) 2007-2012 Alexander Barth |
6702 | 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 | |
7016 | 9 Free Software Foundation; either version 3 of the License, or (at your |
10 option) any later version. | |
6702 | 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 | |
7016 | 18 along with Octave; see the file COPYING. If not, see |
19 <http://www.gnu.org/licenses/>. | |
6702 | 20 |
21 */ | |
22 | |
23 #ifdef HAVE_CONFIG_H | |
24 #include <config.h> | |
25 #endif | |
26 | |
9003
0631d397fbe0
replace lo_ieee_isnan by xisnan, add missing includes
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
27 #include "lo-ieee.h" |
6702 | 28 #include "dNDArray.h" |
8377
25bc2d31e1bf
improve OCTAVE_LOCAL_BUFFER
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
29 #include "oct-locbuf.h" |
6702 | 30 |
15039
e753177cde93
maint: Move non-dynamically linked functions from DLD-FUNCTIONS/ to corefcn/ directory
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
31 #include "defun.h" |
6702 | 32 #include "error.h" |
33 #include "oct-obj.h" | |
34 | |
35 // equivalent to isvector.m | |
36 | |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
37 template <class T> |
6702 | 38 bool |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
39 isvector (const T& array) |
6702 | 40 { |
41 const dim_vector dv = array.dims (); | |
42 return dv.length () == 2 && (dv(0) == 1 || dv(1) == 1); | |
43 } | |
44 | |
45 // lookup a value in a sorted table (lookup.m) | |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
46 template <class T> |
6702 | 47 octave_idx_type |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
48 lookup (const T *x, octave_idx_type n, T y) |
6702 | 49 { |
7561
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
50 octave_idx_type j; |
6702 | 51 |
7561
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
52 if (x[0] < x[n-1]) |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
53 { |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
54 // increasing x |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
55 |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
56 if (y > x[n-1] || y < x[0]) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
57 return -1; |
6702 | 58 |
59 #ifdef EXHAUSTIF | |
7561
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
60 for (j = 0; j < n - 1; j++) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
61 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
62 if (x[j] <= y && y <= x[j+1]) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
63 return j; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
64 } |
6702 | 65 #else |
7561
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
66 octave_idx_type j0 = 0; |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
67 octave_idx_type j1 = n - 1; |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
68 |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
69 while (true) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
70 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
71 j = (j0+j1)/2; |
7561
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
72 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
73 if (y <= x[j+1]) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
74 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
75 if (x[j] <= y) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
76 return j; |
6702 | 77 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
78 j1 = j; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
79 } |
7561
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
80 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
81 if (x[j] <= y) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
82 j0 = j; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
83 } |
7561
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
84 #endif |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
85 } |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
86 else |
6702 | 87 { |
7561
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
88 // decreasing x |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
89 // previous code with x -> -x and y -> -y |
6702 | 90 |
7561
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
91 if (y > x[0] || y < x[n-1]) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
92 return -1; |
7561
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
93 |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
94 #ifdef EXHAUSTIF |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
95 for (j = 0; j < n - 1; j++) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
96 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
97 if (x[j+1] <= y && y <= x[j]) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
98 return j; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
99 } |
7561
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
100 #else |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
101 octave_idx_type j0 = 0; |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
102 octave_idx_type j1 = n - 1; |
6702 | 103 |
7561
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
104 while (true) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
105 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
106 j = (j0+j1)/2; |
7561
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
107 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
108 if (y >= x[j+1]) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
109 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
110 if (x[j] >= y) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
111 return j; |
7561
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
112 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
113 j1 = j; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
114 } |
7561
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
115 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
116 if (x[j] >= y) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
117 j0 = j; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
118 } |
7561
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
119 #endif |
6702 | 120 } |
121 } | |
122 | |
123 // n-dimensional linear interpolation | |
124 | |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
125 template <class T> |
6702 | 126 void |
127 lin_interpn (int n, const octave_idx_type *size, const octave_idx_type *scale, | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
128 octave_idx_type Ni, T extrapval, const T **x, |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
129 const T *v, const T **y, T *vi) |
6702 | 130 { |
131 bool out = false; | |
132 int bit; | |
133 | |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
134 OCTAVE_LOCAL_BUFFER (T, coef, 2*n); |
6702 | 135 OCTAVE_LOCAL_BUFFER (octave_idx_type, index, 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++) | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
142 { |
6702 | 143 index[i] = lookup (x[i], size[i], y[i][m]); |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
144 out = index[i] == -1; |
6702 | 145 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
146 if (out) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
147 break; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
148 else |
6702 | 149 { |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
150 octave_idx_type j = index[i]; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
151 coef[2*i+1] = (y[i][m] - x[i][j])/(x[i][j+1] - x[i][j]); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
152 coef[2*i] = 1 - coef[2*i+1]; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
153 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
154 } |
6702 | 155 |
156 | |
157 if (out) | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
158 vi[m] = extrapval; |
6702 | 159 else |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
160 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
161 vi[m] = 0; |
6702 | 162 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
163 // loop over all corners of hypercube (1<<n = 2^n) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
164 for (int i = 0; i < (1 << n); i++) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
165 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
166 T c = 1; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
167 octave_idx_type l = 0; |
6702 | 168 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
169 // loop over all dimensions |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
170 for (int j = 0; j < n; j++) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
171 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
172 // test if the jth bit in i is set |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
173 bit = i >> j & 1; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
174 l += scale[j] * (index[j] + bit); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
175 c *= coef[2*j+bit]; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
176 } |
6702 | 177 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
178 vi[m] += c * v[l]; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
179 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
180 } |
6702 | 181 } |
182 } | |
183 | |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
184 template <class T, class M> |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
185 octave_value |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
186 lin_interpn (int n, M *X, const M V, M *Y) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
187 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
188 octave_value retval; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
189 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
190 M Vi = M (Y[0].dims ()); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
191 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
192 OCTAVE_LOCAL_BUFFER (const T *, y, n); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
193 OCTAVE_LOCAL_BUFFER (octave_idx_type, size, n); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
194 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
195 for (int i = 0; i < n; i++) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
196 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
197 y[i] = Y[i].data (); |
14846
460a3c6d8bf1
maint: Use Octave coding convention for cuddled parenthis in function calls with empty argument lists.
Rik <octave@nomad.inbox5.com>
parents:
14501
diff
changeset
|
198 size[i] = V.dims ()(i); |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
199 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
200 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
201 OCTAVE_LOCAL_BUFFER (const T *, x, n); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
202 OCTAVE_LOCAL_BUFFER (octave_idx_type, scale, n); |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
203 |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
204 const T *v = V.data (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
205 T *vi = Vi.fortran_vec (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
206 octave_idx_type Ni = Vi.numel (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
207 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
208 T extrapval = octave_NA; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
209 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
210 // offset in memory of each dimension |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
211 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
212 scale[0] = 1; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
213 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
214 for (int i = 1; i < n; i++) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
215 scale[i] = scale[i-1] * size[i-1]; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
216 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
217 // tests if X[0] is a vector, if yes, assume that all elements of X are |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
218 // in the ndgrid format. |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
219 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
220 if (! isvector (X[0])) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
221 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
222 for (int i = 0; i < n; i++) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
223 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
224 if (X[i].dims () != V.dims ()) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
225 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
226 error ("interpn: incompatible size of argument number %d", i+1); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
227 return retval; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
228 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
229 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
230 { |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
231 M tmp = M (dim_vector (size[i], 1)); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
232 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
233 for (octave_idx_type j = 0; j < size[i]; j++) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
234 tmp(j) = X[i](scale[i]*j); |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
235 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
236 X[i] = tmp; |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
237 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
238 } |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
239 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
240 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
241 for (int i = 0; i < n; i++) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
242 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
243 if (! isvector (X[i]) && X[i].numel () != size[i]) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
244 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
245 error ("interpn: incompatible size of argument number %d", i+1); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
246 return retval; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
247 } |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
248 else |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
249 x[i] = X[i].data (); |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
250 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
251 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
252 lin_interpn (n, size, scale, Ni, extrapval, x, v, y, vi); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
253 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
254 retval = Vi; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
255 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
256 return retval; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
257 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
258 |
6945 | 259 // Perform @var{n}-dimensional interpolation. Each element of then |
260 // @var{n}-dimensional array @var{v} represents a value at a location | |
261 // given by the parameters @var{x1}, @var{x2},...,@var{xn}. The parameters | |
262 // @var{x1}, @var{x2}, @dots{}, @var{xn} are either @var{n}-dimensional | |
263 // arrays of the same size as the array @var{v} in the \"ndgrid\" format | |
264 // or vectors. The parameters @var{y1}, @var{y2}, @dots{}, @var{yn} are | |
265 // all @var{n}-dimensional arrays of the same size and represent the | |
266 // points at which the array @var{vi} is interpolated. | |
267 // | |
268 //This function only performs linear interpolation. | |
269 | |
15039
e753177cde93
maint: Move non-dynamically linked functions from DLD-FUNCTIONS/ to corefcn/ directory
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
270 DEFUN (__lin_interpn__, args, , |
6702 | 271 "-*- texinfo -*-\n\ |
15039
e753177cde93
maint: Move non-dynamically linked functions from DLD-FUNCTIONS/ to corefcn/ directory
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
272 @deftypefn {Built-in Function} {@var{vi} =} __lin_interpn__ (@var{x1}, @var{x2}, @dots{}, @var{xn}, @var{v}, @var{y1}, @var{y2}, @dots{}, @var{yn})\n\ |
6945 | 273 Undocumented internal function.\n\ |
6702 | 274 @end deftypefn") |
275 { | |
276 octave_value retval; | |
277 | |
278 int nargin = args.length (); | |
279 | |
280 if (nargin < 2 || nargin % 2 == 0) | |
281 { | |
282 print_usage (); | |
283 return retval; | |
284 } | |
285 | |
286 // dimension of the problem | |
287 int n = (nargin-1)/2; | |
288 | |
14846
460a3c6d8bf1
maint: Use Octave coding convention for cuddled parenthis in function calls with empty argument lists.
Rik <octave@nomad.inbox5.com>
parents:
14501
diff
changeset
|
289 if (args(n).is_single_type ()) |
6702 | 290 { |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
291 OCTAVE_LOCAL_BUFFER (FloatNDArray, X, n); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
292 OCTAVE_LOCAL_BUFFER (FloatNDArray, Y, n); |
6702 | 293 |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
294 const FloatNDArray V = args(n).float_array_value (); |
6702 | 295 |
296 if (error_state) | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
297 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
298 print_usage (); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
299 return retval; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
300 } |
6702 | 301 |
302 for (int i = 0; i < n; i++) | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
303 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
304 X[i] = args(i).float_array_value (); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
305 Y[i] = args(n+i+1).float_array_value (); |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
306 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
307 if (error_state) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
308 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
309 print_usage (); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
310 return retval; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
311 } |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
312 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
313 if (Y[0].dims () != Y[i].dims ()) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
314 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
315 error ("interpn: incompatible size of argument number %d", n+i+2); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
316 return retval; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
317 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
318 } |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
319 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
320 retval = lin_interpn<float, FloatNDArray> (n, X, V, Y); |
6702 | 321 } |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
322 else |
6702 | 323 { |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
324 OCTAVE_LOCAL_BUFFER (NDArray, X, n); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
325 OCTAVE_LOCAL_BUFFER (NDArray, Y, n); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
326 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
327 const NDArray V = args(n).array_value (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
328 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
329 if (error_state) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
330 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
331 print_usage (); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
332 return retval; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
333 } |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
334 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
335 for (int i = 0; i < n; i++) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
336 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
337 X[i] = args(i).array_value (); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
338 Y[i] = args(n+i+1).array_value (); |
6702 | 339 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
340 if (error_state) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
341 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
342 print_usage (); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
343 return retval; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
344 } |
6702 | 345 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
346 if (Y[0].dims () != Y[i].dims ()) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
347 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
348 error ("interpn: incompatible size of argument number %d", n+i+2); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
349 return retval; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
350 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
351 } |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
352 |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
353 retval = lin_interpn<double, NDArray> (n, X, V, Y); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
354 } |
6702 | 355 |
356 return retval; | |
357 } | |
12805
3641167e5b75
codesprint: *.cc helper functions do not need tests
Rik <octave@nomad.inbox5.com>
parents:
11586
diff
changeset
|
358 |
3641167e5b75
codesprint: *.cc helper functions do not need tests
Rik <octave@nomad.inbox5.com>
parents:
11586
diff
changeset
|
359 /* |
3641167e5b75
codesprint: *.cc helper functions do not need tests
Rik <octave@nomad.inbox5.com>
parents:
11586
diff
changeset
|
360 ## No test needed for internal helper function. |
3641167e5b75
codesprint: *.cc helper functions do not need tests
Rik <octave@nomad.inbox5.com>
parents:
11586
diff
changeset
|
361 %!assert (1) |
3641167e5b75
codesprint: *.cc helper functions do not need tests
Rik <octave@nomad.inbox5.com>
parents:
11586
diff
changeset
|
362 */ |