Mercurial > hg > octave-lyh

comparison src/DLD-FUNCTIONS/__lin_interpn__.cc @ 7789:82be108cc558

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
First attempt at single precision tyeps * * * corrections to qrupdate single precision routines * * * prefer demotion to single over promotion to double * * * Add single precision support to log2 function * * * Trivial PROJECT file update * * * Cache optimized hermitian/transpose methods * * * Add tests for tranpose/hermitian and ChangeLog entry for new transpose code
author David Bateman <dbateman@free.fr>
date Sun, 27 Apr 2008 22:34:17 +0200
parents a938cd7869b2
children 25bc2d31e1bf
comparison
equal deleted inserted replaced
7788:45f5faba05a2 7789:82be108cc558
30 #include "error.h" 30 #include "error.h"
31 #include "oct-obj.h" 31 #include "oct-obj.h"
32 32
33 // equivalent to isvector.m 33 // equivalent to isvector.m
34 34
35 template <class T>
35 bool 36 bool
36 isvector (const NDArray& array) 37 isvector (const T& array)
37 { 38 {
38 const dim_vector dv = array.dims (); 39 const dim_vector dv = array.dims ();
39 return dv.length () == 2 && (dv(0) == 1 || dv(1) == 1); 40 return dv.length () == 2 && (dv(0) == 1 || dv(1) == 1);
40 } 41 }
41 42
42 // lookup a value in a sorted table (lookup.m) 43 // lookup a value in a sorted table (lookup.m)
44 template <class T>
43 octave_idx_type 45 octave_idx_type
44 lookup (const double *x, octave_idx_type n, double y) 46 lookup (const T *x, octave_idx_type n, T y)
45 { 47 {
46 octave_idx_type j; 48 octave_idx_type j;
47 49
48 if (x[0] < x[n-1]) 50 if (x[0] < x[n-1])
49 { 51 {
116 } 118 }
117 } 119 }
118 120
119 // n-dimensional linear interpolation 121 // n-dimensional linear interpolation
120 122
123 template <class T>
121 void 124 void
122 lin_interpn (int n, const octave_idx_type *size, const octave_idx_type *scale, 125 lin_interpn (int n, const octave_idx_type *size, const octave_idx_type *scale,
123 octave_idx_type Ni, double extrapval, const double **x, 126 octave_idx_type Ni, T extrapval, const T **x,
124 const double *v, const double **y, double *vi) 127 const T *v, const T **y, T *vi)
125 { 128 {
126 bool out = false; 129 bool out = false;
127 int bit; 130 int bit;
128 131
129 OCTAVE_LOCAL_BUFFER (double, coef, 2*n); 132 OCTAVE_LOCAL_BUFFER (T, coef, 2*n);
130 OCTAVE_LOCAL_BUFFER (octave_idx_type, index, n); 133 OCTAVE_LOCAL_BUFFER (octave_idx_type, index, n);
131 134
132 // loop over all points 135 // loop over all points
133 for (octave_idx_type m = 0; m < Ni; m++) 136 for (octave_idx_type m = 0; m < Ni; m++)
134 { 137 {
156 vi[m] = 0; 159 vi[m] = 0;
157 160
158 // loop over all corners of hypercube (1<<n = 2^n) 161 // loop over all corners of hypercube (1<<n = 2^n)
159 for (int i = 0; i < (1 << n); i++) 162 for (int i = 0; i < (1 << n); i++)
160 { 163 {
161 double c = 1; 164 T c = 1;
162 octave_idx_type l = 0; 165 octave_idx_type l = 0;
163 166
164 // loop over all dimensions 167 // loop over all dimensions
165 for (int j = 0; j < n; j++) 168 for (int j = 0; j < n; j++)
166 { 169 {
172 175
173 vi[m] += c * v[l]; 176 vi[m] += c * v[l];
174 } 177 }
175 } 178 }
176 } 179 }
180 }
181
182 template <class T, class M>
183 octave_value
184 lin_interpn (int n, M *X, const M V, M *Y)
185 {
186 octave_value retval;
187
188 M Vi = M (Y[0].dims ());
189
190 OCTAVE_LOCAL_BUFFER (const T *, y, n);
191 OCTAVE_LOCAL_BUFFER (octave_idx_type, size, n);
192
193 for (int i = 0; i < n; i++)
194 {
195 y[i] = Y[i].data ();
196 size[i] = V.dims()(i);
197 }
198
199 OCTAVE_LOCAL_BUFFER (const T *, x, n);
200 OCTAVE_LOCAL_BUFFER (octave_idx_type, scale, n);
201
202 const T *v = V.data ();
203 T *vi = Vi.fortran_vec ();
204 octave_idx_type Ni = Vi.numel ();
205
206 T extrapval = octave_NA;
207
208 // offset in memory of each dimension
209
210 scale[0] = 1;
211
212 for (int i = 1; i < n; i++)
213 scale[i] = scale[i-1] * size[i-1];
214
215 // tests if X[0] is a vector, if yes, assume that all elements of X are
216 // in the ndgrid format.
217
218 if (! isvector (X[0]))
219 {
220 for (int i = 0; i < n; i++)
221 {
222 if (X[i].dims () != V.dims ())
223 {
224 error ("interpn: incompatible size of argument number %d", i+1);
225 return retval;
226 }
227 else
228 {
229 M tmp = M (dim_vector (size[i], 1));
230
231 for (octave_idx_type j = 0; j < size[i]; j++)
232 tmp(j) = X[i](scale[i]*j);
233
234 X[i] = tmp;
235 }
236 }
237 }
238
239 for (int i = 0; i < n; i++)
240 {
241 if (! isvector (X[i]) && X[i].numel () != size[i])
242 {
243 error ("interpn: incompatible size of argument number %d", i+1);
244 return retval;
245 }
246 else
247 x[i] = X[i].data ();
248 }
249
250 lin_interpn (n, size, scale, Ni, extrapval, x, v, y, vi);
251
252 retval = Vi;
253
254 return retval;
177 } 255 }
178 256
179 // Perform @var{n}-dimensional interpolation. Each element of then 257 // Perform @var{n}-dimensional interpolation. Each element of then
180 // @var{n}-dimensional array @var{v} represents a value at a location 258 // @var{n}-dimensional array @var{v} represents a value at a location
181 // given by the parameters @var{x1}, @var{x2},...,@var{xn}. The parameters 259 // given by the parameters @var{x1}, @var{x2},...,@var{xn}. The parameters
204 } 282 }
205 283
206 // dimension of the problem 284 // dimension of the problem
207 int n = (nargin-1)/2; 285 int n = (nargin-1)/2;
208 286
209 OCTAVE_LOCAL_BUFFER (NDArray, X, n); 287 if (args(n).is_single_type())
210 OCTAVE_LOCAL_BUFFER (NDArray, Y, n); 288 {
211 289 OCTAVE_LOCAL_BUFFER (FloatNDArray, X, n);
212 OCTAVE_LOCAL_BUFFER (const double *, x, n); 290 OCTAVE_LOCAL_BUFFER (FloatNDArray, Y, n);
213 OCTAVE_LOCAL_BUFFER (const double *, y, n); 291
214 OCTAVE_LOCAL_BUFFER (octave_idx_type, scale, n); 292 const FloatNDArray V = args(n).float_array_value ();
215 OCTAVE_LOCAL_BUFFER (octave_idx_type, size, n);
216
217 const NDArray V = args(n).array_value ();
218 NDArray Vi = NDArray (args(n+1).dims ());
219
220 if (error_state)
221 {
222 print_usage ();
223 return retval;
224 }
225
226 const double *v = V.data ();
227 double *vi = Vi.fortran_vec ();
228 octave_idx_type Ni = Vi.numel ();
229
230 double extrapval = octave_NA;
231
232 for (int i = 0; i < n; i++)
233 {
234 X[i] = args(i).array_value ();
235 Y[i] = args(n+i+1).array_value ();
236 293
237 if (error_state) 294 if (error_state)
238 { 295 {
239 print_usage (); 296 print_usage ();
240 return retval; 297 return retval;
241 } 298 }
242 299
243 y[i] = Y[i].data (); 300 for (int i = 0; i < n; i++)
244 size[i] = V.dims()(i); 301 {
245 302 X[i] = args(i).float_array_value ();
246 if (Y[0].dims () != Y[i].dims ()) 303 Y[i] = args(n+i+1).float_array_value ();
247 { 304
248 error ("interpn: incompatible size of argument number %d", n+i+2); 305 if (error_state)
306 {
307 print_usage ();
308 return retval;
309 }
310
311 if (Y[0].dims () != Y[i].dims ())
312 {
313 error ("interpn: incompatible size of argument number %d", n+i+2);
314 return retval;
315 }
316 }
317
318 retval = lin_interpn<float, FloatNDArray> (n, X, V, Y);
319 }
320 else
321 {
322 OCTAVE_LOCAL_BUFFER (NDArray, X, n);
323 OCTAVE_LOCAL_BUFFER (NDArray, Y, n);
324
325 const NDArray V = args(n).array_value ();
326
327 if (error_state)
328 {
329 print_usage ();
249 return retval; 330 return retval;
250 } 331 }
251 } 332
252
253 // offset in memory of each dimension
254
255 scale[0] = 1;
256
257 for (int i = 1; i < n; i++)
258 scale[i] = scale[i-1] * size[i-1];
259
260 // tests if X[0] is a vector, if yes, assume that all elements of X are
261 // in the ndgrid format.
262
263 if (! isvector (X[0]))
264 {
265 for (int i = 0; i < n; i++) 333 for (int i = 0; i < n; i++)
266 { 334 {
267 if (X[i].dims () != V.dims ()) 335 X[i] = args(i).array_value ();
268 { 336 Y[i] = args(n+i+1).array_value ();
269 error ("interpn: incompatible size of argument number %d", i+1); 337
270 return retval; 338 if (error_state)
271 } 339 {
272 else 340 print_usage ();
273 { 341 return retval;
274 NDArray tmp = NDArray (dim_vector (size[i], 1)); 342 }
275 343
276 for (octave_idx_type j = 0; j < size[i]; j++) 344 if (Y[0].dims () != Y[i].dims ())
277 tmp(j) = X[i](scale[i]*j); 345 {
278 346 error ("interpn: incompatible size of argument number %d", n+i+2);
279 X[i] = tmp; 347 return retval;
280 } 348 }
281 } 349 }
282 } 350
283 351 retval = lin_interpn<double, NDArray> (n, X, V, Y);
284 for (int i = 0; i < n; i++) 352 }
285 {
286 if (! isvector (X[i]) && X[i].numel () != size[i])
287 {
288 error ("interpn: incompatible size of argument number %d", i+1);
289 return retval;
290 }
291 else
292 x[i] = X[i].data ();
293 }
294
295 lin_interpn (n, size, scale, Ni, extrapval, x, v, y, vi);
296
297 retval = Vi;
298 353
299 return retval; 354 return retval;
300 } 355 }