Commit 10f70c22 authored by Rafael Laboissière's avatar Rafael Laboissière
Browse files

New upstream version 6.1.03

parent f84e4506

EEG/manual_EEG.cpp

+1 −1
Original line number Diff line number Diff line
@@ -63,7 +63,7 @@ NORMAL (U"##Subtracting a reference signal.# The voltage on the whole scalp may 
	"and use ##Subtract mean channel...# to subtract this from each electrode channel.")

NORMAL (U"##Detrending.# With #Detrend, you subtract from each electrode channel a line in such a way that the first sample and the last sample become zero. "

	"Detrending and reference subtraction can be performed in either order.")

NORMAL (U"##Filtering.# With ##Filter...#, you band-pass filter each electrode channel. Filtering has to be done after detrending, but "

NORMAL (U"##Filtering.# With ##Filter...#, you band-pass filter each electrode channel. Filtering should be done after detrending, but "

	"filtering and reference subtraction can be performed in either order.")

ENTRY (U"4. How to do an ERP analysis")

NORMAL (U"An ERP is an Event-Related Potential. Events are marked somewhere in S1, S2, ... S8. In the above example, "

FFNet/FFNet.cpp

+1 −1
Original line number Diff line number Diff line
@@ -261,8 +261,8 @@ double FFNet_getBias (FFNet me, integer layer, integer unit) { 
		integer bias_unit = my wLast [node];

		return my w [bias_unit];

	} catch (MelderError) {

		return undefined;

		Melder_clearError ();

		return undefined;

	}

}

FFNet/FFNet_Eigen.cpp

+1 −1
Original line number Diff line number Diff line
@@ -26,7 +26,7 @@ 
#include "NUM2.h"



void FFNet_Eigen_drawIntersection (FFNet me, Eigen eigen, Graphics g, integer pcx, integer pcy, double xmin, double xmax, double ymin, double ymax) {

	integer ix = labs (pcx), iy = labs (pcy);

	integer ix = integer_abs (pcx), iy = integer_abs (pcy);

	integer numberOfEigenvalues = eigen -> numberOfEigenvalues;

	integer dimension = eigen -> dimension;

FFNet/manual_FFNet.cpp

+1 −1
Original line number Diff line number Diff line
@@ -329,7 +329,7 @@ DEFINITION (U"the number of different categories that you want the net to learn. 
TAG (U"##Number of units in hidden layer 1#, ##Number of units in hidden layer 2#") \

DEFINITION (U"determine the number of units in the hidden layers. " \

	"If you want a neural net with no hidden layers, both numbers have to be 0. "\

	"If you want a neural net with only 1 hidden layer then one of these numbers has to differ from 0. ")

	"If you want a neural net with only 1 hidden layer then one of these numbers should differ from 0. ")



MAN_BEGIN (U"Create FFNet...", U"djmw", 20040420)

INTRO (U"Create a new feedforward neural net of type @FFNet.")

dwsys/NUM2.cpp

+79 −149
Original line number Diff line number Diff line
@@ -34,7 +34,7 @@ 
 	invStudentQ, invChiSquareQ: modifications for 'undefined' return values.

 djmw 20030830 Corrected a bug in NUMtriangularfilter_amplitude

 djmw 20031111 Added NUMdmatrix_transpose, NUMdmatrix_printMatlabForm

 djmw 20040105 Added NUMmahalanobisDistance_chi

 djmw 20040105 Added NUMmahalanobisDistanceSquared_chi

 djmw 20040303 Added NUMstring_containsPrintableCharacter.

 djmw 20050406 NUMprocrutus->NUMprocrustes

 djmw 20060319 NUMinverse_cholesky: calculation of determinant is made optional
@@ -278,7 +278,7 @@ autoMAT newMATinverse_fromLowerCholeskyInverse (constMAT m) { 
	return result;

}



double NUMmahalanobisDistance (constMAT lowerInverse, constVEC v, constVEC m) {

double NUMmahalanobisDistanceSquared (constMAT lowerInverse, constVEC v, constVEC m) {

	Melder_assert (lowerInverse.ncol == v.size && v.size == m.size);

	longdouble chisq = 0.0;

	if (lowerInverse.nrow == 1) { // diagonal matrix is one row matrix
@@ -1193,66 +1193,28 @@ double NUMlnBeta (double a, double b) { 
	return status == GSL_SUCCESS ? result.val : undefined;

}



double NUMnormalityTest_HenzeZirkler (constMAT data, double *inout_beta, double *out_tnb, double *out_lnmu, double *out_lnvar) {

	const integer n = data.nrow, p = data.ncol;

	if (*inout_beta <= 0)

		*inout_beta = (1.0 / sqrt (2.0)) * pow ((1.0 + 2 * p) / 4.0, 1.0 / (p + 4.0)) * pow (n, 1.0 / (p + 4.0));

	const double p2 = p / 2.0;

	const double beta2 = *inout_beta * *inout_beta, beta4 = beta2 * beta2, beta8 = beta4 * beta4;

	const double gamma = 1.0 + 2.0 * beta2, gamma2 = gamma * gamma, gamma4 = gamma2 * gamma2;

	const double delta = 1.0 + beta2 * (4.0 + 3.0 * beta2), delta2 = delta * delta;

	double prob = undefined;



	double tnb = undefined, lnmu = undefined, lnvar = undefined;



	if (n < 2 || p < 1)

		return prob;



	autoVEC zero = newVECzero (p);

	autoMAT x = newMATcopy (data);



	MATcentreEachColumn_inplace (x.get()); // x - xmean



	autoMAT covar = MATcovarianceFromColumnCentredMatrix (x.get(), 0);



void MATscaledResiduals (MAT const& residuals, constMAT const& data, constMAT const& covariance, constVEC const& means) {

	try {

		MATlowerCholeskyInverse_inplace (covar.get(), nullptr);

		longdouble sumjk = 0.0, sumj = 0.0;

		const double b1 = beta2 / 2.0, b2 = b1 / (1.0 + beta2);

		/* Heinze & Wagner (1997), page 3

			We use d [j] [k] = ||Y [j]-Y [k]||^2 = (Y [j]-Y [k])'S^(-1)(Y [j]-Y [k])

			So d [j] [k]= d [k] [j] and d [j] [j] = 0

		*/

		for (integer j = 1; j <= n; j ++) {

			for (integer k = 1; k < j; k ++) {

				const double djk = NUMmahalanobisDistance (covar.get(), x.row (j), x.row(k));

				sumjk += 2.0 * exp (-b1 * djk); // factor 2 because d [j] [k] == d [k] [j]

		Melder_require (residuals.nrow == data.nrow && residuals.ncol == data.ncol,

			U"The data and the residuals should have the same dimensions.");

		Melder_require (covariance.ncol == means.size && data.ncol == means.size,

			U"The dimensions of the means and the covariance have to conform with the data.");

		autoVEC dif = newVECraw (data.ncol);

		autoMAT lowerInverse = newMATcopy (covariance);

		MATlowerCholeskyInverse_inplace (lowerInverse.get(), nullptr);

		for (integer irow = 1; irow <= data.nrow; irow ++) {

			dif <<= data.row (irow)  -  means;

			residuals.row(irow) <<= 0.0;

			if (lowerInverse.nrow == 1) { // diagonal matrix is one row matrix

				residuals.row(irow) <<= lowerInverse.row(1)  *  dif.get();

			} else {// square matrix

				for (integer icol = 1; icol <= data.ncol; icol ++)

					residuals [irow] [icol] = NUMinner (lowerInverse.row(icol).part (1, icol), dif.part (1, icol));

			}

			sumjk += 1.0; // for k == j

			const double djj = NUMmahalanobisDistance (covar.get(), x.row (j), zero.get());

			sumj += exp (-b2 * djj);

		}

		tnb = (1.0 / n) * (double) sumjk - 2.0 * pow (1.0 + beta2, - p2) * (double) sumj + n * pow (gamma, - p2); // n *

	} catch (MelderError) {

		Melder_clearError ();

		tnb = 4.0 * n;

	}



	double mu = 1.0 - pow (gamma, -p2) * (1.0 + p * beta2 / gamma + p * (p + 2) * beta4 / (2.0 * gamma2));

	double var = 2.0 * pow (1.0 + 4.0 * beta2, -p2)

		 + 2.0 * pow (gamma,  -p) * (1.0 + 2.0 * p * beta4 / gamma2  + 3.0 * p * (p + 2) * beta8 / (4.0 * gamma4))

		 - 4.0 * pow (delta, -p2) * (1.0 + 3.0 * p * beta4 / (2.0 * delta) + p * (p + 2) * beta8 / (2.0 * delta2));

	double mu2 = mu * mu;

	lnmu = log (sqrt (mu2 * mu2 / (mu2 + var)));

	lnvar = sqrt (log ((mu2 + var) / mu2));

	if (out_lnmu)

		*out_lnmu = lnmu;

	if (out_lnvar)

		*out_lnvar = lnvar;

	if (out_tnb)

		*out_tnb = tnb;

	prob = NUMlogNormalQ (tnb, lnmu, lnvar);

	return prob;

		Melder_throw (U"MATscaleResiduals: not performed.");

	}

}



/*************** Hz <--> other freq reps *********************/
@@ -1760,7 +1722,8 @@ bool NUMclipLineWithinRectangle (double xl1, double yl1, double xl2, double yl2, 
	}

	if (ncrossings == 0)

		return false;

	Melder_require (ncrossings <= 2, U"Too many crossings found.");

	Melder_require (ncrossings <= 2,

		U"Too many crossings found.");



	/*

		if start and endpoint of line are outside rectangle and ncrossings == 1, than the line only touches.
@@ -1788,7 +1751,8 @@ end: 
}



void NUMgetEllipseBoundingBox (double a, double b, double cospsi, double *out_width, double *out_height) {

	Melder_assert (cospsi >= -1.0 && cospsi <= 1.0);

	Melder_require (cospsi >= -1.0 && cospsi <= 1.0,

		U"NUMgetEllipseBoundingBox: abs (cospi) should not exceed 1.", cospsi);

	double width, height;

	if (cospsi == 1.0) { // a-axis along x-axis

		width = a;
@@ -1921,33 +1885,24 @@ double NUMminimize_brent (double (*f) (double x, void *closure), double a, doubl 
}



/*

	probs is probability vector, i.e. all 0 <= probs [i] <= 1 and sum(i=1;i=nprobs, probs [i])= 1

	p is a probability

	probs is probability vector, i.e. all 0 <= probs [i] <= 1 and sum(i=1...probs.size, probs [i])= 1

*/

integer NUMgetIndexFromProbability (constVEC probs, double p) { //TODO HMM zero start matrices

integer NUMgetIndexFromProbability (constVEC probs, double p) {

	integer index = 1;

	double psum = probs [index];

	while (p > psum && index < probs.size) {

	while (p > psum && index < probs.size)

		psum += probs [++ index];

	}

	return index;

}

integer NUMgetIndexFromProbability (double *probs, integer nprobs, double p) { //TODO HMM zero start matrices

	integer index = 1;

	double psum = probs [index];

	while (p > psum && index < nprobs) {

		psum += probs [++ index];

	}

	return index;

}



// straight line fitting



void NUMlineFit_theil (constVEC x, constVEC y, double *out_m, double *out_intercept, bool incompleteMethod)

{

void NUMlineFit_theil (constVEC x, constVEC y, double *out_m, double *out_intercept, bool incompleteMethod) {

	try {

		Melder_assert (x.size == y.size);

		/* Theil's incomplete method:

		Melder_require (x.size == y.size,

			U"NUMlineFit_theil: the sizes of the two vectors should be equal.");

		/*

			Theil's incomplete method:

			Split (x [i],y [i]) as

			(x [i],y [i]), (x [N+i],y [N=i], i=1..numberOfPoints/2

			m [i] = (y [N+i]-y [i])/(x [N+i]-x [i])
@@ -1985,15 +1940,18 @@ void NUMlineFit_theil (constVEC x, constVEC y, double *out_m, double *out_interc 
			VECsort_inplace (mbs.part (1, x.size));

			intercept = NUMquantile (mbs.part (1, x.size), 0.5);

		}

		if (out_m) *out_m = m;

		if (out_intercept) *out_intercept = intercept;

		if (out_m)

			*out_m = m;

		if (out_intercept)

			*out_intercept = intercept;

	} catch (MelderError) {

		Melder_throw (U"No line fit (Theil's method).");

	}

}



void NUMlineFit_LS (constVEC x, constVEC y, double *out_m, double *out_intercept) {

	Melder_assert (x.size == y.size);

	Melder_require (x.size == y.size,

		U"NUMlineFit_LS: the sizes of the two vectors should be equal.");

	double sx = NUMsum (x);

	double sy = NUMsum (y);

	double xmean = sx / x.size;
@@ -2005,19 +1963,20 @@ void NUMlineFit_LS (constVEC x, constVEC y, double *out_m, double *out_intercept 
	}

	// y = m*x + b

	m /= st2;

	if (out_intercept) *out_intercept = (sy - m * sx) / x.size;

	if (out_m) *out_m = m;

	if (out_intercept)

		*out_intercept = (sy - m * sx) / x.size;

	if (out_m)

		*out_m = m;

}



void NUMlineFit (constVEC x, constVEC y, double *out_m, double *out_intercept, int method) {

	if (method == 1) {

	if (method == 1)

		NUMlineFit_LS (x, y, out_m, out_intercept);

	} else if (method == 3) {

	else if (method == 3)

		NUMlineFit_theil (x, y, out_m, out_intercept, false);

	} else {

	else

		NUMlineFit_theil (x, y, out_m, out_intercept, true);

}

}



// IEEE: Programs for digital signal processing section 4.3 LPTRN

// lpc [1..n] to rc [1..n]
@@ -2031,11 +1990,10 @@ void VECrc_from_lpc (VEC rc, constVEC lpc) { 
		rc [m] = a [m];

		Melder_require (fabs (rc [m]) <= 1.0, U"Relection coefficient [", m, U"] larger than 1.");

		b.part (1, m) <<= a.part (1, m);

		for (integer i = 1; i < m; i ++) {

		for (integer i = 1; i < m; i ++)

			a [i] = (b [i] - rc [m] * b [m - i]) / (1.0 - rc [m] * rc [m]);

	}

}

}



void VEClpc_from_rc (VEC lpc, constVEC rc) {

	Melder_assert (lpc.size == rc.size);
@@ -2226,7 +2184,8 @@ integer NUMrandomBinomial (double p, integer n) { 
	integer ix;			// return value

	bool flipped = false;



	if (n == 0) return 0;

	if (n == 0)

		return 0;

	if (p > 0.5) {

		p = 1.0 - p;	// work with small p

		flipped = true;
@@ -2257,9 +2216,8 @@ integer NUMrandomBinomial (double p, integer n) { 
			double u = NUMrandomUniform (0.0, 1.0); // djmw gsl_rng_uniform (rng);



			for (ix = 0; ix <= BINV_CUTOFF; ++ ix) {

				if (u < f) {

				if (u < f)

					goto Finish;

				}

				u -= f;

				// Use recursion f(x+1) = f(x)* [(n-x)/(x+1)]* [p/(1-p)]

				f *= s * (n - ix) / (ix + 1.0);
@@ -2352,23 +2310,20 @@ TryAgain: 
			// Parallelogram region

			double x = xl + (u - p1) / c;

			v = v * c + 1.0 - fabs (x - xm) / p1;

			if (v > 1.0 || v <= 0.0) {

			if (v > 1.0 || v <= 0.0)

				goto TryAgain;

			}

			ix = (integer) x;

		} else if (u <= p3) {

			// Left tail

			ix = (integer) (xl + log (v) / lambda_l);

			if (ix < 0) {

			if (ix < 0)

				goto TryAgain;

			}

			v *= ((u - p2) * lambda_l);

		} else {

			// Right tail

			ix = (integer) (xr - log (v) / lambda_r);

			if (ix > (double) n) {

			if (ix > (double) n)

				goto TryAgain;

			}

			v *= ((u - p3) * lambda_r);

		}
@@ -2393,7 +2348,7 @@ TryAgain: 
			More efficient determination of whether v < f(x)/f(M)

		 */



		integer k = labs (ix - m);

		integer k = integer_abs (ix - m);



		if (k <= FAR_FROM_MEAN) {

			/*
@@ -2406,14 +2361,12 @@ TryAgain: 
			var = v;



			if (m < ix) {

				for (integer i = m + 1; i <= ix; i ++) {

				for (integer i = m + 1; i <= ix; i ++)

					f *= (g / i - s);

				}

			} else if (m > ix) {

				for (integer i = ix + 1; i <= m; i ++) {

				for (integer i = ix + 1; i <= m; i ++)

					f /= (g / i - s);

			}

			}



			accept = f;

		} else {
@@ -2429,13 +2382,11 @@ TryAgain: 
				*/

				double amaxp = k / npq * ((k * (k / 3.0 + 0.625) + (1.0 / 6.0)) / npq + 0.5);

				double ynorm = -(k * k / (2.0 * npq));

				if (var < ynorm - amaxp) {

				if (var < ynorm - amaxp)

					goto Finish;

				}

				if (var > ynorm + amaxp) {

				if (var > ynorm + amaxp)

					goto TryAgain;

			}

			}



			/* 

				Now, again: do the test log(v) vs. log f(x)/f(M)
@@ -2486,12 +2437,11 @@ TryAgain: 
		}





		if (var <= accept) {

		if (var <= accept)

			goto Finish;

		} else {

		else

			goto TryAgain;

	}

	}



Finish:
@@ -2499,12 +2449,11 @@ Finish: 
}



double NUMrandomBinomial_real (double p, integer n) {

	if (p < 0.0 || p > 1.0 || n < 0) {

	if (p < 0.0 || p > 1.0 || n < 0)

		return undefined;

	} else {

	else

		return (double) NUMrandomBinomial (p, n);

}

}



double NUMrandomGamma (const double alpha, const double beta) {

	Melder_require (alpha > 0 && beta > 0,
@@ -2591,8 +2540,7 @@ void NUMfixIndicesInRange (integer lowerLimit, integer upperLimit, integer *lowI 
	}

}



void NUMgetEntropies (constMAT m, double *out_h, double *out_hx, 

	double *out_hy,	double *out_hygx, double *out_hxgy, double *out_uygx, double *out_uxgy, double *out_uxy) {

void NUMgetEntropies (constMATVU const& m, double *out_h, double *out_hx, double *out_hy, double *out_hygx, double *out_hxgy, double *out_uygx, double *out_uxgy, double *out_uxy) {

	

	double h = undefined, hx = undefined, hy = undefined;

	double hxgy = undefined, hygx = undefined, uygx = undefined, uxgy = undefined, uxy = undefined;
@@ -2621,8 +2569,7 @@ void NUMgetEntropies (constMAT m, double *out_h, double *out_hx, 
		

		longdouble hx_t = 0.0;

		for (integer j = 1; j <= m.ncol; j ++) {

			longdouble colsum = 0.0;

			for (integer i = 1; i <= m.nrow; i++) colsum += m [i] [j];

			double colsum = NUMsum (m.column (j));

			if (colsum > 0.0) {

				longdouble p = colsum / totalSum;

				hx_t -= p * NUMlog2 (p);
@@ -2667,27 +2614,6 @@ void NUMgetEntropies (constMAT m, double *out_h, double *out_hx, 
}

#undef TINY



double NUMfrobeniusnorm (constMAT x) {

	longdouble scale = 0.0;

	longdouble ssq = 1.0;

	for (integer i = 1; i <= x.nrow; i ++) {

		for (integer j = 1; j <= x.ncol; j ++) {

			if (x [i] [j] != 0.0) {

				longdouble absxi = fabs (x [i] [j]);

				if (scale < absxi) {

					longdouble t = scale / absxi;

					ssq = 1.0 + ssq * t * t;

					scale = absxi;

				} else {

					longdouble t = absxi / scale;

					ssq += t * t;

				}

			}

		}

	}

	return scale * sqrt ((double) ssq);

}



double NUMtrace (const constMATVU& a) {

	Melder_assert (a.nrow == a.ncol);

	longdouble trace = 0.0;
@@ -2762,15 +2688,15 @@ void NUMeigencmp22 (double a, double b, double c, double *out_rt1, double *out_r 
	longdouble acs = fabs (cs), cs1, sn1;

	if (acs > ab) {

		longdouble ct = -tb / cs;

		sn1 = 1 / sqrt (1 + ct * ct);

		sn1 = 1.0 / sqrt (1.0 + ct * ct);

		cs1 = ct * sn1;

	} else {

		if (ab == 0) {

			cs1 = 1;

			sn1 = 0;

			cs1 = 1.0;

			sn1 = 0.0;

		} else {

			tn = -cs / tb;

			cs1 = 1 / sqrt (1 + tn * tn);

			cs1 = 1.0 / sqrt (1.0 + tn * tn);

			sn1 = tn * cs1;

		}

	}
@@ -2779,6 +2705,10 @@ void NUMeigencmp22 (double a, double b, double c, double *out_rt1, double *out_r 
		cs1 = -sn1;

		sn1 = tn;

	}

	if (fabs (cs1) > 1.0)

		cs1 = copysign (1.0, cs1);

	if (fabs (sn1) > 1.0)

		sn1 = copysign (1.0, sn1);

	if (out_rt1)

		*out_rt1 = (double) rt1;

	if (out_rt2)