LPC/LPC_and_LineSpectralFrequencies.cpp

 /* LPC_and_LineSpectralFrequencies.cpp
  *
- * Copyright (C) 2016-2018 David Weenink
+ * Copyright (C) 2016-2019 David Weenink
  *
  * This program is free software; you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
@@ -31,11 +31,11 @@
 	From: Joseph Rothweiler (1999), "On Polynomial Reduction in the Computation of LSP Frequencies." 
 	IEEE Trans. on ASSP 7, 592--594.
 */
-static void cos2x (double *g, integer order) {
-	for (integer i = 2; i <= order; i ++) {
-		for (integer j = order; j > i; j --) {
+
+static void cos2x (VECVU const& g) {
+	for (integer i = 3; i <= g.size; i ++) {
+		for (integer j = g.size; j > i; j --)
 			g [j - 2] -= g [j];
-		}
 		g [i - 2] -= 2.0 * g [i];
 	}
 }
@@ -45,50 +45,48 @@ static void Polynomial_fromLPC_Frame_lspsum (Polynomial me, LPC_Frame lpc) {
 	
 	integer order = lpc -> nCoefficients, g_order = (order + 1) / 2; // orders
 	my coefficients [order + 2] = 1.0;
-	for (integer i = 1; i <= order; i ++) {
+	for (integer i = 1; i <= order; i ++)
 		my coefficients [order + 2 - i] = lpc -> a [i] + lpc -> a [order + 1 - i];
-	}
+
 	my coefficients [1] = 1.0;
 	my numberOfCoefficients = order + 2;
 
-	if (order % 2 == 0) { // order even
+	if (order % 2 == 0) // order even
 		Polynomial_divide_firstOrderFactor (me, -1.0, nullptr);
-	}
 	
 	// transform to cos(w) terms
-	for (integer i = 1; i <= g_order + 1; i ++) {
+	for (integer i = 1; i <= g_order + 1; i ++)
 		my coefficients [i] = my coefficients [g_order + i];
-	}
+
 	my numberOfCoefficients = g_order + 1;
 	
 	// to Chebychev
-	cos2x (& my coefficients [1], g_order);
+	cos2x (my coefficients.part(1, my numberOfCoefficients));
 }
 
 static void Polynomial_fromLPC_Frame_lspdif (Polynomial me, LPC_Frame lpc) {
 	// Fa (z) = A(z) - z^-(p+1)A(1/z)
 	integer order = lpc -> nCoefficients;
 	my coefficients [order + 2] = -1.0;
-	for (integer i = 1; i <= order; i ++) {
+	for (integer i = 1; i <= order; i ++)
 		my coefficients [order + 2 - i] = - lpc -> a [i] + lpc -> a [order + 1 - i];
-	}
+
 	my coefficients [1] = 1.0;
 	my numberOfCoefficients = order + 2;
 
-	if (order % 2 == 0) { // Fa(z)/(z-1)
+	if (order % 2 == 0) // Fa(z)/(z-1)
 		Polynomial_divide_firstOrderFactor (me, 1.0, nullptr);
-	} else { // Fa(z) / (z^2 - 1)
+	else // Fa(z) / (z^2 - 1)
 		Polynomial_divide_secondOrderFactor (me, 1.0);
-	}
 	
 	// transform to cos(w) terms
 	integer g_order = my numberOfCoefficients / 2;
-	for (integer i = 1; i <= g_order + 1; i ++) {
+	for (integer i = 1; i <= g_order + 1; i ++)
 		my coefficients [i] = my coefficients [g_order + i];
-	}
+
 	my numberOfCoefficients = g_order + 1;
 	// to Chebychev
-	cos2x (& my coefficients [1], g_order);
+	cos2x (my coefficients.part(1, my numberOfCoefficients));
 }
 
 #if 0
@@ -205,16 +203,16 @@ static void LineSpectralFrequencies_Frame_initFromLPC_Frame_grid (LineSpectralFr
 	integer numberOfBisections = 0, numberOfRootsFound = 0;
 	while (numberOfRootsFound  < half_order_g1 && numberOfBisections < 10) {
 		numberOfRootsFound = Roots_fromPolynomial_grid (roots, g1, gridSize);
-		gridSize *= 0.5; numberOfBisections++;
+		gridSize *= 0.5;
+		numberOfBisections++;
 	}
 	
 	Melder_require (numberOfBisections < 10, U"Too many bisections.");
 	
 	// [g1-> xmin, g1 -> xmax] <==> [nyquistFrequency, 0] i.e. highest root corresponds to lowest frequency
 	
-	for (integer i = 1; i <= half_order_g1; i ++) {
+	for (integer i = 1; i <= half_order_g1; i ++)
 		my frequencies [2 * i - 1] = acos (roots -> v [half_order_g1 + 1 - i].re / 2.0) / NUMpi * maximumFrequency; 
-	}
 	
 	// the roots of g2 lie inbetween the roots of g1
 	
@@ -222,19 +220,18 @@ static void LineSpectralFrequencies_Frame_initFromLPC_Frame_grid (LineSpectralFr
 		double xmax = roots -> v [half_order_g1 + 1 - i].re;
 		double xmin = i == half_order_g1 ? g1 -> xmin : roots -> v [half_order_g1 - i].re;
 		double root = Polynomial_findOneSimpleRealRoot_ridders (g2, xmin, xmax);
-		if (isdefined (root)) {
+		if (isdefined (root))
 			my frequencies [2 * i] = acos (root / 2.0) / NUMpi * maximumFrequency;
-		} else { 
+		else
 			my numberOfFrequencies --;
 	}
 }
-}
 
 autoLineSpectralFrequencies LPC_to_LineSpectralFrequencies (LPC me, double gridSize) {
 	try {
-		if (gridSize == 0.0) {
+		if (gridSize == 0.0)
 			gridSize = 0.02;
-		}
+
 		double nyquistFrequency = 0.5 / my samplingPeriod;
 		autoLineSpectralFrequencies thee = LineSpectralFrequencies_create (my xmin, my xmax, my nx, my dx, my x1, my maxnCoefficients, nyquistFrequency);
 		autoPolynomial g1 = Polynomial_create (-2.0, 2.0, my maxnCoefficients + 1); // large enough
@@ -272,7 +269,7 @@ static void LPC_Frame_initFromLineSpectralFrequencies_Frame (LPC_Frame me, LineS
 		double omega = thy frequencies [2 * i -1] / maximumFrequency * NUMpi;
 		my a [i] = -2.0 * cos (omega);
 	}
-	Polynomial_initFromProductOfSecondOrderTerms (fs, constVEC (my a.at, numberOfOmegas));
+	Polynomial_initFromProductOfSecondOrderTerms (fs, my a.part (1, numberOfOmegas));
 
 	/*
 		Reconstruct Fa (z)
@@ -282,7 +279,7 @@ static void LPC_Frame_initFromLineSpectralFrequencies_Frame (LPC_Frame me, LineS
 		double omega = thy frequencies [2 * i] / maximumFrequency * NUMpi;
 		my a [i] = -2.0 * cos (omega);
 	}
-	Polynomial_initFromProductOfSecondOrderTerms (fa, constVEC (my a.at, numberOfOmegas));
+	Polynomial_initFromProductOfSecondOrderTerms (fa, my a.part (1, numberOfOmegas));
 	
 	if (thy numberOfFrequencies % 2 == 0) {
 		Polynomial_multiply_firstOrderFactor (fs, -1.0);   // * (z + 1)
@@ -294,10 +291,9 @@ static void LPC_Frame_initFromLineSpectralFrequencies_Frame (LPC_Frame me, LineS
 	/*
 		A(z) = (Fs(z) + Fa(z) / 2
 	*/
-	for (integer i = 1; i <= fs -> numberOfCoefficients - 2; i ++) {
+	for (integer i = 1; i <= fs -> numberOfCoefficients - 2; i ++)
 		my a [thy numberOfFrequencies - i + 1] = 0.5 * (fs -> coefficients [i + 1] + fa -> coefficients [i + 1]);
 }
-}
 
 autoLPC LineSpectralFrequencies_to_LPC (LineSpectralFrequencies me) {
 	try {

LPC/LPC_to_Spectrogram.cpp

 /* LPC_to_Spectrogram.cpp
  *
- * Copyright (C) 1994-2018 David Weenink
+ * Copyright (C) 1994-2019 David Weenink
  *
  * This code is free software; you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
@@ -31,7 +31,8 @@ autoSpectrogram LPC_to_Spectrogram (LPC me, double dfMin, double bandwidthReduct
 			nfft = 512; 
 			dfMin = samplingFrequency / nfft;
 		}
-		while (samplingFrequency / nfft > dfMin || nfft <= my maxnCoefficients) nfft *= 2;
+		while (samplingFrequency / nfft > dfMin || nfft <= my maxnCoefficients)
+			nfft *= 2;
 		double freqStep = samplingFrequency / nfft;
 
 		autoSpectrogram thee = Spectrogram_create (my xmin, my xmax, my nx, my dx, my x1, 0.0, samplingFrequency / 2.0, nfft / 2 + 1, freqStep, 0.0);

LPC/VocalTractTier.cpp

@@ -92,22 +92,16 @@ autoVocalTractTier VocalTractTier_create (double fromTime, double toTime) {
 autoVocalTractTier VocalTract_to_VocalTractTier (VocalTract me, double startTime, double endTime, double time) {
 	try {
 		autoVocalTractTier thee = VocalTractTier_create (startTime, endTime);
-		VocalTractTier_addVocalTract_copy (thee.get(), time, me);
+		VocalTractTier_addVocalTract (thee.get(), time, me);
 		return thee;
 	} catch (MelderError) {
 		Melder_throw (me, U": not converted to VocalTractTier");
 	}
 }
 
-void VocalTractTier_addVocalTract_copy (VocalTractTier me, double time, VocalTract vocaltract) {
+void VocalTractTier_addVocalTract (VocalTractTier me, double time, VocalTract vocaltract) {
 	try {
 		autoVocalTractPoint thee = VocalTract_to_VocalTractPoint (vocaltract, time);
-		if (my d_vocalTracts.size > 0) {
-			VocalTractPoint vtp = my d_vocalTracts.at [1];
-			integer numberOfSections = vtp -> d_vocalTract -> nx;
-			Melder_require (numberOfSections == vocaltract -> nx,
-				U"The number of sections should be equal to ", numberOfSections, U".");
-		}
 		my d_vocalTracts. addItem_move (thee.move());
 	} catch (MelderError) {
 		Melder_throw (me, U": no VocalTract added.");
@@ -137,9 +131,8 @@ autoVocalTract VocalTractTier_to_VocalTract (VocalTractTier me, double time) {
 
 autoLPC VocalTractTier_to_LPC (VocalTractTier me, double timeStep) {
 	try {
-		if (my d_vocalTracts.size == 0) {
-			Melder_throw (U"Empty VocalTractTier");
-		}
+		Melder_require (my d_vocalTracts.size > 0.0,
+			U"The VocalTractTier should not be empty.");
 		integer numberOfFrames = Melder_ifloor ((my xmax - my xmin) / timeStep);
 		VocalTractPoint vtp = my d_vocalTracts.at [1];
 		integer numberOfSections = vtp -> d_vocalTract -> nx;

LPC/VocalTractTier.h

@@ -2,7 +2,7 @@
 #define _VocalTractTier_h_
 /* VocalTractTier.h
  *
- * Copyright (C) 2012, 2015 David Weenink
+ * Copyright (C) 2012-2019 David Weenink
  *
  * This code is free software; you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
@@ -32,7 +32,7 @@ autoVocalTractTier VocalTractTier_create (double fromTime, double toTime);
 
 autoVocalTractTier VocalTract_to_VocalTractTier (VocalTract me, double startTime, double endTime, double time);
 
-void VocalTractTier_addVocalTract_copy (VocalTractTier me, double time, VocalTract thee);
+void VocalTractTier_addVocalTract (VocalTractTier me, double time, VocalTract thee);
 
 autoLPC VocalTractTier_to_LPC (VocalTractTier me, double timeStep);
 

LPC/praat_LPC_init.cpp

 /* praat_LPC_init.cpp
  *
- * Copyright (C) 1994-2016 David Weenink
+ * Copyright (C) 1994-2019 David Weenink
  *
  * This code is free software; you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
@@ -962,7 +962,7 @@ FORM (MODIFY_VocalTractTier_addVocalTract, U"VocalTractTier: Add VocalTract", nu
 	OK
 DO
 	MODIFY_FIRST_OF_TWO (VocalTractTier, VocalTract)
-		VocalTractTier_addVocalTract_copy (me, time, you);
+		VocalTractTier_addVocalTract (me, time, you);
 	MODIFY_FIRST_OF_TWO_END
 }
 

artsynth/Artword_Speaker.cpp

 /* Artword_Speaker.cpp
  *
- * Copyright (C) 1992-2005,2011,2015-2017 Paul Boersma
+ * Copyright (C) 1992-2005,2011,2015-2017,2019 Paul Boersma
  *
  * This code is free software; you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
@@ -30,7 +30,7 @@ void Artword_Speaker_draw (Artword artword, Speaker speaker, Graphics g, int num
 	Graphics_setLineWidth (g, oldLineWidth);
 }
 
-void Artword_Speaker_movie (Artword artword, Speaker speaker, Graphics graphics) {
+void Artword_Speaker_playMovie (Artword artword, Speaker speaker, Graphics graphics) {
 	constexpr double timeStep = 0.03;
 	autoArt art = Art_create ();
 	for (double tim = 0.0; tim < artword -> totalTime; tim += timeStep) {

artsynth/Artword_Speaker.h

 /* Artword_Speaker.h
  *
- * Copyright (C) 1992-2005,2011,2016,2017 Paul Boersma
+ * Copyright (C) 1992-2005,2011,2016,2017,2019 Paul Boersma
  *
  * This code is free software; you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
@@ -22,6 +22,6 @@
 
 void Artword_Speaker_draw (Artword artword, Speaker speaker, Graphics g, int numberOfSteps);
 
-void Artword_Speaker_movie (Artword artword, Speaker speaker, Graphics g);
+void Artword_Speaker_playMovie (Artword artword, Speaker speaker, Graphics g);
 
 /* End of file Artword_Speaker.h */

artsynth/Artword_Speaker_Sound.cpp

 /* Artword_Speaker_Sound.cpp
  *
- * Copyright (C) 1992-2005,2006,2011,2015-2017 Paul Boersma
+ * Copyright (C) 1992-2005,2006,2011,2015-2017,2019 Paul Boersma
  *
  * This code is free software; you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
@@ -40,7 +40,7 @@ static int playCallback (Artword_Speaker_Sound_PlayInfo me, int /* phase */, dou
 	return 1;
 }
 
-void Artword_Speaker_Sound_movie (Artword artword, Speaker speaker, Sound sound, Graphics graphics) {
+void Artword_Speaker_Sound_playMovie (Artword artword, Speaker speaker, Sound sound, Graphics graphics) {
 	try {
 		static autoArtword_Speaker_Sound_PlayInfo info;
 		if (! info)

artsynth/Artword_Speaker_Sound.h

 /* Artword_Speaker_Sound.h
  *
- * Copyright (C) 1992-2005,2011,2016-2017 Paul Boersma
+ * Copyright (C) 1992-2005,2011,2016-2017,2019 Paul Boersma
  *
  * This code is free software; you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
@@ -20,7 +20,7 @@
 #include "Speaker.h"
 #include "Sound.h"
 
-void Artword_Speaker_Sound_movie
+void Artword_Speaker_Sound_playMovie
 	(Artword artword, Speaker speaker, Sound sound, Graphics graphics);
 
 /* End of file Artword_Speaker_Sound.h */

artsynth/praat_Artsynth.cpp

 /* praat_Artsynth.cpp
  *
- * Copyright (C) 1992-2009,2011,2012,2014-2018 Paul Boersma
+ * Copyright (C) 1992-2009,2011,2012,2014-2019 Paul Boersma
  *
  * This code is free software; you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
@@ -198,17 +198,17 @@ DO
 	END
 }
 
-DIRECT (MOVIE_Artword_Speaker_movie) {
+DIRECT (MOVIE_Artword_Speaker_playMovie) {
 	MOVIE_TWO (Artword, Speaker, U"Artword & Speaker movie", 300, 300)
-		Artword_Speaker_movie (me, you, graphics);
+		Artword_Speaker_playMovie (me, you, graphics);
 	MOVIE_TWO_END
 }
 
 // MARK: - ARTWORD & SPEAKER & SOUND
 
-DIRECT (MOVIE_Artword_Speaker_Sound_movie) {
+DIRECT (MOVIE_Artword_Speaker_Sound_playMovie) {
 	MOVIE_THREE (Artword, Speaker, Sound, U"Artword & Speaker & Sound movie", 300, 300)
-		Artword_Speaker_Sound_movie (me, you, him, graphics);
+		Artword_Speaker_Sound_playMovie (me, you, him, graphics);
 	MOVIE_THREE_END
 }
 
@@ -353,13 +353,15 @@ void praat_uvafon_Artsynth_init () {
 	praat_addAction2 (classArt, 1, classSpeaker, 1, U"Synthesize", nullptr, 0, nullptr);
 	praat_addAction2 (classArt, 1, classSpeaker, 1, U"To VocalTract", nullptr, 0, NEW1_Art_Speaker_to_VocalTract);
 
-	praat_addAction2 (classArtword, 1, classSpeaker, 1, U"Movie", nullptr, 0, MOVIE_Artword_Speaker_movie);
+	praat_addAction2 (classArtword, 1, classSpeaker, 1, U"Play movie", nullptr, 0, MOVIE_Artword_Speaker_playMovie);
+	praat_addAction2 (classArtword, 1, classSpeaker, 1, U"Movie", nullptr, praat_HIDDEN, MOVIE_Artword_Speaker_playMovie);
 	praat_addAction2 (classArtword, 1, classSpeaker, 1, U"Draw", nullptr, 0, nullptr);
 	praat_addAction2 (classArtword, 1, classSpeaker, 1, U"Draw...", nullptr, 0, GRAPHICS_Artword_Speaker_draw);
 	praat_addAction2 (classArtword, 1, classSpeaker, 1, U"Synthesize", nullptr, 0, nullptr);
 	praat_addAction2 (classArtword, 1, classSpeaker, 1, U"To Sound...", nullptr, 0, NEW1_Artword_Speaker_to_Sound);
 
-	praat_addAction3 (classArtword, 1, classSpeaker, 1, classSound, 1, U"Movie", nullptr, 0, MOVIE_Artword_Speaker_Sound_movie);
+	praat_addAction3 (classArtword, 1, classSpeaker, 1, classSound, 1, U"Play movie", nullptr, 0, MOVIE_Artword_Speaker_Sound_playMovie);
+	praat_addAction3 (classArtword, 1, classSpeaker, 1, classSound, 1, U"Movie", nullptr, praat_HIDDEN, MOVIE_Artword_Speaker_Sound_playMovie);
 
 	praat_addAction1 (classSpeaker, 0, U"Speaker help", nullptr, 0, HELP_Speaker_help);
 

debian/What_s_new_.html

@@ -7,6 +7,20 @@ What's new?
 <p>
 Latest changes in Praat.</p>
 <p>
+<b>6.1</b> (13 July 2019)</p>
+<p>
+<b>6.0.57</b> (1 July 2019)</p>
+<ul>
+<li>
+&nbsp;Repaired a crash when an editor window zoomed in to a zero time range.
+<li>
+&nbsp;Repaired a crash when playing a matrix movie with one row or a constant value.
+<li>
+&nbsp;Repaired a crash when reconstructing a TableOfReal from a PCA.
+<li>
+&nbsp;Scripting: allow comments after `else`, `endif`, `endfor`...
+</ul>
+<p>
 <b>6.0.56</b> (20 June 2019)</p>
 <ul>
 <li>
@@ -466,208 +480,12 @@ Latest changes in Praat.</p>
 <li>
 &nbsp;Windows: fixed a bug that caused Praat to crash at start-up in some Windows versions.
 </ul>
-<p>
-<b>6.0</b> (28 October 2015)</p>
-<ul>
-<li>
-&nbsp;Linux: use PulseAudio as the server for playing sound.
-<li>
-&nbsp;Windows: removed Praatcon.exe, because Praat.exe can now be used from the command line.
-<li>
-&nbsp;Praat can now open files from the command line on all platforms.
-<li>
-&nbsp;The API change that requires raising the major Praat version to 6: no longer any need for Praatcon.exe, because Praat itself can now run a script from the command line on all platforms.
-<li>
-&nbsp;See <a href="Scripting_6_9__Calling_from_the_command_line.html">Scripting 6.9. Calling from the command line</a> for all command line options.
-<li>
-&nbsp;SpeechSynthesizer: corrected a potentially crashing bug.
-<li>
-&nbsp;Mac Retina graphics: made function drawing (e.g. in the TextGrid window) twice as fast.
-<li>
-&nbsp;LongSound: corrected playing of more than 2 channels.
-</ul>
-<p>
-<b>5.4.22</b> (8 October 2015)</p>
-<ul>
-<li>
-&nbsp;Sound files: added support for MP3 files with ID3 version 2 headers.
-<li>
-&nbsp;Table: Line graph where: removed a bug that caused nothing to be drawn if the number of columns was larger than the number of rows and the column index was larger than the number of rows.
-<li>
-&nbsp;Mac: corrected horizontal spacing in graphical text on El Capit&#225;n.
-<li>
-&nbsp;Mac: corrected a bug that would cause Praat to crash if a PNG file could not be saved (from a script).
-<li>
-&nbsp;Windows: enabled PNG saving from Praatcon.
-</ul>
-<p>
-<b>5.4.21</b> (29 September 2015)</p>
-<ul>
-<li>
-&nbsp;Corrected a bug introduced in 5.4.20 that could cause the Sound window to crash.
-<li>
-&nbsp;Repaired a very old hang with pulses in the Sound window.
-</ul>
-<p>
-<b>5.4.20</b> (26 September 2015)</p>
-<ul>
-<li>
-&nbsp;Corrected the graphical representation of numbers with exponents (wrong since 5.4.10).
-<li>
-&nbsp;Windows: prevented flashing console window in <code>runScript()</code>.
-</ul>
-<p>
-<b>5.4.19</b> (16 September 2015)</p>
-<ul>
-<li>
-&nbsp;LongSound: corrected saving of 8-bit files as 16-bit.
-</ul>
-<p>
-<b>5.4.18</b> (7 September 2015)</p>
-<ul>
-<li>
-&nbsp;Regular expressions: if the pattern is illegal, you'll now get an error message instead of an empty string.
-<li>
-&nbsp;LongSound: removed a bug introduced in 1999 that could play long stereo sounds much too fast if they had unusual sampling frequencies. This bug has become more noticeable on the Mac since January 2015.
-<li>
-&nbsp;Guard better against sound files that contain 0 samples.
-<li>
-&nbsp;Scripting: faster look-up of variables.
-</ul>
-<p>
-<b>5.4.17</b> (20 August 2015)</p>
-<ul>
-<li>
-&nbsp;Scripting: removed a bug that could cause Praat to crash when calling <code>writeFile()</code> with an empty text.
-</ul>
-<p>
-<b>5.4.16</b> (16 August 2015)</p>
-<ul>
-<li>
-&nbsp;Mac 64-bit: removed a bug by which graphical text passed over the right edge in TextGrid intervals.
-<li>
-&nbsp;Mac 64-bit: removed a bug that allowed the general "Zoom" and "Close" commands to be accessible from a script.
-</ul>
-<p>
-<b>5.4.15</b> (1 August 2015)</p>
-<ul>
-<li>
-&nbsp;Introduced several command line options that make it easier to call Praat from a web server.
-<li>
-&nbsp;64-bit Mac: removed a bug introduced in 5.4.11 that caused incorrect spacing of vertical text in picture files.
-</ul>
-<p>
-<b>5.4.14</b> (24 July 2015)</p>
-<ul>
-<li>
-&nbsp;Windows and 32-bit Mac: removed a bug introduced in 5.4.12 that could cause weird behavior (even crashing) when pressing the Tab key in the script window or in the manual.
-</ul>
-<p>
-<b>5.4.13</b> (22 July 2015)</p>
-<ul>
-<li>
-&nbsp;The Demo window can now save itself to a PNG or PDF file.
-<li>
-&nbsp;Windows: removed a bug introduced in 5.4.10 whereby Praatcon would write garbled text to the console.
-</ul>
-<p>
-<b>5.4.12</b> (10 July 2015)</p>
-<ul>
-<li>
-&nbsp;Windows: removed a bug introduced in 5.4.10 whereby SpeechSynthesizer would only pronounce the first letter of the text.
-</ul>
-<p>
-<b>5.4.11</b> (8 July 2015)</p>
-<ul>
-<li>
-&nbsp;Windows: removed a bug introduced in 5.4.10 whereby the file selection window could try to open a file with a garbled name.
-<li>
-&nbsp;Removed a bug whereby the PointProcess window would crash when showing a perfectly silent sound.
-</ul>
-<p>
-<b>5.4.10</b> (27 June 2015)</p>
-<ul>
-<li>
-&nbsp;Removed a bug introduced in version 4.5.09 (January 2007) that caused incorrect upsampling of stereo sounds when the upsampling factor was exactly 2. This bug has caused incorrect playing of stereo sounds with a sampling frequency of 22050 Hz on the Mac since January 2015.
-<li>
-&nbsp;Removed a bug introduced in 2011 that could cause Praat to crash when you scrolled the LongSound window.
-<li>
-&nbsp;TextGrid: <b>Count intervals where...</b> and <b>Count points where...</b>.
-</ul>
-<p>
-<b>5.4.09</b> (1 June 2015)</p>
-<ul>
-<li>
-&nbsp;Linux: the tab key can be used to play sounds on more computers.
-<li>
-&nbsp;Windows: TextGrid files with non-BMP characters are now read correctly.
-<li>
-&nbsp;Windows: files with names that contain non-BMP characters are now saved correctly.
-<li>
-&nbsp;Updated manual.
-</ul>
-<p>
-<b>5.4.08</b> (24 March 2015)</p>
-<ul>
-<li>
-&nbsp;Sound window: removed a bug that caused a crash in "Editor info".
-</ul>
-<p>
-<b>5.4.07</b> (22 March 2015)</p>
-<ul>
-<li>
-&nbsp;TextGrid window: <b>Add interval</b> now does the same thing as <b>Add boundary</b> if a time stretch is selected.
-<li>
-&nbsp;Linux: better redrawing in TextGrid window.
-</ul>
-<p>
-<b>5.4.06</b> (21 February 2015)</p>
-<ul>
-<li>
-&nbsp;Windows: repaired a bug that could cause Praat to crash if a metafile resolution was 200 dpi.
-</ul>
-<p>
-<b>5.4.05</b> (13 February 2015)</p>
-<ul>
-<li>
-&nbsp;Better support for big integer numbers on Windows XP.
-<li>
-&nbsp;Sound window: guarded against empty view ranges for the intensity curve.
-</ul>
-<p>
-<b>5.4.04</b> (28 December 2014)</p>
-<ul>
-<li>
-&nbsp;Windows audio playback: if the sound has more channels than the audio hardware, distribute them evenly.
-</ul>
-<p>
-<b>5.4.03</b> (18 December 2014)</p>
-<ul>
-<li>
-&nbsp;TextGrid reading: Praat now corrects some incomplete TextGrid files created by others.
-<li>
-&nbsp;Better support for text files larger than 2 GB.
-</ul>
-<p>
-<b>5.4.02</b> (26 November 2014)</p>
-<ul>
-<li>
-&nbsp;Mac: repaired a bug by which quote characters typed into the script window could become curly instead of straight.
-</ul>
-<p>
-<b>5.4.01</b> (9 November 2014)</p>
-<ul>
-<li>
-&nbsp;<a href="MelSpectrogram.html">MelSpectrogram</a>, <a href="BarkSpectrogram.html">BarkSpectrogram</a>.
-<li>
-&nbsp;Linux: removed a bug that could cause too many flashes when scrolling a Sound window.
-<li>
-&nbsp;Mac: repaired a bug that could cause a crash in the VowelEditor window.
-</ul>
 <h3>
 What used to be new?</h3>
 <ul>
 <li>
+&nbsp;<a href="What_was_new_in_6_0_.html">What was new in 6.0?</a>
+<li>
 &nbsp;<a href="What_was_new_in_5_4_.html">What was new in 5.4?</a>
 <li>
 &nbsp;<a href="What_was_new_in_5_3_.html">What was new in 5.3?</a>
@@ -710,7 +528,7 @@ What used to be new?</h3>
 </ul>
 <hr>
 <address>
-	<p>&copy; ppgb, June 20, 2019</p>
+	<p>&copy; ppgb, July 13, 2019</p>
 </address>
 </body>
 </html>

debian/changelog

+praat (6.1-1) unstable; urgency=medium
+
+  * New upstream version 6.1
+  * d/watch: Cope with less than three parts in version number
+
+ -- Rafael Laboissiere <rafael@debian.org>  Tue, 16 Jul 2019 09:30:18 -0300
+
 praat (6.0.56-2) unstable; urgency=medium
 
   * Copy the lacking WAV file in d/t/run-tests instead of d/rules.

debian/watch

 version=4
-opts="filenamemangle=s/(?:.*?)?v?(\d+)\.(\d+)\.0?(\d+)\.tar\.gz/praat-$1.$2.$3.tar.gz/"	\
+opts="filenamemangle=s/(?:.*?)?v?([\.\d]+)\.tar\.gz/praat-$1.tar.gz/"	\
     https://github.com/praat/praat/tags (?:.*?/)?v?(\d[\d.]*)\.tar\.gz
 

dwsys/NMF.cpp

@@ -16,6 +16,7 @@
  * along with this work. If not, see <http://www.gnu.org/licenses/>.
  */
 
+#include "Graphics.h"
 #include "NMF.h"
 #include "NUMmachar.h"
 #include "NUM2.h"
@@ -54,6 +55,64 @@ void structNMF :: v_info () {
 Thing_implement (NMF, Daata, 0);
 
 
+double MATgetDivergence_ItakuraSaito (constMATVU const& ref, constMATVU const& x) {
+	Melder_assert (ref.nrow == x.nrow);
+	Melder_assert (ref.ncol == x.ncol);
+	double divergence = 0.0;
+	for (integer irow = 1; irow <= ref.nrow; irow ++)
+		for (integer icol = 1; icol <= ref.ncol; icol ++) {
+			double refval = ref [irow] [icol];
+			if (refval == 0.0)
+				return undefined;
+			divergence += x [irow] [icol] / refval - log (x [irow] [icol] / refval) - 1.0;
+		}
+	return divergence;
+}
+
+void NMF_paintFeatures (NMF me, Graphics g, integer fromFeature, integer toFeature, integer fromRow, integer toRow, double minimum, double maximum, int amplitudeScale, int scaling, bool garnish) {
+	fixUnspecifiedColumnRange (& fromFeature, & toFeature, my features.get());
+	fixUnspecifiedRowRange (& fromRow, & toRow, my features.get());
+	
+	autoMAT part = newMATcopy (my features.part (fromRow, toRow, fromFeature, toFeature));
+	
+	if (minimum == 0.0 && maximum == 0.0) {
+		minimum = NUMmin (part.get());
+		maximum = NUMmax (part.get());
+	}
+	
+	Graphics_setInner (g);
+	Graphics_setWindow (g, fromFeature, toFeature, fromRow, toRow);
+	
+
+	Graphics_cellArray (g, my features.part (fromRow, toRow, fromFeature, toFeature), 
+	fromFeature, toFeature, fromRow, toRow, minimum, maximum);
+	Graphics_unsetInner (g);
+	if (garnish)
+		Graphics_drawInnerBox (g);
+}
+
+void NMF_paintWeights (NMF me, Graphics g, integer fromWeight, integer toWeight, integer fromRow, integer toRow, double minimum, double maximum, int amplitudeScale, int scaling, bool garnish) {
+	fixUnspecifiedColumnRange (& fromWeight, & toWeight, my weights.get());
+	fixUnspecifiedRowRange (& fromRow, & toRow, my weights.get());
+	
+	autoMAT part = newMATcopy (my weights.part (fromRow, toRow, fromWeight, toWeight));
+	
+	if (minimum == 0.0 && maximum == 0.0) {
+		minimum = NUMmin (part.get());
+		maximum = NUMmax (part.get());
+	}
+	
+	Graphics_setInner (g);
+	Graphics_setWindow (g, fromWeight, toWeight, fromRow, toRow);
+	
+
+	Graphics_cellArray (g, my weights.part (fromRow, toRow, fromWeight, toWeight), 
+	fromWeight, toWeight, fromRow, toRow, minimum, maximum);
+	Graphics_unsetInner (g);
+	if (garnish)
+		Graphics_drawInnerBox (g);
+}
+
 autoNMF NMF_create (integer numberOfRows, integer numberOfColumns, integer numberOfFeatures) {
 	try {
 		autoNMF me = Thing_new (NMF);
@@ -117,8 +176,14 @@ double NMF_getEuclideanDistance (NMF me, constMATVU const& data) {
 		U"Dimensions should match.");
 	autoMAT synthesis = NMF_synthesize (me);
 	synthesis.get()  -=  data;
-	double dist = NUMnorm (synthesis.get(), 2.0);
-	return dist;
+	return NUMnorm (synthesis.get(), 2.0);
+}
+
+double NMF_getItakuraSaitoDivergence (NMF me, constMATVU const& data) {
+	Melder_require (data.nrow == my numberOfRows && data.ncol == my numberOfColumns,
+		U"Dimensions should match.");
+	autoMAT synthesis = NMF_synthesize (me);
+	return MATgetDivergence_ItakuraSaito (data, synthesis.get());
 }
 
 static double getMaximumChange (constMATVU const& m, MATVU const& m0, const double sqrteps) {
@@ -272,7 +337,6 @@ void NMF_improveFactorization_als (NMF me, constMATVU const& data, integer maxim
 		const double eps = NUMfpp -> eps;
 		const double sqrteps = sqrt (eps);
 		double dnorm0 = 0.0;
-		double maximum = NUMmax (data);
 		integer iter = 1;
 		bool convergence = false;
 		
@@ -328,6 +392,85 @@ void NMF_improveFactorization_als (NMF me, constMATVU const& data, integer maxim
 	}
 }
 
+static void VECinvertAndScale (VECVU const& target, constVECVU const& source, double scaleFactor) {
+	Melder_assert (target.size == source.size);
+	for (integer i = 1; i <= target.size; i ++)
+		target [i] = scaleFactor / source [i];
+}
+
+void NMF_improveFactorization_is (NMF me, constMATVU const& data, integer maximumNumberOfIterations, double changeTolerance, double approximationTolerance, bool info) {
+	try {
+		Melder_require (my numberOfColumns == data.ncol, U"The number of columns should be equal.");
+		Melder_require (my numberOfRows == data.nrow, U"The number of rows should be equal.");
+		Melder_require (NUMhasZeroElement (data) == false,
+			U"The data matrix should not have cells that are zero.");
+		autoMAT vk = newMATraw (data.nrow, data.ncol);
+		autoMAT fw = newMATraw (data.nrow, data.ncol);
+		autoMAT fcol_x_wrow = newMATraw (data.nrow, data.ncol);
+		autoVEC fcolumn_inv = newVECraw (data.nrow); // feature column
+		autoVEC wrow_inv = newVECraw (data.ncol); // weight row
+		MATmul (fw.get(), my features.get(), my weights.get());
+		double divergence = MATgetDivergence_ItakuraSaito (data, fw.get());
+		double divergence0 = divergence;
+		if (info)
+			MelderInfo_writeLine (U"Iteration: 0", U" divergence: ", divergence, U" delta: ", divergence);
+		integer iter = 1;
+		bool convergence = false;
+		while (iter <= maximumNumberOfIterations && not convergence) {
+			/*
+				algorithm 2, page 806
+				until convergence {
+					for k to numberOfFeateres {
+						G(k) = fcol(k) x wrow(k) / F.H                           (1)
+						V(k) = G(k)^(.2).V+(1-G(k)).(fcol(k) x wrow(k))          (2)
+						wrow(k) <-- (1/fcol(k))' . V(k) / numberOfRows           (3)
+						fcol(k) <-- V(k).(1/wrow(k))' / numberOfColumns          (4)
+						Normalize fcol(k) and wrow(k)                            (5)
+						F.H - old(fcol(k) x wrow(k)) + new(fcol(k) x wrow(k))    (6)
+					}
+				}
+				There is no need to calculate G(k) explicitly as in (1).
+				We can calculate the elements of G(k) while we are doing (2).
+			*/
+			for (integer kf = 1; kf <= my numberOfFeatures; kf ++) {
+				// (1) and (2)
+				MATouter (fcol_x_wrow.get(), my features.column (kf), my weights.row (kf));
+				for (integer irow = 1; irow <= data.nrow; irow ++)
+					for (integer icol = 1; icol <= data.ncol; icol ++) {
+						double gk = fcol_x_wrow [irow] [icol] / fw [irow] [icol];
+						vk [irow] [icol] = gk * gk * data [irow] [icol] + (1.0 - gk) * fcol_x_wrow [irow] [icol];
+					}
+				// (3)
+				VECinvertAndScale (fcolumn_inv.get(), my features.column (kf), 1.0 / my numberOfRows);	
+				VECmul (my weights.row (kf), fcolumn_inv.get(), vk.get());
+				// (4)
+				VECinvertAndScale (wrow_inv.get(), my weights.row (kf), 1.0 / my numberOfColumns);
+				VECmul (my features.column (kf), vk.get(), wrow_inv.get());
+				// (5)
+				double fcolumn_norm = NUMnorm (my features.column (kf), 2.0);
+				my features.column (kf)  /=  fcolumn_norm;
+				my weights.row (kf)  *=  fcolumn_norm;
+				// (6)
+				fw.get()  -=  fcol_x_wrow.get();
+				MATouter (fcol_x_wrow.get(), my features.column (kf), my weights.row (kf));
+				fw.get()  +=  fcol_x_wrow.get();
+			}
+			double divergence_update = MATgetDivergence_ItakuraSaito (data, fw.get());
+			double delta = divergence - divergence_update;
+			convergence = ( iter > 1 && (fabs (delta) < changeTolerance || divergence_update < divergence0 * approximationTolerance) );
+			if (info)
+				MelderInfo_writeLine (U"Iteration: ", iter, U" divergence: ", divergence_update, U" delta: ", delta);
+			++ iter;
+			divergence = divergence_update;
+		}
+		if (info)
+			MelderInfo_drain();
+		
+	} catch (MelderError) {
+		Melder_throw (me, U" IS factorization cannot be improved.");
+	}
+}
+
 autoMAT NMF_synthesize (NMF me) {
 	try {
 		autoMAT result = newMATmul (my features.get(), my weights.get());

dwsys/NMF.h

@@ -30,6 +30,10 @@
 	where H and W are non-negative nrow x k and k x ncol matrices, respectively. 
 */
 
+void NMF_paintFeatures (NMF me, Graphics g, integer fromFeature, integer toFeature, integer fromRow, integer toRow, double minimum, double maximum, int amplitudeScale, int scaling, bool garnish);
+
+void NMF_paintWeights (NMF me, Graphics g, integer fromWeight, integer toWeight, integer fromRow, integer toRow, double minimum, double maximum, int amplitudeScale, int scaling, bool garnish);
+
 autoNMF NMF_create (integer numberOfRows, integer numberOfColumns, integer numberOfFeatures);
 
 autoNMF NMF_createFromGeneralMatrix (constMATVU const& data, integer numberOfFeatures);
@@ -60,8 +64,19 @@ void NMF_improveFactorization_mu (NMF me, constMATVU const& data, integer maximu
 */
 void NMF_improveFactorization_als (NMF me, constMATVU const & data, integer maximumNumberOfIterations, double changeTolerance, double approximationTolerance, bool info);
 
+
+/*
+	Factorize D as F*W, where D, F and W >= 0.
+	
+	C. Févotte, N. Berin & J.-L. Durrieu (2009), Nonnegative matrix facorization with the Itakura-Saito divergene: with 
+	applications to music analysis, Neural Computation 21, 793--830.
+*/
+void NMF_improveFactorization_is (NMF me, constMATVU const& data, integer maximumNumberOfIterations, double changeTolerance, double approximationTolerance, bool info);
+
 autoMAT NMF_synthesize (NMF me); // result = features * weights
 
 double NMF_getEuclideanDistance (NMF me, constMATVU const& data); // sqrt (||data - features*weights||²)
 
+double NMF_getItakuraSaitoDivergence (NMF me, constMATVU const& data);
+
 #endif /* _NMF_h_ */

dwsys/NMF_enums.h

@@ -19,7 +19,6 @@
 enums_begin (kNMF_Initialization, 1)
 	enums_add (kNMF_Initialization, 1, RandomUniform, U"RandomUniform")
 	enums_add (kNMF_Initialization, 2, SVDAbsNegatives, U"SVDAbsNegatives")	
-	
 enums_end (kNMF_Initialization, 2, SVDAbsNegatives)
 
 /* End of file NMF_enums.h */

dwsys/NUM2.h

@@ -172,6 +172,50 @@ inline void VECclip_inplace_inline (VEC x, double min, double max) {
 			x [i] = max;
 }
 
+inline void VECabs (VECVU const& result, constVECVU const& v) {
+	Melder_assert (result.size == v.size);
+	for (integer i = 1; i <= result.size; i ++)
+		result [i] = fabs (v [i]);
+}
+
+inline autoVEC newVECabs (constVECVU const& v) {
+	autoVEC result = newVECraw (v.size);
+	VECabs (result.get(), v);
+	return result;
+}
+
+inline void VECabs_inplace (VECVU const& v) {
+	for (integer i = 1; i <= v.size; i ++)
+		v [i] = fabs (v [i]);
+}
+
+inline void INTVEClinear (INTVEC const& v, integer start, integer step) {
+	for (integer i = 1; i <= v.size; i ++)
+		v [i] = start + (i - 1) * step;
+}
+
+inline autoINTVEC newINTVEClinear (integer size, integer start, integer step) {
+	autoINTVEC result = newINTVECraw (size);
+	INTVEClinear (result, start, step);
+	return result;
+}
+
+inline bool NUMhasZeroElement (constMATVU const m) {
+	for (integer irow = 1; irow <= m.nrow; irow ++)
+		for (integer icol = 1; icol <= m.ncol; icol++)
+			if (m [irow][icol] == 0.0)
+				return true;
+	return false;
+}
+
+inline integer NUMcountNumberOfNonZeroElements (constVECVU const& v) {
+	integer count = 0;
+	for (integer i = 1; i <= v.size; i ++)
+		if (v [i] != 0.0)
+			++ count;
+	return count;
+}
+
 inline double NUMmul (constVECVU const& x, constMATVU const& m, constVECVU const& y) { // x'. M . y
 	Melder_assert (x.size == m.nrow);
 	Melder_assert (y.size == m.ncol);
@@ -256,7 +300,7 @@ void NUMstatistics_huber (constVEC x, double *inout_location, bool wantlocation,
 	k_stdev Winsorizes at `k_stdev' standard deviations.
 */
 
-autoVEC VECmonotoneRegression (constVEC x);
+autoVEC newVECmonotoneRegression (constVEC x);
 /*
 	Find numbers xs[1..n] that have a monotone relationship with
 	the numbers in x[1..n].
@@ -400,13 +444,13 @@ void MATlowerCholeskyInverse_inplace (MAT a, double *out_lnd);
 	and ln(determinant). L^-1 in lower, leave upper part intact.
 */
 
-autoMAT MATinverse_fromLowerCholeskyInverse (constMAT m);
+autoMAT newMATinverse_fromLowerCholeskyInverse (constMAT m);
 /*
 	Return the complete matrix inverse when only the inverse of the lower Cholesky part is given.
 	Input m is a square matrix, in the lower part is the inverse of the lower Cholesky part as calculated by NUMlowerCholeskyInverse.
 */
 
-double MATdeterminant_fromSymmetricMatrix (constMAT m);
+double NUMdeterminant_fromSymmetricMatrix (constMAT m);
 /*
 	ln(determinant) of a symmetric p.s.d. matrix
 */
@@ -462,21 +506,33 @@ integer NUMsolveQuadraticEquation (double a, double b, double c, double *x1, dou
 	If no roots found then x1 and x2 will not be changed.
 */
 
-autoVEC NUMsolveEquation (constMATVU const& a, constVECVU const& b, double tol);
+autoVEC newVECsolve (constMATVU const& a, constVECVU const& b, double tol);
 /*
 	Solve the equation: A.x = b for x;
 	a[1..nr][1..nc], b[1..nr] and the unknown x[1..nc]
 	Algorithm: s.v.d.
 */
 
-autoMAT NUMsolveEquations (constMATVU const& a, constMATVU const& b, double tol);
+autoMAT newMATsolve (constMATVU const& a, constMATVU const& b, double tol);
 /*
 	Solve the equations: A.X = B;
 	a[1..nr][1..nc], b[1..nr][1..nc2] and the unknown x[1..nc][1..nc2]
 	Algorithm: s.v.d.
 */
 
-autoVEC NUMsolveNonNegativeLeastSquaresRegression (constMAT m, constVEC d, double tol, integer itermax);
+
+/*
+	Solve y = D.x + e for x, where x is sparse and e is observation noise.
+	Minimize the 2-norm (y - D.x), where maximally K elements of x may be non-zero, by an iterative hard thresholding algorithm.
+	D is a MxN real matrix with (many) more columns than rows, i.e. N > M. We need to find a vector x
+	with maximally K non-zero elements (sparse).
+	The algorithm is described in T. Blumensath & M.E. Davies, "Normalised iterative hard thresholding;
+	guaranteed stability and performance", IEEE Journal of Selected Topics in Signal Processing #4, 298-309.
+	x in/out: the start value (you typically would start the iteration with all zeros).
+*/
+void newVECsolveSparse_IHT (VECVU const& x, constMATVU const& p, constVECVU const& y, integer numberOfNonZeros, integer maximumNumberOfIterations, double tolerance, bool info);
+
+autoVEC newVECsolveNonNegativeLeastSquaresRegression (constMAT m, constVEC d, double tol, integer itermax);
 /*
 	Solve the equation: M.b = d for b under the constraint: all b[i] >= 0;
 	m[1..nr][1..nc], d[1..nr] and b[1..nc].
@@ -497,7 +553,7 @@ void NUMsolveConstrainedLSQuadraticRegression (constMAT o, constVEC y, double *o
 	Psychometrika 48, 631-638.
 */
 
-autoVEC NUMsolveWeaklyConstrainedLinearRegression (constMAT f, constVEC phi, double alpha, double delta);
+autoVEC newVECsolveWeaklyConstrainedLinearRegression (constMAT f, constVEC phi, double alpha, double delta);
 /*
 	Solve g(t) = ||Ft - phi||^2 + alpha (t't - delta)^2 for t[1..m],
 	where F[1..n][1..m] is a matrix, phi[1..n] a given vector, and alpha
@@ -850,7 +906,7 @@ double NUMcubicSplineInterpolation (constVEC xa, constVEC ya, constVEC y2a, doub
 	a value of x, this routine returns an interpolated value y.
 */
 
-autoVEC NUMbiharmonic2DSplineInterpolation_getWeights (constVECVU const& x, constVECVU const& y, constVECVU const& w);
+autoVEC newVECbiharmonic2DSplineInterpolation_getWeights (constVECVU const& x, constVECVU const& y, constVECVU const& w);
 /*
 	Input: x[1..numberOfPoints], y[1..numberOfPoints], (xp,yp)
 	Output: interpolated result
@@ -860,7 +916,7 @@ double NUMbiharmonic2DSplineInterpolation (constVECVU const& x, constVECVU const
 /* Biharmonic spline interpolation based on Green's function.
 	. Given z[i] values at points (x[i],y[i]) for i=1..n, 
 	Get value at new point (px,py).
-	1. Calculate weights w once: NUMbiharmonic2DSplineInterpolation_getWeights
+	1. Calculate weights w once: newVECbiharmonic2DSplineInterpolation_getWeights
 	2. Interpolate at (xp,yp): NUMbiharmonic2DSplineInterpolation
 	Input: x[1..numberOfPoints], y[1..numberOfPoints], z[1..numberOfPoints], weights[1..numberOfPoints]
 	Output: weights[1..numberOfPoints]
@@ -1161,6 +1217,21 @@ void NUMlineFit_LS (constVEC x, constVEC y, double *out_m, double *out_intercept
 integer NUMrandomBinomial (double p, integer n);
 double NUMrandomBinomial_real (double p, integer n);
 
+/*
+	Generates random numbers according to a Gamma distribution with shape parameter "alpha"
+	and rate parameter "beta".
+	
+	The Gamma distribution of order (shape) parameter alpha and rate (beta) is defined as:
+
+		f(x; alpha, beta) = (1 / Gamma (alpha)) beta^alpha x^(alpha-1) e^(-beta.x),
+		for x > 0, alpha > 0 && beta > 0.
+
+	The method is described in
+		G. Marsaglia & W. Tsang (2000): A simple method for generating gamma variables. ACM Transactions on Mathematical Software, 26(3):363-372.
+	Preconditions: alpha > 0 && beta > 0.
+*/
+double NUMrandomGamma (const double alpha, const double beta);
+
 // IEEE: Programs for digital signal processing section 4.3 LPTRN (modfied)
 // lpc[1..n] to rc[1..n]
 void VECrc_from_lpc (VEC rc, constVEC lpc);

dwsys/NUMstring.cpp

@@ -31,7 +31,7 @@ autoVEC VEC_createFromString (conststring32 s) {
 		Melder_throw (U"Empty string.");
 	autoVEC numbers = newVECraw (tokens.size);
 	for (integer inum = 1; inum <= tokens.size; inum ++)
-		Interpreter_numericExpression (0, tokens [inum].get(), & numbers [inum]);
+		Interpreter_numericExpression (nullptr, tokens [inum].get(), & numbers [inum]);
 	return numbers;
 }
 

dwsys/Permutation.cpp

@@ -28,6 +28,7 @@
 
 #include <time.h>
 #include "Permutation.h"
+#include "NUM2.h"
 
 #include "oo_DESTROY.h"
 #include "Permutation_def.h"
@@ -82,8 +83,7 @@ void structPermutation :: v_readText (MelderReadText text, int /*formatVersion*/
 
 void Permutation_init (Permutation me, integer numberOfElements) {
 	my numberOfElements = numberOfElements;
-	my p = newINTVECraw (numberOfElements);
-	Permutation_sort (me);  // performs: p [i] = i
+	my p = newINTVEClinear (numberOfElements, 1, 1);
 }
 
 void Permutation_tableJump_inline (Permutation me, integer jumpSize, integer first) {
@@ -116,8 +116,7 @@ autoPermutation Permutation_create (integer numberOfElements) {
 }
 
 void Permutation_sort (Permutation me) {
-	for (integer i = 1; i <= my numberOfElements; i ++)
-		my p [i] = i;
+	INTVEClinear (my p.get(), 1, 1);
 }
 
 void Permutation_swapPositions (Permutation me, integer i1, integer i2) {
@@ -419,7 +418,8 @@ autoPermutation Permutations_multiply2 (Permutation me, Permutation thee) {
 
 autoPermutation Permutations_multiply (OrderedOf<structPermutation>* me) {
 	try {
-	Melder_require (my size > 1, U"There should be at least two Permutations to multiply.");
+	Melder_require (my size > 1,
+		U"There should be at least two Permutations to multiply.");
 		autoPermutation thee = Permutations_multiply2 (my at [1], my at [2]);
 		for (integer i = 3; i <= my size; i ++) {
 			autoPermutation him = Permutations_multiply2 (thee.get(), my at [i]);