EEG/manual_EEG.cpp

@@ -63,7 +63,7 @@ NORMAL (U"##Subtracting a reference signal.# The voltage on the whole scalp may
 	"and use ##Subtract mean channel...# to subtract this from each electrode channel.")
 NORMAL (U"##Detrending.# With #Detrend, you subtract from each electrode channel a line in such a way that the first sample and the last sample become zero. "
 	"Detrending and reference subtraction can be performed in either order.")
-NORMAL (U"##Filtering.# With ##Filter...#, you band-pass filter each electrode channel. Filtering has to be done after detrending, but "
+NORMAL (U"##Filtering.# With ##Filter...#, you band-pass filter each electrode channel. Filtering should be done after detrending, but "
 	"filtering and reference subtraction can be performed in either order.")
 ENTRY (U"4. How to do an ERP analysis")
 NORMAL (U"An ERP is an Event-Related Potential. Events are marked somewhere in S1, S2, ... S8. In the above example, "

FFNet/FFNet.cpp

@@ -261,8 +261,8 @@ double FFNet_getBias (FFNet me, integer layer, integer unit) {
 		integer bias_unit = my wLast [node];
 		return my w [bias_unit];
 	} catch (MelderError) {
-		return undefined;
 		Melder_clearError ();
+		return undefined;
 	}
 }
 

FFNet/FFNet_Eigen.cpp

@@ -26,7 +26,7 @@
 #include "NUM2.h"
 
 void FFNet_Eigen_drawIntersection (FFNet me, Eigen eigen, Graphics g, integer pcx, integer pcy, double xmin, double xmax, double ymin, double ymax) {
-	integer ix = labs (pcx), iy = labs (pcy);
+	integer ix = integer_abs (pcx), iy = integer_abs (pcy);
 	integer numberOfEigenvalues = eigen -> numberOfEigenvalues;
 	integer dimension = eigen -> dimension;
 

FFNet/manual_FFNet.cpp

@@ -329,7 +329,7 @@ DEFINITION (U"the number of different categories that you want the net to learn.
 TAG (U"##Number of units in hidden layer 1#, ##Number of units in hidden layer 2#") \
 DEFINITION (U"determine the number of units in the hidden layers. " \
 	"If you want a neural net with no hidden layers, both numbers have to be 0. "\
-	"If you want a neural net with only 1 hidden layer then one of these numbers has to differ from 0. ")
+	"If you want a neural net with only 1 hidden layer then one of these numbers should differ from 0. ")
 
 MAN_BEGIN (U"Create FFNet...", U"djmw", 20040420)
 INTRO (U"Create a new feedforward neural net of type @FFNet.")

debian/changelog

+praat (6.1.03-1) unstable; urgency=medium
+
+  * New upstream version 6.1.03
+  * d/p/fix-test-electroglottogram.patch: New patch
+  * d/t/run-tests: Exclude unit test MacRoman.praat
+
+ -- Rafael Laboissiere <rafael@debian.org>  Tue, 24 Sep 2019 16:21:24 -0300
+
 praat (6.1-1) unstable; urgency=medium
 
   * New upstream version 6.1

debian/patches/fix-test-electroglottogram.patch

+Description: Fix Electroglottogram unit test
+Author: Rafael Laboissiere <rafael@debian.org>
+Origin: https://github.com/rlaboiss/praat/commit/87b0067a86b73dadd35cf119d8082d3790f0594d
+Last-Update: 2019-09-24
+
+--- praat-6.1.03.orig/dwtest/test_Electroglottogram.praat
++++ praat-6.1.03/dwtest/test_Electroglottogram.praat
+@@ -8,18 +8,16 @@ pitchCeiling = 500
+ closingThreshold = 0.30
+ silenceThreshold = 0.03
+ sound = Read from file: "s_egg_test.wav"
+-egg = Extract Electroglottogram: 1
+-intervalTier = To IntervalTier: pitchFloor, pitchCeiling, closingThreshold, silenceThreshold
++egg = Extract Electroglottogram: 1, "no"
++intervalTier = Get closed glottis intervals: pitchFloor, pitchCeiling, closingThreshold, silenceThreshold
+ selectObject: egg
+-degg = To Electroglottogram (derivative): 5000, 50
++degg = Derivative: 5000, 50, "yes"
+ selectObject: egg
+ To AmplitudeTier (levels): pitchFloor, pitchCeiling, closingThreshold, "yes", "yes"
+ peaks = selected ("AmplitudeTier", 1)
+ valleys = selected ("AmplitudeTier", 2)
+ levels = selected ("AmplitudeTier", 3)
+-selectObject: sound, egg, degg
+-combined = Combine to multi-channel
+ 
+-removeObject: levels, valleys, peaks, degg, intervalTier, egg, sound, combined
++removeObject: levels, valleys, peaks, degg, intervalTier, egg, sound
+ 
+ appendInfoLine : "test_Electroglottogram.praat OK"

debian/patches/series

@@ -2,3 +2,4 @@ use-ldflags.patch
 remove-time-date-macros.patch
 fix-procrustes-unit-test.patch
 cross-build.patch
+fix-test-electroglottogram.patch

debian/tests/run-tests

@@ -4,6 +4,7 @@
 ### distributed in the upstream sources.  For now, just the *.praat
 ### scripts in the test/* directories are executed.
 
+pwd
 ### Allow specification of the Praat executable
 PRAAT=${PRAAT:-praat}
 
@@ -56,6 +57,7 @@ for d in $(ls -d */) ; do
 		 -a $f != fon/resample48_8.praat			\
 		 -a $f != fon/resample96_8.praat			\
 		 -a $f != fon/Spectrum_draw.praat			\
+		 -a $f != kar/MacRoman.praat				\
 		 -a $f != sys/graphics.praat				\
 		 -a $f != sys/graphicsText.praat			\
 		 -a $f != sys/graphicsTextSpeed.praat ] ; then

dwsys/NUM2.cpp

@@ -34,7 +34,7 @@
  	invStudentQ, invChiSquareQ: modifications for 'undefined' return values.
  djmw 20030830 Corrected a bug in NUMtriangularfilter_amplitude
  djmw 20031111 Added NUMdmatrix_transpose, NUMdmatrix_printMatlabForm
- djmw 20040105 Added NUMmahalanobisDistance_chi
+ djmw 20040105 Added NUMmahalanobisDistanceSquared_chi
  djmw 20040303 Added NUMstring_containsPrintableCharacter.
  djmw 20050406 NUMprocrutus->NUMprocrustes
  djmw 20060319 NUMinverse_cholesky: calculation of determinant is made optional
@@ -278,7 +278,7 @@ autoMAT newMATinverse_fromLowerCholeskyInverse (constMAT m) {
 	return result;
 }
 
-double NUMmahalanobisDistance (constMAT lowerInverse, constVEC v, constVEC m) {
+double NUMmahalanobisDistanceSquared (constMAT lowerInverse, constVEC v, constVEC m) {
 	Melder_assert (lowerInverse.ncol == v.size && v.size == m.size);
 	longdouble chisq = 0.0;
 	if (lowerInverse.nrow == 1) { // diagonal matrix is one row matrix
@@ -1193,66 +1193,28 @@ double NUMlnBeta (double a, double b) {
 	return status == GSL_SUCCESS ? result.val : undefined;
 }
 
-double NUMnormalityTest_HenzeZirkler (constMAT data, double *inout_beta, double *out_tnb, double *out_lnmu, double *out_lnvar) {
-	const integer n = data.nrow, p = data.ncol;
-	if (*inout_beta <= 0)
-		*inout_beta = (1.0 / sqrt (2.0)) * pow ((1.0 + 2 * p) / 4.0, 1.0 / (p + 4.0)) * pow (n, 1.0 / (p + 4.0));
-	const double p2 = p / 2.0;
-	const double beta2 = *inout_beta * *inout_beta, beta4 = beta2 * beta2, beta8 = beta4 * beta4;
-	const double gamma = 1.0 + 2.0 * beta2, gamma2 = gamma * gamma, gamma4 = gamma2 * gamma2;
-	const double delta = 1.0 + beta2 * (4.0 + 3.0 * beta2), delta2 = delta * delta;
-	double prob = undefined;
-
-	double tnb = undefined, lnmu = undefined, lnvar = undefined;
-
-	if (n < 2 || p < 1)
-		return prob;
-
-	autoVEC zero = newVECzero (p);
-	autoMAT x = newMATcopy (data);
-
-	MATcentreEachColumn_inplace (x.get()); // x - xmean
-
-	autoMAT covar = MATcovarianceFromColumnCentredMatrix (x.get(), 0);
-
+void MATscaledResiduals (MAT const& residuals, constMAT const& data, constMAT const& covariance, constVEC const& means) {
 	try {
-		MATlowerCholeskyInverse_inplace (covar.get(), nullptr);
-		longdouble sumjk = 0.0, sumj = 0.0;
-		const double b1 = beta2 / 2.0, b2 = b1 / (1.0 + beta2);
-		/* Heinze & Wagner (1997), page 3
-			We use d [j] [k] = ||Y [j]-Y [k]||^2 = (Y [j]-Y [k])'S^(-1)(Y [j]-Y [k])
-			So d [j] [k]= d [k] [j] and d [j] [j] = 0
-		*/
-		for (integer j = 1; j <= n; j ++) {
-			for (integer k = 1; k < j; k ++) {
-				const double djk = NUMmahalanobisDistance (covar.get(), x.row (j), x.row(k));
-				sumjk += 2.0 * exp (-b1 * djk); // factor 2 because d [j] [k] == d [k] [j]
+		Melder_require (residuals.nrow == data.nrow && residuals.ncol == data.ncol,
+			U"The data and the residuals should have the same dimensions.");
+		Melder_require (covariance.ncol == means.size && data.ncol == means.size,
+			U"The dimensions of the means and the covariance have to conform with the data.");
+		autoVEC dif = newVECraw (data.ncol);
+		autoMAT lowerInverse = newMATcopy (covariance);
+		MATlowerCholeskyInverse_inplace (lowerInverse.get(), nullptr);
+		for (integer irow = 1; irow <= data.nrow; irow ++) {
+			dif <<= data.row (irow)  -  means;
+			residuals.row(irow) <<= 0.0;
+			if (lowerInverse.nrow == 1) { // diagonal matrix is one row matrix
+				residuals.row(irow) <<= lowerInverse.row(1)  *  dif.get();
+			} else {// square matrix
+				for (integer icol = 1; icol <= data.ncol; icol ++)
+					residuals [irow] [icol] = NUMinner (lowerInverse.row(icol).part (1, icol), dif.part (1, icol));
 			}
-			sumjk += 1.0; // for k == j
-			const double djj = NUMmahalanobisDistance (covar.get(), x.row (j), zero.get());
-			sumj += exp (-b2 * djj);
 		}
-		tnb = (1.0 / n) * (double) sumjk - 2.0 * pow (1.0 + beta2, - p2) * (double) sumj + n * pow (gamma, - p2); // n *
 	} catch (MelderError) {
-		Melder_clearError ();
-		tnb = 4.0 * n;
-	}
-
-	double mu = 1.0 - pow (gamma, -p2) * (1.0 + p * beta2 / gamma + p * (p + 2) * beta4 / (2.0 * gamma2));
-	double var = 2.0 * pow (1.0 + 4.0 * beta2, -p2)
-		 + 2.0 * pow (gamma,  -p) * (1.0 + 2.0 * p * beta4 / gamma2  + 3.0 * p * (p + 2) * beta8 / (4.0 * gamma4))
-		 - 4.0 * pow (delta, -p2) * (1.0 + 3.0 * p * beta4 / (2.0 * delta) + p * (p + 2) * beta8 / (2.0 * delta2));
-	double mu2 = mu * mu;
-	lnmu = log (sqrt (mu2 * mu2 / (mu2 + var)));
-	lnvar = sqrt (log ((mu2 + var) / mu2));
-	if (out_lnmu)
-		*out_lnmu = lnmu;
-	if (out_lnvar)
-		*out_lnvar = lnvar;
-	if (out_tnb)
-		*out_tnb = tnb;
-	prob = NUMlogNormalQ (tnb, lnmu, lnvar);
-	return prob;
+		Melder_throw (U"MATscaleResiduals: not performed.");
+	}
 }
 
 /*************** Hz <--> other freq reps *********************/
@@ -1760,7 +1722,8 @@ bool NUMclipLineWithinRectangle (double xl1, double yl1, double xl2, double yl2,
 	}
 	if (ncrossings == 0)
 		return false;
-	Melder_require (ncrossings <= 2, U"Too many crossings found.");
+	Melder_require (ncrossings <= 2,
+		U"Too many crossings found.");
 
 	/*
 		if start and endpoint of line are outside rectangle and ncrossings == 1, than the line only touches.
@@ -1788,7 +1751,8 @@ end:
 }
 
 void NUMgetEllipseBoundingBox (double a, double b, double cospsi, double *out_width, double *out_height) {
-	Melder_assert (cospsi >= -1.0 && cospsi <= 1.0);
+	Melder_require (cospsi >= -1.0 && cospsi <= 1.0,
+		U"NUMgetEllipseBoundingBox: abs (cospi) should not exceed 1.", cospsi);
 	double width, height;
 	if (cospsi == 1.0) { // a-axis along x-axis
 		width = a;
@@ -1921,33 +1885,24 @@ double NUMminimize_brent (double (*f) (double x, void *closure), double a, doubl
 }
 
 /*
-	probs is probability vector, i.e. all 0 <= probs [i] <= 1 and sum(i=1;i=nprobs, probs [i])= 1
-	p is a probability
+	probs is probability vector, i.e. all 0 <= probs [i] <= 1 and sum(i=1...probs.size, probs [i])= 1
 */
-integer NUMgetIndexFromProbability (constVEC probs, double p) { //TODO HMM zero start matrices
+integer NUMgetIndexFromProbability (constVEC probs, double p) {
 	integer index = 1;
 	double psum = probs [index];
-	while (p > psum && index < probs.size) {
+	while (p > psum && index < probs.size)
 		psum += probs [++ index];
-	}
-	return index;
-}
-integer NUMgetIndexFromProbability (double *probs, integer nprobs, double p) { //TODO HMM zero start matrices
-	integer index = 1;
-	double psum = probs [index];
-	while (p > psum && index < nprobs) {
-		psum += probs [++ index];
-	}
 	return index;
 }
 
 // straight line fitting
 
-void NUMlineFit_theil (constVEC x, constVEC y, double *out_m, double *out_intercept, bool incompleteMethod)
-{
+void NUMlineFit_theil (constVEC x, constVEC y, double *out_m, double *out_intercept, bool incompleteMethod) {
 	try {
-		Melder_assert (x.size == y.size);
-		/* Theil's incomplete method:
+		Melder_require (x.size == y.size,
+			U"NUMlineFit_theil: the sizes of the two vectors should be equal.");
+		/*
+			Theil's incomplete method:
 			Split (x [i],y [i]) as
 			(x [i],y [i]), (x [N+i],y [N=i], i=1..numberOfPoints/2
 			m [i] = (y [N+i]-y [i])/(x [N+i]-x [i])
@@ -1985,15 +1940,18 @@ void NUMlineFit_theil (constVEC x, constVEC y, double *out_m, double *out_interc
 			VECsort_inplace (mbs.part (1, x.size));
 			intercept = NUMquantile (mbs.part (1, x.size), 0.5);
 		}
-		if (out_m) *out_m = m;
-		if (out_intercept) *out_intercept = intercept;
+		if (out_m)
+			*out_m = m;
+		if (out_intercept)
+			*out_intercept = intercept;
 	} catch (MelderError) {
 		Melder_throw (U"No line fit (Theil's method).");
 	}
 }
 
 void NUMlineFit_LS (constVEC x, constVEC y, double *out_m, double *out_intercept) {
-	Melder_assert (x.size == y.size);
+	Melder_require (x.size == y.size,
+		U"NUMlineFit_LS: the sizes of the two vectors should be equal.");
 	double sx = NUMsum (x);
 	double sy = NUMsum (y);
 	double xmean = sx / x.size;
@@ -2005,19 +1963,20 @@ void NUMlineFit_LS (constVEC x, constVEC y, double *out_m, double *out_intercept
 	}
 	// y = m*x + b
 	m /= st2;
-	if (out_intercept) *out_intercept = (sy - m * sx) / x.size;
-	if (out_m) *out_m = m;
+	if (out_intercept)
+		*out_intercept = (sy - m * sx) / x.size;
+	if (out_m)
+		*out_m = m;
 }
 
 void NUMlineFit (constVEC x, constVEC y, double *out_m, double *out_intercept, int method) {
-	if (method == 1) {
+	if (method == 1)
 		NUMlineFit_LS (x, y, out_m, out_intercept);
-	} else if (method == 3) {
+	else if (method == 3)
 		NUMlineFit_theil (x, y, out_m, out_intercept, false);
-	} else {
+	else
 		NUMlineFit_theil (x, y, out_m, out_intercept, true);
 }
-}
 
 // IEEE: Programs for digital signal processing section 4.3 LPTRN
 // lpc [1..n] to rc [1..n]
@@ -2031,11 +1990,10 @@ void VECrc_from_lpc (VEC rc, constVEC lpc) {
 		rc [m] = a [m];
 		Melder_require (fabs (rc [m]) <= 1.0, U"Relection coefficient [", m, U"] larger than 1.");
 		b.part (1, m) <<= a.part (1, m);
-		for (integer i = 1; i < m; i ++) {
+		for (integer i = 1; i < m; i ++)
 			a [i] = (b [i] - rc [m] * b [m - i]) / (1.0 - rc [m] * rc [m]);
 	}
 }
-}
 
 void VEClpc_from_rc (VEC lpc, constVEC rc) {
 	Melder_assert (lpc.size == rc.size);
@@ -2226,7 +2184,8 @@ integer NUMrandomBinomial (double p, integer n) {
 	integer ix;			// return value
 	bool flipped = false;
 
-	if (n == 0) return 0;
+	if (n == 0)
+		return 0;
 	if (p > 0.5) {
 		p = 1.0 - p;	// work with small p
 		flipped = true;
@@ -2257,9 +2216,8 @@ integer NUMrandomBinomial (double p, integer n) {
 			double u = NUMrandomUniform (0.0, 1.0); // djmw gsl_rng_uniform (rng);
 
 			for (ix = 0; ix <= BINV_CUTOFF; ++ ix) {
-				if (u < f) {
+				if (u < f)
 					goto Finish;
-				}
 				u -= f;
 				// Use recursion f(x+1) = f(x)* [(n-x)/(x+1)]* [p/(1-p)]
 				f *= s * (n - ix) / (ix + 1.0);
@@ -2352,23 +2310,20 @@ TryAgain:
 			// Parallelogram region
 			double x = xl + (u - p1) / c;
 			v = v * c + 1.0 - fabs (x - xm) / p1;
-			if (v > 1.0 || v <= 0.0) {
+			if (v > 1.0 || v <= 0.0)
 				goto TryAgain;
-			}
 			ix = (integer) x;
 		} else if (u <= p3) {
 			// Left tail
 			ix = (integer) (xl + log (v) / lambda_l);
-			if (ix < 0) {
+			if (ix < 0)
 				goto TryAgain;
-			}
 			v *= ((u - p2) * lambda_l);
 		} else {
 			// Right tail
 			ix = (integer) (xr - log (v) / lambda_r);
-			if (ix > (double) n) {
+			if (ix > (double) n)
 				goto TryAgain;
-			}
 			v *= ((u - p3) * lambda_r);
 		}
 
@@ -2393,7 +2348,7 @@ TryAgain:
 			More efficient determination of whether v < f(x)/f(M)
 		 */
 
-		integer k = labs (ix - m);
+		integer k = integer_abs (ix - m);
 
 		if (k <= FAR_FROM_MEAN) {
 			/*
@@ -2406,14 +2361,12 @@ TryAgain:
 			var = v;
 
 			if (m < ix) {
-				for (integer i = m + 1; i <= ix; i ++) {
+				for (integer i = m + 1; i <= ix; i ++)
 					f *= (g / i - s);
-				}
 			} else if (m > ix) {
-				for (integer i = ix + 1; i <= m; i ++) {
+				for (integer i = ix + 1; i <= m; i ++)
 					f /= (g / i - s);
 			}
-			}
 
 			accept = f;
 		} else {
@@ -2429,13 +2382,11 @@ TryAgain:
 				*/
 				double amaxp = k / npq * ((k * (k / 3.0 + 0.625) + (1.0 / 6.0)) / npq + 0.5);
 				double ynorm = -(k * k / (2.0 * npq));
-				if (var < ynorm - amaxp) {
+				if (var < ynorm - amaxp)
 					goto Finish;
-				}
-				if (var > ynorm + amaxp) {
+				if (var > ynorm + amaxp)
 					goto TryAgain;
 			}
-			}
 
 			/* 
 				Now, again: do the test log(v) vs. log f(x)/f(M)
@@ -2486,12 +2437,11 @@ TryAgain:
 		}
 
 
-		if (var <= accept) {
+		if (var <= accept)
 			goto Finish;
-		} else {
+		else
 			goto TryAgain;
 	}
-	}
 
 Finish:
 
@@ -2499,12 +2449,11 @@ Finish:
 }
 
 double NUMrandomBinomial_real (double p, integer n) {
-	if (p < 0.0 || p > 1.0 || n < 0) {
+	if (p < 0.0 || p > 1.0 || n < 0)
 		return undefined;
-	} else {
+	else
 		return (double) NUMrandomBinomial (p, n);
 }
-}
 
 double NUMrandomGamma (const double alpha, const double beta) {
 	Melder_require (alpha > 0 && beta > 0,
@@ -2591,8 +2540,7 @@ void NUMfixIndicesInRange (integer lowerLimit, integer upperLimit, integer *lowI
 	}
 }
 
-void NUMgetEntropies (constMAT m, double *out_h, double *out_hx, 
-	double *out_hy,	double *out_hygx, double *out_hxgy, double *out_uygx, double *out_uxgy, double *out_uxy) {
+void NUMgetEntropies (constMATVU const& m, double *out_h, double *out_hx, double *out_hy, double *out_hygx, double *out_hxgy, double *out_uygx, double *out_uxgy, double *out_uxy) {
 	
 	double h = undefined, hx = undefined, hy = undefined;
 	double hxgy = undefined, hygx = undefined, uygx = undefined, uxgy = undefined, uxy = undefined;
@@ -2621,8 +2569,7 @@ void NUMgetEntropies (constMAT m, double *out_h, double *out_hx,
 		
 		longdouble hx_t = 0.0;
 		for (integer j = 1; j <= m.ncol; j ++) {
-			longdouble colsum = 0.0;
-			for (integer i = 1; i <= m.nrow; i++) colsum += m [i] [j];
+			double colsum = NUMsum (m.column (j));
 			if (colsum > 0.0) {
 				longdouble p = colsum / totalSum;
 				hx_t -= p * NUMlog2 (p);
@@ -2667,27 +2614,6 @@ void NUMgetEntropies (constMAT m, double *out_h, double *out_hx,
 }
 #undef TINY
 
-double NUMfrobeniusnorm (constMAT x) {
-	longdouble scale = 0.0;
-	longdouble ssq = 1.0;
-	for (integer i = 1; i <= x.nrow; i ++) {
-		for (integer j = 1; j <= x.ncol; j ++) {
-			if (x [i] [j] != 0.0) {
-				longdouble absxi = fabs (x [i] [j]);
-				if (scale < absxi) {
-					longdouble t = scale / absxi;
-					ssq = 1.0 + ssq * t * t;
-					scale = absxi;
-				} else {
-					longdouble t = absxi / scale;
-					ssq += t * t;
-				}
-			}
-		}
-	}
-	return scale * sqrt ((double) ssq);
-}
-
 double NUMtrace (const constMATVU& a) {
 	Melder_assert (a.nrow == a.ncol);
 	longdouble trace = 0.0;
@@ -2762,15 +2688,15 @@ void NUMeigencmp22 (double a, double b, double c, double *out_rt1, double *out_r
 	longdouble acs = fabs (cs), cs1, sn1;
 	if (acs > ab) {
 		longdouble ct = -tb / cs;
-		sn1 = 1 / sqrt (1 + ct * ct);
+		sn1 = 1.0 / sqrt (1.0 + ct * ct);
 		cs1 = ct * sn1;
 	} else {
 		if (ab == 0) {
-			cs1 = 1;
-			sn1 = 0;
+			cs1 = 1.0;
+			sn1 = 0.0;
 		} else {
 			tn = -cs / tb;
-			cs1 = 1 / sqrt (1 + tn * tn);
+			cs1 = 1.0 / sqrt (1.0 + tn * tn);
 			sn1 = tn * cs1;
 		}
 	}
@@ -2779,6 +2705,10 @@ void NUMeigencmp22 (double a, double b, double c, double *out_rt1, double *out_r
 		cs1 = -sn1;
 		sn1 = tn;
 	}
+	if (fabs (cs1) > 1.0)
+		cs1 = copysign (1.0, cs1);
+	if (fabs (sn1) > 1.0)
+		sn1 = copysign (1.0, sn1);
 	if (out_rt1)
 		*out_rt1 = (double) rt1;
 	if (out_rt2)

dwsys/NUM2.h

@@ -439,6 +439,11 @@ void NUMrank (vector<T> a) {
 }
 
 void MATlowerCholeskyInverse_inplace (MAT a, double *out_lnd);
+inline autoMAT newMATlowerCholeskyInverse (constMAT const& a) {
+	autoMAT result = newMATcopy (a);
+	MATlowerCholeskyInverse_inplace (result.get(), nullptr);
+	return result;
+}
 /*
 	Calculates L^-1, where A = L.L' is a symmetric positive definite matrix
 	and ln(determinant). L^-1 in lower, leave upper part intact.
@@ -455,7 +460,7 @@ double NUMdeterminant_fromSymmetricMatrix (constMAT m);
 	ln(determinant) of a symmetric p.s.d. matrix
 */
 
-double NUMmahalanobisDistance (constMAT lowerInverse, constVEC v, constVEC m);
+double NUMmahalanobisDistanceSquared (constMAT lowerInverse, constVEC v, constVEC m);
 /*
 	Calculates squared Mahalanobis distance: (v-m)'S^-1(v-m).
 	Input matrix (li) is the inverse L^-1 of the Cholesky decomposition S = L.L'
@@ -730,11 +735,6 @@ double NUMinvTukeyQ (double p, double cc, double df, double rr);
  *  df = degrees of freedom of error term
  */
 
-double NUMnormalityTest_HenzeZirkler (constMAT data, double *inout_beta, double *out_tnb, double *out_lnmu, double *out_lnvar);
-/*
-	Multivariate normality test of nxp data matrix according to the method described in Henze & Wagner (1997).
-	The test statistic is returned in tnb, together with the lognormal mean 'lnmu' and the lognormal variance 'lnvar'.
-*/
 
 /******  Frequency in Hz to other frequency reps ****/
 
@@ -1252,14 +1252,9 @@ void VECarea_from_lpc (VEC area, constVEC lpc);
 */
 void NUMfixIndicesInRange (integer lowerLimit, integer upperLimit, integer *lowIndex, integer *highIndex);
 
-void NUMgetEntropies (constMAT m, double *out_h, double *out_hx, 
+void NUMgetEntropies (constMATVU const& m, double *out_h, double *out_hx, 
 	double *out_hy,	double *out_hygx, double *out_hxgy, double *out_uygx, double *out_uxgy, double *out_uxy);
 
-double NUMfrobeniusnorm (constMAT x);
-/*
-	Returns frobenius norm of matrix sqrt (sum (i=1:nrow, j=1:ncol, x[i][j]^2))
-*/
-
 inline double NUMmean_weighted (constVEC x, constVEC w) {
 	Melder_assert (x.size == w.size);
 	double inproduct = NUMinner (x, w);

dwsys/NUMcblas.cpp

@@ -1111,7 +1111,7 @@ L10:
 			}
 
 		} else if (ngpmin == gpmin && ngnmin == gnmin) {
-			if ( (i__1 = ngpmin - ngnmin, labs (i__1)) == 1) {
+			if ( (i__1 = ngpmin - ngnmin, integer_abs (i__1)) == 1) {
 				lemin = MAX (ngpmin, ngnmin);
 				/* ( Twos-complement machines, no gradual underflow; e.g.,
 				   CYBER 205 ) */
@@ -1121,7 +1121,7 @@ L10:
 				iwarn = true;
 			}
 
-		} else if ( (i__1 = ngpmin - ngnmin, labs (i__1)) == 1 && gpmin == gnmin) {
+		} else if ( (i__1 = ngpmin - ngnmin, integer_abs (i__1)) == 1 && gpmin == gnmin) {
 			if (gpmin - MIN (ngpmin, ngnmin) == 3) {
 				lemin = MAX (ngpmin, ngnmin) - 1 + lt;
 				/* ( Twos-complement machines with gradual underflow; no

dwtest/test_Electroglottogram.praat

+# test_Electroglottogram.praat
+# djmw 20190402
+
+appendInfoLine : "test_Electroglottogram.praat"
+
+pitchFloor = 50
+pitchCeiling = 500
+closingThreshold = 0.30
+silenceThreshold = 0.03
+sound = Read from file: "s_egg_test.wav"
+egg = Extract Electroglottogram: 1
+intervalTier = To IntervalTier: pitchFloor, pitchCeiling, closingThreshold, silenceThreshold
+selectObject: egg
+degg = To Electroglottogram (derivative): 5000, 50
+selectObject: egg
+To AmplitudeTier (levels): pitchFloor, pitchCeiling, closingThreshold, "yes", "yes"
+peaks = selected ("AmplitudeTier", 1)
+valleys = selected ("AmplitudeTier", 2)
+levels = selected ("AmplitudeTier", 3)
+selectObject: sound, egg, degg
+combined = Combine to multi-channel
+
+removeObject: levels, valleys, peaks, degg, intervalTier, egg, sound, combined
+
+appendInfoLine : "test_Electroglottogram.praat OK"

dwtest/test_GaussianMixture.praat

@@ -4,7 +4,7 @@
 appendInfoLine: "test_GaussianMixture.praat"
 appendInfoLine: tab$, "Interface"
 pols = Create TableOfReal (Pols 1973): "no"
-gm = To GaussianMixture (row labels): "Complete"
+gm = To GaussianMixture (row labels): "Diagonals"
 numberOfComponents = Get number of components
 assert numberOfComponents == 12
 dimension = Get dimension of component
@@ -45,20 +45,20 @@ Draw marginal pdf: 1, 0, 0, 0, 0, 500, 0, "yes"
 Select outer viewport: 3, 6, 3, 6
 Draw marginal pdf: 2, 0, 0, 0, 0, 500, 0, "yes"
 
-appendInfoLine: tab$, "GaussinaMixture and TableOfReal"
+appendInfoLine: tab$, "GaussianMixture and TableOfReal"
 
 selectObject: pols, gm
 p = Get likelihood value: "Likelihood"
-cemm = To GaussianMixture (CEMM): 1, 0.001, 200, 0.001, "MessageLength"
+cemm = To GaussianMixture (CEMM): 1, 0.001, 200, "yes"
 selectObject: pols, gm
 ct = To ClassificationTable
 selectObject: pols, gm
 correlation = To Correlation
 selectObject: pols, gm
-tor = To TableOfReal (BHEP normality tests): 1.41
-removeObject:  tor, correlation, ct, cemm
+table = To Table (BHEP normality tests): 1.41
+removeObject:  table, correlation, ct, cemm
 
-;removeObject: pols, gm
+removeObject: pols, gm
 
 appendInfoLine: "test_GaussianMixture.praat OK"
 

dwtest/test_Matrix_solve.praat

@@ -59,7 +59,7 @@ procedure solve_undetermined: .nrow, .ncol
 		.c = Get value in cell: 1, .irow
 		assert .c > 1-1e-7 and .c < 1+1e-7
 	endfor
-	;removeObject: .m, .ms
+	removeObject: .m, .ms
 endproc
 
 # test for several dimensions

dwtest/test_bss_twoSoundsMixed.praat

 # bss_twoSoundsMixed.praat
 # test_bss_example_manual.praat
 # djmw 20151030
+# ppgb 20190811
 
 Erase all
 
@@ -8,6 +9,8 @@ synth = Create SpeechSynthesizer: "English (Great Britain)", "Female1"
 s1 = To Sound: "This is some text.", "no"
 selectObject: synth
 Set speech output settings: 44100, 0.01, 80, 50, 145, "no", "IPA"
+Speech output settings: 44100, 0.01, 1.2, 1.0, 145, "IPA"
+Estimate speech rate from speech: "no"
 
 s2 =To Sound: "abracadabra, abra", "no"
 plusObject: s1

dwtest/test_henzeZirklerMultivariateNormalityTest.praat

@@ -16,7 +16,7 @@ lvariance = extractNumber (report$, " variance:")
 assert abs (lvariance - 0.13790691675095398)<tol
 beta = extractNumber (report$, "Smoothing:")
 printline 'beta'
-assert abs (beta -0.5541226910830553)<tol
+assert abs (beta -1.27608342442082)<tol
 sampleSize = extractNumber (report$, " size:")
 assert sampleSize = 50
 nvars = extractNumber (report$, " variables:")

dwtools/CCA.cpp

@@ -110,8 +110,8 @@ autoCCA TableOfReal_to_CCA (TableOfReal me, integer numberOfDependents) {
 		autoSVD svdx = SVD_create (numberOfObservations, numberOfIndependents);	 // numberOfObservations >= numberOfIndependents, hence no transposition
 		svdy -> u.all() <<= my data.verticalBand (1, numberOfDependents);
 		svdx -> u.all() <<= my data.verticalBand (numberOfDependents + 1, my numberOfColumns);
-		double fnormy = NUMfrobeniusnorm (svdy -> u.get());
-		double fnormx = NUMfrobeniusnorm (svdx -> u.get());
+		double fnormy = NUMnorm (svdy -> u.get(), 2.0);
+		double fnormx = NUMnorm (svdx -> u.get(), 2.0);
 		
 		Melder_require (fnormy > 0.0 && fnormx > 0.0,
 			U"One of the parts of the table contains only zeros.");

dwtools/Discriminant.cpp

@@ -590,7 +590,7 @@ autoClassificationTable Discriminant_TableOfReal_to_ClassificationTable (Discrim
 			double norm = 0.0, pt_max = -1e308;
 			for (integer j = 1; j <= numberOfGroups; j ++) {
 				SSCP t = groups->at [j];
-				double md = NUMmahalanobisDistance (sscpvec [j] -> data.get(), thy data.row (i), t -> centroid.get());
+				double md = NUMmahalanobisDistanceSquared (sscpvec [j] -> data.get(), thy data.row (i), t -> centroid.get());
 				double pt = log_apriori [j] - 0.5 * (ln_determinant [j] + md);
 				if (pt > pt_max)
 					pt_max = pt;
@@ -725,7 +725,7 @@ autoClassificationTable Discriminant_TableOfReal_to_ClassificationTable_dw (Disc
 				x [k] = thy data [i] [k] + displacement [k];
 			for (integer j = 1; j <= g; j ++) {
 				SSCP t = groups->at [j];
-				double md = NUMmahalanobisDistance (sscpvec [j] -> data.get(), x.get(), t -> centroid.get());
+				double md = NUMmahalanobisDistanceSquared (sscpvec [j] -> data.get(), x.get(), t -> centroid.get());
 				double pt = log_apriori [j] - 0.5 * (ln_determinant [j] + md);
 				if (pt > pt_max) {
 					pt_max = pt;