EEG/EEG.cpp

 /* EEG.cpp
  *
- * Copyright (C) 2011-2018 Paul Boersma
+ * Copyright (C) 2011-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
@@ -97,28 +97,38 @@ autoEEG EEG_readFromBdfFile (MelderFile file) {
 	try {
 		autofile f = Melder_fopen (file, "rb");
 		char buffer [81];
-		fread (buffer, 1, 8, f); buffer [8] = '\0';
+		fread (buffer, 1, 8, f);
+		buffer [8] = '\0';
 		bool is24bit = buffer [0] == (char) 255;
-		fread (buffer, 1, 80, f); buffer [80] = '\0';
+		fread (buffer, 1, 80, f);
+		buffer [80] = '\0';
 		trace (U"Local subject identification: \"", Melder_peek8to32 (buffer), U"\"");
-		fread (buffer, 1, 80, f); buffer [80] = '\0';
+		fread (buffer, 1, 80, f);
+		buffer [80] = '\0';
 		trace (U"Local recording identification: \"", Melder_peek8to32 (buffer), U"\"");
-		fread (buffer, 1, 8, f); buffer [8] = '\0';
+		fread (buffer, 1, 8, f);
+		buffer [8] = '\0';
 		trace (U"Start date of recording: \"", Melder_peek8to32 (buffer), U"\"");
-		fread (buffer, 1, 8, f); buffer [8] = '\0';
+		fread (buffer, 1, 8, f);
+		buffer [8] = '\0';
 		trace (U"Start time of recording: \"", Melder_peek8to32 (buffer), U"\"");
-		fread (buffer, 1, 8, f); buffer [8] = '\0';
+		fread (buffer, 1, 8, f);
+		buffer [8] = '\0';
 		integer numberOfBytesInHeaderRecord = atol (buffer);
 		trace (U"Number of bytes in header record: ", numberOfBytesInHeaderRecord);
-		fread (buffer, 1, 44, f); buffer [44] = '\0';
+		fread (buffer, 1, 44, f);
+		buffer [44] = '\0';
 		trace (U"Version of data format: \"", Melder_peek8to32 (buffer), U"\"");
-		fread (buffer, 1, 8, f); buffer [8] = '\0';
+		fread (buffer, 1, 8, f);
+		buffer [8] = '\0';
 		integer numberOfDataRecords = strtol (buffer, nullptr, 10);
 		trace (U"Number of data records: ", numberOfDataRecords);
-		fread (buffer, 1, 8, f); buffer [8] = '\0';
+		fread (buffer, 1, 8, f);
+		buffer [8] = '\0';
 		double durationOfDataRecord = atof (buffer);
 		trace (U"Duration of a data record: ", durationOfDataRecord);
-		fread (buffer, 1, 4, f); buffer [4] = '\0';
+		fread (buffer, 1, 4, f);
+		buffer [4] = '\0';
 		integer numberOfChannels = atol (buffer);
 		trace (U"Number of channels in data record: ", numberOfChannels);
 		if (numberOfBytesInHeaderRecord != (numberOfChannels + 1) * 256)
@@ -126,7 +136,8 @@ autoEEG EEG_readFromBdfFile (MelderFile file) {
 				U") doesn't match number of channels (", numberOfChannels, U").");
 		autostring32vector channelNames (numberOfChannels);
 		for (integer ichannel = 1; ichannel <= numberOfChannels; ichannel ++) {
-			fread (buffer, 1, 16, f); buffer [16] = '\0';   // labels of the channels
+			fread (buffer, 1, 16, f);
+			buffer [16] = '\0';   // labels of the channels
 			/*
 			 * Strip all final spaces.
 			 */
