1 #ifndef _NUM2_h_
2 #define _NUM2_h_
3 /* NUM2.h
4 *
5 * Copyright (C) 1997-2021 David Weenink
6 *
7 * This code is free software; you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation; either version 2 of the License, or (at
10 * your option) any later version.
11 *
12 * This code is distributed in the hope that it will be useful, but
13 * WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15 * General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with this work. If not, see <http://www.gnu.org/licenses/>.
19 */
20
21 /*
22 djmw 20020815 GPL header
23 */
24
25 #include <algorithm>
26 #include <limits.h>
27 #include "melder.h"
28 #include "MAT_numerics.h"
29
30 /*
31 A tipical iteration goes like
32 do {
33 previous is ...
34 get value of current
35 while (fabs (current - previous) > fabs (NUMeps * current);
36 Because of floating point arithmatic we have a rounding differences which are large or equal to the
37 machine precision eps = 2.2...e-16. To have an escape if fabs (current - previous) e<= eps we choose
38 NUMeps just a little larger than eps.
39 */
40 #define NUMeps 2.3e-16
41
42 /*
43 only used to append the time info to the name of an object.
44 1.2345678 -> 1_2345678
45 1.2 -> 1_2
46 2 -> 2
47 */
48 char32 * NUMnumber_as_stringWithDotReplacedByUnderscore (double time);
49
50 regexp *NUMregexp_compile (conststring32 regexp);
51 /* Compiles a regular expression to a datastructure used by the regexp engine */
52
53 char32 *strstr_regexp (conststring32 string, conststring32 search_regexp);
54 /*
55 Returns a pointer to the first occurrence in 'string' of the
56 regular expression 'searchRE'. It returns a null pointer if
57 no match is found.
58 */
59
60 autoSTRVEC string32vector_searchAndReplace (constSTRVEC me,
61 conststring32 search, conststring32 replace, integer maximumNumberOfReplaces,
62 integer *nmatches, integer *nstringmatches, bool use_regexp);
63 /*
64 Searches and replaces in string array of strings.
65 If use_regexp != 0, 'search' and 'replace' will be interpreted
66 as regular expressions. Else these strings are taken literally.
67
68 'maximumNumberOfReplaces' is the maximum number of replaces in EACH string
69 in the array of strings (you can replace ALL occurrences by making this
70 number <= 0).
71 The totalnumber of matches found is returned in 'nmatches'.
72 The number of strings with at least one match is returned in
73 'nstringmatches'.
74 */
75
76 void MATprintMatlabForm (constMATVU const& m, conststring32 name);
77 /*
78 Print a matrix in a form that can be used as input for octave/matlab.
79 1 2 3
80 Let A be the matrix: 4 5 6
81 7 8 9
82 The output from MATprintMatlabForm (A, "M") will be
83 M= [1, 2, 3;
84 4, 5, 6;
85 7, 8, 9];
86 */
87
NUMisNonNegative(constVECVU const & vec)88 inline bool NUMisNonNegative (constVECVU const& vec) {
89 for (integer i = 1; i <= vec.size; i ++)
90 if (vec [i] < 0.0)
91 return false;
92 return true;
93 }
NUMisNonNegative(constMATVU const & mat)94 inline bool NUMisNonNegative (constMATVU const& mat) {
95 for (integer irow = 1; irow <= mat.nrow; irow ++)
96 for (integer icol = 1; icol <= mat.ncol; icol ++)
97 if (mat [irow] [icol] < 0.0)
98 return false;
99 return true;
100 }
101
NUMmaxPos(constVECVU const & v)102 inline integer NUMmaxPos (constVECVU const& v) {
103 if (NUMisEmpty (v)) return 0;
104 integer index = 1;
105 double maximum = v [1];
106 for (integer i = 2; i <= v.size; i ++) {
107 if (v [i] > maximum) {
108 maximum = v [i];
109 index = i;
110 }
111 }
112 return index;
113 }
114
NUMmaxPos(constINTVECVU const & v)115 inline integer NUMmaxPos (constINTVECVU const& v) {
116 if (NUMisEmpty (v)) return 0;
117 integer index = 1;
118 integer maximum = v [1];
119 for (integer i = 2; i <= v.size; i ++) {
120 if (v [i] > maximum) {
121 maximum = v [i];
122 index = i;
123 }
124 }
125 return index;
126 }
127
NUMminPos(constVECVU const & v)128 inline integer NUMminPos (constVECVU const& v) {
129 if (NUMisEmpty (v)) return 0;
130 integer index = 1;
131 double minimum = v [1];
132 for (integer i = 2; i <= v.size; i ++) {
133 if (v [i] < minimum) {
134 minimum = v [i];
135 index = i;
136 }
137 }
138 return index;
139 }
140
NUMminPos(constINTVECVU const & v)141 inline integer NUMminPos (constINTVECVU const& v) {
142 if (NUMisEmpty (v))
143 return 0;
144 integer index = 1;
145 integer minimum = v [1];
146 for (integer i = 2; i <= v.size; i ++) {
147 if (v [i] < minimum) {
148 minimum = v [i];
149 index = i;
150 }
151 }
152 return index;
153 }
154
NUMextrema(constVECVU const & x,double * out_minimum,double * out_maximum)155 inline void NUMextrema (constVECVU const& x, double *out_minimum, double *out_maximum) {
156 if (out_minimum)
157 *out_minimum = NUMmin (x);
158 if (out_maximum)
159 *out_maximum = NUMmax (x);
160 }
161
162 /*
163 Clip array values.
164 c [i] = c [i] < min ? min : (c [i] > max ? max : c [i])
165 */
VECclip_inplace(double min,VECVU const & x,double max)166 inline void VECclip_inplace (double min, VECVU const& x, double max) {
167 for (integer i = 1; i <= x.size; i ++)
168 Melder_clip (min, & x [i], max);
169 }
170
VECabs(VECVU const & result,constVECVU const & v)171 inline void VECabs (VECVU const& result, constVECVU const& v) {
172 Melder_assert (result.size == v.size);
173 for (integer i = 1; i <= result.size; i ++)
174 result [i] = fabs (v [i]);
175 }
176
newVECabs(constVECVU const & v)177 inline autoVEC newVECabs (constVECVU const& v) {
178 autoVEC result = raw_VEC (v.size);
179 VECabs (result.get(), v);
180 return result;
181 }
182
VECabs_inplace(VECVU const & v)183 inline void VECabs_inplace (VECVU const& v) {
184 for (integer i = 1; i <= v.size; i ++)
185 v [i] = fabs (v [i]);
186 }
187
NUMhasZeroElement(constMATVU const m)188 inline bool NUMhasZeroElement (constMATVU const m) {
189 for (integer irow = 1; irow <= m.nrow; irow ++)
190 for (integer icol = 1; icol <= m.ncol; icol++)
191 if (m [irow] [icol] == 0.0)
192 return true;
193 return false;
194 }
195
NUMcountNumberOfNonZeroElements(constVECVU const & v)196 inline integer NUMcountNumberOfNonZeroElements (constVECVU const& v) {
197 integer count = 0;
198 for (integer i = 1; i <= v.size; i ++)
199 if (v [i] != 0.0)
200 ++ count;
201 return count;
202 }
203
NUMmul(constVECVU const & x,constMATVU const & m,constVECVU const & y)204 inline double NUMmul (constVECVU const& x, constMATVU const& m, constVECVU const& y) { // x'. M . y
205 Melder_assert (x.size == m.nrow);
206 Melder_assert (y.size == m.ncol);
207 longdouble result = 0.0;
208 for (integer k = 1; k <= x.size; k ++)
209 result += x [k] * NUMinner (m.row (k), y);
210 return (double) result;
211 }
212
VECnorm_rows(constMATVU const & x,double power)213 inline autoVEC VECnorm_rows (constMATVU const& x, double power) {
214 autoVEC norm = raw_VEC (x.nrow);
215 for (integer irow = 1; irow <= norm.size; irow ++)
216 norm [irow] = NUMnorm (x.row (irow), power);
217 return norm;
218 }
219
VECnormalize_inplace(VECVU const & vec,double power,double newNorm)220 inline void VECnormalize_inplace (VECVU const& vec, double power, double newNorm) {
221 Melder_assert (newNorm > 0.0);
222 double oldnorm = NUMnorm (vec, power);
223 if (oldnorm > 0.0)
224 vec *= newNorm / oldnorm;
225 }
226
MATnormalize_inplace(MATVU const & mat,double power,double newNorm)227 inline void MATnormalize_inplace (MATVU const& mat, double power, double newNorm) {
228 Melder_assert (newNorm > 0.0);
229 double oldnorm = NUMnorm (mat, power);
230 if (oldnorm > 0.0)
231 mat *= newNorm / oldnorm;
232 }
233
MATnormalizeRows_inplace(MATVU const & a,double power,double norm)234 inline void MATnormalizeRows_inplace (MATVU const& a, double power, double norm) {
235 Melder_assert (norm > 0.0);
236 for (integer irow = 1; irow <= a.nrow; irow ++)
237 VECnormalize_inplace (a.row (irow), power, norm);
238 }
239
MATnormalizeColumns_inplace(MATVU const & a,double power,double norm)240 inline void MATnormalizeColumns_inplace (MATVU const& a, double power, double norm) {
241 MATnormalizeRows_inplace (a.transpose(), power, norm);
242 }
243 /*
244 Scale a [.] [j] such that sqrt (Sum(a [i] [j]^2, i=1..nPoints)) = norm.
245 */
246
247 void VECsmoothByMovingAverage_preallocated (VECVU const& out, constVECVU const& in, integer window);
248
249 autoMAT MATcovarianceFromColumnCentredMatrix (constMATVU const& x, integer ndf);
250 /*
251 Calculate covariance matrix(ncols x ncols) from data matrix (nrows x ncols);
252 The matrix x must be column centered.
253 covar [i] [j] = sum (k=1..nrows, x [i] [k] * x [k] [j]) / (nrows - ndf)
254 */
255
256 void MATmtm_weighRows (MATVU const& result, constMATVU const& data, constVECVU const& rowWeights);
257
newMATmtm_weighRows(constMATVU const & data,constVECVU const & rowWeights)258 inline autoMAT newMATmtm_weighRows (constMATVU const& data, constVECVU const& rowWeights) {
259 autoMAT result = raw_MAT (data.ncol, data.ncol);
260 MATmtm_weighRows (result.get(), data, rowWeights);
261 return result;
262 }
263
264 double NUMmultivariateKurtosis (constMATVU const& x, integer method);
265 /*
266 calculate multivariate kurtosis.
267 method = 1 : Schott (2001), J. of Statistical planning and Inference 94, 25-36.
268 */
269
270 void NUMmad (constVEC x, double *inout_location, bool wantlocation, double *out_mad, VEC const& workSpace);
271 /*
272 Computes the median absolute deviation, i.e., the median of the
273 absolute deviations from the median, and adjust by a factor for
274 asymptotically normal consistency, i.e. the returned value is 1.4826*mad which
275 makes the returned value "equal" to the standard deviation if the data is normally distributed.
276 You either GIVE the median location (if wantlocation = false) or it
277 will be calculated (if wantlocation = true);
278
279 Precondition: workSpace.size >= x.size
280 */
281
282 void NUMstatistics_huber (constVEC x, double *inout_location, bool wantlocation, double *inout_scale, bool wantscale, double k_stdev, double tol, integer maximumNumberOfiterations, VEC const& workSpace);
283 /*
284 Finds the Huber M-estimator for location with scale specified,
285 scale with location specified, or both if neither is specified.
286 k_stdev Winsorizes at `k_stdev' standard deviations.
287
288 Precondition: workSpace.size >= x.size
289 */
290
291 autoVEC newVECmonotoneRegression (constVEC x);
292 /*
293 Find numbers xs [1..n] that have a monotone relationship with
294 the numbers in x [1..n].
295 The xs [i] will be ascending.
296 */
297
298
299 /* NUMsort2:
300 NUMsort2 uses heapsort to sort the second array in parallel with the first one.
301
302 Algorithm follows p. 145 and 642 in:
303 Donald E. Knuth (1998): The art of computer programming. Third edition. Vol. 3: sorting and searching.
304 Boston: Addison-Wesley, printed may 2002.
305 Modification: there is no distinction between record and key and
306 Floyd's optimization (page 642) is used.
307 Sorts (inplace) an array a [1..n] into ascending order using the Heapsort algorithm,
308 while making the corresponding rearrangement of the companion
309 array b [1..n]. A characteristic of heapsort is that it does not conserve
310 the order of equals: e.g., the array 3,1,1,2 will be sorted as 1,1,2,3 and
311 it may occur that the first 1 after sorting came from position 3 and the second
312 1 came from position 2.
313 */
314 template<typename T1, typename T2>
NUMsortTogether(vector<T1> a,vector<T2> b)315 void NUMsortTogether (vector<T1> a, vector<T2> b) {
316 Melder_assert (a.size == b.size);
317 T1 k, min;
318 T2 kb, min2;
319 if (a.size < 2) return; /* Already sorted. */
320 if (a.size == 2) {
321 if (a [1] > a [2]) {
322 min = a [2];
323 a [2] = a [1];
324 a [1] = min;
325 min2 = b [2];
326 b [2] = b [1];
327 b [1] = min2;
328 }
329 return;
330 }
331 if (a.size <= 12) {
332 for (integer i = 1; i < a.size; i ++) {
333 min = a [i];
334 integer imin = i;
335 for (integer j = i + 1; j <= a.size; j ++)
336 if (a [j] < min) {
337 min = a [j];
338 imin = j;
339 }
340 a [imin] = a [i];
341 a [i] = min;
342 min2 = b [imin];
343 b [imin] = b [i];
344 b [i] = min2;
345 }
346 return;
347 }
348 /* H1 */
349 integer l = (a.size >> 1) + 1;
350 integer r = a.size;
351 for (;;) {
352 if (l > 1) {
353 l --;
354 k = a [l];
355 kb = b [l];
356 } else /* l == 1 */ {
357 k = a [r];
358 kb = b [r];
359 a [r] = a [1];
360 b [r] = b [1];
361 r --;
362 if (r == 1) {
363 a [1] = k;
364 b [1] = kb;
365 return;
366 }
367 }
368 /* H3 */
369 integer i, j = l;
370 for (;;) { /* H4 */
371 i = j;
372 j = j << 1;
373 if (j > r) break;
374 if (j < r && a [j] < a [j + 1]) j ++; /* H5 */
375 /* if (k >= a [j]) break; H6 */
376 a [i] = a [j];
377 b [i] = b [j]; /* H7 */
378 }
379 /* a [i] = k; b [i] = kb; H8 */
380 for (;;) { /*H8' */
381 j = i;
382 i = j >> 1;
383 /* H9' */
384 if (j == l || k <= a [i]) {
385 a [j] = k;
386 b [j] = kb;
387 break;
388 }
389 a [j] = a [i];
390 b [j] = b [i];
391 }
392 }
393 }
394
395 void VECsort3_inplace (VEC const& a, INTVEC const& iv1, INTVEC const& iv2, bool descending); // TODO template
396 /* Sort a together with iv1 and iv2 */
397
398 void INTVECindex (INTVEC const& target, constVEC const& a);
399 void INTVECindex (INTVEC const& target, constSTRVEC const& s);
400
newINTVECindex(constVEC const & a)401 inline autoINTVEC newINTVECindex (constVEC const& a) {
402 autoINTVEC result = raw_INTVEC (a.size);
403 INTVECindex (result.get(), a);
404 return result;
405 }
406
newINTVECindex(constSTRVEC const & s)407 inline autoINTVEC newINTVECindex (constSTRVEC const& s) {
408 autoINTVEC result = raw_INTVEC (s.size);
409 INTVECindex (result.get(), s);
410 return result;
411 }
412
413 void MATrankColumns (MAT m, integer cb, integer ce);
414
415 /* rank:
416 * Replace content of sorted array by rank number, including midranking of ties.
417 * E.g. The elements {10, 20.1, 20.1, 20.1, 20.1, 30} in array a will be replaced
418 * by {1, 3.5, 3.5, 3.5, 3.5, 4}, respectively. *
419 */
420
VECrankSorted(VECVU const & a)421 inline void VECrankSorted (VECVU const& a) {
422 integer jt, j = 1;
423 while (j < a.size) {
424 for (jt = j + 1; jt <= a.size && a [jt] == a [j]; jt ++) {}
425 double rank = (j + jt - 1) * 0.5;
426 for (integer i = j; i <= jt - 1; i ++)
427 a [i] = rank;
428 j = jt;
429 }
430 if (j == a.size)
431 a [a.size] = a.size;
432 }
433
434 autoMAT newMATlowerCholeslyInverse_fromLowerCholesky (constMAT const& m);
435
436 void MATlowerCholesky_inplace (MAT a, double *out_lnd);
437
438 autoMAT newMATlowerCholesky (constMATVU const& a, double *out_lnd);
439
440 void MATlowerCholeskyInverse_inplace (MAT a, double *out_lnd);
441
newMATlowerCholeskyInverse(constMAT const & a)442 inline autoMAT newMATlowerCholeskyInverse (constMAT const& a) {
443 autoMAT result = copy_MAT (a);
444 MATlowerCholeskyInverse_inplace (result.get(), nullptr);
445 return result;
446 }
447 /*
448 Calculates L^-1, where A = L.L' is a symmetric positive definite matrix
449 and ln(determinant). L^-1 in lower, leave upper part intact.
450 */
451
452 autoMAT newMATinverse_fromLowerCholeskyInverse (constMAT m);
453 /*
454 Return the complete matrix inverse when only the inverse of the lower Cholesky part is given.
455 Input m is a square matrix, in the lower part is the inverse of the lower Cholesky part as calculated by NUMlowerCholeskyInverse.
456 */
457
458 double NUMdeterminant_fromSymmetricMatrix (constMAT m);
459 /*
460 ln(determinant) of a symmetric p.s.d. matrix
461 */
462
463 double NUMmahalanobisDistanceSquared (constMAT lowerInverse, constVEC v, constVEC m);
464 /*
465 Calculates squared Mahalanobis distance: (v-m)'S^-1(v-m).
466 Input matrix (li) is the inverse L^-1 of the Cholesky decomposition S = L.L'
467 as calculated by NUMlowerCholeskyInverse or 1-row for a diagonal matrix (nr =1)
468 Mahalanobis distance calculation. S = L.L' -> S**-1 = L**-1' . L**-1
469 (x-m)'S**-1 (x-m) = (x-m)'L**-1' . L**-1. (x-m) =
470 (L**-1.(x-m))' . (L**-1.(x-m))
471 */
472
473 double NUMtrace (constMATVU const& a);
474 double NUMtrace2 (constMATVU const& x, constMATVU const& y);
475 /*
476 Calculates the trace from a product matrix x*y
477 */
478
479 void MATprojectColumnsOnEigenspace_preallocated (MAT projection, constMATVU const& data, constMATVU const& eigenvectors);
480 /* Input:
481 data [dimension, numberOfColumns]
482 contains the column vectors to be projected on the eigenspace.
483 eigenvectors [numberOfEigenvectors] [dimension]
484 the eigenvectors stored as rowvectors
485 Input/Output
486 projection [numberOfEigenvectors, numberOfColumns]
487 the projected vectors from 'data'
488
489 Project the columnvectors in matrix 'data' along the 'numberOfEigenvectors' eigenvectors into the matrix 'projection'.
490 */
491
492
493 double VECdominantEigenvector_inplace (VEC inout_q, constMAT m, double tolerance);
494 /*
495 Determines the first dominant eigenvector from a square GENERAL matrix m.
496 Besides the matrix m, a first guess for the eigenvector q must
497 be supplied (e.g. 1,0,...,0) and a value for tolerance (iteration
498 stops when fabs(lamda [k] - lambda [k-1]) <= tolerance, where lamda [k] is
499 the eigenvalue at the k-th iteration step.
500 The methos is described in:
501 G. Golub & C. van Loan (1996), Matrix computations, third edition,
502 The Johns Hopkins University Press Ltd.,
503 London, (Par. 7.3.1 The Power Method)
504 */
505
506 integer NUMsolveQuadraticEquation (double a, double b, double c, double *x1, double *x2);
507 /*
508 Finds the real roots of ax^2 + bx + c = 0.
509 The number of real roots is returned and their locations in x1 and x2.
510 If only one root found it is stored in x1.
511 If no roots found then x1 and x2 will not be changed.
512 */
513
514 autoVEC newVECsolve (constMATVU const& a, constVECVU const& b, double tol);
515 /*
516 Solve the equation: A.x = b for x;
517 a [1..nr] [1..nc], b [1..nr] and the unknown x [1..nc]
518 Algorithm: s.v.d.
519 */
520
521 autoMAT newMATsolve (constMATVU const& a, constMATVU const& b, double tol);
522 /*
523 Solve the equations: A.X = B;
524 a [1..nr] [1..nc], b [1..nr] [1..nc2] and the unknown x [1..nc] [1..nc2]
525 Algorithm: s.v.d.
526 */
527
528
529 /*
530 Solve y = D.x + e for x, where x is sparse and e is observation noise.
531 Minimize the 2-norm (y - D.x), where maximally K elements of x may be non-zero, by an iterative hard thresholding algorithm.
532 D is a MxN real matrix with (many) more columns than rows, i.e. N > M. We need to find a vector x
533 with maximally K non-zero elements (sparse).
534 The algorithm is described in T. Blumensath & M.E. Davies (2010): "Normalised iterative hard thresholding;
535 guaranteed stability and performance", IEEE Journal of Selected Topics in Signal Processing #4: 298-309.
536 x in/out: the start value (you typically would start the iteration with all zeros).
537 */
538 void VECsolveSparse_IHT (VECVU const& x, constMATVU const& d, constVECVU const& y, integer numberOfNonZeros, integer maximumNumberOfIterations, double tolerance, integer infoLevel);
539 autoVEC newVECsolveSparse_IHT (constMATVU const& d, constVECVU const& y, integer numberOfNonZeros, integer maximumNumberOfIterations, double tolerance, integer infoLevel);
540
541 void VECsolveNonnegativeLeastSquaresRegression (VECVU const& result, constMATVU const& m, constVECVU const& y, integer itermax, double tol, integer infoLevel);
542
newVECsolveNonnegativeLeastSquaresRegression(constMATVU const & a,constVECVU const & y,integer itermax,double tol,integer infoLevel)543 inline autoVEC newVECsolveNonnegativeLeastSquaresRegression (constMATVU const& a, constVECVU const& y, integer itermax, double tol, integer infoLevel) {
544 autoVEC result = zero_VEC (a.ncol);
545 VECsolveNonnegativeLeastSquaresRegression (result.get(), a, y, itermax, tol, infoLevel);
546 return result;
547 }
548 /*
549 Solve the equation: A.x = y for x under the constraint: all x [i] >= 0;
550 a [1..nr][1..nc], y [1..nr] and x [1..nc].
551 Algorithm: Alternating least squares.
552 Borg & Groenen (1997), Modern multidimensional scaling, Springer, page 180.
553 */
554
555 void NUMsolveConstrainedLSQuadraticRegression (constMAT const& x, constVEC y, double *out_alpha, double *out_gamma);
556 /*
557 Solve y [i] = alpha + beta * x [i] + gamma * x [i]^2, with i = 1..n,
558 subject to the constraint beta^2 = 4 * alpha * gamma, for alpha and
559 gamma (Least Squares).
560 The input Vandermonde-matrix o [1..n,1..3] has columns with 1, x [i] and
561 x [i]^2, respectively.
562 The algorithm is according to:
563 Jos M.F. Ten Berge (1983), A generalization of Verhelst's solution for
564 a constrained regression problem in ALSCAL and related MDS-algorithms,
565 Psychometrika 48, 631-638.
566 */
567
568 autoVEC newVECsolveWeaklyConstrainedLinearRegression (constMAT const& a, constVEC const& y, double alpha, double delta);
569 /*
570 Solve g(x) = ||A*x - y||^2 + alpha (x'*x - delta)^2 for x [1..m],
571 where A [1..n] [1..m] is a matrix, y [1..n] a given vector, and alpha
572 and delta are fixed numbers.
573 This class of functions is composed of a linear regression function and
574 a penalty function for the sum of squared regression weights. It is weakly
575 constrained because the penalty function prohibits a relatively large
576 departure of x'x from delta.
577 The solution is due to:
578 Jos M.F. ten Berge (1991), A general solution for a class of weakly
579 constrained linear regression problems, Psychometrika 56, 601-609.
580 Preconditions:
581 a.nrow >= a.ncol
582 alpha >= 0
583 */
584
585 void NUMprocrustes (constMATVU const& x, constMATVU const& y, autoMAT *out_rotation, autoVEC *out_translation, double *out_scale);
586 /*
587 Given two configurations x and y (nPoints x nDimensions), find the
588 the Procrustes rotation/reflection matrix T, the translation vector v and the scaling
589 factor s such that Y = sXT+1v' (1 is the nPoints vector with ones).
590 Solution: see Borg and Groenen (1997), Modern Multidimensional Scaling, pp 340-346.
591 When out_translation == nullptr or out_scale == nullptr, only the matrix T will be solved for:
592 the orthogonal Procrustes transform.
593 */
594
595 double NUMnrbis (double (*f)(double x, double *dfx, void *closure), double xmin, double xmax, void *closure);
596 /*
597 Find the root of a function between xmin and xmax.
598 Method: Newton-Raphson with bisection (i.e., derivative is known!).
599 Error condition:
600 return undefined if root not bracketed.
601 */
602
603 double NUMridders (double (*f) (double x, void *closure), double xmin, double xmax, void *closure);
604 /*
605 Return the root of a function f bracketed in [xmin, xmax].
606 Error condition:
607 root not bracketed.
608 */
609
610 double NUMmspline (constVEC const & knot, integer order, integer i, double x);
611 /*
612 Calculates an M-spline for a knot sequence.
613 After Ramsay (1988), Monotone splines in action, Statistical Science 4.
614
615 M-splines of order k have degree k-1.
616 M-splines are zero outside interval [ knot [i], knot [i+order] ).
617 First and last 'order' knots are equal, i.e.,
618 knot [1] = ... = knot [order] && knot [nKnots-order+1] = ... knot [nKnots].
619 Error condition: no memory.
620 */
621
622 double NUMispline (constVEC const & aknot, integer order, integer i, double x);
623 /*
624 Calculates an I-spline for simple knot sequences: only one knot at each
625 interior boundary.
626 After Ramsay (1988), Monotone splines in action, Statistical Science 4.
627
628 I-splines of order k have degree k (because they Integrate an M-spline
629 of degree k-1).
630 In the calculation of the integral of M(x|k,t), M-splines are used that
631 have two more knots, i.e., M(x|k+1,t). For reasons of efficiency we
632 demand that these extra knots are given, i.e., the 'aknot []' argument
633 contains the knot positions as if the spline to be integrated were an
634 M(x|k+1,t) spline.
635 knot [1] = ... = knot [order+1] && knot [nKnots-order] = ... knot [nKnots]
636 Error condition: no memory.
637 */
638
639 double NUMwilksLambda (constVEC const& lambda, integer from, integer to);
640 /*
641 Calculate: Product (i=from..to; 1/(1+lambda [i]))
642 Preconditions: to >= from
643 */
644
645 double NUMlnBeta (double a, double b);
646 /*
647 Computes the logarithm of the beta function log(B(a,b) subject to
648 a and b not being negative integers.
649 */
650
651 double NUMbeta2 (double z, double w);//temporarily
652
653 double NUMbetaContinuedFraction(double a, double b, double x);
654
655 double NUMfactln (integer n);
656 /* Returns ln (n!) */
657
658 void NUMlngamma_complex (double zr, double zi, double *out_lnr, double *out_arg);
659 /* Log [Gamma(z)] for z complex, z not a negative integer
660 * Uses complex Lanczos method. Note that the phase part (arg)
661 * is not well-determined when |z| is very large, due
662 * to inevitable roundoff in restricting to (-pi, pi].
663 * The absolute value part (lnr), however, never suffers.
664 *
665 * Calculates:
666 * lnr = log|Gamma(z)|
667 * arg = arg(Gamma(z)) in (-Pi, Pi]
668 */
669
670 /***** STATISTICS: PROBABILITY DENSITY FUNCTIONS ********************/
671
672 double NUMlogNormalP (double x, double zeta, double sigma);
673 /* Area under log normal from 0 to x */
674
675 double NUMlogNormalQ (double x, double zeta, double sigma);
676 /* Area under log normal from x to +infinity */
677
678 double NUMstudentP (double t, double df);
679 /*
680 The area under the student T-distribution from -infinity to t.
681 Precondition: df > 0
682 */
683
684 double NUMstudentQ (double t, double df);
685 /*
686 The area under the student T distribution from t to +infinity.
687 Precondition: df > 0
688 */
689
690 double NUMfisherP (double f, double df1, double df2);
691 /*
692 The area under Fisher's F-distribution from 0 to f
693 Preconditions: f >= 0, df1 > 0, df2 > 0
694 */
695
696 double NUMfisherQ (double f, double df1, double df2);
697 /*
698 The area under Fisher's F-distribution from f to +infinity
699 Preconditions: f >= 0, df1 > 0, df2 > 0
700 */
701
702 double NUMinvGaussQ (double p);
703 /*
704 Solves NUMgaussQ (x) == p for x, given p.
705 Precondition: 0 < p < 1
706 Method: Abramovitz & Stegun 26.2.23
707 Precision: |eps(p)| < 4.5 10^-4
708 */
709
710 double NUMinvChiSquareQ (double p, double df);
711 /*
712 Solves NUMchiSquareQ (chiSquare, df) == p for chiSquare, given p, df.
713 Preconditions: 0 < p < 1, df > 0
714 */
715
716 double NUMinvStudentQ (double p, double df);
717 /*
718 Solves NUMstudentQ (t, df) == p for t, given p, df.
719 Preconditions: 0 < p < 1, df > 0
720 */
721
722 double NUMinvFisherQ (double p, double df1, double df2);
723 /*
724 Solves NUMfisherQ (f, df1, df2) == p for f, given p, df1, df2
725 Precondition: 0 < p < 1
726 */
727
728 double NUMtukeyQ (double q, double cc, double df, double rr);
729 /* Computes the probability that the maximum of rr studentized
730 * ranges, each based on cc means and with df degrees of freedom
731 * for the standard error, is larger than q.
732 */
733
734 double NUMinvTukeyQ (double p, double cc, double df, double rr);
735 /* Solves NUMtukeyQ (q, rr, cc, df) == p for q given p, rr, cc and df.
736 * Computes the quantiles of the maximum of rr studentized
737 * ranges, each based on cc means and with df degrees of freedom
738 * for the standard error, is larger than q.
739 * p = probability (alpha)
740 * rr = no. of rows or groups
741 * cc = no. of columns or treatments
742 * df = degrees of freedom of error term
743 */
744
745
746 /****** Frequency in Hz to other frequency reps ****/
747
748 double NUMmelToHertz2 (double mel);
749 /*
750 Return 700 * (pow (10.0, mel / 2595.0) - 1)
751 */
752
753 double NUMhertzToMel2 (double f);
754 /*
755 Return 2595 * log10 (1 + f/700)
756 */
757
758 double NUMmelToHertz3 (double mel);
759 /*
760 Return mel < 1000 ? mel : 1000 * (exp (mel * log10(2) / 1000) - 1)
761 */
762
763 double NUMhertzToMel3 (double hz);
764 /*
765 Return hz < 1000 ? hz : 1000 * log10 (1 + hz / 1000) / log10 (2)
766 */
767
768 double NUMbarkToHertz2 (double bark);
769 /*
770 Return 650 * sinh (bark / 7)
771 */
772
773 double NUMhertzToBark2 (double hz);
774 /*
775 Return 7 * ln (hz / 650 + sqrt(1 + (hz / 650)^2))
776 */
777
778 double NUMhertzToBark_traunmueller (double hz);
779 /*
780 return 26.81 * hz /(1960 + hz) -0.53;
781 */
782
783 double NUMbarkToHertz_traunmueller (double bark);
784 /*
785 return 1960* (bark + 0.53) / (26.28 - bark);
786 */
787
788 double NUMbarkToHertz_zwickerterhardt (double hz);
789 /*
790 return 13 * atan (0.00076 * hz) + 3.5 * atan (hz / 7500);
791 */
792
793 double NUMbladonlindblomfilter_amplitude (double zc, double z);
794 /*
795 Amplitude of filter at dz (barks) from centre frequency.
796 dz may be positive and negative.
797
798 The bladonlindblomfilter function is:
799
800 z' = zc - z + 0.474
801 10 log10 F(z') = 15.81 + 7.5 z' - 17.5 sqrt( 1 + z'^2 )
802
803 Reference: Bladon, R.A.W & Lindblom, B., (1980),
804 "Modeling the judgment of vowel quality differences", JASA 69, 1414-1422.
805 The filter has a bandwidth of 1.43 Bark, the maximum occurs at z = zc,
806 and the slopes are -10 dB/Bark and +25 dB/Bark.
807 */
808
809 double NUMsekeyhansonfilter_amplitude (double zc, double z);
810 /*
811 Amplitude of filter at dz (barks) from centre frequency.
812 dz may be positive and negative.
813
814 The sekeyhansonfilter function is:
815 z' = zc - z - 0.215
816 10 log10 F(z') = 7 - 7.5 * z' - 17.5 * sqrt( 0.196 + z'^2 )
817
818 Reference: Sekey, A. & Hanson, B.A. (1984),
819 "Improved 1-Bark bandwidth auditory filter", JASA 75, 1902-1904.
820 The filter function has a bandwidth of 1 Bark, the maximum response
821 occurs at z=zc, and the slopes are +10 dB/Bark and -25 dB/Bark.
822 It is an improved version of bladonlindblomfilter.
823 */
824
825 double NUMtriangularfilter_amplitude (double fl, double fc, double fh,
826 double f);
827 /*
828 Filterfunction that intermediates in Mel frequency cepstral coefficients
829 calculation.
830 The filter function is
831
832 (f-fl)/(fc-fl) fl < f < fc
833 H(z) = (fh-f)/(fh-fc) fc < f < fh
834 0 otherwise
835 Preconditions:
836 0 < fl < fh
837 */
838
839 double NUMformantfilter_amplitude (double fc, double bw, double f);
840 /*
841 Filterfunction with a formant-like shape on a linear freq. scale.
842
843 H(f) = 1.0 / (dq * dq + 1.0), where
844 dq = (fc * fc - f * f) / (bw * f)
845 Preconditions: f > 0 && bw > 0
846 */
847
848 double VECburg (VEC const& a, constVEC const& x);
849 /*
850 Calculates linear prediction coefficients according to the algorithm
851 from J.P. Burg as described by N.Anderson in Childers, D. (ed), Modern
852 Spectrum Analysis, IEEE Press, 1978, 252-255.
853 Returns the sum of squared sample values or 0.0 if failure
854 */
855
856 autoVEC newVECburg (constVEC const& x, integer numberOfPredictionCoefficients, double *out_xms);
857
858 void VECfilterInverse_inplace (VEC const& s, constVEC const& filter, VEC const& filterMemory);
859
860 void NUMdmatrix_to_dBs (MAT const& m, double ref, double factor, double floor);
861 /*
862 Transforms the values in the matrix m [rb..re] [cb..ce] to dB's
863
864 m [i] [j] = factor * 10 * log10 (m [i] [j] / ref)
865 if (m [i] [j] < floor) m [i] [j] = floor;
866
867 Preconditions:
868 rb <= re
869 cb <= ce
870 ref > 0
871 factor > 0
872 Errors:
873 Matrix elements < 0;
874 */
875
876 autoMAT MATcosinesTable (integer n);
877 /*
878 Generate table with cosines.
879
880 result [i] [j] = cos (i * pi * (j - 1/2) / npoints)
881 */
882 void VECcosineTransform_preallocated (VEC const& target, constVEC const& x, constMAT const& cosinesTable);
883 void VECinverseCosineTransform_preallocated (VEC const& target, constVEC const& x, constMAT const& cosinesTable);
884
885 /****** Interpolation ****/
886
887 void NUMcubicSplineInterpolation_getSecondDerivatives (VEC const& out_y, constVEC const& x, constVEC const& y, double yp1, double ypn);
888 /*
889 Given arrays x [1..n] and y [1..n] containing a tabulated function, i.e.,
890 y [i] = f(x [i]), with x [1] < x [2] < ... < x [n], and given values yp1 and
891 ypn for the first derivative of the interpolating function at point
892 1 and n, respectively, this routine returns an array out_y [1..n] that
893 contains the second derivative of the interpolating function at the
894 tabulated point x.
895 If yp1 and/or ypn are >= 10^30, the routine is signaled to
896 set the corresponding boundary condition for a natural spline, with
897 zero second derivative on that boundary.
898 */
899
900 double NUMcubicSplineInterpolation (constVEC const& xa, constVEC const& ya, constVEC const& y2a, double x);
901 /*
902 Given arrays xa [1..n] and ya [1..n] containing a tabulated function,
903 i.e., y [i] = f(x [i]), with x [1] < x [2] < ... < x [n], and given the
904 array y2a [1..n] which is the output of NUMcubicSplineInterpolation_getSecondDerivatives above, and given
905 a value of x, this routine returns an interpolated value y.
906 */
907
908 autoVEC newVECbiharmonic2DSplineInterpolation_getWeights (constVECVU const& x, constVECVU const& y, constVECVU const& w);
909 /*
910 Input: x [1..numberOfPoints], y [1..numberOfPoints], (xp,yp)
911 Output: interpolated result
912 */
913
914 double NUMbiharmonic2DSplineInterpolation (constVECVU const& x, constVECVU const& y, constVECVU const& w, double xp, double yp);
915 /* Biharmonic spline interpolation based on Green's function.
916 . Given z [i] values at points (x [i],y [i]) for i=1..n,
917 Get value at new point (px,py).
918 1. Calculate weights w once: newVECbiharmonic2DSplineInterpolation_getWeights
919 2. Interpolate at (xp,yp): NUMbiharmonic2DSplineInterpolation
920 Input: x [1..numberOfPoints], y [1..numberOfPoints], z [1..numberOfPoints], weights [1..numberOfPoints]
921 Output: weights [1..numberOfPoints]
922
923 Preconditions: all x [i] are different and all y [i] are different.
924
925 This routine inializes the numberOfPoints weigts by inverting a numberOfPoints x numberOfPoints matrix.
926 D. Sandwell (1987), Biharmonic spline interpolation of GEOS-3 and SEASAT altimetr data, Geophysical Research Letters 14, 139--142
927 X. Deng & Z. Tang (2011), Moving surface spline interpolation based on Green's function, Math. Geosci 43: 663--680
928 */
929
930
931 double NUMsincpi (const double x);
932 /* Calculates sin(pi*x)/(pi*x) */
933 double NUMsinc (const double x);
934 /* Calculates sin(x)/(x) */
935
936 /*********************** Geometry *************************************/
937
938 integer NUMgetOrientationOfPoints (double x1, double y1, double x2, double y2, double x3, double y3);
939 /* Traverse points 1, 2 and 3. If we travel counter-clockwise the result will be 1,
940 if we travel clockwise the result will be -1 and the result will be 0 if 3 is on the line segment between 1 and 2.
941 J. O'Rourke: Computational Geometry, 2nd Edition, Code 1.5
942 */
943
944 bool NUMdoLineSegmentsIntersect (double x1, double y1, double x2, double y2, double x3, double y3,
945 double x4, double y4);
946 /* Does the line segment from (x1,y1) to (x2,y2) intersect with the line segment from (x3,y3) to (x4,y4)? */
947
948 integer NUMgetIntersectionsWithRectangle (double x1, double y1, double x2, double y2,
949 double xmin, double ymin, double xmax, double ymax, double *xi, double *yi);
950 /* Get the intersection points of the line through the points (x1,y1) and (x2,y2) with the
951 rectangle with corners (xmin, ymin) and (xmax,ymax).
952 The returned value is the number of intersections found and is either 0 or 1 or 2.
953 */
954
955 bool NUMclipLineWithinRectangle (double line_x1, double line_y1, double line_x2, double line_y2, double rect_x1, double rect_y1, double rect_x2, double rect_y2, double *out_line_x1, double *out_line_y1, double *out_line_x2, double *out_line_y2);
956 /*
957 If true, then returns in (out_line_x1, out_line_y1) and (out_line_x2, out_line_y2) the coordinates of start and end points of the line (line_x1, line_y1)..(line_x2, line_y2) that can be drawn within the rectangle with lowerleft corner (rect_x1, rect_y1) and upperright (rect_x2, rect_y2).
958 Returns false if there is nothing to be drawn inside.
959 */
960
961 void NUMgetEllipseBoundingBox (double a, double b, double cospsi,
962 double *out_width, double *out_height);
963 /*
964 Get the width and the height of the bonding box around an ellipse.
965 a and b are the lengths of the long axes.
966 cospsi is the cosine of the angle between the a-axis and the horizontal
967 x-axis (cs == 0 when a-axis and x-axis are perpendicular).
968
969 Parametrisation of the ellipse:
970 x(phi) = a cos(psi) cos(phi) - b sin (psi) sin(phi)
971 y(phi) = a sin(psi) cos(phi) + b cos(psi) sin(phi) 0 <= phi <= 2 pi
972 Extrema:
973 d x(phi) / dphi == 0 and d y(phi) / dphi == 0
974 Solution:
975 x(phi1) = a cos(psi) cos(phi1) - b sin (psi) sin(phi1)
976 y(phi2) = a sin(psi) cos(phi2) + b cos(psi) sin(phi2),
977 where
978 phi1 = arctg ( -b/a tg(psi))
979 phi2 = arctg ( b/a cotg(psi))
980 Special cases are psi = 0 and pi /2
981 */
982
983 double NUMminimize_brent (double (*f) (double x, void *closure), double a, double b,
984 void *closure, double tol, double *fx);
985 /*
986 The function returns an estimate for the minimum location with accuracy
987 3 * SQRT_EPSILON * abs(x) + tol.
988 The function always obtains a local minimum which coincides with
989 the global one only if a function under investigation being unimodular.
990 If a function being examined possesses no local minimum within
991 the given range, the function returns 'a' (if f(a) < f(b)), otherwise
992 it returns the right range boundary value b.
993
994 Algorithm
995
996 The function makes use of the golden section procedure combined with
997 parabolic interpolation.
998 At every step, the program operates at three abscissae - x, v, and w.
999 x - the last and the best approximation to the minimum location,
1000 i.e. f(x) <= f(a) or/and f(x) <= f(b)
1001 (if the function f has a local minimum in (a,b), then both
1002 conditions are fulfiled after one or two steps).
1003 v, w are previous approximations to the minimum location. They may
1004 coincide with a, b, or x (although the algorithm tries to make all
1005 u, v, and w distinct). Points x, v, and w are used to construct
1006 interpolating parabola whose minimum will be treated as a new
1007 approximation to the minimum location if the former falls within
1008 [a,b] and reduces the range enveloping minimum more efficient than
1009 the golden section procedure.
1010 When f(x) has a second derivative positive at the minimum location
1011 (not coinciding with a or b) the procedure converges superlinearly
1012 at a rate order about 1.324
1013 */
1014
1015 /********************** fft ******************************************/
1016
1017 struct structNUMfft_Table
1018 {
1019 integer n;
1020 autoVEC trigcache;
1021 autoINTVEC splitcache;
1022 };
1023
1024 typedef struct structNUMfft_Table *NUMfft_Table;
1025
1026 void NUMfft_Table_init (NUMfft_Table table, integer n);
1027 /*
1028 n : data size
1029 */
1030
1031 struct autoNUMfft_Table : public structNUMfft_Table {
throwautoNUMfft_Table1032 autoNUMfft_Table () throw () {
1033 n = 0;
1034 }
~autoNUMfft_TableautoNUMfft_Table1035 ~autoNUMfft_Table () { }
1036 };
1037
1038 void NUMfft_forward (NUMfft_Table table, VEC data);
1039 /*
1040 Function:
1041 Calculates the Fourier Transform of a set of n real-valued data points.
1042 Replaces this data in array data [1...n] by the positive frequency half
1043 of its complex Fourier Transform, with a minus sign in the exponent.
1044 Preconditions:
1045 data != NULL;
1046 table must have been initialised with NUMfft_Table_init
1047 Postconditions:
1048 data [1] contains real valued first component (Direct Current)
1049 data [2..n-1] even index : real part; odd index: imaginary part of DFT.
1050 data [n] contains real valued last component (Nyquist frequency)
1051
1052 Output parameters:
1053
1054 data r(1) = the sum from i=1 to i=n of r(i)
1055
1056 If l =(int) (n+1)/2
1057
1058 then for k = 2,...,l
1059
1060 r(2*k-2) = the sum from i = 1 to i = n of r(i)*cos((k-1)*(i-1)*2*pi/n)
1061
1062 r(2*k-1) = the sum from i = 1 to i = n of -r(i)*sin((k-1)*(i-1)*2*pi/n)
1063
1064 if n is even
1065
1066 r(n) = the sum from i = 1 to i = n of (-1)**(i-1)*r(i)
1067
1068 i.e., the ordering of the output array will be for n even
1069 r(1),(r(2),i(2)),(r(3),i(3)),...,(r(l-1),i(l-1)),r(l).
1070 Or ...., (r(l),i(l)) for n uneven.
1071
1072 ***** note
1073 this transform is unnormalized since a call of NUMfft_forward
1074 followed by a call of NUMfft_backward will multiply the input sequence by n.
1075 */
1076
1077 void NUMfft_backward (NUMfft_Table table, VEC data);
1078 /*
1079 Function:
1080 Calculates the inverse transform of a complex array if it is the transform of real data.
1081 (Result in this case should be multiplied by 1/n.)
1082 Preconditions:
1083 n is an integer power of 2.
1084 data [1] contains real valued first component (Direct Current)
1085 data [2..n-1] even index : real part; odd index: imaginary part of DFT.
1086 data [n] contains real valued last component (Nyquist frequency)
1087
1088 table must have been initialised with NUMfft_Table_init
1089
1090 Output parameters
1091
1092 data for n even and for i = 1,...,n
1093
1094 r(i) = r(1)+(-1)**(i-1)*r(n)
1095
1096 plus the sum from k=2 to k=n/2 of
1097
1098 2.0*r(2*k-2)*cos((k-1)*(i-1)*2*pi/n) -2.0*r(2*k-1)*sin((k-1)*(i-1)*2*pi/n)
1099
1100 for n odd and for i = 1,...,n
1101
1102 r(i) = r(1) plus the sum from k=2 to k=(n+1)/2 of
1103
1104 2.0*r(2*k-2)*cos((k-1)*(i-1)*2*pi/n) -2.0*r(2*k-1)*sin((k-1)*(i-1)*2*pi/n)
1105
1106 ***** note
1107 this transform is unnormalized since a call of NUMfft_forward
1108 followed by a call of NUMfft_backward will multiply the input
1109 sequence by n.
1110 */
1111
1112 /**** Compatibility with NR fft's */
1113
1114 void NUMforwardRealFastFourierTransform (VEC data);
1115 /*
1116 Function:
1117 Calculates the Fourier Transform of a set of n real-valued data points.
1118 Replaces this data in array data [1...n] by the positive frequency half
1119 of its complex Fourier Transform, with a minus sign in the exponent.
1120 Preconditions:
1121 n is an integer power of 2.
1122 data != NULL;
1123 Postconditions:
1124 data [1] contains real valued first component (Direct Current)
1125 data [2] contains real valued last component (Nyquist frequency)
1126 data [3..n] odd index : real part; even index: imaginary part of DFT.
1127 */
1128 void NUMreverseRealFastFourierTransform (VEC data);
1129 /*
1130 Function:
1131 Calculates the inverse transform of a complex array if it is the transform of real data.
1132 (Result in this case should be multiplied by 1/n.)
1133 Preconditions:
1134 n is an integer power of 2.
1135 data [1] contains real valued first component (Direct Current)
1136 data [2] contains real valued last component (Nyquist frequency)
1137 data [3..n] odd index : real part; even index: imaginary part of DFT.
1138 */
1139 void NUMrealft (VEC data, integer direction);
1140
1141 void VECsmooth_gaussian_inplace (VECVU const& in_out, double sigma);
1142 void VECsmooth_gaussian_inplace (VECVU const& in_out, double sigma, NUMfft_Table fftTable);
1143 void VECsmooth_gaussian (VECVU const& out, constVECVU const& in, double sigma, NUMfft_Table fftTable);
1144 /*
1145 Smooth the vector 'in/in_out' by convolving with a Gaussian, i.e. convolve with gaussian by
1146 using the Fourier Transform. Normally an FFT is used unless otherwise specified in 'fftTable"
1147 If fftTable == nullptr the FFT of size 2^k is used, where 2^(k-1) < n <= 2^k.
1148 For a given fftTable we require that fftTable->n >= n.
1149 */
1150
1151 integer NUMgetIndexFromProbability (constVEC probs, double p); //TODO HMM zero start matrices
1152 integer NUMgetIndexFromProbability (double *probs, integer nprobs, double p);
1153
1154 // Fit the line y= ax+b
1155 void NUMlineFit (constVEC x, constVEC y, double *out_m, double *out_intercept, integer method);
1156 /* method
1157 * 1 least squares
1158 * 2 rubust incomplete Theil O(N/2)
1159 * 3 robust complete Theil (very slow for large N, O(N^2))
1160 */
1161
1162 void NUMfitLine_theil_preallocated (VEC const& out_lineParameters, constVEC const& x, constVEC const& y, bool wantIntercept, double oneTailedUnconfidence, bool completeMethod);
1163
1164 void NUMfitLine_theil (constVEC const& x, constVEC const& y, double *out_m, double *out_intercept, bool completeMethod);
1165 /*
1166 * Preconditions:
1167 * all x [i] should be different, i.e. x [i] != x [j] for all i = 1..(numberOfPoints - 1), j = (i+1) ..numberOfPoints
1168 * Algorithm:
1169 * Theils robust line fit method:
1170 * 1. Use all combination of pairs (x [i],y [i]), (x [j],y [j]) to calculate an intercept m [k] as
1171 * m [k] = (y [j] - y [i]) / (x [j] - x [i]).
1172 * There will be (numberOfPoints - 1) * numberOfPoints / 2 numbers m [k].
1173 * 2. Take the median value m of all the m [k].
1174 * 3. Calculate the numberOfPoints intercepts b [i] as b [i] = y [i] - m * x [i]
1175 * 4. Take the median value b of all the b [i] values
1176 *
1177 * If incompleteMethod we use Theil's incomplete method to reduce the number of combinations.
1178 * I.e. split the data in two equal parts at n2 = numberOfPoints / 2 and then calculate the numberOfPoints/2 intercepts m [i] as
1179 * m [i] = (y [n2+i] - y [i]) / (x [n2 + i] - x [i]).
1180 * The rest proceeds as outlined above
1181 */
1182
1183
1184 void NUMfitLine_LS (constVEC const& x, constVEC const& y, double *out_m, double *out_intercept);
1185
1186 /*
1187 scale_lambda, shape_k > 0
1188 */
1189 double NUMrandomWeibull (double scale_lambda, double shape_k);
1190
1191 /* The binomial distribution has the form,
1192
1193 f(x) = n!/(x!(n-x)!) * p^x (1-p)^(n-x) for integer 0 <= x <= n
1194 = 0 otherwise
1195
1196 This implementation follows the public domain ranlib function
1197 "ignbin", the bulk of which is the BTPE (Binomial Triangle
1198 Parallelogram Exponential) algorithm introduced in
1199 Kachitvichyanukul and Schmeiser [1]. It has been translated to use
1200 modern C coding standards.
1201
1202 If n is small and/or p is near 0 or near 1 (specifically, if
1203 n*min(p,1-p) < SMALL_MEAN), then a different algorithm, called
1204 BINV, is used which has an average runtime that scales linearly
1205 with n*min(p,1-p).
1206
1207 But for larger problems, the BTPE algorithm takes the form of two
1208 functions b(x) and t(x) -- "bottom" and "top" -- for which b(x) <
1209 f(x)/f(M) < t(x), with M = floor(n*p+p). b(x) defines a triangular
1210 region, and t(x) includes a parallelogram and two tails. Details
1211 (including a nice drawing) are in the paper.
1212
1213 [1] Kachitvichyanukul, V. and Schmeiser, B. W. Binomial Random
1214 Variate Generation. Communications of the ACM, 31, 2 (February,
1215 1988) 216.
1216
1217 Note, Bruce Schmeiser (personal communication) points out that if
1218 you want very fast binomial deviates, and you are happy with
1219 approximate results, and/or n and n*p are both large, then you can
1220 just use gaussian estimates: mean=n*p, variance=n*p*(1-p).
1221
1222 This implementation by James Theiler, April 2003, after obtaining
1223 permission -- and some good advice -- from Drs. Kachitvichyanukul
1224 and Schmeiser to use their code as a starting point, and then doing
1225 a little bit of tweaking.
1226
1227 Additional polishing for GSL coding standards by Brian Gough.
1228 */
1229 integer NUMrandomBinomial (double p, integer n);
1230 double NUMrandomBinomial_real (double p, integer n);
1231
1232 /*
1233 Generates random numbers according to a Gamma distribution with shape parameter "alpha"
1234 and rate parameter "beta".
1235
1236 The Gamma distribution of order (shape) parameter alpha and rate (beta) is defined as:
1237
1238 f(x; alpha, beta) = (1 / Gamma (alpha)) beta^alpha x^(alpha-1) e^(-beta.x),
1239 for x > 0, alpha > 0 && beta > 0.
1240
1241 The method is described in
1242 G. Marsaglia & W. Tsang (2000): A simple method for generating gamma variables. ACM Transactions on Mathematical Software, 26(3):363-372.
1243 Preconditions: alpha > 0 && beta > 0.
1244 */
1245 double NUMrandomGamma (const double alpha, const double beta);
1246
1247 // IEEE: Programs for digital signal processing section 4.3 LPTRN (modfied)
1248 // lpc [1..n] to rc [1..n]
1249 void VECrc_from_lpc (VEC rc, constVEC lpc);
1250
1251 // rc [1..n] to area [1..n], implicit: area [n+1] = 0.0001; (1 cm^2)
1252 void VECarea_from_rc (VEC area, constVEC rc);
1253
1254 // area [1..n] to rc [1..n-1] (modification: LPTRN assumes area [n+1])
1255 void VECrc_from_area (VEC rc, constVEC area);
1256
1257 // area [1..n] to lpc [1..n-1]! (modification: lptrn gives lpc [1] = 1 we don't)
1258 void VEClpc_from_area (VEC lpc, constVEC area);
1259
1260 // lpc [1..n] to area [1..n+1], area [m+1] = 0.0001; (1 cm^2)
1261 void VECarea_from_lpc (VEC area, constVEC lpc);
1262 /*
1263 Fix indices to be in the range [lowerLimit, upperLimit].
1264 */
1265 void NUMfixIndicesInRange (integer lowerLimit, integer upperLimit, integer *lowIndex, integer *highIndex);
1266
1267 void NUMgetEntropies (constMATVU const& m, double *out_h, double *out_hx,
1268 double *out_hy, double *out_hygx, double *out_hxgy, double *out_uygx, double *out_uxgy, double *out_uxy);
1269
NUMmean_weighted(constVEC x,constVEC w)1270 inline double NUMmean_weighted (constVEC x, constVEC w) {
1271 Melder_assert (x.size == w.size);
1272 double inproduct = NUMinner (x, w);
1273 double wsum = NUMsum (w);
1274 return inproduct / wsum;
1275 }
1276
VECchainRows_preallocated(VECVU const & v,constMATVU const & m)1277 inline void VECchainRows_preallocated (VECVU const& v, constMATVU const& m) {
1278 Melder_assert (m.nrow * m.ncol == v.size);
1279 integer k = 1;
1280 for (integer irow = 1; irow <= m.nrow; irow ++)
1281 for (integer icol = 1; icol <= m.ncol; icol ++)
1282 v [k ++] = m [irow] [icol];
1283 }
1284
VECchainRows(constMATVU const & m)1285 inline autoVEC VECchainRows (constMATVU const& m) {
1286 autoVEC result = raw_VEC (m.nrow * m.ncol);
1287 VECchainRows_preallocated (result.get(), m);
1288 return result;
1289 }
1290
VECchainColumns_preallocated(VEC const & v,constMATVU const & m)1291 inline void VECchainColumns_preallocated (VEC const& v, constMATVU const& m) {
1292 Melder_assert (m.nrow * m.ncol == v.size);
1293 integer k = 1;
1294 for (integer icol = 1; icol <= m.ncol; icol ++)
1295 for (integer irow = 1; irow <= m.nrow; irow ++)
1296 v [k ++] = m [irow] [icol];
1297 }
1298
VECchainColumns(constMATVU const & m)1299 inline autoVEC VECchainColumns (constMATVU const& m) {
1300 autoVEC result = raw_VEC (m.nrow * m.ncol);
1301 VECchainColumns_preallocated (result.get(), m);
1302 return result;
1303 }
1304
1305 /* R = X.Y.Z */
1306 void MATmul3 (MATVU const & target, constMATVU const& X, constMATVU const& Y, constMATVU const& Z);
1307
1308 /* Z = X.Y.X' */
1309 void MATmul3_XYXt (MATVU const& target, constMATVU const& X, constMATVU const& Y);
1310
1311 /* Z = X.Y.X where Y is a symmetric matrix */
1312 void MATmul3_XYsXt (MATVU const& target, constMATVU const& X, constMATVU const& Y);
1313
1314 /*
1315 First row (n elements) is at v [1]..v [n],
1316 second row (n-1 elements) is at v [n+1],..,v [n+n-1],
1317 third row (n-2 elements) is at v [n+n],..,v [n+n+n-2]
1318 last row (1 element) is at v [n(n+1)/2].
1319 */
MATfromUpperTriangularVector_preallocated(MAT m,constVEC v)1320 inline void MATfromUpperTriangularVector_preallocated (MAT m, constVEC v) {
1321 Melder_assert (v.size == m.ncol * (m.ncol + 1) / 2);
1322 integer irow = 1;
1323 for (integer inum = 1; inum <= v.size; inum ++) {
1324 integer nskipped = (irow - 1) * irow / 2;
1325 integer inumc = inum + nskipped;
1326 irow = (inumc - 1) / m.ncol + 1;
1327 integer icol = ( (inumc - 1) % m.ncol) + 1;
1328 m [irow] [icol] = m [icol] [irow] = v [inum];
1329 if (icol == m.ncol) irow ++;
1330 }
1331 }
1332
1333 void NUMeigencmp22 (double a, double b, double c, double *out_rt1, double *out_rt2, double *out_cs1, double *out_sn1 );
1334 /*
1335 This routine is copied from LAPACK.
1336 Computes the eigendecomposition of a 2-by-2 symmetric matrix
1337 [ a b ]
1338 [ b c ].
1339 on return, rt1 is the eigenvalue of larger absolute value, rt2 is the
1340 eigenvalue of smaller absolute value, and (cs1,sn1) is the unit right
1341 eigenvector for rt1, giving the decomposition
1342
1343 [ cs1 sn1 ] [ a b ] [ cs1 -sn1 ] = [ rt1 0 ]
1344 [-sn1 cs1 ] [ b c ] [ sn1 cs1 ] [ 0 rt2 ].
1345
1346
1347
1348 rt1 is accurate to a few ulps barring over/underflow.
1349
1350 rt2 may be inaccurate if there is massive cancellation in the
1351 determinant a*c-b*b; higher precision or correctly rounded or
1352 correctly truncated arithmetic would be needed to compute rt2
1353 accurately in all cases.
1354
1355 cs1 and sn1 are accurate to a few ulps barring over/underflow.
1356
1357 overflow is possible only if rt1 is within a factor of 5 of overflow.
1358 underflow is harmless if the input data is 0 or exceeds
1359 underflow_threshold / macheps.
1360 */
1361
1362 void NUMpolynomial_recurrence (VEC const& pn, double a, double b, double c, constVEC const& pnm1, constVEC const& pnm2);
1363
1364
1365 /* 20200405 djmw This functions resides here temporarily until MelderThread.h copes with lambda's */
NUMgetThreadingInfo(integer numberOfFrames,integer maximumNumberOfThreads,integer * inout_numberOfFramesPerThread,integer * out_numberOfThreads)1366 static inline void NUMgetThreadingInfo (integer numberOfFrames, integer maximumNumberOfThreads, integer *inout_numberOfFramesPerThread, integer *out_numberOfThreads) {
1367 if (*inout_numberOfFramesPerThread <= 0)
1368 *inout_numberOfFramesPerThread = 25;
1369 integer numberOfThreads = (numberOfFrames - 1) / *inout_numberOfFramesPerThread + 1;
1370 Melder_clip (1_integer, & numberOfThreads, maximumNumberOfThreads);
1371 *inout_numberOfFramesPerThread = (numberOfFrames - 1) / numberOfThreads + 1;
1372 if (out_numberOfThreads)
1373 *out_numberOfThreads = numberOfThreads;
1374 }
1375
1376 #endif // _NUM2_h_
1377