.gitignore

@@ -15,3 +15,8 @@
 /praat_nogui
 /sendpraat
 /*.1
+/dwtest/kanweg.FFNet
+/test/fon/kanweg.txt
+/test/kar/kanweg.TextGrid
+/test/kar/kanweg_abc*.txt
+/test/fon/kanweg.Object

EEG/EEG.cpp

@@ -636,8 +636,7 @@ autoEEG EEGs_concatenate (OrderedOf<structEEG>* me) {
 			Melder_throw (U"Cannot concatenate zero EEG objects.");
 		EEG first = my at [1];
 		integer numberOfChannels = first -> numberOfChannels;
-		autostring32vector channelNames;
-		channelNames. copyFrom (first -> channelNames);
+		autostring32vector channelNames = STRVECclone (first -> channelNames.get());
 		for (integer ieeg = 2; ieeg <= my size; ieeg ++) {
 			EEG other = my at [ieeg];
 			if (other -> numberOfChannels != numberOfChannels)
@@ -672,7 +671,7 @@ autoEEG EEG_extractPart (EEG me, double tmin, double tmax, bool preserveTimes) {
 	try {
 		autoEEG thee = Thing_new (EEG);
 		thy numberOfChannels = my numberOfChannels;
-		thy channelNames. copyFrom (my channelNames);
+		thy channelNames = STRVECclone (my channelNames.get());
 		thy sound = Sound_extractPart (my sound.get(), tmin, tmax, kSound_windowShape::RECTANGULAR, 1.0, preserveTimes);
 		thy textgrid = TextGrid_extractPart (my textgrid.get(), tmin, tmax, preserveTimes);
 		thy xmin = thy textgrid -> xmin;
@@ -724,7 +723,7 @@ autoEEG EEG_MixingMatrix_to_EEG_unmix (EEG me, MixingMatrix you) {
 	his sound = Sound_MixingMatrix_unmix (my sound.get(), you);
 	his textgrid = Data_copy (my textgrid.get());
 	his numberOfChannels = your numberOfColumns;
-	his channelNames. copyFrom (your columnLabels);
+	his channelNames = STRVECclone (your columnLabels.get());
 	return him;
 }
 
@@ -742,7 +741,7 @@ autoEEG EEG_MixingMatrix_to_EEG_mix (EEG me, MixingMatrix you) {
 	his sound = Sound_MixingMatrix_mix (my sound.get(), you);
 	his textgrid = Data_copy (my textgrid.get());
 	his numberOfChannels = your numberOfRows;
-	his channelNames. copyFrom (your rowLabels);
+	his channelNames = STRVECclone (your rowLabels.get());
 	return him;
 }
 

EEG/ERPTier.cpp

@@ -47,10 +47,9 @@ Thing_implement (ERPTier, AnyTier, 0);
 
 integer ERPTier_getChannelNumber (ERPTier me, conststring32 channelName) {
 	for (integer ichan = 1; ichan <= my numberOfChannels; ichan ++) {
-		if (Melder_equ (my channelNames [ichan].get(), channelName)) {
+		if (Melder_equ (my channelNames [ichan].get(), channelName))
 			return ichan;
 	}
-	}
 	return 0;
 }
 
@@ -71,7 +70,7 @@ static autoERPTier EEG_PointProcess_to_ERPTier (EEG me, PointProcess events, dou
 		Function_init (thee.get(), fromTime, toTime);
 		thy numberOfChannels = my numberOfChannels - EEG_getNumberOfExtraSensors (me);
 		Melder_assert (thy numberOfChannels > 0);
-		thy channelNames. copyFrom (my channelNames);
+		thy channelNames = STRVECclone (my channelNames.get());
 		integer numberOfEvents = events -> nt;
 		double soundDuration = toTime - fromTime;
 		double samplingPeriod = my sound -> dx;
@@ -246,7 +245,7 @@ autoERP ERPTier_extractERP (ERPTier me, integer eventNumber) {
 				newChannel [isample] = oldChannel [isample];
 			}
 		}
-		thy channelNames. copyFrom (my channelNames);
+		thy channelNames = STRVECclone (my channelNames.get());
 		return thee;
 	} catch (MelderError) {
 		Melder_throw (me, U": ERP not extracted.");
@@ -281,7 +280,7 @@ autoERP ERPTier_to_ERP_mean (ERPTier me) {
 			}
 		}
 		Melder_assert (mean -> ny == my numberOfChannels);
-		mean -> channelNames. copyFrom (my channelNames);
+		mean -> channelNames = STRVECclone (my channelNames.get());
 		return mean;
 	} catch (MelderError) {
 		Melder_throw (me, U": mean not computed.");
@@ -298,7 +297,7 @@ autoERPTier ERPTier_extractEventsWhereColumn_number (ERPTier me, Table table, in
 		autoERPTier thee = Thing_new (ERPTier);
 		Function_init (thee.get(), my xmin, my xmax);
 		thy numberOfChannels = my numberOfChannels;
-		thy channelNames. copyFrom (my channelNames);
+		thy channelNames = STRVECclone (my channelNames.get());
 		for (integer ievent = 1; ievent <= my points.size; ievent ++) {
 			ERPPoint oldEvent = my points.at [ievent];
 			TableRow row = table -> rows.at [ievent];
@@ -327,7 +326,7 @@ autoERPTier ERPTier_extractEventsWhereColumn_string (ERPTier me, Table table,
 		autoERPTier thee = Thing_new (ERPTier);
 		Function_init (thee.get(), my xmin, my xmax);
 		thy numberOfChannels = my numberOfChannels;
-		thy channelNames. copyFrom (my channelNames);
+		thy channelNames = STRVECclone (my channelNames.get());
 		for (integer ievent = 1; ievent <= my points.size; ievent ++) {
 			ERPPoint oldEvent = my points.at [ievent];
 			TableRow row = table -> rows.at [ievent];

LPC/Sound_and_LPC.cpp

@@ -246,8 +246,10 @@ static int Sound_into_LPC_Frame_marple (Sound me, LPC_Frame thee, double tol1, d
 	}
 	e0 *= 2.0;
 	if (e0 == 0.0) {
-		m = 0; thy gain *= 0.5; /* because e0 is twice the energy */
-		thy nCoefficients = m; return 0; // warning no signal
+		m = 0;
+		thy gain *= 0.5; /* because e0 is twice the energy */
+		thy nCoefficients = m;
+		return 0; // warning no signal
 	}
 	double q1 = 1.0 / e0;
 	double q2 = q1 * x [1], q = q1 * x [1] * x [1], w = q1 * x [n] * x [n];
