Mercurial > octave-nkf

comparison src/DLD-FUNCTIONS/__lin_interpn__.cc @ 10154:40dfc0c99116

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
DLD-FUNCTIONS/*.cc: untabify
author John W. Eaton <jwe@octave.org>
date Wed, 20 Jan 2010 17:33:41 -0500
parents 0631d397fbe0
children fd0a3ac60b0e
comparison
equal deleted inserted replaced
10153:2c28f9d0360f 10154:40dfc0c99116
52 if (x[0] < x[n-1]) 52 if (x[0] < x[n-1])
53 { 53 {
54 // increasing x 54 // increasing x
55 55
56 if (y > x[n-1] || y < x[0]) 56 if (y > x[n-1] || y < x[0])
57 return -1; 57 return -1;
58 58
59 #ifdef EXHAUSTIF 59 #ifdef EXHAUSTIF
60 for (j = 0; j < n - 1; j++) 60 for (j = 0; j < n - 1; j++)
61 { 61 {
62 if (x[j] <= y && y <= x[j+1]) 62 if (x[j] <= y && y <= x[j+1])
63 return j; 63 return j;
64 } 64 }
65 #else 65 #else
66 octave_idx_type j0 = 0; 66 octave_idx_type j0 = 0;
67 octave_idx_type j1 = n - 1; 67 octave_idx_type j1 = n - 1;
68 68
69 while (true) 69 while (true)
70 { 70 {
71 j = (j0+j1)/2; 71 j = (j0+j1)/2;
72 72
73 if (y <= x[j+1]) 73 if (y <= x[j+1])
74 { 74 {
75 if (x[j] <= y) 75 if (x[j] <= y)
76 return j; 76 return j;
77 77
78 j1 = j; 78 j1 = j;
79 } 79 }
80 80
81 if (x[j] <= y) 81 if (x[j] <= y)
82 j0 = j; 82 j0 = j;
83 } 83 }
84 #endif 84 #endif
85 } 85 }
86 else 86 else
87 { 87 {
88 // decreasing x 88 // decreasing x
89 // previous code with x -> -x and y -> -y 89 // previous code with x -> -x and y -> -y
90 90
91 if (y > x[0] || y < x[n-1]) 91 if (y > x[0] || y < x[n-1])
92 return -1; 92 return -1;
93 93
94 #ifdef EXHAUSTIF 94 #ifdef EXHAUSTIF
95 for (j = 0; j < n - 1; j++) 95 for (j = 0; j < n - 1; j++)
96 { 96 {
97 if (x[j+1] <= y && y <= x[j]) 97 if (x[j+1] <= y && y <= x[j])
98 return j; 98 return j;
99 } 99 }
100 #else 100 #else
101 octave_idx_type j0 = 0; 101 octave_idx_type j0 = 0;
102 octave_idx_type j1 = n - 1; 102 octave_idx_type j1 = n - 1;
103 103
104 while (true) 104 while (true)
105 { 105 {
106 j = (j0+j1)/2; 106 j = (j0+j1)/2;
107 107
108 if (y >= x[j+1]) 108 if (y >= x[j+1])
109 { 109 {
110 if (x[j] >= y) 110 if (x[j] >= y)
111 return j; 111 return j;
112 112
113 j1 = j; 113 j1 = j;
114 } 114 }
115 115
116 if (x[j] >= y) 116 if (x[j] >= y)
117 j0 = j; 117 j0 = j;
118 } 118 }
119 #endif 119 #endif
120 } 120 }
121 } 121 }
122 122
123 // n-dimensional linear interpolation 123 // n-dimensional linear interpolation
124 124
125 template <class T> 125 template <class T>
126 void 126 void
127 lin_interpn (int n, const octave_idx_type *size, const octave_idx_type *scale, 127 lin_interpn (int n, const octave_idx_type *size, const octave_idx_type *scale,
128 octave_idx_type Ni, T extrapval, const T **x, 128 octave_idx_type Ni, T extrapval, const T **x,
129 const T *v, const T **y, T *vi) 129 const T *v, const T **y, T *vi)
130 { 130 {
131 bool out = false; 131 bool out = false;
132 int bit; 132 int bit;
133 133
134 OCTAVE_LOCAL_BUFFER (T, coef, 2*n); 134 OCTAVE_LOCAL_BUFFER (T, coef, 2*n);
137 // loop over all points 137 // loop over all points
138 for (octave_idx_type m = 0; m < Ni; m++) 138 for (octave_idx_type m = 0; m < Ni; m++)
139 { 139 {
140 // loop over all dimensions 140 // loop over all dimensions
141 for (int i = 0; i < n; i++) 141 for (int i = 0; i < n; i++)
142 { 142 {
143 index[i] = lookup (x[i], size[i], y[i][m]); 143 index[i] = lookup (x[i], size[i], y[i][m]);
144 out = index[i] == -1; 144 out = index[i] == -1;
145 145
146 if (out) 146 if (out)
147 break; 147 break;
148 else 148 else
149 { 149 {
150 octave_idx_type j = index[i]; 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]); 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]; 152 coef[2*i] = 1 - coef[2*i+1];
153 } 153 }
154 } 154 }
155 155
156 156
157 if (out) 157 if (out)
158 vi[m] = extrapval; 158 vi[m] = extrapval;
159 else 159 else
160 { 160 {
161 vi[m] = 0; 161 vi[m] = 0;
162 162
163 // loop over all corners of hypercube (1<<n = 2^n) 163 // loop over all corners of hypercube (1<<n = 2^n)
164 for (int i = 0; i < (1 << n); i++) 164 for (int i = 0; i < (1 << n); i++)
165 { 165 {
166 T c = 1; 166 T c = 1;
167 octave_idx_type l = 0; 167 octave_idx_type l = 0;
168 168
169 // loop over all dimensions 169 // loop over all dimensions
170 for (int j = 0; j < n; j++) 170 for (int j = 0; j < n; j++)
171 { 171 {
172 // test if the jth bit in i is set 172 // test if the jth bit in i is set
173 bit = i >> j & 1; 173 bit = i >> j & 1;
174 l += scale[j] * (index[j] + bit); 174 l += scale[j] * (index[j] + bit);
175 c *= coef[2*j+bit]; 175 c *= coef[2*j+bit];
176 } 176 }
177 177
178 vi[m] += c * v[l]; 178 vi[m] += c * v[l];
179 } 179 }
180 } 180 }
181 } 181 }
182 } 182 }
183 183
184 template <class T, class M> 184 template <class T, class M>
185 octave_value 185 octave_value
218 // in the ndgrid format. 218 // in the ndgrid format.
219 219
220 if (! isvector (X[0])) 220 if (! isvector (X[0]))
221 { 221 {
222 for (int i = 0; i < n; i++) 222 for (int i = 0; i < n; i++)
223 { 223 {
224 if (X[i].dims () != V.dims ()) 224 if (X[i].dims () != V.dims ())
225 { 225 {
226 error ("interpn: incompatible size of argument number %d", i+1); 226 error ("interpn: incompatible size of argument number %d", i+1);
227 return retval; 227 return retval;
228 } 228 }
229 else 229 else
230 { 230 {
231 M tmp = M (dim_vector (size[i], 1)); 231 M tmp = M (dim_vector (size[i], 1));
232 232
233 for (octave_idx_type j = 0; j < size[i]; j++) 233 for (octave_idx_type j = 0; j < size[i]; j++)
234 tmp(j) = X[i](scale[i]*j); 234 tmp(j) = X[i](scale[i]*j);
235 235
236 X[i] = tmp; 236 X[i] = tmp;
237 } 237 }
238 } 238 }
239 } 239 }
240 240
241 for (int i = 0; i < n; i++) 241 for (int i = 0; i < n; i++)
242 { 242 {
243 if (! isvector (X[i]) && X[i].numel () != size[i]) 243 if (! isvector (X[i]) && X[i].numel () != size[i])
244 { 244 {
245 error ("interpn: incompatible size of argument number %d", i+1); 245 error ("interpn: incompatible size of argument number %d", i+1);
246 return retval; 246 return retval;
247 } 247 }
248 else 248 else
249 x[i] = X[i].data (); 249 x[i] = X[i].data ();
250 } 250 }
251 251
252 lin_interpn (n, size, scale, Ni, extrapval, x, v, y, vi); 252 lin_interpn (n, size, scale, Ni, extrapval, x, v, y, vi);
253 253
254 retval = Vi; 254 retval = Vi;
292 OCTAVE_LOCAL_BUFFER (FloatNDArray, Y, n); 292 OCTAVE_LOCAL_BUFFER (FloatNDArray, Y, n);
293 293
294 const FloatNDArray V = args(n).float_array_value (); 294 const FloatNDArray V = args(n).float_array_value ();
295 295
296 if (error_state) 296 if (error_state)
297 { 297 {
298 print_usage (); 298 print_usage ();
299 return retval; 299 return retval;
300 } 300 }
301 301
302 for (int i = 0; i < n; i++) 302 for (int i = 0; i < n; i++)
303 { 303 {
304 X[i] = args(i).float_array_value (); 304 X[i] = args(i).float_array_value ();
305 Y[i] = args(n+i+1).float_array_value (); 305 Y[i] = args(n+i+1).float_array_value ();
306 306
307 if (error_state) 307 if (error_state)
308 { 308 {
309 print_usage (); 309 print_usage ();
310 return retval; 310 return retval;
311 } 311 }
312 312
313 if (Y[0].dims () != Y[i].dims ()) 313 if (Y[0].dims () != Y[i].dims ())
314 { 314 {
315 error ("interpn: incompatible size of argument number %d", n+i+2); 315 error ("interpn: incompatible size of argument number %d", n+i+2);
316 return retval; 316 return retval;
317 } 317 }
318 } 318 }
319 319
320 retval = lin_interpn<float, FloatNDArray> (n, X, V, Y); 320 retval = lin_interpn<float, FloatNDArray> (n, X, V, Y);
321 } 321 }
322 else 322 else
323 { 323 {
325 OCTAVE_LOCAL_BUFFER (NDArray, Y, n); 325 OCTAVE_LOCAL_BUFFER (NDArray, Y, n);
326 326
327 const NDArray V = args(n).array_value (); 327 const NDArray V = args(n).array_value ();
328 328
329 if (error_state) 329 if (error_state)
330 { 330 {
331 print_usage (); 331 print_usage ();
332 return retval; 332 return retval;
333 } 333 }
334 334
335 for (int i = 0; i < n; i++) 335 for (int i = 0; i < n; i++)
336 { 336 {
337 X[i] = args(i).array_value (); 337 X[i] = args(i).array_value ();
338 Y[i] = args(n+i+1).array_value (); 338 Y[i] = args(n+i+1).array_value ();
339 339
340 if (error_state) 340 if (error_state)
341 { 341 {
342 print_usage (); 342 print_usage ();
343 return retval; 343 return retval;
344 } 344 }
345 345
346 if (Y[0].dims () != Y[i].dims ()) 346 if (Y[0].dims () != Y[i].dims ())
347 { 347 {
348 error ("interpn: incompatible size of argument number %d", n+i+2); 348 error ("interpn: incompatible size of argument number %d", n+i+2);
349 return retval; 349 return retval;
350 } 350 }
351 } 351 }
352 352
353 retval = lin_interpn<double, NDArray> (n, X, V, Y); 353 retval = lin_interpn<double, NDArray> (n, X, V, Y);
354 } 354 }
355 355
356 return retval; 356 return retval;