Mercurial > octave
annotate libinterp/corefcn/__lin_interpn__.cc @ 27919:1891570abac8
update Octave Project Developers copyright for the new year
In files that have the "Octave Project Developers" copyright notice,
update for 2020.
author | John W. Eaton <jwe@octave.org> |
---|---|
date | Mon, 06 Jan 2020 22:29:51 -0500 |
parents | b442ec6dda5c |
children | bd51beb6205e |
rev | line source |
---|---|
6702 | 1 /* |
2 | |
27919
1891570abac8
update Octave Project Developers copyright for the new year
John W. Eaton <jwe@octave.org>
parents:
27918
diff
changeset
|
3 Copyright (C) 2007-2020 The Octave Project Developers |
27918
b442ec6dda5c
use centralized file for copyright info for individual contributors
John W. Eaton <jwe@octave.org>
parents:
26376
diff
changeset
|
4 |
b442ec6dda5c
use centralized file for copyright info for individual contributors
John W. Eaton <jwe@octave.org>
parents:
26376
diff
changeset
|
5 See the file COPYRIGHT.md in the top-level directory of this distribution |
b442ec6dda5c
use centralized file for copyright info for individual contributors
John W. Eaton <jwe@octave.org>
parents:
26376
diff
changeset
|
6 or <https://octave.org/COPYRIGHT.html/>. |
b442ec6dda5c
use centralized file for copyright info for individual contributors
John W. Eaton <jwe@octave.org>
parents:
26376
diff
changeset
|
7 |
6702 | 8 |
9 This file is part of Octave. | |
10 | |
24534
194eb4bd202b
maint: Update punctuation for GPL v3 license text.
Rik <rik@octave.org>
parents:
23803
diff
changeset
|
11 Octave is free software: you can redistribute it and/or modify it |
22755
3a2b891d0b33
maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents:
22407
diff
changeset
|
12 under the terms of the GNU General Public License as published by |
24534
194eb4bd202b
maint: Update punctuation for GPL v3 license text.
Rik <rik@octave.org>
parents:
23803
diff
changeset
|
13 the Free Software Foundation, either version 3 of the License, or |
22755
3a2b891d0b33
maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents:
22407
diff
changeset
|
14 (at your option) any later version. |
6702 | 15 |
22755
3a2b891d0b33
maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents:
22407
diff
changeset
|
16 Octave is distributed in the hope that it will be useful, but |
3a2b891d0b33
maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents:
22407
diff
changeset
|
17 WITHOUT ANY WARRANTY; without even the implied warranty of |
3a2b891d0b33
maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents:
22407
diff
changeset
|
18 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
3a2b891d0b33
maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents:
22407
diff
changeset
|
19 GNU General Public License for more details. |
6702 | 20 |
21 You should have received a copy of the GNU General Public License | |
7016 | 22 along with Octave; see the file COPYING. If not, see |
24534
194eb4bd202b
maint: Update punctuation for GPL v3 license text.
Rik <rik@octave.org>
parents:
23803
diff
changeset
|
23 <https://www.gnu.org/licenses/>. |
6702 | 24 |
25 */ | |
26 | |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21660
diff
changeset
|
27 #if defined (HAVE_CONFIG_H) |
21301
40de9f8f23a6
Use '#include "config.h"' rather than <config.h>.
Rik <rik@octave.org>
parents:
21200
diff
changeset
|
28 # include "config.h" |
6702 | 29 #endif |
30 | |
9003
0631d397fbe0
replace lo_ieee_isnan by xisnan, add missing includes
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
31 #include "lo-ieee.h" |
6702 | 32 #include "dNDArray.h" |
8377
25bc2d31e1bf
improve OCTAVE_LOCAL_BUFFER
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
33 #include "oct-locbuf.h" |
6702 | 34 |
15039
e753177cde93
maint: Move non-dynamically linked functions from DLD-FUNCTIONS/ to corefcn/ directory
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
35 #include "defun.h" |
6702 | 36 #include "error.h" |
20940
48b2ad5ee801
maint: Rename oct-obj.[cc|h] to ovl.[cc|h] for clarity.
Rik <rik@octave.org>
parents:
20918
diff
changeset
|
37 #include "ovl.h" |
6702 | 38 |
39 // equivalent to isvector.m | |
40 | |
21139
538b57866b90
consistently use "typename" intead of "class" in template declarations
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
41 template <typename T> |
6702 | 42 bool |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
43 isvector (const T& array) |
6702 | 44 { |
45 const dim_vector dv = array.dims (); | |
21660
53728df3e4c9
maint: for clarity, call ndims() rather than length() on dimension vectors.
Rik <rik@octave.org>
parents:
21568
diff
changeset
|
46 return dv.ndims () == 2 && (dv(0) == 1 || dv(1) == 1); |
6702 | 47 } |
48 | |
49 // lookup a value in a sorted table (lookup.m) | |
21139
538b57866b90
consistently use "typename" intead of "class" in template declarations
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
50 template <typename T> |
6702 | 51 octave_idx_type |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
52 lookup (const T *x, octave_idx_type n, T y) |
6702 | 53 { |
7561
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
54 octave_idx_type j; |
6702 | 55 |
7561
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
56 if (x[0] < x[n-1]) |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
57 { |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
58 // increasing x |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
59 |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
60 if (y > x[n-1] || y < x[0]) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
61 return -1; |
6702 | 62 |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21660
diff
changeset
|
63 #if defined (EXHAUSTIF) |
7561
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
64 for (j = 0; j < n - 1; j++) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
65 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
66 if (x[j] <= y && y <= x[j+1]) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
67 return j; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
68 } |
6702 | 69 #else |
7561
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
70 octave_idx_type j0 = 0; |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
71 octave_idx_type j1 = n - 1; |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
72 |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
73 while (true) |
10154
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 j = (j0+j1)/2; |
7561
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
76 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
77 if (y <= x[j+1]) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
78 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
79 if (x[j] <= y) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
80 return j; |
6702 | 81 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
82 j1 = 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 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
85 if (x[j] <= y) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
86 j0 = j; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
87 } |
7561
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
88 #endif |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
89 } |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
90 else |
6702 | 91 { |
7561
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
92 // decreasing x |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
93 // previous code with x -> -x and y -> -y |
6702 | 94 |
7561
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
95 if (y > x[0] || y < x[n-1]) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
96 return -1; |
7561
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
97 |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21660
diff
changeset
|
98 #if defined (EXHAUSTIF) |
7561
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
99 for (j = 0; j < n - 1; j++) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
100 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
101 if (x[j+1] <= y && y <= x[j]) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
102 return j; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
103 } |
7561
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
104 #else |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
105 octave_idx_type j0 = 0; |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
106 octave_idx_type j1 = n - 1; |
6702 | 107 |
7561
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
108 while (true) |
10154
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 j = (j0+j1)/2; |
7561
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
111 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
112 if (y >= x[j+1]) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
113 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
114 if (x[j] >= y) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
115 return j; |
7561
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
116 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
117 j1 = 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 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
120 if (x[j] >= y) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
121 j0 = j; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
122 } |
7561
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
123 #endif |
6702 | 124 } |
125 } | |
126 | |
127 // n-dimensional linear interpolation | |
128 | |
21139
538b57866b90
consistently use "typename" intead of "class" in template declarations
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
129 template <typename T> |
6702 | 130 void |
131 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
|
132 octave_idx_type Ni, T extrapval, const T **x, |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
133 const T *v, const T **y, T *vi) |
6702 | 134 { |
135 bool out = false; | |
136 int bit; | |
137 | |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
138 OCTAVE_LOCAL_BUFFER (T, coef, 2*n); |
6702 | 139 OCTAVE_LOCAL_BUFFER (octave_idx_type, index, n); |
140 | |
141 // loop over all points | |
142 for (octave_idx_type m = 0; m < Ni; m++) | |
143 { | |
144 // loop over all dimensions | |
145 for (int i = 0; i < n; i++) | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
146 { |
6702 | 147 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
|
148 out = index[i] == -1; |
6702 | 149 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
150 if (out) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
151 break; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
152 else |
6702 | 153 { |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
154 octave_idx_type j = index[i]; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
155 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
|
156 coef[2*i] = 1 - coef[2*i+1]; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
157 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
158 } |
6702 | 159 |
160 if (out) | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
161 vi[m] = extrapval; |
6702 | 162 else |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
163 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
164 vi[m] = 0; |
6702 | 165 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
166 // 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
|
167 for (int i = 0; i < (1 << n); i++) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
168 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
169 T c = 1; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
170 octave_idx_type l = 0; |
6702 | 171 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
172 // loop over all dimensions |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
173 for (int j = 0; j < n; j++) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
174 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
175 // test if the jth bit in i is set |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
176 bit = i >> j & 1; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
177 l += scale[j] * (index[j] + bit); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
178 c *= coef[2*j+bit]; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
179 } |
6702 | 180 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
181 vi[m] += c * v[l]; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
182 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
183 } |
6702 | 184 } |
185 } | |
186 | |
21139
538b57866b90
consistently use "typename" intead of "class" in template declarations
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
187 template <typename T, typename M> |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
188 octave_value |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
189 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
|
190 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
191 octave_value retval; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
192 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
193 M Vi = M (Y[0].dims ()); |
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 OCTAVE_LOCAL_BUFFER (const T *, y, n); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
196 OCTAVE_LOCAL_BUFFER (octave_idx_type, size, n); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
197 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
198 for (int i = 0; i < n; i++) |
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 y[i] = Y[i].data (); |
21568
3d60ed163b70
maint: Eliminate bad spacing around '='.
Rik <rik@octave.org>
parents:
21301
diff
changeset
|
201 size[i] = V.dims ()(i); |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
202 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
203 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
204 OCTAVE_LOCAL_BUFFER (const T *, x, n); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
205 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
|
206 |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
207 const T *v = V.data (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
208 T *vi = Vi.fortran_vec (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
209 octave_idx_type Ni = Vi.numel (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
210 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
211 T extrapval = octave_NA; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
212 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
213 // offset in memory of each dimension |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
214 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
215 scale[0] = 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 for (int i = 1; i < n; i++) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
218 scale[i] = scale[i-1] * size[i-1]; |
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 // 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
|
221 // in the ndgrid format. |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
222 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
223 if (! isvector (X[0])) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
224 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
225 for (int i = 0; i < n; i++) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
226 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
227 if (X[i].dims () != V.dims ()) |
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20812
diff
changeset
|
228 error ("interpn: incompatible size of argument number %d", i+1); |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20812
diff
changeset
|
229 |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20812
diff
changeset
|
230 M tmp = M (dim_vector (size[i], 1)); |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
231 |
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20812
diff
changeset
|
232 for (octave_idx_type j = 0; j < size[i]; j++) |
21568
3d60ed163b70
maint: Eliminate bad spacing around '='.
Rik <rik@octave.org>
parents:
21301
diff
changeset
|
233 tmp(j) = X[i](scale[i]*j); |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
234 |
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20812
diff
changeset
|
235 X[i] = tmp; |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
236 } |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
237 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
238 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
239 for (int i = 0; i < n; i++) |
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 if (! isvector (X[i]) && X[i].numel () != size[i]) |
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20812
diff
changeset
|
242 error ("interpn: incompatible size of argument number %d", i+1); |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20812
diff
changeset
|
243 |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20812
diff
changeset
|
244 x[i] = X[i].data (); |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
245 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
246 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
247 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
|
248 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
249 retval = Vi; |
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 return retval; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
252 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
253 |
6945 | 254 // Perform @var{n}-dimensional interpolation. Each element of then |
255 // @var{n}-dimensional array @var{v} represents a value at a location | |
21751
b571fc85953f
maint: Use two spaces after period to indicate sentence break.
Rik <rik@octave.org>
parents:
21724
diff
changeset
|
256 // given by the parameters @var{x1}, @var{x2},...,@var{xn}. The parameters |
6945 | 257 // @var{x1}, @var{x2}, @dots{}, @var{xn} are either @var{n}-dimensional |
23803
90689bdbe048
Use C++11 raw string literals to avoid escaping double quotes.
Rik <rik@octave.org>
parents:
23220
diff
changeset
|
258 // arrays of the same size as the array @var{v} in the "ndgrid" format |
6945 | 259 // or vectors. The parameters @var{y1}, @var{y2}, @dots{}, @var{yn} are |
260 // all @var{n}-dimensional arrays of the same size and represent the | |
261 // points at which the array @var{vi} is interpolated. | |
262 // | |
263 //This function only performs linear interpolation. | |
264 | |
15039
e753177cde93
maint: Move non-dynamically linked functions from DLD-FUNCTIONS/ to corefcn/ directory
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
265 DEFUN (__lin_interpn__, args, , |
21966
112b20240c87
move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents:
21751
diff
changeset
|
266 doc: /* -*- texinfo -*- |
112b20240c87
move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents:
21751
diff
changeset
|
267 @deftypefn {} {@var{vi} =} __lin_interpn__ (@var{x1}, @var{x2}, @dots{}, @var{xn}, @var{v}, @var{y1}, @var{y2}, @dots{}, @var{yn}) |
112b20240c87
move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents:
21751
diff
changeset
|
268 Undocumented internal function. |
112b20240c87
move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents:
21751
diff
changeset
|
269 @end deftypefn */) |
6702 | 270 { |
271 int nargin = args.length (); | |
272 | |
20812
d9ca869ca124
maint: Clean-up more instances of print_usage().
Rik <rik@octave.org>
parents:
20802
diff
changeset
|
273 if (nargin < 2 || nargin % 2 == 0) |
20802
8bb38ba1bad6
eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents:
20678
diff
changeset
|
274 print_usage (); |
6702 | 275 |
20918
6f0bd96f93c0
maint: Use new C++ archetype in more files.
Rik <rik@octave.org>
parents:
20853
diff
changeset
|
276 octave_value retval; |
6f0bd96f93c0
maint: Use new C++ archetype in more files.
Rik <rik@octave.org>
parents:
20853
diff
changeset
|
277 |
6702 | 278 // dimension of the problem |
279 int n = (nargin-1)/2; | |
280 | |
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
|
281 if (args(n).is_single_type ()) |
6702 | 282 { |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
283 OCTAVE_LOCAL_BUFFER (FloatNDArray, X, n); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
284 OCTAVE_LOCAL_BUFFER (FloatNDArray, Y, n); |
6702 | 285 |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
286 const FloatNDArray V = args(n).float_array_value (); |
6702 | 287 |
288 for (int i = 0; i < n; i++) | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
289 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
290 X[i] = args(i).float_array_value (); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
291 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
|
292 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
293 if (Y[0].dims () != Y[i].dims ()) |
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20812
diff
changeset
|
294 error ("interpn: incompatible size of argument number %d", n+i+2); |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
295 } |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
296 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
297 retval = lin_interpn<float, FloatNDArray> (n, X, V, Y); |
6702 | 298 } |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
299 else |
6702 | 300 { |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
301 OCTAVE_LOCAL_BUFFER (NDArray, X, n); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
302 OCTAVE_LOCAL_BUFFER (NDArray, Y, n); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
303 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
304 const NDArray V = args(n).array_value (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
305 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
306 for (int i = 0; i < n; i++) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
307 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
308 X[i] = args(i).array_value (); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
309 Y[i] = args(n+i+1).array_value (); |
6702 | 310 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
311 if (Y[0].dims () != Y[i].dims ()) |
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20812
diff
changeset
|
312 error ("interpn: incompatible size of argument number %d", n+i+2); |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9003
diff
changeset
|
313 } |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
314 |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
315 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
|
316 } |
6702 | 317 |
318 return retval; | |
319 } | |
12805
3641167e5b75
codesprint: *.cc helper functions do not need tests
Rik <octave@nomad.inbox5.com>
parents:
11586
diff
changeset
|
320 |
3641167e5b75
codesprint: *.cc helper functions do not need tests
Rik <octave@nomad.inbox5.com>
parents:
11586
diff
changeset
|
321 /* |
3641167e5b75
codesprint: *.cc helper functions do not need tests
Rik <octave@nomad.inbox5.com>
parents:
11586
diff
changeset
|
322 ## 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
|
323 %!assert (1) |
3641167e5b75
codesprint: *.cc helper functions do not need tests
Rik <octave@nomad.inbox5.com>
parents:
11586
diff
changeset
|
324 */ |