@@ -143,37 +154,45 @@ autoEEG EEG_readFromBdfFile (MelderFile file) {
 		bool hasLetters = str32equ (channelNames [numberOfChannels].get(), U"EDF Annotations");
 		double samplingFrequency = undefined;
 		for (integer channel = 1; channel <= numberOfChannels; channel ++) {
-			fread (buffer, 1, 80, f); buffer [80] = '\0';   // transducer type
+			fread (buffer, 1, 80, f);
+			buffer [80] = '\0';   // transducer type
 		}
 		for (integer channel = 1; channel <= numberOfChannels; channel ++) {
-			fread (buffer, 1, 8, f); buffer [8] = '\0';   // physical dimension of channels
+			fread (buffer, 1, 8, f);
+			buffer [8] = '\0';   // physical dimension of channels
 		}
 		autoVEC physicalMinimum (numberOfChannels, kTensorInitializationType::RAW);
 		for (integer ichannel = 1; ichannel <= numberOfChannels; ichannel ++) {
-			fread (buffer, 1, 8, f); buffer [8] = '\0';
+			fread (buffer, 1, 8, f);
+			buffer [8] = '\0';
 			physicalMinimum [ichannel] = atof (buffer);
 		}
 		autoVEC physicalMaximum (numberOfChannels, kTensorInitializationType::RAW);
 		for (integer ichannel = 1; ichannel <= numberOfChannels; ichannel ++) {
-			fread (buffer, 1, 8, f); buffer [8] = '\0';
+			fread (buffer, 1, 8, f);
+			buffer [8] = '\0';
 			physicalMaximum [ichannel] = atof (buffer);
 		}
 		autoVEC digitalMinimum (numberOfChannels, kTensorInitializationType::RAW);
 		for (integer ichannel = 1; ichannel <= numberOfChannels; ichannel ++) {
-			fread (buffer, 1, 8, f); buffer [8] = '\0';
+			fread (buffer, 1, 8, f);
+			buffer [8] = '\0';
 			digitalMinimum [ichannel] = atof (buffer);
 		}
 		autoVEC digitalMaximum (numberOfChannels, kTensorInitializationType::RAW);
 		for (integer ichannel = 1; ichannel <= numberOfChannels; ichannel ++) {
-			fread (buffer, 1, 8, f); buffer [8] = '\0';
+			fread (buffer, 1, 8, f);
+			buffer [8] = '\0';
 			digitalMaximum [ichannel] = atof (buffer);
 		}
 		for (integer channel = 1; channel <= numberOfChannels; channel ++) {
-			fread (buffer, 1, 80, f); buffer [80] = '\0';   // prefiltering
+			fread (buffer, 1, 80, f);
+			buffer [80] = '\0';   // prefiltering
 		}
 		integer numberOfSamplesPerDataRecord = 0;
 		for (integer channel = 1; channel <= numberOfChannels; channel ++) {
-			fread (buffer, 1, 8, f); buffer [8] = '\0';   // number of samples in each data record
+			fread (buffer, 1, 8, f);
+			buffer [8] = '\0';   // number of samples in each data record
 			integer numberOfSamplesInThisDataRecord = atol (buffer);
 			if (isundef (samplingFrequency)) {
 				numberOfSamplesPerDataRecord = numberOfSamplesInThisDataRecord;
@@ -185,7 +204,8 @@ autoEEG EEG_readFromBdfFile (MelderFile file) {
 					U") doesn't match sampling frequency of channel 1 (", samplingFrequency, U").");
 		}
 		for (integer channel = 1; channel <= numberOfChannels; channel ++) {
-			fread (buffer, 1, 32, f); buffer [32] = '\0';   // reserved
+			fread (buffer, 1, 32, f);
+			buffer [32] = '\0';   // reserved
 		}
 		double duration = numberOfDataRecords * durationOfDataRecord;
 		autoEEG him = EEG_create (0, duration);
@@ -282,8 +302,10 @@ autoEEG EEG_readFromBdfFile (MelderFile file) {
 				time = undefined;   // defensive
 			}
 		} else {
-			thee = TextGrid_create (0, duration,
-				numberOfStatusBits == 8 ? U"S1 S2 S3 S4 S5 S6 S7 S8" : U"S1 S2 S3 S4 S5 S6 S7 S8 S9 S10 S11 S12 S13 S14 S15 S16", U"");
+			thee = TextGrid_create (0.0, duration,
+				numberOfStatusBits == 8 ? U"S1 S2 S3 S4 S5 S6 S7 S8" : U"S1 S2 S3 S4 S5 S6 S7 S8 S9 S10 S11 S12 S13 S14 S15 S16",
+				U""
+			);
 			for (int bit = 1; bit <= numberOfStatusBits; bit ++) {
 				uint32 bitValue = 1 << (bit - 1);
 				IntervalTier tier = (IntervalTier) thy tiers->at [bit];
@@ -435,11 +457,10 @@ void EEG_filter (EEG me, double lowFrequency, double lowWidth, double highFreque
 			autoSpectrum spec = Sound_to_Spectrum (channel.get(), true);
 			Spectrum_passHannBand (spec.get(), lowFrequency, 0.0, lowWidth);
 			Spectrum_passHannBand (spec.get(), 0.0, highFrequency, highWidth);
-			if (doNotch50Hz) {
+			if (doNotch50Hz)
 				Spectrum_stopHannBand (spec.get(), 48.0, 52.0, 1.0);
-			}
 			autoSound him = Spectrum_to_Sound (spec.get());
-			NUMvector_copyElements (& his z [1] [0], & my sound -> z [ichan] [0], 1, my sound -> nx);
+			my sound -> z.row (ichan) <<= his z.row (1).part (1, my sound -> nx);
 		}
 	} catch (MelderError) {
 		Melder_throw (me, U": not filtered.");
@@ -472,15 +493,13 @@ void EEG_subtractReference (EEG me, conststring32 channelName1, conststring32 ch
 	if (channelNumber1 == 0)
 		Melder_throw (me, U": no channel named \"", channelName1, U"\".");
 	integer channelNumber2 = EEG_getChannelNumber (me, channelName2);
-	if (channelNumber2 == 0 && channelName2 [0] != '\0')
+	if (channelNumber2 == 0 && channelName2 [0] != U'\0')
 		Melder_throw (me, U": no channel named \"", channelName2, U"\".");
 	const integer numberOfElectrodeChannels = my numberOfChannels - EEG_getNumberOfExtraSensors (me);
 	for (integer isamp = 1; isamp <= my sound -> nx; isamp ++) {
-		double referenceValue = channelNumber2 == 0 ? my sound -> z [channelNumber1] [isamp] :
-			0.5 * (my sound -> z [channelNumber1] [isamp] + my sound -> z [channelNumber2] [isamp]);
-		for (integer ichan = 1; ichan <= numberOfElectrodeChannels; ichan ++) {
-			my sound -> z [ichan] [isamp] -= referenceValue;
-		}
+		const double referenceValue = ( channelNumber2 == 0 ? my sound -> z [channelNumber1] [isamp] :
+			0.5 * (my sound -> z [channelNumber1] [isamp] + my sound -> z [channelNumber2] [isamp]) );
+		my sound -> z.column (isamp).part (1, numberOfElectrodeChannels)  -=  referenceValue;
 	}
 }
 
@@ -493,12 +512,8 @@ void EEG_subtractMeanChannel (EEG me, integer fromChannel, integer toChannel) {
 		Melder_throw (U"Channel range cannot run from ", fromChannel, U" to ", toChannel, U". Please reverse.");
 	const integer numberOfElectrodeChannels = my numberOfChannels - EEG_getNumberOfExtraSensors (me);
 	for (integer isamp = 1; isamp <= my sound -> nx; isamp ++) {
-		double referenceValue = 0.0;
-		for (integer ichan = fromChannel; ichan <= toChannel; ichan ++)
-			referenceValue += my sound -> z [ichan] [isamp];
-		referenceValue /= (toChannel - fromChannel + 1);
-		for (integer ichan = 1; ichan <= numberOfElectrodeChannels; ichan ++)
-			my sound -> z [ichan] [isamp] -= referenceValue;
+		const double referenceValue = NUMmean (my sound -> z.column (isamp).part (fromChannel, toChannel));
+		my sound -> z.column (isamp).part (1, numberOfElectrodeChannels)  -=  referenceValue;
 	}
 }
 
@@ -506,10 +521,7 @@ void EEG_setChannelToZero (EEG me, integer channelNumber) {
 	try {
 		if (channelNumber < 1 || channelNumber > my numberOfChannels)
 			Melder_throw (U"No channel ", channelNumber, U".");
-		integer numberOfSamples = my sound -> nx;
-		VEC channel = my sound -> z.row (channelNumber);
-		for (integer isample = 1; isample <= numberOfSamples; isample ++)
-			channel [isample] = 0.0;
+		my sound -> z.row (channelNumber) <<= 0.0;
 	} catch (MelderError) {
 		Melder_throw (me, U": channel ", channelNumber, U" not set to zero.");
 	}
@@ -563,7 +575,7 @@ autoEEG EEG_extractChannel (EEG me, conststring32 channelName) {
 	}
 }
 
-autoEEG EEG_extractChannels (EEG me, constVEC channelNumbers) {
+autoEEG EEG_extractChannels (EEG me, constVECVU const& channelNumbers) {
 	try {
 		integer numberOfChannels = channelNumbers.size;
 		Melder_require (numberOfChannels > 0,
@@ -590,7 +602,7 @@ static void Sound_removeChannel (Sound me, integer channelNumber) {
 		Melder_require (my ny > 1,
 			U"Cannot remove last remaining channel.");
 		for (integer ichan = channelNumber; ichan < my ny; ichan ++)
-			NUMvector_copyElements (& my z [ichan + 1] [0], & my z [ichan] [0], 1, my nx);
+			my z.row (ichan) <<= my z.row (ichan + 1);
 		my ymax -= 1.0;
 		my ny -= 1;
 	} catch (MelderError) {

EEG/EEG.h

@@ -53,7 +53,7 @@ void EEG_setChannelToZero (EEG me, conststring32 channelName);
 void EEG_removeTriggers (EEG me, kMelder_string which, conststring32 criterion);
 autoEEG EEG_extractChannel (EEG me, integer channelNumber);
 autoEEG EEG_extractChannel (EEG me, conststring32 channelName);
-autoEEG EEG_extractChannels (EEG me, constVEC channelNumbers);
+autoEEG EEG_extractChannels (EEG me, constVECVU const& channelNumbers);
 void EEG_removeChannel (EEG me, integer channelNumber);
 void EEG_removeChannel (EEG me, conststring32 channelName);
 static inline autoSound EEG_extractSound (EEG me) { return Data_copy (my sound.get()); }

LPC/Cepstrogram.cpp

@@ -184,22 +184,18 @@ autoPowerCepstrogram PowerCepstrogram_smooth (PowerCepstrogram me, double timeAv
 		integer numberOfFrames = Melder_ifloor (timeAveragingWindow / my dx);
 		if (numberOfFrames > 1) {
 			autoVEC qin = newVECraw (my nx);
-			autoVEC qout = newVECraw (my nx);
 			for (integer iq = 1; iq <= my ny; iq ++) {
 				qin.all() <<= thy z.row (iq);   // ppgb: why this extra copying?
-				VECsmoothByMovingAverage_preallocated (qout.get(), qin.get(), numberOfFrames);
-				thy z.row (iq) <<= qout.all();
+				VECsmoothByMovingAverage_preallocated (thy z.row (iq), qin.get(), numberOfFrames);
 			}
 		}
 		// 2. average across quefrencies
 		integer numberOfQuefrencyBins = Melder_ifloor (quefrencyAveragingWindow / my dy);
 		if (numberOfQuefrencyBins > 1) {
 			autoVEC qin = newVECraw (thy ny);
-			autoVEC qout = newVECraw (thy ny);
 			for (integer iframe = 1; iframe <= my nx; iframe ++) {
 				qin.get() <<= thy z.column (iframe);
-				VECsmoothByMovingAverage_preallocated (qout.get(), qin.get(), numberOfQuefrencyBins);
-				thy z.column (iframe) <<= qout.get();
+				VECsmoothByMovingAverage_preallocated (thy z.column (iframe), qin.get(), numberOfQuefrencyBins);
 			}
 		}
 		return thee;

LPC/Cepstrum.cpp

@@ -110,7 +110,7 @@ static void _Cepstrum_draw (Cepstrum me, Graphics g, double qmin, double qmax, d
 	if (autoscaling)
 		NUMextrema (y.get(), & minimum, & maximum);
 	else
-		VECclip_inplace (y.get(), minimum, maximum);
+		VECclip_inplace_inline (y.get(), minimum, maximum);
 
 	Graphics_setWindow (g, qmin, qmax, minimum, maximum);
 	Graphics_function (g, y.at, 1, numberOfSelected, Matrix_columnToX (me, imin), Matrix_columnToX (me, imax));

LPC/LPC.cpp

@@ -125,7 +125,7 @@ void LPC_drawGain (LPC me, Graphics g, double tmin, double tmax, double gmin, do
 void LPC_drawPoles (LPC me, Graphics g, double time, bool garnish) {
 	autoPolynomial p = LPC_to_Polynomial (me, time);
 	autoRoots r = Polynomial_to_Roots (p.get());
-	Roots_draw (r.get(), g, -1.0, 1.0, -1.0, 1.0, U"+", 12, garnish);
+	Roots_draw (r.get(), g, -1.0, 1.0, -1.0, 1.0, U"+", 12.0, garnish);
 }
 
 autoMatrix LPC_downto_Matrix_lpc (LPC me) {

LPC/Sound_and_LPC.cpp

@@ -210,7 +210,7 @@ end:
 }
 
 static int Sound_into_LPC_Frame_burg (Sound me, LPC_Frame thee) {
-	thy gain = NUMburg_preallocated (thy a.get(), my z.row(1));
+	thy gain = VECburg (thy a.get(), my z.row(1));
 	thy gain *= my nx;
 	for (integer i = 1; i <= thy nCoefficients; i ++) {
 		thy a [i] = -thy a [i];

LPC/manual_LPC.cpp

@@ -529,7 +529,7 @@ NORMAL (U"The settings for ##Time step (s)#, ##Maximum number of formants#, ##Ma
 	"##Window length (s)# and ##Pre emphasis from (Hz)# are as in @@Sound: To Formant (burg)...@. "
 	" The following settings determine aspects of the iterative formant frequency refinement.")
 TAG (U"%%Number of std. dev.%,")
-DEFINITION (U"determines the number of standard deviation from where selective weighing of samples starts. ")
+DEFINITION (U"determines the number of standard deviations from where selective weighing of samples starts. ")
 TAG (U"%%Maximum number of iterations%,")
 DEFINITION (U"determines the maximum number of iterations allowed in the refinement step.")
 TAG (U"%%Tolerance%,")

debian/What_s_new_.html

@@ -7,6 +7,30 @@ What's new?
 <p>
 Latest changes in Praat.</p>
 <p>
+<b>6.0.52</b> (2 May 2019)</p>
+<ul>
+<li>
+&nbsp;Removed a bug introduced in 6.0.51 that could cause turning a Discriminant into a Configuration to crash.
+<li>
+&nbsp;Removed a bug introduced in contour grey drawing in August 2017.
+</ul>
+<p>
+<b>6.0.51</b> (29 April 2019)</p>
+<ul>
+<li>
+&nbsp;Script window: Use Selection for Find.
+<li>
+&nbsp;Removed a bug introduced in 6.0.41 that could cause Praat to crash after removing an element from a Strings or a row from a TableOfReal.
+</ul>
+<p>
+<b>6.0.50</b> (31 March 2019)</p>
+<ul>
+<li>
+&nbsp;Manual updates, speed, more tests.
+<li>
+&nbsp;Scripting: rowSums#, columnSums#; randomGauss## finally with four arguments.
+</ul>
+<p>
 <b>6.0.49</b> (2 March 2019)</p>
 <ul>
 <li>
@@ -656,7 +680,7 @@ What used to be new?</h3>
 </ul>
 <hr>
 <address>
-	<p>&copy; ppgb, March 2, 2019</p>
+	<p>&copy; ppgb, May 2, 2019</p>
 </address>
 </body>
 </html>

debian/changelog

+praat (6.0.52-1) experimental; urgency=medium
+
+  * New upstream version 6.0.52
+
+ -- Rafael Laboissiere <rafael@debian.org>  Mon, 06 May 2019 04:15:47 -0300
+
 praat (6.0.49-2) unstable; urgency=medium
 
   * t/run-tests: Drop test_Discriminant from the list of excluded unit tests

debian/clean

@@ -5,6 +5,7 @@ test/fon/kanweg.nsp
 test/fon/kanweg.Object
 test/kar/kanweg*.txt
 test/kar/kanweg.TextGrid
+test/stat/kanweg.txt
 test/sys/kanweg.txt
 test/sys/kanweg.wav
 test/sys/kanweg2.txt

dwsys/Eigen.cpp

 /* Eigen.cpp
  *
- * Copyright (C) 1993-2018 David Weenink
+ * Copyright (C) 1993-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
@@ -110,8 +110,9 @@ void Eigen_init (Eigen me, integer numberOfEigenvalues, integer dimension) {
 	Eigenvectors: the columns of the matrix V
 	Eigenvalues: D_i^2
 */
-void Eigen_initFromSquareRoot (Eigen me, constMAT a) {
-	Melder_assert (a.nrow >= 1);
+void Eigen_initFromSquareRoot (Eigen me, constMATVU const& a) {
+	Melder_require (a.nrow >= 1,
+		U"The matrix must at least have one row.");
 	integer nsv = std::min (a.nrow, a.ncol);
 	my dimension = a.ncol;
 	autoSVD svd = SVD_createFromGeneralMatrix (a);
@@ -141,7 +142,6 @@ void Eigen_initFromSquareRoot (Eigen me, constMAT a) {
 
 
 void Eigen_initFromSquareRootPair (Eigen me, constMAT a, constMAT b) {
-	//Melder_assert (a.nrow >= a.ncol && b.nrow >= b.ncol);   // pggb: this cannot be an assert, and seems too strict anyway
 	Melder_require (a.ncol == b.ncol,
 		U"The numbers of columns should be equal, not ", a.ncol, U" and ", b.ncol, U".");
 	// Eigen has not been inited yet.
@@ -166,7 +166,8 @@ void Eigen_initFromSquareRootPair (Eigen me, constMAT a, constMAT b) {
 	(void) NUMlapack_dggsvd (& jobu, & jobv, & jobq, & m, & n, & p, & k, & ll,
 		& ac [1][1], & lda, & bc [1][1], & ldb, alpha.begin(), beta.begin(), u, & ldu,
 		v, & ldv, & q [1][1], & ldq, work.begin(), iwork.begin(), & info);
-	Melder_require (info == 0, U"dggsvd fails with code ", info, U".");
+	Melder_require (info == 0,
+		U"dggsvd fails with code ", info, U".");
 
 	// Calculate the eigenvalues (alpha[i]/beta[i])^2 and store in alpha[i].
 
@@ -188,7 +189,8 @@ void Eigen_initFromSquareRootPair (Eigen me, constMAT a, constMAT b) {
 		}
 	}
 
-	Melder_require (ll - numberOfDeselected > 0, U"No eigenvectors can be found. Matrix too singular.");
+	Melder_require (ll - numberOfDeselected > 0,
+		U"No eigenvectors can be found. Matrix too singular.");
 
 	Eigen_init (me, ll - numberOfDeselected, a.ncol);
 
@@ -205,7 +207,7 @@ void Eigen_initFromSquareRootPair (Eigen me, constMAT a, constMAT b) {
 	MATnormalizeRows_inplace (my eigenvectors.get(), 2.0, 1.0);
 }
 
-void Eigen_initFromSymmetricMatrix (Eigen me, constMAT a) {
+void Eigen_initFromSymmetricMatrix (Eigen me, constMATVU const& a) {
 	Melder_assert (a.ncol == a.nrow);
 	if (NUMisEmpty (my eigenvectors))   // ppgb: BUG dubious logic
 		Eigen_init (me, a.ncol, a.ncol);
@@ -293,10 +295,8 @@ void Eigen_sort (Eigen me) {
 }
 
 void Eigen_invertEigenvector (Eigen me, integer ivec) {
-
-	if (ivec < 1 || ivec > my numberOfEigenvalues) {
-		return;
-	}
+	Melder_require (ivec >= 1 and ivec <= my numberOfEigenvalues,
+		U"The eigenvector number should be in the interval from 1 to ", my numberOfEigenvalues, U".");
 	
 	for (integer j = 1; j <= my dimension; j ++) {
 		my eigenvectors [ivec] [j] = - my eigenvectors [ivec] [j];
@@ -366,7 +366,7 @@ void Eigen_drawEigenvector (Eigen me, Graphics g, integer ivec, integer first, i
 	// If ymax < ymin the eigenvector will automatically be drawn inverted.
 
 	if (ymax == ymin) {
-		NUMextrema (vec, first, last, & ymin, & ymax);
+		NUMextrema (vec.part (first, last), & ymin, & ymax);
 		ymax *= w;
 		ymin *= w;
 	}
@@ -423,8 +423,10 @@ static autoVEC Eigens_getAnglesBetweenSubspaces (Eigen me, Eigen thee, integer i
 
 	integer nmin = std::min (my numberOfEigenvalues, thy numberOfEigenvalues);
 
-	Melder_require (my dimension == thy dimension, U"The eigenvectors should have equal dimensions.");
-	Melder_require (ivec_from > 0 && ivec_from <= ivec_to && ivec_to <= nmin, U"Eigenvector range too large.");
+	Melder_require (my dimension == thy dimension,
+		U"The eigenvectors should have equal dimensions.");
+	Melder_require (ivec_from > 0 && ivec_from <= ivec_to && ivec_to <= nmin,
+		U"Eigenvector range too large.");
 
 	/*
 		Algorithm 12.4.3 Golub & van Loan
@@ -435,10 +437,8 @@ static autoVEC Eigens_getAnglesBetweenSubspaces (Eigen me, Eigen thee, integer i
 
 	autoVEC angles_degrees = newVECraw (numberOfVectors);
 
-	autoMAT c = newMATmul (
-			my eigenvectors.horizontalBand (ivec_from, ivec_to),
-			thy eigenvectors. horizontalBand (ivec_from, ivec_to). transpose()
-	);
+	autoMAT c = newMATmul (my eigenvectors.horizontalBand (ivec_from, ivec_to),
+			thy eigenvectors. horizontalBand (ivec_from, ivec_to). transpose());
 	autoSVD svd = SVD_createFromGeneralMatrix (c.get());
 	for (integer i = 1; i <= numberOfVectors; i ++) {
 		angles_degrees [i] = acos (svd -> d [i]) * (180.0 / NUMpi);

dwsys/Eigen.h

@@ -2,7 +2,7 @@
 #define _Eigen_h_
 /* Eigen.h
  *
- * Copyright (C) 1993-2018 David Weenink
+ * Copyright (C) 1993-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
@@ -28,9 +28,9 @@ autoEigen Eigen_create (integer numberOfEigenvalues, integer dimension);
 
 void Eigen_init (Eigen me, integer numberOfEigenvalues, integer dimension);
 
-void Eigen_initFromSymmetricMatrix (Eigen me, constMAT a);
+void Eigen_initFromSymmetricMatrix (Eigen me, constMATVU const& a);
 
-void Eigen_initFromSquareRoot (Eigen me, constMAT a);
+void Eigen_initFromSquareRoot (Eigen me, constMATVU const& a);
 /*
 	Calculate eigenstructure for symmetric matrix A'A (e.g. covariance matrix),
 	when only A is given.

dwsys/FileInMemory.cpp

@@ -85,7 +85,7 @@ autoFileInMemory FileInMemory_createWithData (integer numberOfBytes, const char
 		} else {
 			my _dontOwnData = false;
 			my d_data = NUMvector <unsigned char> ((integer) 0, numberOfBytes);
-			NUMvector_copyElements <unsigned char> (reinterpret_cast<unsigned char *> (const_cast<char *> (data)), my d_data, 0, numberOfBytes);
+			memcpy (my d_data, data, (size_t) numberOfBytes + 1);
 		} 
 		return me;
 	} catch (MelderError) {

dwsys/Graphics_extensions.cpp

@@ -163,9 +163,9 @@ void Graphics_boxAndWhiskerPlot (Graphics g, constVEC data, double x, double r,
 }
 
 void Graphics_quantileQuantilePlot (Graphics g, integer numberOfQuantiles, constVEC x, constVEC y, 
-	double xmin, double xmax, double ymin, double ymax, int labelSize, conststring32 plotLabel)
+	double xmin, double xmax, double ymin, double ymax, double labelSize, conststring32 plotLabel)
 {
-	int fontSize = Graphics_inqFontSize (g);
+	const double fontSize = Graphics_inqFontSize (g);
 
 	Graphics_setTextAlignment (g, Graphics_CENTRE, Graphics_HALF);
 	Graphics_setFontSize (g, labelSize);
@@ -200,11 +200,11 @@ void Graphics_quantileQuantilePlot (Graphics g, integer numberOfQuantiles, const
 	Graphics_setFontSize (g, fontSize);
 }
 
-void Graphics_lagPlot (Graphics g, constVEC data, double xmin, double xmax, integer lag, int labelSize, conststring32 plotLabel)
+void Graphics_lagPlot (Graphics g, constVEC data, double xmin, double xmax, integer lag, double labelSize, conststring32 plotLabel)
 {
 	if (lag < 0 || lag >= data.size)
 		return;
-	int fontSize = Graphics_inqFontSize (g);
+	const double fontSize = Graphics_inqFontSize (g);
 	Graphics_setFontSize (g, labelSize);
 	Graphics_setTextAlignment (g, Graphics_CENTRE, Graphics_HALF);
 	// plot x[i] vertically and x[i-lag] horizontally

dwsys/Graphics_extensions.h

@@ -28,8 +28,8 @@
 void Graphics_boxAndWhiskerPlot (Graphics g, constVEC data, double x, double r, double w, double ymin, double ymax);
 
 void Graphics_quantileQuantilePlot (Graphics g, integer numberOfQuantiles, constVEC xdata, constVEC ydata,
-	double xmin, double xmax, double ymin, double ymax, int labelSize, conststring32 plotLabel);
+	double xmin, double xmax, double ymin, double ymax, double labelSize, conststring32 plotLabel);
 
-void Graphics_lagPlot (Graphics g, constVEC x, double xmin, double xmax, integer lag, int labelSize, conststring32 plotLabel);
+void Graphics_lagPlot (Graphics g, constVEC x, double xmin, double xmax, integer lag, double labelSize, conststring32 plotLabel);
 
 #endif /* _Graphics_extensions_h_ */

dwsys/MAT_numerics.cpp

 /* MAT_numerics.cpp
  *
- * Copyright (C) 2018 David Weenink
+ * Copyright (C) 2018-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
@@ -19,22 +19,23 @@
 #include "NUMclapack.h"
 #include "MAT_numerics.h"
 #include "SVD.h"
+#include "PAIRWISE_SUM.h"
 
-void MAT_getEigenSystemFromSymmetricMatrix_preallocated (MAT eigenvectors, VEC eigenvalues, constMAT m, bool sortAscending) {
+void MAT_getEigenSystemFromSymmetricMatrix_preallocated (MAT eigenvectors, VEC eigenvalues, constMATVU const& m, bool sortAscending) {
 	Melder_assert (m.nrow == m.ncol);
 	Melder_assert (eigenvalues.size == m.ncol);
 	Melder_assert (eigenvectors.nrow == eigenvectors.ncol);
 	Melder_assert (m.nrow == eigenvectors.nrow);
 
 	char jobz = 'V', uplo = 'U';
-	integer workSize = -1, info;
+	integer workSize = -1, info, ncol = m.ncol;
 	double wt [1];
 	
 	eigenvectors <<= m;
 	// 0. No need to transpose a because it is a symmetric matrix
 	// 1. Query for the size of the work array
 	
-	(void) NUMlapack_dsyev (& jobz, & uplo, & m.ncol, & eigenvectors [1] [1], & m.ncol, eigenvalues.begin(), wt, & workSize, & info);
+	(void) NUMlapack_dsyev (& jobz, & uplo, & ncol, & eigenvectors [1] [1], & ncol, eigenvalues.begin(), wt, & workSize, & info);
 	Melder_require (info == 0, U"dsyev initialization code = ", info, U").");
 	
 	workSize = Melder_iceiling (wt [0]);
@@ -42,7 +43,7 @@ void MAT_getEigenSystemFromSymmetricMatrix_preallocated (MAT eigenvectors, VEC e
 
 	// 2. Calculate the eigenvalues and eigenvectors (row-wise)
 	
-	(void) NUMlapack_dsyev (& jobz, & uplo, & m.ncol, & eigenvectors [1] [1], & m.ncol, eigenvalues.begin(), work.begin(), & workSize, & info);
+	(void) NUMlapack_dsyev (& jobz, & uplo, & ncol, & eigenvectors [1] [1], & ncol, eigenvalues.begin(), work.begin(), & workSize, & info);
 	Melder_require (info == 0, U"dsyev code = ", info, U").");
 
 	// 3. Eigenvalues are returned in ascending order
@@ -105,8 +106,7 @@ void MAT_getEigenSystemFromGeneralMatrix (constMAT a, autoMAT *out_lefteigenvect
 	integer right_nvecs = out_righteigenvectors ? n : 1;
 	double wt;
 	
-	autoMAT data = newMATraw (n, n);
-	MATtranspose_preallocated (data.get(), a); // lapack is fortran storage
+	autoMAT data = newMATtranspose (a);   // lapack is fortran storage
 	autoVEC eigenvalues_re = newVECraw (n);
 	autoVEC eigenvalues_im = newVECraw (n);
 	autoMAT righteigenvectors = newMATraw (right_nvecs, right_nvecs); // 1x1 if not needed
@@ -129,39 +129,41 @@ void MAT_getEigenSystemFromGeneralMatrix (constMAT a, autoMAT *out_lefteigenvect
 	if (out_eigenvalues_im) *out_eigenvalues_im = eigenvalues_im.move();
 }
 
-void MAT_asPrincipalComponents_preallocated (MAT pc, constMAT m) {
-	Melder_assert (pc.nrow == m.nrow && pc.ncol <= m.ncol);
+void MAT_asPrincipalComponents_preallocated (MATVU result, constMATVU const& m, integer numberOfComponents) {
+	Melder_assert (numberOfComponents  > 0 && numberOfComponents <= m.ncol);
+	Melder_assert (result.nrow == m.nrow && result.ncol == numberOfComponents);
 	autoSVD svd = SVD_createFromGeneralMatrix (m);
-	MATVUmul (pc, m, svd -> v.verticalBand (1, pc.ncol));
+	MATmul (result, m, svd -> v.verticalBand (1, result.ncol));
 }
 
-autoMAT MAT_asPrincipalComponents (constMAT m, integer numberOfComponents) {
+autoMAT MAT_asPrincipalComponents (constMATVU m, integer numberOfComponents) {
 	Melder_assert (numberOfComponents  > 0 && numberOfComponents <= m.ncol);
 	autoMAT result = newMATraw (m.nrow, numberOfComponents);
-	MAT_asPrincipalComponents_preallocated (result.get(), m);
+	MAT_asPrincipalComponents_preallocated (result.get(), m, numberOfComponents);
 	return result;
 }
 
-void MATpseudoInverse_preallocated (MAT target, constMAT m, double tolerance) {
-	Melder_assert (target.nrow == m.ncol && target.ncol == m.nrow);
-	autoSVD me = SVD_createFromGeneralMatrix (m);
+void MATpseudoInverse (MATVU const& target, constMATVU const& mat, double tolerance) {
+	Melder_assert (target.nrow == mat.ncol && target.ncol == mat.nrow);
+	autoSVD me = SVD_createFromGeneralMatrix (mat);
 	(void) SVD_zeroSmallSingularValues (me.get(), tolerance);
-	for (integer i = 1; i <= m.ncol; i ++) {
-		for (integer j = 1; j <= m.nrow; j ++) {
-			longdouble s = 0.0;
-			for (integer k = 1; k <= m.ncol; k ++) {
-				if (my d [k] != 0.0) {
-					s += my v [i] [k] * my u [j] [k] / my d [k];
-				}
-			}
-			target [i] [j] = (double) s;
+	for (integer irow = 1; irow <= target.nrow; irow ++) {
+		for (integer icol = 1; icol <= target.ncol; icol ++) {
+			PAIRWISE_SUM (
+				longdouble, sum,
+				integer, mat.ncol,
+				integer k = 1,
+				my d [k] == 0.0 ? 0.0 : my v [irow] [k] * my u [icol] [k] / my d [k],
+				k += 1
+			)
+			target [irow] [icol] = double (sum);
 		}
 	}
 }
 
-autoMAT MATpseudoInverse (constMAT m, double tolerance) {
-	autoMAT result = newMATraw (m.ncol, m.nrow);
-	MATpseudoInverse_preallocated (result.get(), m, tolerance);
+autoMAT newMATpseudoInverse (constMATVU const& mat, double tolerance) {
+	autoMAT result = newMATraw (mat.ncol, mat.nrow);
+	MATpseudoInverse (result.all(), mat, tolerance);
 	return result;
 }
 

dwsys/MAT_numerics.h

 #pragma once
 /* MAT_numerics.h
  *
- * Copyright (C) 2018 David Weenink
+ * Copyright (C) 2018-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,7 @@ void MAT_getEigenSystemFromSymmetricMatrix (constMAT a, autoMAT *out_eigenvector
 		if (out_eigenvectors) eigenvectors corresponding to the eigenvalues, stored as row-wise vectors.
 */
 
-void MAT_getEigenSystemFromSymmetricMatrix_preallocated (MAT eigenvectors, VEC eigenvalues, constMAT a, bool sortAscending);
+void MAT_getEigenSystemFromSymmetricMatrix_preallocated (MAT eigenvectors, VEC eigenvalues, constMATVU const& a, bool sortAscending);
 /*
 	Input:
 		a, a symmetric a.ncol x a.ncol matrix
@@ -66,11 +66,11 @@ void MAT_eigenvectors_decompress (constMAT eigenvectors, constVEC eigenvalues_re
 	Decompresses each eigenvector row into two consecutive columns (real and imaginary part)
 */
 
-void MAT_asPrincipalComponents_preallocated (MAT pc, constMAT m);
-autoMAT MAT_asPrincipalComponents (constMAT m, integer numberOfComponents);
+void MAT_asPrincipalComponents_preallocated (MATVU pc, constMATVU const& m, integer numberOfComponents);
+autoMAT MAT_asPrincipalComponents (constMATVU m, integer numberOfComponents);
 
-void MATpseudoInverse_preallocated (MAT target, constMAT m, double tolerance);
-autoMAT MATpseudoInverse (constMAT m, double tolerance);
+void MATpseudoInverse (MATVU const& target, constMATVU const& mat, double tolerance);
+autoMAT newMATpseudoInverse (constMATVU const& mat, double tolerance);
 /*
 	Determines the pseudo-inverse Y^-1 of Y[1..nrow][1..ncol] via s.v.d.
 	Alternative notation for pseudo-inverse: (Y'.Y)^-1.Y'

dwsys/Makefile

@@ -13,7 +13,7 @@ OBJECTS = Collection_extensions.o Command.o \
 	FileInMemory.o FileInMemorySet.o FileInMemoryManager.o\
 	Graphics_extensions.o Index.o \
 	MAT_numerics.o \
-	NUM2.o NUMhuber.o NUMmachar.o \
+	NMF.o NUM2.o NUMhuber.o NUMmachar.o \
 	NUMf2c.o NUMcblas.o NUMclapack.o NUMcomplex.o NUMfft_d.o NUMsort2.o \
 	NUMmathlib.o NUMstring.o \
 	Permutation.o Permutation_and_Index.o \

dwsys/NMF.cpp

+/* NMF.cpp
+ *
+ * Copyright (C) 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
+ * the Free Software Foundation; either version 2 of the License, or (at
+ * your option) any later version.
+ *
+ * This code is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this work. If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include "NMF.h"
+#include "NUMmachar.h"
+#include "NUM2.h"
+#include "SVD.h"
+
+#include "oo_DESTROY.h"
+#include "NMF_def.h"
+#include "oo_COPY.h"
+#include "NMF_def.h"
+#include "oo_EQUAL.h"
+#include "NMF_def.h"
+#include "oo_CAN_WRITE_AS_ENCODING.h"
+#include "NMF_def.h"
+#include "oo_WRITE_TEXT.h"
+#include "NMF_def.h"
+#include "oo_WRITE_BINARY.h"
+#include "NMF_def.h"
+#include "oo_READ_TEXT.h"
+#include "NMF_def.h"
+#include "oo_READ_BINARY.h"
+#include "NMF_def.h"
+#include "oo_DESCRIPTION.h"
+#include "NMF_def.h"
+
+#include "enums_getText.h"
+#include "NMF_enums.h"
+#include "enums_getValue.h"
+#include "NMF_enums.h"
+
+void structNMF :: v_info () {
+	MelderInfo_writeLine (U"Number of rows: ", numberOfRows);
+	MelderInfo_writeLine (U"Number of columns: ", numberOfColumns);
+	MelderInfo_writeLine (U"Number of features: ", numberOfFeatures);
+}
+
+Thing_implement (NMF, Daata, 0);
+
+
+autoNMF NMF_create (integer numberOfRows, integer numberOfColumns, integer numberOfFeatures) {
+	try {
+		autoNMF me = Thing_new (NMF);
+		my numberOfRows = numberOfRows;
+		my numberOfColumns = numberOfColumns;
+		my numberOfFeatures = numberOfFeatures;
+		my features = newMATzero (numberOfRows, numberOfFeatures);
+		my weights = newMATzero (numberOfFeatures, numberOfColumns);
+		return me;
+	} catch (MelderError) {
+		Melder_throw (U"NMF not created.");
+	}
+}
+
+static void MATmakeElementsNonNegative (MATVU const& m, int strategy) {
+	for (integer irow = 1; irow <= m.nrow; irow ++)
+		for (integer icol = 1; icol <= m.ncol; icol ++)
+			if (m [irow][icol] < 0.0)
+				m [irow] [icol] = strategy == 0 ? 0.0 : fabs (m [irow] [icol]);	
+}
+
+static void NMF_initializeFactorization_svd (NMF me, constMATVU const& data, kNMF_Initialization initializationMethod) {
+	try {
+		autoSVD thee = SVD_createFromGeneralMatrix (data);
+		MATmakeElementsNonNegative (thy u.get(), 1);
+		MATmakeElementsNonNegative (thy v.get(), 1);
+		my features.get() <<= thy u.verticalBand (1, my numberOfFeatures);
+		for (integer irow = 1; irow <= my numberOfFeatures; irow ++)
+			my weights.row (irow) <<= thy d [irow]  *  thy v.row (irow);
+		
+	} catch (MelderError) {
+		Melder_throw (me, U": cpuld not initialize by svd method.");
+	}	
+}
+
+void NMF_initializeFactorization (NMF me, constMATVU const& data, kNMF_Initialization initializationMethod) {
+	if (initializationMethod == kNMF_Initialization::RandomUniform) {
+		double rmin = 0.0, rmax = 1.0;
+		MATrandomUniform (my features.all(), rmin, rmax);
+		MATrandomUniform (my weights.all(), rmin, rmax);
+	} else {
+		NMF_initializeFactorization_svd (me, data, initializationMethod);
+	}
+}
+
+autoNMF NMF_createFromGeneralMatrix (constMATVU const& m, integer numberOfFeatures) {
+	try {
+		Melder_require (NUMisNonNegative (m),
+			U"No matrix elements should be negative.");
+		Melder_require (numberOfFeatures <= m.ncol,
+			U"The number of features should not exceed the number of columns.");
+		autoNMF me = NMF_create (m.nrow, m.ncol, numberOfFeatures);
+		return me;
+	} catch (MelderError) {
+		Melder_throw (U"NMF cannot be created.");
+	}
+}
+
+double NMF_getEuclideanDistance (NMF me, constMATVU const& data) {
+	Melder_require (data.nrow == my numberOfRows && data.ncol == my numberOfColumns,
+		U"Dimensions should match.");
+	autoMAT synthesis = NMF_synthesize (me);
+	synthesis.get()  -=  data;
+	double dist = NUMnorm (synthesis.get(), 2.0);
+	return dist;
+}
+
+static double getMaximumChange (constMATVU const& m, MATVU const& m0, const double sqrteps) {
+	double min = NUMmin (m0);
+	double max = NUMmax (m0);
+	double extremum1 = std::max (fabs (min), fabs (max));
+	m0  -=  m;
+	min = NUMmin (m0);
+	max = NUMmax (m0);
+	double extremum2 = std::max (fabs (min), fabs (max));
+	double dmat = extremum2 / (sqrteps + extremum1);
+	return dmat;
+}
+
+/*
+	Calculating elementwise matrix multiplication, division and addition m = m0 .* (numerm ./(denom + eps))
+	Set elements < zero_threshold to zero
+*/
+
+static const void update (MATVU const& m, constMATVU const& m0, constMATVU const& numer, constMATVU const& denom, double zeroThreshold, double maximum) {
+	Melder_assert (m.nrow == m0.nrow && m.ncol == m0.ncol);
+	Melder_assert (m.nrow == numer.nrow && m.ncol == numer.ncol);
+	Melder_assert (m.nrow == denom.nrow && m.ncol == denom.ncol);
+	/*
+		The value 1e-9 is OK for matrices with values that are larger than 1.
+		For matrices with very small values we have to scale the divByZeroAvoidance value
+		otherwise the precision would suffer. 
+		A scaling with the maximum value seems reasonable.
+	*/
+	const double divByZeroAvoidance = 1e-09 * ( maximum < 1.0 ? maximum : 1.0 );
+	for (integer irow = 1; irow <= m.nrow; irow ++) 
+		for (integer icol = 1; icol <= m.ncol; icol++) {
+			if (m0 [irow] [icol] == 0.0 || numer [irow] [icol]  == 0.0) {
+				m [irow] [icol] = 0.0;
+			} else {
+				double update = m0 [irow] [icol] * (numer [irow] [icol] / (denom [irow] [icol] + divByZeroAvoidance));
+				m [irow] [icol] = ( update < zeroThreshold ? 0.0 : update );
+			}
+		}
+}
+
+/*
+	Algorithm for Non-negative Matrix Factorization by multiplicative updates.
+	The algoritm is inspired by the nmf_mu.c algorithm in libNMF by 
+	A. Janecek, S. Schulze Grotthoff & W.N. Gangsterer (2011): 
+		"LIBNMF - A library for nonnegative matrix factorization."
+		Computing and informatics% #30: 205--224.
+
+*/
+void NMF_improveFactorization_mu (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.");
+		
+		autoMAT productFtD = newMATzero (my numberOfFeatures, my numberOfColumns); // calculations of F'D
+		autoMAT productFtFW = newMATzero (my numberOfFeatures, my numberOfColumns); // calculations of F'F W
+		autoMAT weights0 = newMATzero (my numberOfFeatures, my numberOfColumns);
+		
+		autoMAT productDWt = newMATzero (my numberOfRows, my numberOfFeatures); // calculations of DW'
+		autoMAT productFWWt = newMATzero (my numberOfRows, my numberOfFeatures); // calculations of FWW'
+		autoMAT features0 = newMATzero (my numberOfRows, my numberOfFeatures);
+		
+		autoMAT productWWt = newMATzero (my numberOfFeatures, my numberOfFeatures); // calculations of WW'
+		autoMAT productFtF = newMATzero (my numberOfFeatures, my numberOfFeatures); // calculations of F'F
+		
+		double traceDtD = NUMtrace2 (data.transpose(), data); // for distance calculation
+		features0.get() <<= my features.get();
+		weights0.get() <<= my weights.get();
+		
+		if (! NUMfpp) NUMmachar ();
+		const double eps = NUMfpp -> eps;
+		const double sqrteps = sqrt (eps);
+		double dnorm0 = 0.0;
+		double maximum = NUMmax (data);
+		integer iter = 1;
+		bool convergence = false;		
+		
+		while (iter <= maximumNumberOfIterations && not convergence) {
+			/*
+				while iter < maxinter and not convergence
+					(1) W = W .* (F'*D) ./ (F'*F*W + 10^^−9^)
+					(2) F = F .* (D*W') ./ (F*W*W' + 10^^−9^)
+					(3) test for convergence
+				endwhile
+			*/
+			
+			// 1. Update W matrix
+			features0.get() <<= my features.get();
+			weights0.get() <<= my weights.get();
+			MATmul (productFtD.get(), features0.transpose(), data);
+			MATmul (productFtF.get(), features0.transpose(), features0.get());
+			MATmul (productFtFW.get(), productFtF.get(), weights0.get());
+			update (my weights.get(), weights0.get(), productFtD.get(), productFtFW.get(), eps, maximum);
+
+			// 2. Update F matrix
+			MATmul (productDWt.get(), data, my weights.transpose()); // productDWt = data*weights'
+			MATmul (productWWt.get(), my weights.get(), my weights.transpose()); // work1 = weights*weights'
+			MATmul (productFWWt.get(), features0.get(), productWWt.get()); // productFWWt = features0 * work1
+			update (my features.get(), features0.get(), productDWt.get(), productFWWt.get(), eps, maximum);
+			
+			/* 3. Convergence test:
+				The Frobenius norm ||D-FW|| of a matrix can be written as
+				||D-FW||=trace(D'D) − 2trace(W'F'D) + trace(W'F'FW)
+						=trace(D'D) - 2trace(W'(F'D))+trace((F'F)(WW'))
+				This saves us from explicitly calculating the reconstruction FW because we already have performed most of
+				the needed matrix multiplications in the update step.
+			*/
+			
+			double traceWtFtD  = NUMtrace2 (my weights.transpose(), productFtD.get());
+			double traceWtFtFW = NUMtrace2 (productFtF.get(), productWWt.get());
+			double distance = sqrt (std::max (traceDtD - 2.0 * traceWtFtD + traceWtFtFW, 0.0)); // just in case
+			double dnorm = distance / (my numberOfRows * my numberOfColumns);
+			double df = getMaximumChange (my features.get(), features0.get(), sqrteps);
+			double dw = getMaximumChange (my weights.get(), weights0.get(), sqrteps);
+			double delta = std::max (df, dw);
+			
+			convergence = ( iter > 1 && (delta < changeTolerance || dnorm < dnorm0 * approximationTolerance) );
+			if (info)
+				MelderInfo_writeLine (U"Iteration: ", iter, U", dnorm: ", dnorm, U", delta: ", delta);
+			
+			dnorm0 = dnorm;
+			++ iter;
+		}
+		if (info)
+			MelderInfo_drain();
+	} catch (MelderError) {
+		Melder_throw (me, U" factorization cannot be improved.");
+	}
+}
+
+void NMF_improveFactorization_als (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.");
+		
+		autoMAT productFtD = newMATzero (my numberOfFeatures, my numberOfColumns); // calculations of F'D
+		autoMAT productWDt = newMATzero (my numberOfFeatures, my numberOfRows); // calculations of WD'
+		
+		autoMAT weights0 = newMATzero (my numberOfFeatures, my numberOfColumns);
+		autoMAT features0 = newMATzero (my numberOfRows, my numberOfFeatures);
+
+		autoMAT productFtF = newMATzero (my numberOfFeatures, my numberOfFeatures); // calculations of F'F
+		autoMAT productWWt = newMATzero (my numberOfFeatures, my numberOfFeatures); // calculations of WW'
+		
+		autoSVD svd_WWt = SVD_create (my numberOfFeatures, my numberOfFeatures); // solving W*W'*F' = W*D'
+		autoSVD svd_FtF = SVD_create (my numberOfFeatures, my numberOfFeatures); // solving F´*F*W = F'*D
+				
+		double traceDtD = NUMtrace2 (data.transpose(), data); // for distance calculation
+		
+		if (! NUMfpp) NUMmachar ();
+		const double eps = NUMfpp -> eps;
+		const double sqrteps = sqrt (eps);
+		double dnorm0 = 0.0;
+		double maximum = NUMmax (data);
+		integer iter = 1;
+		bool convergence = false;
+		
+		while (iter <= maximumNumberOfIterations && not convergence) {
+			/*
+				for iter to maxiter
+					(1) W is solution of F´*F*W = F'*D.
+					    Set all negative elements in W to 0.
+					(2) F is solution of W*W'*F' = W*D' .
+					    Set all negative elements in F to 0.
+					(3) test for convergence
+				endfor
+			*/
+			
+			// 1. Solve equations for new W:  F´*F*W = F'*D.
+			weights0.get() <<= my weights.get(); // save previous weigts for convergence test
+			MATmul (productFtD.get(), my features.transpose(), data);
+			MATmul (productFtF.get(), my features.transpose(), my features.get());
+
+			svd_FtF -> u.get() <<= productFtF.get();
+			SVD_compute (svd_FtF.get());
+			SVD_solve_preallocated (svd_FtF.get(), productFtD.get(), my weights.get());
+			MATmakeElementsNonNegative (my weights.get(), 0);
+			
+			// 2. Solve equations for new F:  W*W'*F' = W*D'
+			features0.get() <<= my features.get(); // save previous features for convergence test
+			MATmul  (productWDt.get(), my weights.get(), data.transpose());
+			MATmul (productWWt.get(), my weights.get(), my weights.transpose());
+
+			svd_WWt -> u.get() <<= productWWt.get();
+			SVD_compute (svd_WWt.get());
+			SVD_solve_preallocated (svd_WWt.get(), productWDt.get(), my features.transpose());
+
+			// 3. Convergence test
+			double traceWtFtD  = NUMtrace2 (my weights.transpose(), productFtD.get());
+			double traceWtFtFW = NUMtrace2 (productFtF.get(), productWWt.get());
+			double distance = sqrt (std::max (traceDtD - 2.0 * traceWtFtD + traceWtFtFW, 0.0)); // just in case
+			double dnorm = distance / (my numberOfRows * my numberOfColumns);
+			double df = getMaximumChange (my features.get(), features0.get(), sqrteps);
+			double dw = getMaximumChange (my weights.get(), weights0.get(), sqrteps);
+			double delta = std::max (df, dw);
+			
+			convergence = ( iter > 1 && (delta < changeTolerance || dnorm < dnorm0 * approximationTolerance) );
+			if (info)
+				MelderInfo_writeLine (U"Iteration: ", iter, U", dnorm: ", dnorm, U", delta: ", delta);
+			dnorm0 = dnorm;
+			++ iter;
+		}
+		if (info)
+			MelderInfo_drain();
+	} catch (MelderError) {
+		Melder_throw (me, U" ALS factorization cannot be improved.");
+	}
+}
+
+autoMAT NMF_synthesize (NMF me) {
+	try {
+		autoMAT result = newMATmul (my features.get(), my weights.get());
+		return result;
+	} catch (MelderError) {
+		Melder_throw (me, U": No matrix created.");
+	}
+}
+
+/* End of file NMF.cpp */

dwsys/NMF.h

+/* NMF.h
+ *
+ * Copyright (C) 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
+ * the Free Software Foundation; either version 2 of the License, or (at
+ * your option) any later version.
+ *
+ * This code is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this work. If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef _NMF_h_
+#define _NMF_h_
+
+#include "NUM2.h"
+#include "Data.h"
+#include "NMF_enums.h"
+
+#include "NMF_def.h"
+
+/*
+	Non-negative Matrix Factorization. I.e. we try to factorize a non-negative nrow x ncol matrix M as H*W,
+	where H and W are non-negative nrow x k and k x ncol matrices, respectively. 
+*/
+
+autoNMF NMF_create (integer numberOfRows, integer numberOfColumns, integer numberOfFeatures);
+
+autoNMF NMF_createFromGeneralMatrix (constMATVU const& data, integer numberOfFeatures);
+
+void NMF_initializeFactorization (NMF me, constMATVU const& data, kNMF_Initialization initializationMethod);
+
+/*
+	Factorize D as F*W, where D, F and W >= 0
+	
+	initiaze F and W;
+	for iter to maxiter
+		W(n+1) = W(n).(F(n)'*D) / (F(n)'*F(n))*W(n) + eps)
+		F(n+1) = F(n).(D*W(n+1)') / (F(n)*(W(n+1)*W(n+1)') + eps)
+	endfor
+*/
+void NMF_improveFactorization_mu (NMF me, constMATVU const& data, integer maximumNumberOfIterations, double changeTolerance, double approximationTolerance, bool info);
+
+/*
+	Factorize D as F*W, where D, F and W >= 0
+	
+	inialize F;
+	for iter to maxiter
+		W is solution of F´*F*W = F'*D.
+		Set all negative elements in W to 0.
+		F is solution of W*W'*F' = W*D' .
+		Set all negative elements in F to 0.
+	endfor
+*/
+void NMF_improveFactorization_als (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||²)
+
+#endif /* _NMF_h_ */