@@ -370,7 +372,7 @@ end:
 
 static autoLPC _Sound_to_LPC (Sound me, int predictionOrder, double analysisWidth, double dt, double preEmphasisFrequency, int method, double tol1, double tol2) {
 	double t1, samplingFrequency = 1.0 / my dx;
-	double windowDuration = 2 * analysisWidth; /* gaussian window */
+	double windowDuration = 2.0 * analysisWidth; /* gaussian window */
 	integer numberOfFrames, frameErrorCount = 0;
 	Melder_require (Melder_roundDown (windowDuration / my dx) > predictionOrder, 
 		U"Analysis window duration too short.\n For a prediction order of ", predictionOrder,
@@ -388,12 +390,12 @@ static autoLPC _Sound_to_LPC (Sound me, int predictionOrder, double analysisWidt
 
 	autoMelderProgress progress (U"LPC analysis");
 
-	if (preEmphasisFrequency < samplingFrequency / 2) {
+	if (preEmphasisFrequency < samplingFrequency / 2.0) {
 		Sound_preEmphasis (sound.get(), preEmphasisFrequency);
 	}
 
 	for (integer i = 1; i <= numberOfFrames; i ++) {
-		LPC_Frame lpcframe = (LPC_Frame) & thy d_frames [i];
+		LPC_Frame lpcframe = & thy d_frames [i];
 		double t = Sampled_indexToX (thee.get(), i);
 		LPC_Frame_init (lpcframe, predictionOrder);
 		Sound_into_Sound (sound.get(), sframe.get(), t - windowDuration / 2);
@@ -416,10 +418,9 @@ static autoLPC _Sound_to_LPC (Sound me, int predictionOrder, double analysisWidt
 				frameErrorCount ++;
 			}
 		}
-		if ((i % 10) == 1) {
+		if (i % 10 == 1)
 			Melder_progress ( (double) i / numberOfFrames, U"LPC analysis of frame ", i, U" out of ", numberOfFrames, U".");
 	}
-	}
 	return thee;
 }
 

artsynth/Artword.cpp

 /* Artword.cpp
  *
- * Copyright (C) 1992-2009,2011,2015-2017 Paul Boersma
+ * Copyright (C) 1992-2009,2011,2015-2018 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
@@ -49,10 +49,8 @@ autoArtword Artword_create (double totalTime) {
 
 void Artword_setDefault (Artword me, kArt_muscle muscle) {
 	ArtwordData f = & my data [(int) muscle];
-	NUMvector_free <double> (f -> times, 1);
-	NUMvector_free <double> (f -> targets, 1);
-	f -> times = NUMvector <double> (1, 2);
-	f -> targets = NUMvector <double> (1, 2);
+	f -> times = VECzero (2);
+	f -> targets = VECzero (2);
 	f -> numberOfTargets = 2;
 	f -> times [1] = 0.0;
 	f -> targets [1] = 0.0;
@@ -61,49 +59,52 @@ void Artword_setDefault (Artword me, kArt_muscle muscle) {
 	f -> _iTarget = 1;
 }
 
+static void ArtwordData_setTarget (ArtwordData me, double time, double target) {
+	Melder_assert (my numberOfTargets >= 2);
+	int32 insertionPosition = 1;   // should be able to go up to 32768
+	while (insertionPosition <= my numberOfTargets && my times [insertionPosition] < time)
+		insertionPosition ++;
+	Melder_assert (insertionPosition <= my numberOfTargets);   // can never insert past totalTime
+	if (my times [insertionPosition] != time) {
+		if (my numberOfTargets == INT16_MAX)
+			Melder_throw (U"An Artword cannot have more than ", INT16_MAX, U" targets.");
+		my times.insert (insertionPosition);
+		my targets.insert (insertionPosition);
+		my numberOfTargets ++;   // maintain invariant
+	}
+	my targets [insertionPosition] = target;
+	my times [insertionPosition] = time;
+}
+
 void Artword_setTarget (Artword me, kArt_muscle muscle, double time, double target) {
 	try {
 		Melder_assert ((int) muscle >= 1);
 		Melder_assert ((int) muscle <= (int) kArt_muscle::MAX);
 		Melder_assert (muscle <= kArt_muscle::MAX);
-		ArtwordData f = & my data [(int) muscle];
-		Melder_assert (f -> numberOfTargets >= 2);
-		int32 insertionPosition = 1;   // should be able to go up to 32768
 		if (time < 0.0) time = 0.0;
 		if (time > my totalTime) time = my totalTime;
-		while (insertionPosition <= f -> numberOfTargets && f -> times [insertionPosition] < time)
-			insertionPosition ++;
-		Melder_assert (insertionPosition <= f -> numberOfTargets);   // can never insert past totalTime
-		if (f -> times [insertionPosition] != time) {
-			if (f -> numberOfTargets == INT16_MAX)
-				Melder_throw (U"An Artword cannot have more than ", INT16_MAX, U" targets.");
-			integer numberOfTargets = f -> numberOfTargets;
-			NUMvector_insert <double> (& f -> times, 1, & numberOfTargets, insertionPosition);
-			numberOfTargets = f -> numberOfTargets;
-			NUMvector_insert <double> (& f -> targets, 1, & numberOfTargets, insertionPosition);
-			f -> numberOfTargets ++;
-		}
-		f -> targets [insertionPosition] = target;
-		f -> times [insertionPosition] = time;
+		ArtwordData_setTarget (& my data [(int) muscle], time, target);
 	} catch (MelderError) {
 		Melder_throw (me, U": target not set.");
 	}
 }
 
-double Artword_getTarget (Artword me, kArt_muscle muscle, double time) {
-	ArtwordData f = & my data [(int) muscle];
-	double *times = f -> times, *targets = f -> targets;
-	int16 targetNumber = f -> _iTarget;
+static double ArtwordData_getTarget (ArtwordData me, double time) {
+	int16 targetNumber = my _iTarget;
 	if (! targetNumber) targetNumber = 1;
-	while (time > times [targetNumber + 1] && targetNumber < f -> numberOfTargets - 1)
+	while (time > my times [targetNumber + 1] && targetNumber < my numberOfTargets - 1)
 		targetNumber ++;
-	while (time < times [targetNumber] && targetNumber > 1)
+	while (time < my times [targetNumber] && targetNumber > 1)
 		targetNumber --;
-	f -> _iTarget = targetNumber;
-	Melder_assert (targetNumber > 0 && targetNumber < f -> numberOfTargets);
-	return targets [targetNumber] + (time - times [targetNumber]) *
-		(targets [targetNumber + 1] - targets [targetNumber]) /
-		(times [targetNumber + 1] - times [targetNumber]);
+	my _iTarget = targetNumber;
+	Melder_assert (targetNumber > 0 && targetNumber < my numberOfTargets);
+	return my targets [targetNumber] + (time - my times [targetNumber]) *
+		(my targets [targetNumber + 1] - my targets [targetNumber]) /
+		(my times [targetNumber + 1] - my times [targetNumber]);
+}
+
+double Artword_getTarget (Artword me, kArt_muscle muscle, double time) {
+	return ArtwordData_getTarget (& my data [(int) muscle], time);
 }
 
 void Artword_removeTarget (Artword me, kArt_muscle muscle, int16 targetNumber) {
@@ -133,8 +134,7 @@ void Artword_intoArt (Artword me, Art art, double time) {
 void Artword_draw (Artword me, Graphics g, kArt_muscle muscle, bool garnish) {
 	int16 numberOfTargets = my data [(int) muscle]. numberOfTargets;
 	if (numberOfTargets > 0) {
-		autoNUMvector <double> x (1, numberOfTargets);
-		autoNUMvector <double> y (1, numberOfTargets);
+		auto x = VECraw (numberOfTargets), y = VECraw (numberOfTargets);
 		Graphics_setInner (g);
 		Graphics_setWindow (g, 0, my totalTime, -1.0, 1.0);
 		for (int16 i = 1; i <= numberOfTargets; i ++) {

artsynth/Artword_def.h

 /* Artword_def.h
  *
- * Copyright (C) 1992-2005,2008,2009,2011,2015-2017 Paul Boersma
+ * Copyright (C) 1992-2005,2008,2009,2011,2015-2018 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
@@ -21,8 +21,8 @@
 oo_DEFINE_STRUCT (ArtwordData)
 
 	oo_INT16 (numberOfTargets)
-	oo_DOUBLE_VECTOR (targets, numberOfTargets)
-	oo_DOUBLE_VECTOR (times, numberOfTargets)
+	oo_VEC (targets, numberOfTargets)
+	oo_VEC (times, numberOfTargets)
 
 	#if oo_DECLARING
 		oo_INT16 (_iTarget)

artsynth/Speaker.cpp

@@ -112,7 +112,7 @@ autoSpeaker Speaker_create (conststring32 kindOfSpeaker, int16 numberOfVocalCord
 
 	my nose.Dx = 0.007 * scaling;
 	my nose.Dz = 0.014 * scaling;
-	my nose.weq = NUMvector <double> (1, 14);
+	my nose.weq = VECraw (14);
 	my nose.weq [1] = 0.018 * scaling;
 	my nose.weq [2] = 0.016 * scaling;
 	my nose.weq [3] = 0.014 * scaling;

artsynth/Speaker_def.h

 /* Speaker_def.h
  *
- * Copyright (C) 1992-2005,2011,2015-2017 Paul Boersma
+ * Copyright (C) 1992-2005,2011,2015-2018 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
@@ -135,7 +135,7 @@ oo_DEFINE_STRUCT (Speaker_Nose)
 
 	oo_DOUBLE (Dx)
 	oo_DOUBLE (Dz)
-	oo_DOUBLE_VECTOR (weq, 14)
+	oo_VEC (weq, 14)
 
 oo_END_STRUCT (Speaker_Nose)
 #undef ooSTRUCT

debian/changelog

+praat (6.0.43-1) unstable; urgency=medium
+
+  * New upstream version 6.0.43
+  * d/copyright: Reflect upstream changes
+  * d/p/melder-tensor-resize.patch: New patch
+  * d/rules, d/clean: Do file removal in the appropriate place
+
+ -- Rafael Laboissiere <rafael@debian.org>  Tue, 11 Sep 2018 07:35:19 -0300
+
 praat (6.0.42-2) unstable; urgency=medium
 
   * d/t/run-tests: Exclude two problematic unit tests

debian/clean

@@ -2,8 +2,17 @@ dwtest/kanweg.FFNet
 dwtest/kanweg.SpeechSynthesizer
 test/fon/kanweg*.txt
 test/fon/kanweg.nsp
+test/fon/kanweg.Object
 test/kar/kanweg*.txt
 test/kar/kanweg.TextGrid
 test/sys/kanweg.txt
 test/sys/kanweg.wav
 test/sys/kanweg2.txt
+makefile.defs
+debian/praat.dirs
+debian/praat.install
+praat_nogui sendpraat
+praat.1
+praat_nogui.1
+sendpraat.1
+praat-open-files.1

debian/copyright

@@ -11,6 +11,14 @@ Copyright: 1990-2018 Paul Boersma
            2008, 2017 Stefan de Konink
            2010 Franz Brausse
            2013 Tom Naughton
+           2007-2009 Ola Söder
+           2000-2004 Underbit Technologies, Inc
+           2007 Erez Volk
+           2007 Timothy A. Davis, Patrick R. Amestoy, and Iain S. Duff
+           2005 Rafael Laboissière
+           1998, Ross Ihaka
+           2000-2007, The R Core Team
+           2002 The NEdit Developers
 License: GPL-2+
 
 Files: external/glpk/*
@@ -51,11 +59,6 @@ Files: external/flac/flac_share_alloc.h
 Copyright: 2007 Josh Coalson
 License: LGPL-2.1+
 
-Files: external/mp3/*
-Copyright: 2000-2004 Underbit Technologies, Inc
-           2007 Erez Volk
-License: GPL-2+
-
 Files: external/portaudio*/*
 Copyright: 1999-2011 Ross Bencina
            1999-2011 Phil Burk
@@ -89,32 +92,15 @@ License: MIT
  CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
  WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
 
-Files: external/glpk/amd*
-Copyright: 2007 Timothy A. Davis, Patrick R. Amestoy, and Iain S. Duff
-License: GPL-2+
-
 Files: external/glpk/colamd.c
 Copyright: 1998-2007 Timothy A. Davis
 License: LGPL-2.1+
 
 Files: external/espeak/*
-Copyright: 2005-2014 Jonathan Duddington
-           2010 Nicolas Pitre
-           2010 Bill Cox
-           2007 Gilles Casse
-           2008 Sun Microsystems, Inc.
-           1993-1994 Jon Iles and Nick Ing-Simmons
+Copyright: 2005-2015 Jonathan Duddington
+           2008-2016 Reece H. Dunn
 License: GPL-3+
 
-Files: external/espeak/espeakdata_FileInMemory.*
-Copyright: 2012 David Weenink
-License: GPL-2+
-
-Files: contrib/ola/*
-Copyright: 2007-2009 Ola Söder
-           2010-2011, 2015 Paul Boersma
-License: GPL-2+
-
 Files: external/gsl/*
 Copyright: 1996-2004, 2006-2007 Brian Gough
            1996-2006 Gerard Jungman
@@ -163,19 +149,6 @@ Copyright: 1996-2004, 2006-2007 Brian Gough
            1996-2000 Tim Mooney
 License: GPL-3+
 
-Files: kar/ipaSerifRegularPS.cpp
-Copyright: 2005 Rafael Laboissière
-License: GPL-2+
-
-Files: melder/regularExp.h
-Copyright: 2002 The NEdit Developers
-License: GPL-2+
-
-Files: dwsys/NUMmathlib.cpp
-Copyright: 1998, Ross Ihaka
-           2000-2007, The R Core Team
-License: GPL-2+
-
 Files: debian/*
 Copyright: 2005-2009, 2012-2018 Rafael Laboissière <rafael@debian.org>
            2009-2011 Andreas Tille <tille@debian.org>

debian/patches/melder-tensor-resize.patch

+Description: Fix melder tensor resize
+ This is taken from the upstream commit 95af7dd292b3b09c8d6b27ee0af24cc83fdf908d
+ and fixes a bug introduced in version 6.0.43 that makes tthe unit
+ test fon/PointProcess.praat hang forever.
+Origin: upstream, https://github.com/praat/praat/commit/95af7dd292b3b09c8d6b27ee0af24cc83fdf908d
+Forwarded: not-needed
+Last-Update: 2018-09-10
+
+--- praat-6.0.43.orig/melder/melder_tensor.h
++++ praat-6.0.43/melder/melder_tensor.h
+@@ -470,10 +470,17 @@ public:
+ 		elsewhere than at the head of the vector,
+ 		you should shift the elements after resizing.
+ 	*/
+-	bool resize (integer newSize, integer *inout_capacity,
++	void resize (integer newSize, integer *inout_capacity = nullptr,
+ 		kTensorInitializationType initializationType = kTensorInitializationType::ZERO)
+ 	{
+-		if (newSize > *inout_capacity) {
++		const integer currentCapacity = ( inout_capacity ? *inout_capacity : our size );
++		if (newSize > currentCapacity) {
++			/*
++				The new capacity is at least twice the old capacity.
++				When starting at a capacity of 0, and continually upsizing by one,
++				the capacity sequence will be: 0, 11, 33, 77, 165, 341, 693, 1397,
++				2805, 5621, 11253, 22517, 45045, 90101, 180213, 360437, 720885...
++			*/
+ 			integer newCapacity = newSize + our size + 10;   // this is at least a doubling!
+ 			/*
+ 				Create without change.
+@@ -486,9 +493,10 @@ public:
+ 				newAt [i] = our at [i];
+ 			if (our at) NUMvector_free (our at, 1);
+ 			our at = newAt;
+-			*inout_capacity = newCapacity;
++			if (inout_capacity)
++				*inout_capacity = newCapacity;
+ 		}
+-		our size = newSize;   // shrink below capacity
++		our size = newSize;
+ 	}
+ };
+ 

debian/patches/series

 use-ldflags.patch
 remove-time-date-macros.patch
 fix-procrustes-unit-test.patch
+melder-tensor-resize.patch

debian/rules

@@ -35,11 +35,6 @@ clean:
 	dh clean
 	# Fix the lacking cleaning command for the main/ directory
 	[ ! -f makefile.defs ] || make -C main clean
-	rm -f makefile.defs $(DEBDIR)/praat.dirs	\
-		 $(DEBDIR)/praat.install		\
-		 praat_nogui sendpraat			\
-		 praat.1 praat_nogui.1			\
-		 sendpraat.1 praat-open-files.1
 
 override_dh_auto_clean:
 	touch makefile.defs

dwsys/Collection_extensions.cpp

@@ -175,8 +175,8 @@ void OrderedOfString_changeStrings (StringList me, char32 *search, char32 *repla
 			SimpleString ss = my at [i];
 			integer nmatches_sub;
 			autostring32 r = use_regexp ?
-				replace_regexStr (ss -> string.get(), compiled_search, replace, maximumNumberOfReplaces, & nmatches_sub) :
-				replaceStr (ss -> string.get(), search, replace, maximumNumberOfReplaces, & nmatches_sub);
+				STRreplace_regex (ss -> string.get(), compiled_search, replace, maximumNumberOfReplaces, & nmatches_sub) :
+				STRreplace (ss -> string.get(), search, replace, maximumNumberOfReplaces, & nmatches_sub);
 
 			/*
 				Change without error.

dwsys/Eigen.cpp

@@ -32,11 +32,12 @@
  djmw 20040622 Less horizontal labels in Eigen_drawEigenvector.
  djmw 20050706 Shortened horizontal offsets in Eigen_drawEigenvalues from 1 to 0.5
  djmw 20051204 Eigen_initFromSquareRoot adapted for nrows < ncols
- djmw 20071012 Added: o_CAN_WRITE_AS_ENCODING.h
+ djmw 20071012 Added: oo_CAN_WRITE_AS_ENCODING.h
  djmw 20110304 Thing_new
 */
 
 #include "Eigen.h"
+#include "MAT_numerics.h"
 #include "NUMmachar.h"
 #include "NUMlapack.h"
 #include "NUMclapack.h"
@@ -213,41 +214,15 @@ void Eigen_initFromSquareRootPair (Eigen me, double **a, integer numberOfRows, i
 	NUMnormalizeRows (my eigenvectors, my numberOfEigenvalues, numberOfColumns, 1);
 }
 
-void Eigen_initFromSymmetricMatrix (Eigen me, double **a, integer n) {
-	double wt[1], temp;
-	char jobz = 'V', uplo = 'U';
-	integer lwork = -1, info;
-
-	my dimension = my numberOfEigenvalues = n;
-
+void Eigen_initFromSymmetricMatrix (Eigen me, constMAT a) {
+	Melder_assert (a.ncol == a.nrow);
 	if (! my eigenvectors) {
-		Eigen_init (me, n, n);
-	}
-
-	NUMmatrix_copyElements (a, my eigenvectors, 1, n, 1, n);
-
-	// Get size of work array
-
-	(void) NUMlapack_dsyev (& jobz, & uplo, & n, & my eigenvectors [1] [1], & n, & my eigenvalues [1], wt, & lwork, & info);
-	Melder_require (info == 0, U"dsyev initialization fails");
-	
-
-	lwork = Melder_ifloor (wt [0]);
-	autoNUMvector <double> work ((integer) 0, lwork);
-
-	(void) NUMlapack_dsyev (&jobz, &uplo, &n, &my eigenvectors[1][1], &n, &my eigenvalues[1], work.peek(), & lwork, & info);
-	Melder_require (info == 0, U"dsyev fails");
-
-	// We want descending order instead of ascending.
-
-	for (integer i = 1; i <= n / 2; i ++) {
-		integer ilast = n - i + 1;
-
-		SWAP (my eigenvalues [i], my eigenvalues [ilast])
-		for (integer j = 1; j <= n; j ++) {
-			SWAP (my eigenvectors [i] [j], my eigenvectors [ilast] [j])
-		}
+		Eigen_init (me, a.ncol, a.ncol);
 	}	
+	MAT eigenvectors (my eigenvectors, my numberOfEigenvalues, my dimension);
+	VEC eigenvalues (my eigenvalues, my numberOfEigenvalues);
+	MATtranspose_preallocated (eigenvectors, a);
+	MAT_getEigenSystemFromSymmetricMatrix_inplace (eigenvectors, true, eigenvalues, false);
 }
 
 autoEigen Eigen_create (integer numberOfEigenvalues, integer dimension) {

dwsys/Eigen.h

@@ -28,7 +28,7 @@ autoEigen Eigen_create (integer numberOfEigenvalues, integer dimension);
 
 void Eigen_init (Eigen me, integer numberOfEigenvalues, integer dimension);
 
-void Eigen_initFromSymmetricMatrix (Eigen me, double **a, integer n);
+void Eigen_initFromSymmetricMatrix (Eigen me, constMAT a);
 
 void Eigen_initFromSquareRoot (Eigen me, double **a, integer numberOfRows, integer numberOfColumns);
 /*

dwsys/MAT_numerics.cpp

+/* MAT_numerics.cpp
+ *
+ * Copyright (C) 2018 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 "melder.h"
+#include "NUMclapack.h"
+#include "MAT_numerics.h"
+
+#define SWAP(x,y) { tmp = x; x = y; y = tmp; }
+
+void MAT_getEigenSystemFromSymmetricMatrix_inplace (MAT a, bool wantEigenvectors, VEC eigenvalues, bool sortAscending) {
+	Melder_assert (a.nrow == a.ncol);
+	Melder_assert (eigenvalues.size == a.ncol);
+
+	char jobz = wantEigenvectors ? 'V' : 'N', uplo = 'U';
+	integer workSize = -1, n = a.ncol, info;
+	double wt [1], tmp;
+	
+	// 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, & n, & a [1] [1], & n, & eigenvalues [1], wt, & workSize, & info);
+	Melder_require (info == 0, U"dsyev initialization code = ", info, U").");
+	
+	workSize = Melder_iceiling (wt [0]);
+	autoVEC work = VECraw (workSize);
+
+	// 2. Calculate the eigenvalues and eigenvectors (row-wise)
+	
+	(void) NUMlapack_dsyev (& jobz, & uplo, & n, & a [1] [1], & n, & eigenvalues [1], & work [1], & workSize, & info);
+	Melder_require (info == 0, U"dsyev code = ", info, U").");
+
+	// 3. Eigenvalues are returned in ascending order
+
+	if (! sortAscending) {
+		for (integer i = 1; i <= n / 2; i ++) {
+			integer ilast = n - i + 1;
+			SWAP (eigenvalues [i], eigenvalues [ilast])
+			if (wantEigenvectors) {
+				for (integer j = 1; j <= n; j ++) {
+					SWAP (a [i] [j], a [ilast] [j])
+				}
+			}
+		}
+	}
+}
+
+void MAT_getEigenSystemFromSymmetricMatrix (constMAT a, autoMAT *out_eigenvectors, autoVEC *out_eigenvalues, bool sortAscending) {
+	Melder_assert (a.nrow == a.ncol);	
+	autoVEC eigenvalues = VECraw (a.nrow);	
+	autoMAT eigenvectors = MATtranspose (a);	
+	
+	bool wantEigenvectors = out_eigenvectors != nullptr;
+	MAT_getEigenSystemFromSymmetricMatrix_inplace (eigenvectors.get(), wantEigenvectors, eigenvalues.get(), sortAscending);
+	
+	if (out_eigenvectors) *out_eigenvectors = eigenvectors.move ();
+	if (out_eigenvalues) *out_eigenvalues = eigenvalues.move ();
+}
+
+void MAT_eigenvectors_decompress (constMAT eigenvectors, constVEC eigenvalues_re, constVEC eigenvalues_im, autoMAT *out_eigenvectors_reim) {
+	integer n = eigenvalues_re.size;	
+	Melder_assert (eigenvalues_im.size == n);
+	Melder_assert (eigenvectors.nrow == n && eigenvectors.ncol == eigenvectors.nrow);
+		
+	autoMAT eigenvectors_reim = MATzero (n, 2 * n);
+	integer pair_index = 0;
+	for (integer ivec = 1; ivec <= eigenvalues_re.size; ivec ++) {
+		// eigenvalues of a real matrix are either real or occur in complex conjugate pairs
+		
+		if (eigenvalues_im [ivec] == 0.0) { // real eigenvalue
+			for (integer j = 1; j <= n; j ++)
+				eigenvectors_reim [j] [2 * ivec - 1] = eigenvectors [ivec] [j];
+		} else if (ivec > 1 && eigenvalues_re [ivec] == eigenvalues_re [ivec - 1] && 
+			eigenvalues_im [ivec] == -eigenvalues_im [ivec - 1] && ivec - pair_index > 1) {
+			for (integer j = 1; j <= n; j ++) {
+				eigenvectors_reim [j] [2 * (ivec - 1) - 1]= eigenvectors [ivec - 1] [j];
+				eigenvectors_reim [j] [2 * (ivec - 1)] = eigenvectors [ivec] [j];
+				eigenvectors_reim [j] [2 * ivec - 1]  = eigenvectors [ivec - 1] [j];
+				eigenvectors_reim [j] [2 * ivec] = eigenvectors [ivec] [j] != 0.0 ? -eigenvectors [ivec] [j] : 0.0; // avoid -0
+			}
+			pair_index = ivec; // guard for multiple equal complex conjugate pairs
+		}
+	}
+	if (out_eigenvectors_reim) *out_eigenvectors_reim = eigenvectors_reim.move();
+}
+
+void MAT_getEigenSystemFromGeneralMatrix (constMAT a, autoMAT *out_lefteigenvectors, autoMAT *out_righteigenvectors, autoVEC *out_eigenvalues_re, autoVEC *out_eigenvalues_im) {
+	Melder_assert (a.nrow == a.ncol);
+	integer n = a.nrow, info, workSize = -1;
+	char jobvl = out_lefteigenvectors ? 'V' : 'N';
+	integer left_nvecs = out_lefteigenvectors ? n : 1; // 1x1 if not wanted
+	char jobvr = out_righteigenvectors ? 'V' : 'N';
+	integer right_nvecs = out_righteigenvectors ? n : 1;
+	double wt;
+	
+	autoMAT data = MATraw (n, n);
+	MATtranspose_preallocated (data.get(), a); // lapack is fortran storage
+	autoVEC eigenvalues_re = VECraw (n);
+	autoVEC eigenvalues_im = VECraw(n);
+	autoMAT righteigenvectors = MATraw (right_nvecs, right_nvecs); // 1x1 if not needed
+	autoMAT lefteigenvectors = MATraw (left_nvecs, left_nvecs); // 1x1 if not needed
+	
+	(void) NUMlapack_dgeev (& jobvl, & jobvr, & n, & data [1] [1], & n,	& eigenvalues_re [1], & eigenvalues_im [1],	& lefteigenvectors [1] [1],
+		& n, & righteigenvectors [1] [1], & n, & wt, & workSize, & info);
+	Melder_require (info == 0, U"dgeev initialization code = ", info, U").");
+	
+	workSize = Melder_iceiling (wt);
+	autoVEC work = VECraw (workSize);
+	
+	(void) NUMlapack_dgeev (& jobvl, & jobvr, & n, & data [1] [1], & n, & eigenvalues_re [1], & eigenvalues_im [1],	& lefteigenvectors [1] [1],
+		& n, & righteigenvectors [1] [1], & n, & work [1], & workSize, & info);
+	Melder_require (info == 0, U"dgeev code = ", info, U").");
+	
+	if (out_righteigenvectors) *out_righteigenvectors = righteigenvectors.move();
+	if (out_lefteigenvectors) *out_lefteigenvectors = lefteigenvectors.move();
+	if (out_eigenvalues_re) *out_eigenvalues_re = eigenvalues_re.move();
+	if (out_eigenvalues_im) *out_eigenvalues_im = eigenvalues_im.move();
+}
+
+void MAT_getPrincipalComponentsOfSymmetricMatrix_inplace (constMAT a, integer nComponents, MAT pc) {
+	Melder_assert (a.nrow == a.ncol);
+	Melder_assert (nComponents > 0 && nComponents <= a.ncol);
+	Melder_assert (pc.nrow == a.nrow && pc.ncol == nComponents);
+	
+	autoMAT eigenvectors = MATraw (a.nrow, a.nrow);
+	MAT_getEigenSystemFromSymmetricMatrix (a, & eigenvectors, nullptr, false);
+
+	for (integer i = 1; i <= a.nrow; i ++) {
+		for (integer j = 1; j <= nComponents; j ++) {
+			longdouble s = 0.0;
+			for (integer k = 1; k <= a.nrow; k ++) {
+				s += a [k] [i] * eigenvectors [k] [j];
+			}
+			pc [i] [j] = (double) s;
+		}
+	}
+}
+
+bool MAT_isSymmetric (constMAT x) {
+	bool isSymmetric = false;
+	if (x.nrow == x.ncol) {
+		for (integer i = 1; i <= x.nrow; i ++)
+			for (integer j = i + 1; j < x.ncol; j ++)
+				if (x [i] [j] != x [j] [i])
+					return false;
+		isSymmetric = true;
+	}
+	return isSymmetric;
+}
+
+#undef SWAP
+
+/* End of file MAT_numerics.cpp */

dwsys/MAT_numerics.h

+#pragma once
+/* MAT_numerics.h
+ *
+ * Copyright (C) 2018 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/>.
+ */
+
+bool MAT_isSymmetric (constMAT x);
+/*
+	Returns true if the matrix is symmetric else false.
+*/
+
+void MAT_getEigenSystemFromSymmetricMatrix (constMAT a, autoMAT *out_eigenvectors, autoVEC *out_eigenvalues, bool sortAscending);
+/*
+	Calculates the eigenvectors and eigenvalues of a symmetric matrix;
+	Input:
+		a, a symmetric matrix of size nrow x nrow (a.ncol == a.nrow)
+		   only the upper-half of the matrix is used in the calculation
+		sortAscending if true eigenvalues (and corresponding eigenvectors) are sorted ascending
+	Output:
+		if(out_eigenvalues) eigenvalues sorted
+		if(out_eigenvectors) eigenvectors corresponding to the eigenvalues, stored as row-wise vectors.
+*/
+
+void MAT_getEigenSystemFromSymmetricMatrix_inplace (MAT inout_a, bool wantEigenvectors, VEC inout_eigenvalues, bool sortAscending);
+/*
+	Input:
+		inout_a, a symmetric a.ncol x a.ncol matrix
+		   only the upper-half of the matrix is used in the calculation
+		inout_eigenvalues, a vector of size ncol
+		wantEigenvectors, true if you want eigenvectors calculated
+		sortAscending if true eigenvalues (and corresponding eigenvectors) are sorted ascending
+	Output:
+		inout_a, if (wantEigenvectors) eigenvectors, stored row-wise
+		inout_eigenvalues, eigenvalues sorted from large to small
+*/
+
+void MAT_getPrincipalComponentsOfSymmetricMatrix_inplace (constMAT a, integer nComponents, MAT inout_pc);
+/*
+	Input:
+		a, a symmetric nrow x nrow matrix
+		   only the upper-half of the matrix is used in the calculation
+		nComponents, the number of components to determine (1 <= nComponents <= a.nrow)
+		inout_pc, a a.nrow x nComponents matrix
+	Output:
+		inout_pc, the principal components stored column-wise
+*/
+
+void MAT_getEigenSystemFromGeneralMatrix (constMAT a, autoMAT *out_lefteigenvectors, autoMAT *out_righteigenvectors, autoVEC *out_eigenvalues_re, autoVEC *out_eigenvalues_im);
+/* no standard sorting with complex numbers.
+	Compute eigenvalues of general nxn matrix with optionally the left/right eigenvectors.
+	There is no standard sorting with complex numbers.
+	Input:
+	Output:
+		out_eigenvalues_re, out_eigenvalues_im
+			the real and imaginary parts of the eigenvalues.
+			Complex conjugate pairs of eigenvalues appear consecutively
+            with the eigenvalue having the positive imaginary part first.
+		out_righteigenvectors, out_lefteigenvectors:
+			the left and right eigenvectors (stored row-wise!, compressed)
+			if the j-th eigenvalue is real the eigenvector is real and in row j.
+			if j and j+1 form a complex conjugate pair, the two eigenvectors are 
+			complex conjugates, whose real part is in j and its imaginary part in j+1.
+*/
+
+void MAT_eigenvectors_decompress (constMAT eigenvectors, constVEC eigenvalues_re, constVEC eigenvalues_im, autoMAT *out_eigenvectors_reim);
+/*
+	Decompresses each eigenvector row into two consecutive columns (real and imaginary part)
+*/
+
+/* End of file MAT_numerics.h */

dwsys/Makefile

 # makefile for library "dwsys".
-# David Weenink 20171019
+# David Weenink 20180828
 # Paul Boersma 2018-08-10
 
 include ../makefile.defs
@@ -12,6 +12,7 @@ OBJECTS = Collection_extensions.o Command.o \
 	DoublyLinkedList.o Eigen.o \
 	FileInMemory.o FileInMemorySet.o FileInMemoryManager.o\
 	Graphics_extensions.o Index.o \
+	MAT_numerics.o \
 	NUM2.o NUMhuber.o NUMlapack.o NUMmachar.o \
 	NUMf2c.o NUMcblas.o NUMclapack.o NUMcomplex.o NUMfft_d.o NUMsort2.o \
 	NUMmathlib.o NUMstring.o \