Skip to content
Commits on Source (7)
/.pc/
/debian/files
/debian/praat.dirs
/debian/praat.install
/debian/praat.menu
/debian/praat*debhelper*
/debian/praat.substvars
/debian/praat/
/debian/debhelper-build-stamp
/debian/.debhelper/
......@@ -591,6 +591,48 @@ autoEEG EEG_extractChannels (EEG me, numvec channelNumbers) {
}
}
static void Sound_removeChannel (Sound me, integer channelNumber) {
try {
Melder_require (channelNumber >= 1 && channelNumber <= my ny,
U"No channel ", channelNumber, U".");
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], my z [ichan], 1, my nx);
}
my ymax -= 1.0;
my ny -= 1;
} catch (MelderError) {
Melder_throw (me, U": channel ", channelNumber, U" not removed.");
}
}
void EEG_removeChannel (EEG me, integer channelNumber) {
try {
if (channelNumber < 1 || channelNumber > my numberOfChannels)
Melder_throw (U"No channel ", channelNumber, U".");
Melder_free (my channelNames [channelNumber]);
for (integer ichan = channelNumber; ichan < my numberOfChannels; ichan ++) {
my channelNames [ichan] = my channelNames [ichan + 1];
}
my numberOfChannels -= 1;
Sound_removeChannel (my sound.get(), channelNumber);
} catch (MelderError) {
Melder_throw (me, U": channel ", channelNumber, U" not removed.");
}
}
void EEG_removeChannel (EEG me, const char32 *channelName) {
try {
integer channelNumber = EEG_getChannelNumber (me, channelName);
if (channelNumber == 0)
Melder_throw (U"No channel named \"", channelName, U"\".");
EEG_removeChannel (me, channelNumber);
} catch (MelderError) {
Melder_throw (me, U": channel ", channelName, U" not removed.");
}
}
autoEEG EEGs_concatenate (OrderedOf<structEEG>* me) {
try {
if (my size < 1)
......@@ -678,7 +720,7 @@ autoMixingMatrix EEG_to_MixingMatrix (EEG me,
autoEEG EEG_MixingMatrix_to_EEG_unmix (EEG me, MixingMatrix you) {
Melder_require (my numberOfChannels == your numberOfRows,
U"To be able to unmix, the number of channels in ", me, U" (", my numberOfChannels, U")",
U" should be equal to the number of rows in ", you, U"(", your numberOfRows, U").");
U" should be equal to the number of rows in ", you, U" (", your numberOfRows, U").");
for (integer ichan = 1; ichan <= your numberOfRows; ichan ++) {
Melder_require (Melder_equ (my channelNames [ichan], your rowLabels [ichan]),
U"To be able to unmix, the name of channel ", ichan,
......@@ -700,7 +742,7 @@ autoEEG EEG_MixingMatrix_to_EEG_unmix (EEG me, MixingMatrix you) {
autoEEG EEG_MixingMatrix_to_EEG_mix (EEG me, MixingMatrix you) {
Melder_require (my numberOfChannels == your numberOfColumns,
U"To be able to mix, the number of channels in ", me, U" (", my numberOfChannels, U")",
U" should be equal to the number of columns in ", you, U"(", your numberOfColumns, U").");
U" should be equal to the number of columns in ", you, U" (", your numberOfColumns, U").");
for (integer ichan = 1; ichan <= your numberOfColumns; ichan ++) {
Melder_require (Melder_equ (my channelNames [ichan], your columnLabels [ichan]),
U"To be able to mix, the name of channel ", ichan,
......
......@@ -54,6 +54,8 @@ void EEG_removeTriggers (EEG me, kMelder_string which, const char32 *criterion);
autoEEG EEG_extractChannel (EEG me, integer channelNumber);
autoEEG EEG_extractChannel (EEG me, const char32 *channelName);
autoEEG EEG_extractChannels (EEG me, numvec channelNumbers);
void EEG_removeChannel (EEG me, integer channelNumber);
void EEG_removeChannel (EEG me, const char32 *channelName);
static inline autoSound EEG_extractSound (EEG me) { return Data_copy (my sound.get()); }
static inline autoTextGrid EEG_extractTextGrid (EEG me) { return Data_copy (my textgrid.get()); }
autoEEG EEG_extractPart (EEG me, double tmin, double tmax, bool preserveTimes);
......
/* ERPTier.cpp
*
* Copyright (C) 2011-2012,2014,2015,2016,2017 Paul Boersma
* Copyright (C) 2011-2012,2014-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
......@@ -344,7 +344,7 @@ autoERPTier ERPTier_extractEventsWhereColumn_string (ERPTier me, Table table,
for (integer ievent = 1; ievent <= my points.size; ievent ++) {
ERPPoint oldEvent = my points.at [ievent];
TableRow row = table -> rows.at [ievent];
if (Melder_stringMatchesCriterion (row -> cells [columnNumber]. string, which, criterion)) {
if (Melder_stringMatchesCriterion (row -> cells [columnNumber]. string, which, criterion, true)) {
autoERPPoint newEvent = Data_copy (oldEvent);
thy points. addItem_move (newEvent.move());
}
......
......@@ -202,7 +202,7 @@ void ERP_drawScalp (ERP me, Graphics graphics, double tmin, double tmax, double
double x = -1.0 + (icol - 1) * d;
if (x * x + y * y <= 1.0) {
double value = undefined;
real80 sum = 0.0, weight = 0.0;
longdouble sum = 0.0, weight = 0.0;
for (integer ichan = 1; ichan <= numberOfDrawableChannels; ichan ++) {
double dx = x - biosemiLocationData [ichan]. topX;
double dy = y - biosemiLocationData [ichan]. topY;
......@@ -216,7 +216,7 @@ void ERP_drawScalp (ERP me, Graphics graphics, double tmin, double tmax, double
weight += 1.0 / distance;
}
if (isundef (value))
value = ( sum == 0.0 ? 0.0 : real (sum / weight) );
value = ( sum == 0.0 ? 0.0 : double (sum / weight) );
image [irow] [icol] = value;
}
}
......@@ -306,7 +306,7 @@ void structERPWindow :: v_drawSelectionViewer () {
double x = -1.0 + (icol - 1) * d;
if (x * x + y * y <= 1.0) {
double value = undefined;
real80 sum = 0.0, weight = 0.0;
longdouble sum = 0.0, weight = 0.0;
for (integer ichan = 1; ichan <= numberOfDrawableChannels; ichan ++) {
double dx = x - biosemiLocationData [ichan]. topX;
double dy = y - biosemiLocationData [ichan]. topY;
......@@ -320,7 +320,7 @@ void structERPWindow :: v_drawSelectionViewer () {
weight += 1.0 / distance;
}
if (isundef (value))
value = ( sum == 0.0 ? 0.0 : real (sum / weight) );
value = ( sum == 0.0 ? 0.0 : double (sum / weight) );
image [irow] [icol] = value;
}
}
......
......@@ -90,7 +90,7 @@ ENTRY (U"See also")
LIST_ITEM (U"\\bu @@Independent Component Analysis on EEG@")
MAN_END
MAN_BEGIN (U"Independent Component Analysis on EEG", U"ppgb", 20180329)
MAN_BEGIN (U"Independent Component Analysis on EEG", U"ppgb", 20180502)
INTRO (U"Independent Component Analysis (ICA) is often used to improve @EEG signals. "
"See @@blind source separation@ for the algorithm.")
ENTRY (U"1. Selecting your channels")
......@@ -102,10 +102,10 @@ NORMAL (U"Once you have your reduced EEG, you can start to do ICA on it. "
"This starts by creating a @MixingMatrix: select your EEG object and choose ##To MixingMatrix...#. "
"The resulting MixingMatrix has one row for each of your 64 EEG channels, and columns called \"ic1\" through \"ic64\".")
ENTRY (U"3. How to see the independent components")
NORMAL (U"Select you EEG and your MixingMatrix together and choose ##Unmix#. "
NORMAL (U"Select your EEG and your MixingMatrix together and choose ##Unmix#. "
"The resulting ICA-EEG will have 64 channels called \"ic1\" through \"ic64\".")
ENTRY (U"4. Just checking back")
NORMAL (U"If you select your ICA-EEG together with your MixingMatrix choose ##Mix#, "
NORMAL (U"If you select your ICA-EEG together with your MixingMatrix and choose ##Mix#, "
"the resulting EEG should be very similar to your original 64-channel EEG.")
MAN_END
......
......@@ -153,6 +153,15 @@ DO
MODIFY_EACH_END
}
FORM (MODIFY_EEG_removeChannel, U"Remove channel", nullptr) {
SENTENCE (channel, U"Channel", U"Iz")
OK
DO
MODIFY_EACH (EEG)
EEG_removeChannel (me, channel);
MODIFY_EACH_END
}
FORM (MODIFY_EEG_subtractMeanChannel, U"Subtract mean channel", nullptr) {
LABEL (U"Range of reference channels:")
NATURAL (fromChannel, U"From channel", U"1")
......@@ -199,7 +208,7 @@ DO
}
FORM (NEW_EEG_extractChannels, U"EEG: Extract channels", nullptr) {
NUMVEC (channels, U"Channel numbers", U"to# (64)")
NUMVEC (channels, U"Channel numbers:", U"to# (64)")
OK
DO
CONVERT_EACH (EEG)
......@@ -721,10 +730,8 @@ DO
static autoDaata bdfFileRecognizer (integer nread, const char * /* header */, MelderFile file) {
const char32 *fileName = MelderFile_name (file);
bool isBdfFile = Melder_stringMatchesCriterion (fileName, kMelder_string::ENDS_WITH, U".bdf") ||
Melder_stringMatchesCriterion (fileName, kMelder_string::ENDS_WITH, U".BDF");
bool isEdfFile = Melder_stringMatchesCriterion (fileName, kMelder_string::ENDS_WITH, U".edf") ||
Melder_stringMatchesCriterion (fileName, kMelder_string::ENDS_WITH, U".EDF");
bool isBdfFile = Melder_stringMatchesCriterion (fileName, kMelder_string::ENDS_WITH, U".bdf", false);
bool isEdfFile = Melder_stringMatchesCriterion (fileName, kMelder_string::ENDS_WITH, U".edf", false);
if (nread < 512 || (! isBdfFile && ! isEdfFile)) return autoDaata ();
return EEG_readFromBdfFile (file);
}
......@@ -753,6 +760,7 @@ void praat_EEG_init () {
praat_addAction1 (classEEG, 0, U"Filter...", nullptr, 1, MODIFY_EEG_filter);
praat_addAction1 (classEEG, 0, U"Remove triggers...", nullptr, 1, MODIFY_EEG_removeTriggers);
praat_addAction1 (classEEG, 0, U"Set channel to zero...", nullptr, 1, MODIFY_EEG_setChannelToZero);
praat_addAction1 (classEEG, 0, U"Remove channel...", nullptr, 1, MODIFY_EEG_removeChannel);
praat_addAction1 (classEEG, 0, U"Analyse", nullptr, 0, nullptr);
praat_addAction1 (classEEG, 0, U"Extract channel...", nullptr, 0, NEW_EEG_extractChannel);
praat_addAction1 (classEEG, 0, U"Extract channels...", nullptr, 0, NEW_EEG_extractChannels);
......
......@@ -133,7 +133,7 @@ autoDTW Cepstrumc_to_DTW (Cepstrumc me, Cepstrumc thee, double wc, double wle, d
regression (me, i, ri.peek(), nr);
for (integer j = 1; j <= thy nx; j ++) {
Cepstrumc_Frame fj = & thy frame [j];
real80 d, dist = 0, distr = 0;
longdouble d, dist = 0.0, distr = 0.0;
if (wc != 0) { /* cepstral distance */
for (integer k = 1; k <= fj -> nCoefficients; k ++) {
d = fi -> c [k] - fj -> c [k];
......@@ -160,7 +160,7 @@ autoDTW Cepstrumc_to_DTW (Cepstrumc me, Cepstrumc thee, double wc, double wle, d
dist += wer * d * d;
}
dist /= wc + wle + wr + wer;
his z [i] [j] = sqrt ((real) dist); // prototype along y-direction
his z [i] [j] = sqrt ((double) dist); // prototype along y-direction
}
Melder_progress ( (double) i / my nx, U"Calculate distances: frame ",
i, U" from ", my nx, U".");
......
......@@ -817,15 +817,14 @@ DO
CONVERT_EACH_END (my name)
}
static void Sound_to_LPC_addWarning (UiForm dia) {
LABEL (U"Warning 1: for formant analysis, use \"To Formant\" instead.")
LABEL (U"Warning 2: if you do use \"To LPC\", you may want to resample first.")
LABEL (U"Click Help for more details.")
#define Sound_to_LPC_addWarning \
LABEL (U"Warning 1: for formant analysis, use \"To Formant\" instead.") \
LABEL (U"Warning 2: if you do use \"To LPC\", you may want to resample first.") \
LABEL (U"Click Help for more details.") \
LABEL (U"")
}
FORM (NEW_Sound_to_LPC_autocorrelation, U"Sound: To LPC (autocorrelation)", U"Sound: To LPC (autocorrelation)...") {
Sound_to_LPC_addWarning (dia);
Sound_to_LPC_addWarning
NATURAL (predictionOrder, U"Prediction order", U"16")
POSITIVE (windowLength, U"Window length (s)", U"0.025")
POSITIVE (timeStep, U"Time step (s)", U"0.005")
......@@ -839,7 +838,7 @@ DO
}
FORM (NEW_Sound_to_LPC_covariance, U"Sound: To LPC (covariance)", U"Sound: To LPC (covariance)...") {
Sound_to_LPC_addWarning (dia);
Sound_to_LPC_addWarning
NATURAL (predictionOrder, U"Prediction order", U"16")
POSITIVE (windowLength, U"Window length (s)", U"0.025")
POSITIVE (timeStep, U"Time step (s)", U"0.005")
......@@ -853,7 +852,7 @@ DO
}
FORM (NEW_Sound_to_LPC_burg, U"Sound: To LPC (burg)", U"Sound: To LPC (burg)...") {
Sound_to_LPC_addWarning (dia);
Sound_to_LPC_addWarning
NATURAL (predictionOrder, U"Prediction order", U"16")
POSITIVE (windowLength, U"Window length (s)", U"0.025")
POSITIVE (timeStep, U"Time step (s)", U"0.005")
......@@ -867,7 +866,7 @@ DO
}
FORM (NEW_Sound_to_LPC_marple, U"Sound: To LPC (marple)", U"Sound: To LPC (marple)...") {
Sound_to_LPC_addWarning (dia);
Sound_to_LPC_addWarning
NATURAL (predictionOrder, U"Prediction order", U"16")
POSITIVE (windowLength, U"Window length (s)", U"0.025")
POSITIVE (timeStep, U"Time step (s)", U"0.005")
......
......@@ -41,7 +41,7 @@ DO
FORM (WINDOW_Art_viewAndEdit, U"View & Edit Articulation", nullptr) {
static double muscles [1 + (int) kArt_muscle::MAX];
for (kArt_muscle muscle = (kArt_muscle) 1; muscle <= kArt_muscle::MAX; ++ muscle)
UiForm_addReal (dia, & muscles [(int) muscle], nullptr /* GUI-only */, kArt_muscle_getText (muscle), U"0.0");
UiForm_addReal (_dia_, & muscles [(int) muscle], nullptr /* GUI-only */, kArt_muscle_getText (muscle), U"0.0");
OK
FIND_ONE (Art)
for (int i = 1; i <= (int) kArt_muscle::MAX; i ++)
......
......@@ -490,7 +490,7 @@ void * KNN_classifyToTableOfRealAux
switch (input -> dist) {
case kOla_DISTANCE_WEIGHTED_VOTING:
for (integer y = input -> istart; y <= input -> istop; y ++) {
real80 sum = 0.0;
longdouble sum = 0.0;
for (integer c = 1; c <= ncategories; c ++) {
input -> output -> data [y] [c] *= 1.0 / OlaMAX (distances [c], kOla_MINFLOAT);
sum += input -> output -> data [y] [c];
......@@ -502,7 +502,7 @@ void * KNN_classifyToTableOfRealAux
break;
case kOla_SQUARED_DISTANCE_WEIGHTED_VOTING:
for (integer y = input -> istart; y <= input -> istop; y ++) {
real80 sum = 0.0;
longdouble sum = 0.0;
for (integer c = 1; c <= ncategories; c ++) {
input -> output -> data [y] [c] *= 1.0 / OlaMAX (OlaSQUARE (distances [c]), kOla_MINFLOAT);
sum += input -> output -> data [y] [c];
......@@ -514,7 +514,7 @@ void * KNN_classifyToTableOfRealAux
break;
case kOla_FLAT_VOTING:
for (integer y = input -> istart; y <= input -> istop; y ++) {
real80 sum = 0.0;
longdouble sum = 0.0;
for (integer c = 1; c <= ncategories; c ++) {
sum += input -> output -> data [y] [c];
}
......@@ -833,10 +833,10 @@ double KNN_distanceManhattan
)
{
real80 distance = 0.0;
longdouble distance = 0.0;
for (integer x = 1; x <= ps->nx; x ++)
distance += fabs (ps->z[rows][x] - pt->z[rowt][x]);
return (real) distance;
return (double) distance;
}
/////////////////////////////////////////////////////////////////////////////////////////////
......@@ -1414,7 +1414,7 @@ void KNN_normalizeFloatArray
{
integer c = 0;
real80 sum = 0.0;
longdouble sum = 0.0;
while (c < n)
sum += array [c ++]; // this sums over array [0 .. n-1]
......
/files
/praat.dirs
/praat.install
/praat.menu
/praat*debhelper*
/praat.substvars
/praat/
/debhelper-build-stamp
/.debhelper/
praat (6.0.40-1) unstable; urgency=medium
* New upstream version 6.0.40
* d/clean: Remove file test/fon/kanweg.nsp
* d/t/run-tests: Exclude unit test test_TextGrid_extensions.praat
-- Rafael Laboissiere <rafael@debian.org> Thu, 17 May 2018 03:55:36 -0300
praat (6.0.38-2) unstable; urgency=medium
* d/control: Use Salsa.d.o URLs in Vcs-* fields
......
dwtest/kanweg.FFNet
dwtest/kanweg.SpeechSynthesizer
test/fon/kanweg*.txt
test/fon/kanweg.nsp
test/kar/kanweg*.txt
test/kar/kanweg.TextGrid
test/sys/kanweg.txt
......
......@@ -79,7 +79,8 @@ for f in $(ls *.praat) ; do
-a $f != test_Sound_draw_where.praat \
-a $f != test_SpeechSynthesizer.praat \
-a $f != test_Sound_paint_where.praat \
-a $f != test_spectrogramTypes.praat ] ; then
-a $f != test_spectrogramTypes.praat \
-a $f != test_TextGrid_extensions.praat ] ; then
run $f
fi
done
......
......@@ -285,15 +285,15 @@ double Eigen_getSumOfEigenvalues (Eigen me, integer from, integer to) {
if (to > my numberOfEigenvalues || from > to) {
return undefined;
}
double sum = 0.0;
longdouble sum = 0.0;
for (integer i = from; i <= to; i++) {
sum += my eigenvalues[i];
}
return sum;
return (double) sum;
}
double Eigen_getCumulativeContributionOfComponents (Eigen me, integer from, integer to) {
double partial = 0.0, sum = 0.0;
longdouble partial = 0.0, sum = 0.0;
if (to == 0) {
to = my numberOfEigenvalues;
......@@ -318,7 +318,7 @@ integer Eigen_getDimensionOfFraction (Eigen me, double fraction) {
}
integer n = 1;
double p = my eigenvalues[1];
longdouble p = my eigenvalues [1];
while (p / sum < fraction && n < my numberOfEigenvalues) {
p += my eigenvalues [++ n];
}
......@@ -546,12 +546,12 @@ void Eigen_matrix_into_matrix_principalComponents (Eigen me, double **from, inte
for (integer irow = 1; irow <= numberOfRows; irow ++) {
for (integer icol = 1; icol <= numberOfDimensionsToKeep; icol ++) {
real80 r = 0.0;
longdouble r = 0.0;
for (integer k = 1; k <= my dimension; k ++) {
// eigenvector[icol] is in row[icol] of my eigenvectors
r += my eigenvectors [icol] [k] * from [irow] [from_colbegin + k - 1];
}
to [irow] [to_colbegin + icol - 1] = (real) r;
to [irow] [to_colbegin + icol - 1] = (double) r;
}
}
}
......
......@@ -116,7 +116,7 @@ autoFileInMemorySet FileInMemorySet_extractFiles (FileInMemorySet me, kMelder_st
autoFileInMemorySet thee = Thing_new (FileInMemorySet);
for (integer ifile = 1; ifile <= my size; ifile ++) {
FileInMemory fim = static_cast <FileInMemory> (my at [ifile]);
if (Melder_stringMatchesCriterion (fim -> d_path, which, criterion)) {
if (Melder_stringMatchesCriterion (fim -> d_path, which, criterion, true)) {
autoFileInMemory item = Data_copy (fim);
thy addItem_move (item.move());
}
......@@ -132,7 +132,7 @@ autoFileInMemorySet FileInMemorySet_listFiles (FileInMemorySet me, kMelder_strin
autoFileInMemorySet thee = Thing_new (FileInMemorySet);
for (integer ifile = 1; ifile <= my size; ifile ++) {
FileInMemory fim = static_cast<FileInMemory> (my at [ifile]);
if (Melder_stringMatchesCriterion (fim -> d_path, which, criterion)) {
if (Melder_stringMatchesCriterion (fim -> d_path, which, criterion, true)) {
thy addItem_ref (fim);
}
}
......@@ -209,7 +209,7 @@ integer FileInMemorySet_findNumberOfMatches_path (FileInMemorySet me, kMelder_st
integer numberOfMatches = 0;
for (integer ifile = 1; ifile <= my size; ifile ++) {
FileInMemory fim = static_cast <FileInMemory> (my at [ifile]);
if (Melder_stringMatchesCriterion (fim -> d_path, which, criterion)) {
if (Melder_stringMatchesCriterion (fim -> d_path, which, criterion, true)) {
numberOfMatches ++;
}
}
......@@ -222,7 +222,7 @@ bool FileInMemorySet_hasDirectory (FileInMemorySet me, const char32 *name) {
for (integer i = 1; i <= my size; i ++) {
FileInMemory fim = (FileInMemory) my at [i];
MelderString_append (& regex, U".*/", name, U"/.*");
if (Melder_stringMatchesCriterion (fim -> d_path, kMelder_string :: MATCH_REGEXP, regex.string)) {
if (Melder_stringMatchesCriterion (fim -> d_path, kMelder_string :: MATCH_REGEXP, regex.string, true)) {
match = true;
break;
}
......
......@@ -132,8 +132,8 @@ void NUMcentreRows_old (double **a, integer rb, integer re, integer cb, integer
}
void numvec_centre_inplace (numvec x, real *p_mean) {
real xmean;
void numvec_centre_inplace (numvec x, double *p_mean) {
double xmean;
sum_mean_scalar (x, nullptr, & xmean);
for (integer i = 1; i <= x.size; i ++) {
x [i] -= xmean;
......@@ -155,7 +155,7 @@ void NUMcentreColumns (double **a, integer rb, integer re, integer cb, integer c
for (integer i = rb; i <= re; i ++) {
colvec [i - rb + 1] = a [i] [j];
}
real colmean;
double colmean;
numvec_centre_inplace (colvec.get(), & colmean);
for (integer i = rb; i <= re; i ++) {
a [i] [j] = colvec [i - rb + 1];
......@@ -174,14 +174,14 @@ void NUMdoubleCentre (double **a, integer rb, integer re, integer cb, integer ce
void NUMnormalizeColumns (double **a, integer nr, integer nc, double norm) {
Melder_assert (norm > 0.0);
for (integer j = 1; j <= nc; j ++) {
real80 s = 0.0;
longdouble s = 0.0;
for (integer i = 1; i <= nr; i ++) {
s += a [i] [j] * a [i] [j];
}
if (s <= 0.0) {
continue;
}
s = sqrt (norm / (real) s);
s = sqrt (norm / (double) s);
for (integer i = 1; i <= nr; i ++) {
a [i] [j] *= s;
}
......@@ -191,14 +191,14 @@ void NUMnormalizeColumns (double **a, integer nr, integer nc, double norm) {
void NUMnormalizeRows (double **a, integer nr, integer nc, double norm) {
Melder_assert (norm > 0);
for (integer i = 1; i <= nr; i ++) {
real80 s = 0.0;
longdouble s = 0.0;
for (integer j = 1; j <= nc; j ++) {
s += a [i] [j] * a [i] [j];
}
if (s <= 0.0) {
continue;
}
s = sqrt (norm / (real) s);
s = sqrt (norm / (double) s);
for (integer j = 1; j <= nc; j ++) {
a [i] [j] *= s;
}
......@@ -207,7 +207,7 @@ void NUMnormalizeRows (double **a, integer nr, integer nc, double norm) {
void NUMnormalize (double **a, integer nr, integer nc, double norm) {
Melder_assert (norm > 0);
real80 sq = 0.0;
longdouble sq = 0.0;
for (integer i = 1; i <= nr; i ++) {
for (integer j = 1; j <= nc; j ++) {
sq += a [i] [j] * a [i] [j];
......@@ -216,7 +216,7 @@ void NUMnormalize (double **a, integer nr, integer nc, double norm) {
if (sq <= 0.0) {
return;
}
norm = sqrt (norm / (real) sq);
norm = sqrt (norm / (double) sq);
for (integer i = 1; i <= nr; i ++) {
for (integer j = 1; j <= nc; j ++) {
a [i] [j] *= norm;
......@@ -263,11 +263,11 @@ void NUMcovarianceFromColumnCentredMatrix (double **x, integer nrows, integer nc
for (integer i = 1; i <= ncols; i ++) {
for (integer j = i; j <= ncols; j ++) {
real80 sum = 0.0;
longdouble sum = 0.0;
for (integer k = 1; k <= nrows; k ++) {
sum += x [k] [i] * x [k] [j];
}
covar [i] [j] = covar [j] [i] = (real) sum / (nrows - ndf);
covar [i] [j] = covar [j] [i] = (double) sum / (nrows - ndf);
}
}
}
......@@ -367,57 +367,57 @@ void NUMmonotoneRegression (const double x [], integer n, double xs []) {
}
double NUMvector_getNorm1 (const double v [], integer n) {
real80 norm = 0.0;
longdouble norm = 0.0;
for (integer i = 1; i <= n; i ++) {
norm += fabs (v [i]);
}
return (real) norm;
return (double) norm;
}
double NUMvector_getNorm2 (const double v [], integer n) {
real80 norm = 0.0;
longdouble norm = 0.0;
for (integer i = 1; i <= n; i ++) {
norm += v [i] * v [i];
}
return sqrt ((real) norm);
return sqrt ((double) norm);
}
double NUMvector_normalize1 (double v [], integer n) {
real80 norm = NUMvector_getNorm1 (v, n);
longdouble norm = NUMvector_getNorm1 (v, n);
if (norm > 0.0) {
for (integer i = 1; i <= n; i ++) {
v [i] /= norm;
}
}
return (real) norm;
return (double) norm;
}
double NUMvector_normalize2 (double v [], integer n) {
real80 norm = NUMvector_getNorm2 (v, n);
longdouble norm = NUMvector_getNorm2 (v, n);
if (norm > 0) {
for (integer i = 1; i <= n; i ++) {
v [i] /= norm;
}
}
return (real) norm;
return (double) norm;
}
#undef TINY
void NUMcholeskySolve (double **a, integer n, double d [], double b [], double x []) {
for (integer i = 1; i <= n; i++) { /* Solve L.y=b */
real80 sum = b [i];
longdouble sum = b [i];
for (integer k = i - 1; k >= 1; k--) {
sum -= a [i] [k] * x [k];
}
x [i] = (real) sum / d [i];
x [i] = (double) sum / d [i];
}
for (integer i = n; i >= 1; i --) { /* Solve L^T.x=y */
real80 sum = x [i];
longdouble sum = x [i];
for (integer k = i + 1; k <= n; k ++) {
sum -= a [k] [i] * x [k];
}
x [i] = (real) sum / d [i];
x [i] = (double) sum / d [i];
}
}
......@@ -437,7 +437,7 @@ double NUMdeterminant_cholesky (double **a, integer n) {
// Determinant from diagonal, restore diagonal
real80 lnd = 0.0;
longdouble lnd = 0.0;
for (integer i = 1; i <= n; i ++) {
lnd += log (a [i] [i]);
a [i] [i] = d [i];
......@@ -451,7 +451,7 @@ double NUMdeterminant_cholesky (double **a, integer n) {
a [j] [i] = a [i] [j];
}
}
return (real) lnd;
return (double) lnd;
}
void NUMlowerCholeskyInverse (double **a, integer n, double *p_lnd) {
......@@ -485,18 +485,18 @@ double **NUMinverseFromLowerCholesky (double **m, integer n) {
autoNUMmatrix<double> r (1, n, 1, n);
for (integer i = 1; i <= n; i ++) {
for (integer j = 1; j <= i; j ++) {
real80 sum = 0.0;
longdouble sum = 0.0;
for (integer k = i; k <= n; k ++) {
sum += m [k] [i] * m [k] [j];
}
r [i] [j] = r [j] [i] = (real) sum;
r [i] [j] = r [j] [i] = (double) sum;
}
}
return r.transfer();
}
double NUMmahalanobisDistance_chi (double **linv, double *v, double *m, integer nr, integer n) {
real80 chisq = 0.0;
longdouble chisq = 0.0;
if (nr == 1) { // 1xn matrix
for (integer j = 1; j <= n; j ++) {
double t = linv [1] [j] * (v [j] - m [j]);
......@@ -511,25 +511,25 @@ double NUMmahalanobisDistance_chi (double **linv, double *v, double *m, integer
chisq += t * t;
}
}
return (real) chisq;
return (double) chisq;
}
double NUMtrace (double **a, integer n) {
real80 trace = 0.0;
longdouble trace = 0.0;
for (integer i = 1; i <= n; i ++) {
trace += a [i] [i];
}
return (real) trace;
return (double) trace;
}
double NUMtrace2 (double **a1, double **a2, integer n) {
real80 trace = 0.0;
longdouble trace = 0.0;
for (integer i = 1; i <= n; i ++) {
for (integer k = 1; k <= n; k ++) {
trace += a1 [i] [k] * a2 [k] [i];
}
}
return (real) trace;
return (double) trace;
}
void NUMeigensystem (double **a, integer n, double **evec, double eval []) {
......@@ -547,7 +547,7 @@ void NUMdominantEigenvector (double **mns, integer n, double *q, double *p_lambd
autoNUMvector<double> z (1, n);
double lambda0;
real80 lambda = 0.0;
longdouble lambda = 0.0;
for (integer k = 1; k <= n; k ++) {
for (integer l = 1; l <= n; l ++) {
lambda += q [k] * mns [k] [l] * q [l];
......@@ -574,7 +574,7 @@ void NUMdominantEigenvector (double **mns, integer n, double *q, double *p_lambd
q [k] = z [k] / znorm2;
}
lambda0 = (real) lambda;
lambda0 = (double) lambda;
lambda = 0.0;
for (integer k = 1; k <= n; k ++) {
......@@ -583,9 +583,9 @@ void NUMdominantEigenvector (double **mns, integer n, double *q, double *p_lambd
}
}
} while (fabs ((real) lambda - lambda0) > tolerance || ++ iter < 30);
} while (fabs ((double) lambda - lambda0) > tolerance || ++ iter < 30);
if (p_lambda) {
*p_lambda = (real) lambda;
*p_lambda = (double) lambda;
}
}
......@@ -594,11 +594,11 @@ void NUMprincipalComponents (double **a, integer n, integer nComponents, double
NUMeigensystem (a, n, evec.peek(), NULL);
for (integer i = 1; i <= n; i ++) {
for (integer j = 1; j <= nComponents; j ++) {
real80 s = 0.0;
longdouble s = 0.0;
for (integer k = 1; k <= n; k ++) {
s += a [k] [i] * evec [k] [j]; /* times sqrt(eigenvalue) ?? */
}
pc [i] [j] = (real) s;
pc [i] [j] = (double) s;
}
}
}
......@@ -620,11 +620,11 @@ void NUMdmatrix_projectRowsOnEigenspace (double **data, integer numberOfRows, in
for (integer irow = 1; irow <= numberOfRows; irow ++) {
for (integer icol = 1; icol <= numberOfEigenvectors; icol ++) {
real80 r = 0.0;
longdouble r = 0.0;
for (integer k = 1; k <= dimension; k ++) {
r += eigenvectors [icol] [k] * data [irow] [from_col + k - 1];
}
projection [irow] [to_col + icol - 1] = (real) r;
projection [irow] [to_col + icol - 1] = (double) r;
}
}
}
......@@ -644,11 +644,11 @@ void NUMdmatrix_projectColumnsOnEigenspace (double **data, integer numberOfColum
for (integer icol = 1; icol <= numberOfColumns; icol ++) {
for (integer irow = 1; irow <= numberOfEigenvectors; irow ++) {
real80 r = 0.0;
longdouble r = 0.0;
for (integer k = 1; k <= dimension; k ++) {
r += eigenvectors [irow] [k] * data [k] [icol];
}
projection [irow] [icol] = (real) r;
projection [irow] [icol] = (double) r;
}
}
}
......@@ -661,11 +661,11 @@ void NUMdmatrix_into_principalComponents (double **m, integer nrows, integer nco
autoSVD svd = SVD_create_d (mc.peek(), nrows, ncols);
for (integer i = 1; i <= nrows; i ++) {
for (integer j = 1; j <= numberOfComponents; j ++) {
real80 sum = 0.0;
longdouble sum = 0.0;
for (integer k = 1; k <= ncols; k ++) {
sum += svd -> v [k] [j] * m [i] [k];
}
pc [i] [j] = (real) sum;
pc [i] [j] = (double) sum;
}
}
}
......@@ -676,7 +676,7 @@ void NUMpseudoInverse (double **y, integer nr, integer nc, double **yinv, double
(void) SVD_zeroSmallSingularValues (me.get(), tolerance);
for (integer i = 1; i <= nc; i ++) {
for (integer j = 1; j <= nr; j ++) {
real80 s = 0.0;
longdouble s = 0.0;
for (integer k = 1; k <= nc; k ++) {
if (my d [k] != 0.0) {
s += my v [i] [k] * my u [j] [k] / my d [k];
......@@ -2979,7 +2979,7 @@ void NUMbiharmonic2DSplineInterpolation_getWeights (double *x, double *y, double
}
double NUMbiharmonic2DSplineInterpolation (double *x, double *y, integer n, double *w, double xp, double yp) {
real80 result = 0.0;
longdouble result = 0.0;
for (integer i = 1; i <= n; i ++) {
double dx = xp - x [i], dy = yp - y [i];
double d = dx * dx + dy * dy;
......
......@@ -34,7 +34,7 @@ int NUMstring_containsPrintableCharacter (const char32 *s);
void NUMstring_chopWhiteSpaceAtExtremes_inplace (char32 *string);
real * NUMstring_to_numbers (const char32 *s, integer *numbers_found);
double * NUMstring_to_numbers (const char32 *s, integer *numbers_found);
/* return array with the number of numbers found */
/*
......
......@@ -12,6 +12,7 @@
djmw 20020813 GPL header
djmw 20071201 Latest modification
pb 20100120 dlamc3_: declare volatile double ret_val to prevent optimization!
pb 2018 logicals -> bools
*/
......@@ -24,12 +25,12 @@
#define MAX(m,n) ((m) > (n) ? (m) : (n))
#define MIN(m,n) ((m) < (n) ? (m) : (n))
static int dlamc1_ (integer *beta, integer *t, integer *rnd, integer *ieee1);
static int dlamc2_ (integer *beta, integer *t, integer *rnd, double *eps, integer *emin, double *rmin, integer *emax,
static int dlamc1_ (integer *beta, integer *t, bool *rnd, bool *ieee1);
static int dlamc2_ (integer *beta, integer *t, bool *rnd, double *eps, integer *emin, double *rmin, integer *emax,
double *rmax);
static double dlamc3_ (double *, double *);
static int dlamc4_ (integer *emin, double *start, integer *base);
static int dlamc5_ (integer *beta, integer *p, integer *emin, integer *ieee, integer *emax, double *rmax);
static int dlamc5_ (integer *beta, integer *p, integer *emin, bool *ieee, integer *emax, double *rmax);
int NUMblas_daxpy (integer *n, double *da, double *dx, integer *incx, double *dy, integer *incy) {
/* System generated locals */
......@@ -698,7 +699,7 @@ int NUMblas_dgemv (const char *trans, integer *m, integer *n, double *alpha, dou
double NUMblas_dlamch (const char *cmach) {
/* Initialized data */
static integer first = TRUE;
static bool first = true;
/* System generated locals */
integer i__1;
......@@ -710,14 +711,14 @@ double NUMblas_dlamch (const char *cmach) {
static integer beta;
static double emin, prec, emax;
static integer imin, imax;
static integer lrnd;
static bool lrnd;
static double rmin, rmax, t, rmach;
static double smal, sfmin;
static integer it;
static double rnd, eps;
if (first) {
first = FALSE;
first = false;
dlamc2_ (&beta, &it, &lrnd, &eps, &imin, &rmin, &imax, &rmax);
base = (double) beta;
t = (double) it;
......@@ -770,7 +771,7 @@ double NUMblas_dlamch (const char *cmach) {
return ret_val;
} /* NUMblas_dlamch */
static int dlamc1_ (integer *beta, integer *t, integer *rnd, integer *ieee1) {
static int dlamc1_ (integer *beta, integer *t, bool *rnd, bool *ieee1) {
/* -- LAPACK auxiliary routine (version 3.0) -- Univ. of Tennessee, Univ.
of California Berkeley, NAG Ltd., Courant Institute, Argonne National
Lab, and Rice University October 31, 1992
......@@ -810,23 +811,23 @@ static int dlamc1_ (integer *beta, integer *t, integer *rnd, integer *ieee1) {
===================================================================== */
/* Initialized data */
static integer first = TRUE;
static bool first = true;
/* System generated locals */
double d__1, d__2;
/* Local variables */
static integer lrnd;
static bool lrnd;
static double a, b, c, f;
static integer lbeta;
static double savec;
static integer lieee1;
static bool lieee1;
static double t1, t2;
static integer lt;
static double one, qtr;
if (first) {
first = FALSE;
first = false;
one = 1.;
/* LBETA, LIEEE1, LT and LRND are the local values of BETA, IEEE1, T
......@@ -886,16 +887,16 @@ L20:
f = dlamc3_ (&d__1, &d__2);
c = dlamc3_ (&f, &a);
if (c == a) {
lrnd = TRUE;
lrnd = true;
} else {
lrnd = FALSE;
lrnd = false;
}
d__1 = b / 2;
d__2 = b / 100;
f = dlamc3_ (&d__1, &d__2);
c = dlamc3_ (&f, &a);
if (lrnd && c == a) {
lrnd = FALSE;
lrnd = false;
}
/* Try and decide whether rounding is done in the IEEE 'round to
......@@ -940,7 +941,7 @@ L30:
return 0;
} /* dlamc1_ */
static int dlamc2_ (integer *beta, integer *t, integer *rnd, double *eps, integer *emin, double *rmin, integer *emax,
static int dlamc2_ (integer *beta, integer *t, bool *rnd, double *eps, integer *emin, double *rmin, integer *emax,
double *rmax) {
/* -- LAPACK auxiliary routine (version 3.0) -- Univ. of Tennessee, Univ.
of California Berkeley, NAG Ltd., Courant Institute, Argonne National
......@@ -986,8 +987,8 @@ static int dlamc2_ (integer *beta, integer *t, integer *rnd, double *eps, intege
===================================================================== */
/* Table of constant values */
/* Initialized data */
static integer first = TRUE;
static integer iwarn = FALSE;
static bool first = true;
static bool iwarn = false;
/* System generated locals */
integer i__1;
......@@ -995,9 +996,9 @@ static int dlamc2_ (integer *beta, integer *t, integer *rnd, double *eps, intege
/* Builtin functions */
/* Local variables */
static integer ieee;
static bool ieee;
static double half;
static integer lrnd;
static bool lrnd;
static double leps, zero, a, b, c;
static integer i, lbeta;
static double rbase;
......@@ -1005,12 +1006,12 @@ static int dlamc2_ (integer *beta, integer *t, integer *rnd, double *eps, intege
static double smal;
static integer gpmin;
static double third, lrmin, lrmax, sixth;
static integer lieee1;
static bool lieee1;
static integer lt, ngnmin, ngpmin;
static double one, two;
if (first) {
first = FALSE;
first = false;
zero = 0.;
one = 1.;
two = 2.;
......@@ -1091,7 +1092,7 @@ L10:
dlamc4_ (&gpmin, &a, &lbeta);
d__1 = -a;
dlamc4_ (&gnmin, &d__1, &lbeta);
ieee = FALSE;
ieee = false;
if (ngpmin == ngnmin && gpmin == gnmin) {
if (ngpmin == gpmin) {
......@@ -1100,13 +1101,13 @@ L10:
e.g., VAX ) */
} else if (gpmin - ngpmin == 3) {
lemin = ngpmin - 1 + lt;
ieee = TRUE;
ieee = true;
/* ( Non twos-complement machines, with gradual underflow;
e.g., IEEE standard followers ) */
} else {
lemin = MIN (ngpmin, gpmin);
/* ( A guess; no known machine ) */
iwarn = TRUE;
iwarn = true;
}
} else if (ngpmin == gpmin && ngnmin == gnmin) {
......@@ -1117,7 +1118,7 @@ L10:
} else {
lemin = MIN (ngpmin, ngnmin);
/* ( A guess; no known machine ) */
iwarn = TRUE;
iwarn = true;
}
} else if ( (i__1 = ngpmin - ngnmin, labs (i__1)) == 1 && gpmin == gnmin) {
......@@ -1128,7 +1129,7 @@ L10:
} else {
lemin = MIN (ngpmin, ngnmin);
/* ( A guess; no known machine ) */
iwarn = TRUE;
iwarn = true;
}
} else {
......@@ -1136,11 +1137,11 @@ L10:
i__1 = MIN (ngpmin, ngnmin), i__1 = MIN (i__1, gpmin);
lemin = MIN (i__1, gnmin);
/* ( A guess; no known machine ) */
iwarn = TRUE;
iwarn = true;
}
/* Comment out this if block if EMIN is ok */
if (iwarn) {
first = TRUE;
first = true;
Melder_warning (U"\n\n WARNING. The value EMIN may be incorrect:- " "EMIN = ", lemin,
U"\nIf, after inspection, the value EMIN looks acceptable"
"please comment out \n the IF block as marked within the"
......@@ -1278,7 +1279,7 @@ L10:
return 0;
} /* dlamc4_ */
static int dlamc5_ (integer *beta, integer *p, integer *emin, integer *ieee, integer *emax, double *rmax) {
static int dlamc5_ (integer *beta, integer *p, integer *emin, bool *ieee, integer *emax, double *rmax) {
/*
First compute LEXP and UEXP, two powers of 2 that bound abs(EMIN). We
then assume that EMAX + abs(EMIN) will sum approximately to the bound
......