|
57 | 57 | (rarray) += (results)->itemsize;\ |
58 | 58 | }) |
59 | 59 |
|
60 | | -// The mean could be calculated by simply dividing the sum by |
61 | | -// the number of elements, but that method is numerically unstable |
62 | | -#define RUN_MEAN1(type, array, rarray, ss)\ |
63 | | -({\ |
64 | | - mp_float_t M = 0.0;\ |
65 | | - for(size_t i=0; i < (ss).shape[0]; i++) {\ |
66 | | - mp_float_t value = (mp_float_t)(*(type *)(array));\ |
67 | | - M = M + (value - M) / (mp_float_t)(i+1);\ |
68 | | - (array) += (ss).strides[0];\ |
69 | | - }\ |
70 | | - *(rarray)++ = M;\ |
71 | | -}) |
72 | | - |
73 | 60 | // Instead of the straightforward implementation of the definition, |
74 | 61 | // we take the numerically stable Welford algorithm here |
75 | 62 | // https://www.johndcook.com/blog/2008/09/26/comparing-three-methods-of-computing-standard-deviation/ |
76 | | -#define RUN_STD1(type, array, rarray, ss, div)\ |
77 | | -({\ |
78 | | - mp_float_t M = 0.0, m = 0.0, S = 0.0;\ |
79 | | - for(size_t i=0; i < (ss).shape[0]; i++) {\ |
80 | | - mp_float_t value = (mp_float_t)(*(type *)(array));\ |
81 | | - m = M + (value - M) / (mp_float_t)(i+1);\ |
82 | | - S = S + (value - M) * (value - m);\ |
83 | | - M = m;\ |
84 | | - (array) += (ss).strides[0];\ |
85 | | - }\ |
86 | | - *(rarray)++ = MICROPY_FLOAT_C_FUN(sqrt)(S / (div));\ |
87 | | -}) |
88 | | - |
89 | 63 | #define RUN_MEAN_STD1(type, array, rarray, ss, div, isStd)\ |
90 | 64 | ({\ |
91 | 65 | mp_float_t M = 0.0, m = 0.0, S = 0.0;\ |
|
193 | 167 | RUN_SUM1(type, (array), (results), (rarray), (ss));\ |
194 | 168 | } while(0) |
195 | 169 |
|
196 | | -#define RUN_MEAN(type, array, rarray, ss) do {\ |
197 | | - RUN_MEAN1(type, (array), (rarray), (ss));\ |
198 | | -} while(0) |
199 | | - |
200 | | -#define RUN_STD(type, array, rarray, ss, div) do {\ |
201 | | - RUN_STD1(type, (array), (results), (rarray), (ss), (div));\ |
202 | | -} while(0) |
203 | | - |
204 | 170 | #define RUN_MEAN_STD(type, array, rarray, ss, div, isStd) do {\ |
205 | 171 | RUN_MEAN_STD1(type, (array), (rarray), (ss), (div), (isStd));\ |
206 | 172 | } while(0) |
|
234 | 200 | } while(l < (ss).shape[ULAB_MAX_DIMS - 1]);\ |
235 | 201 | } while(0) |
236 | 202 |
|
237 | | -#define RUN_MEAN(type, array, rarray, ss) do {\ |
238 | | - size_t l = 0;\ |
239 | | - do {\ |
240 | | - RUN_MEAN1(type, (array), (rarray), (ss));\ |
241 | | - (array) -= (ss).strides[0] * (ss).shape[0];\ |
242 | | - (array) += (ss).strides[ULAB_MAX_DIMS - 1];\ |
243 | | - l++;\ |
244 | | - } while(l < (ss).shape[ULAB_MAX_DIMS - 1]);\ |
245 | | -} while(0) |
246 | | - |
247 | | -#define RUN_STD(type, array, rarray, ss, div) do {\ |
248 | | - size_t l = 0;\ |
249 | | - do {\ |
250 | | - RUN_STD1(type, (array), (rarray), (ss), (div));\ |
251 | | - (array) -= (ss).strides[0] * (ss).shape[0];\ |
252 | | - (array) += (ss).strides[ULAB_MAX_DIMS - 1];\ |
253 | | - l++;\ |
254 | | - } while(l < (ss).shape[ULAB_MAX_DIMS - 1]);\ |
255 | | -} while(0) |
256 | | - |
257 | 203 | #define RUN_MEAN_STD(type, array, rarray, ss, div, isStd) do {\ |
258 | 204 | size_t l = 0;\ |
259 | 205 | do {\ |
|
325 | 271 | } while(k < (ss).shape[ULAB_MAX_DIMS - 2]);\ |
326 | 272 | } while(0) |
327 | 273 |
|
328 | | -#define RUN_MEAN(type, array, rarray, ss) do {\ |
329 | | - size_t k = 0;\ |
330 | | - do {\ |
331 | | - size_t l = 0;\ |
332 | | - do {\ |
333 | | - RUN_MEAN1(type, (array), (rarray), (ss));\ |
334 | | - (array) -= (ss).strides[0] * (ss).shape[0];\ |
335 | | - (array) += (ss).strides[ULAB_MAX_DIMS - 1];\ |
336 | | - l++;\ |
337 | | - } while(l < (ss).shape[ULAB_MAX_DIMS - 1]);\ |
338 | | - (array) -= (ss).strides[ULAB_MAX_DIMS - 1] * (ss).shape[ULAB_MAX_DIMS - 1];\ |
339 | | - (array) += (ss).strides[ULAB_MAX_DIMS - 2];\ |
340 | | - k++;\ |
341 | | - } while(k < (ss).shape[ULAB_MAX_DIMS - 2]);\ |
342 | | -} while(0) |
343 | | - |
344 | | -#define RUN_STD(type, array, rarray, ss, div) do {\ |
345 | | - size_t k = 0;\ |
346 | | - do {\ |
347 | | - size_t l = 0;\ |
348 | | - do {\ |
349 | | - RUN_STD1(type, (array), (rarray), (ss), (div));\ |
350 | | - (array) -= (ss).strides[0] * (ss).shape[0];\ |
351 | | - (array) += (ss).strides[ULAB_MAX_DIMS - 1];\ |
352 | | - l++;\ |
353 | | - } while(l < (ss).shape[ULAB_MAX_DIMS - 1]);\ |
354 | | - (array) -= (ss).strides[ULAB_MAX_DIMS - 1] * (ss).shape[ULAB_MAX_DIMS - 1];\ |
355 | | - (array) += (ss).strides[ULAB_MAX_DIMS - 2];\ |
356 | | - k++;\ |
357 | | - } while(k < (ss).shape[ULAB_MAX_DIMS - 2]);\ |
358 | | -} while(0) |
359 | | - |
360 | 274 | #define RUN_MEAN_STD(type, array, rarray, ss, div, isStd) do {\ |
361 | 275 | size_t k = 0;\ |
362 | 276 | do {\ |
|
467 | 381 | } while(j < (ss).shape[ULAB_MAX_DIMS - 3]);\ |
468 | 382 | } while(0) |
469 | 383 |
|
470 | | -#define RUN_MEAN(type, array, rarray, ss) do {\ |
471 | | - size_t j = 0;\ |
472 | | - do {\ |
473 | | - size_t k = 0;\ |
474 | | - do {\ |
475 | | - size_t l = 0;\ |
476 | | - do {\ |
477 | | - RUN_MEAN1(type, (array), (rarray), (ss));\ |
478 | | - (array) -= (ss).strides[0] * (ss).shape[0];\ |
479 | | - (array) += (ss).strides[ULAB_MAX_DIMS - 1];\ |
480 | | - l++;\ |
481 | | - } while(l < (ss).shape[ULAB_MAX_DIMS - 1]);\ |
482 | | - (array) -= (ss).strides[ULAB_MAX_DIMS - 1] * (ss).shape[ULAB_MAX_DIMS - 1];\ |
483 | | - (array) += (ss).strides[ULAB_MAX_DIMS - 2];\ |
484 | | - k++;\ |
485 | | - } while(k < (ss).shape[ULAB_MAX_DIMS - 2]);\ |
486 | | - (array) -= (ss).strides[ULAB_MAX_DIMS - 2] * (ss).shape[ULAB_MAX_DIMS - 2];\ |
487 | | - (array) += (ss).strides[ULAB_MAX_DIMS - 3];\ |
488 | | - j++;\ |
489 | | - } while(j < (ss).shape[ULAB_MAX_DIMS - 3]);\ |
490 | | -} while(0) |
491 | | - |
492 | | -#define RUN_STD(type, array, rarray, ss, div) do {\ |
493 | | - size_t j = 0;\ |
494 | | - do {\ |
495 | | - size_t k = 0;\ |
496 | | - do {\ |
497 | | - size_t l = 0;\ |
498 | | - do {\ |
499 | | - RUN_STD1(type, (array), (rarray), (ss), (div));\ |
500 | | - (array) -= (ss).strides[0] * (ss).shape[0];\ |
501 | | - (array) += (ss).strides[ULAB_MAX_DIMS - 1];\ |
502 | | - l++;\ |
503 | | - } while(l < (ss).shape[ULAB_MAX_DIMS - 1]);\ |
504 | | - (array) -= (ss).strides[ULAB_MAX_DIMS - 1] * (ss).shape[ULAB_MAX_DIMS - 1];\ |
505 | | - (array) += (ss).strides[ULAB_MAX_DIMS - 2];\ |
506 | | - k++;\ |
507 | | - } while(k < (ss).shape[ULAB_MAX_DIMS - 2]);\ |
508 | | - (array) -= (ss).strides[ULAB_MAX_DIMS - 2] * (ss).shape[ULAB_MAX_DIMS - 2];\ |
509 | | - (array) += (ss).strides[ULAB_MAX_DIMS - 3];\ |
510 | | - j++;\ |
511 | | - } while(j < (ss).shape[ULAB_MAX_DIMS - 3]);\ |
512 | | -} while(0) |
513 | | - |
514 | 384 | #define RUN_MEAN_STD(type, array, rarray, ss, div, isStd) do {\ |
515 | 385 | size_t j = 0;\ |
516 | 386 | do {\ |
|
0 commit comments