/* NUM.c
 *
 * Copyright (C) 1992-2007 Paul Boersma
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 *
 * This program 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 program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

/*
 * pb 2002/03/07 GPL
 * pb 2003/06/19 ridders3 replaced with ridders
 * pb 2003/07/09 gsl
 * pb 2003/08/27 NUMfisherQ: underflow and iteration excess should not return NUMundefined
 * pb 2005/07/08 NUMpow
 * pb 2006/08/02 NUMinvSigmoid
 * pb 2007/01/27 use #defines for value interpolation
 */

#include "NUM.h"
#include "NUM2.h"
#include <stdlib.h>
#include "melder.h"
#define SIGN(x,s) ((s) < 0 ? -fabs (x) : fabs(x))
#define my  me ->
double NUMpow (double base, double exponent) { return base <= 0.0 ? 0.0 : pow (base, exponent); }
/*
	GSL is more accurate than the other routines, but makes
	"Sound: To Intensity..." 10 times slower...
*/
#include "gsl_errno.h"
#include "gsl_sf_bessel.h"
#include "gsl_sf_gamma.h"
#include "gsl_sf_erf.h"

void NUMfbtoa (double formant, double bandwidth, double dt, double *a1, double *a2) {
	*a1 = 2 * exp (- NUMpi * bandwidth * dt) * cos (2 * NUMpi * formant * dt);
	*a2 = exp (- 2 * NUMpi * bandwidth * dt);
}

void NUMfilterSecondOrderSection_a (float x [], long n, double a1, double a2) {
	long i;
	x [2] += a1 * x [1];
	for (i = 3; i <= n; i ++)
		x [i] += a1 * x [i - 1] - a2 * x [i - 2];
}

void NUMfilterSecondOrderSection_fb (float x [], long n, double dt, double formant, double bandwidth) {
	double a1, a2;
	NUMfbtoa (formant, bandwidth, dt, & a1, & a2);
	NUMfilterSecondOrderSection_a (x, n, a1, a2);
}

double NUMftopreemphasis (double f, double dt) {
	return exp (- 2.0 * NUMpi * f * dt);
}

void NUMpreemphasize_a (float x [], long n, double preemphasis) {
	long i;
	for (i = n; i >= 2; i --)
		x [i] -= preemphasis * x [i - 1];
}

void NUMdeemphasize_a (float x [], long n, double preemphasis) {
	long i;
	for (i = 2; i <= n; i ++)
		x [i] += preemphasis * x [i - 1];
}

void NUMpreemphasize_f (float x [], long n, double dt, double frequency) {
	NUMpreemphasize_a (x, n, NUMftopreemphasis (frequency, dt));
}

void NUMdeemphasize_f (float x [], long n, double dt, double frequency) {
	NUMdeemphasize_a (x, n, NUMftopreemphasis (frequency, dt));
}

void NUMautoscale (float x [], long n, double scale) {
	long i;
	double maximum = 0.0;
	for (i = 1; i <= n; i ++)
		if (fabs (x [i]) > maximum) maximum = fabs (x [i]);
	if (maximum > 0.0) {
		double factor = scale / maximum;
		for (i = 1; i <= n; i ++)
			x [i] *= factor;
	}
}

double NUMlnGamma (double x) {
	gsl_sf_result result;
	int status = gsl_sf_lngamma_e (x, & result);
	return status == GSL_SUCCESS ? result. val : NUMundefined;
}

double NUMbeta (double z, double w) {
	if (z <= 0.0 || w <= 0.0) return NUMundefined;
	return exp (NUMlnGamma (z) + NUMlnGamma (w) - NUMlnGamma (z + w));
}

double NUMincompleteBeta (double a, double b, double x) {
	gsl_sf_result result;
	int status = gsl_sf_beta_inc_e (a, b, x, & result);
	if (status != GSL_SUCCESS && status != GSL_EUNDRFLW && status != GSL_EMAXITER) {
		Melder_fatal ("NUMincompleteBeta status %d", status);
		return NUMundefined;
	}
	return result. val;
}

double NUMbinomialP (double p, double k, double n) {
	double binomialQ;
	if (p < 0.0 || p > 1.0 || n <= 0.0 || k < 0.0 || k > n) return NUMundefined;
	if (k == n) return 1.0;
	binomialQ = NUMincompleteBeta (k + 1, n - k, p);
	if (binomialQ == NUMundefined) return NUMundefined;
	return 1.0 - binomialQ;
}

double NUMbinomialQ (double p, double k, double n) {
	if (p < 0.0 || p > 1.0 || n <= 0.0 || k < 0.0 || k > n) return NUMundefined;
	if (k == 0.0) return 1.0;
	return NUMincompleteBeta (k, n - k + 1, p);
}

struct binomial { double p, k, n; };

static double binomialP (double p, void *binomial_void) {
	struct binomial *binomial = (struct binomial *) binomial_void;
	return NUMbinomialP (p, binomial -> k, binomial -> n) - binomial -> p;
}

static double binomialQ (double p, void *binomial_void) {
	struct binomial *binomial = (struct binomial *) binomial_void;
	return NUMbinomialQ (p, binomial -> k, binomial -> n) - binomial -> p;
}

double NUMinvBinomialP (double p, double k, double n) {
	static struct binomial binomial;
	if (p < 0 || p > 1 || n <= 0 || k < 0 || k > n) return NUMundefined;
	if (k == n) return 1.0;
	binomial. p = p;
	binomial. k = k;
	binomial. n = n;
	return NUMridders (binomialP, 0.0, 1.0, & binomial);
}

double NUMinvBinomialQ (double p, double k, double n) {
	static struct binomial binomial;
	if (p < 0 || p > 1 || n <= 0 || k < 0 || k > n) return NUMundefined;
	if (k == 0) return 0.0;
	binomial. p = p;
	binomial. k = k;
	binomial. n = n;
	return NUMridders (binomialQ, 0.0, 1.0, & binomial);
}

/* Modified Bessel function I0. Abramowicz & Stegun, p. 378.*/
double NUMbessel_i0_f (double x) {
	double t;
	if (x < 0.0) return NUMbessel_i0_f (- x);
	if (x < 3.75) {
		/* Formula 9.8.1. Accuracy 1.6e-7. */
		t = x / 3.75;
		t *= t;
		return 1.0 + t * (3.5156229 + t * (3.0899424 + t * (1.2067492
			+ t * (0.2659732 + t * (0.0360768 + t * 0.0045813)))));
	}
	/*
		otherwise: x >= 3.75
	*/
	/* Formula 9.8.2. Accuracy of the polynomial factor 1.9e-7. */
	t = 3.75 / x;   /* <= 1.0 */
	return exp (x) / sqrt (x) * (0.39894228 + t * (0.01328592
		+ t * (0.00225319 + t * (-0.00157565 + t * (0.00916281
		+ t * (-0.02057706 + t * (0.02635537 + t * (-0.01647633
		+ t * 0.00392377))))))));
}

/* Modified Bessel function I1. Abramowicz & Stegun, p. 378. */
double NUMbessel_i1_f (double x) {
	double t;
	if (x < 0.0) return - NUMbessel_i1_f (- x);
	if (x < 3.75) {
		/* Formula 9.8.3. Accuracy of the polynomial factor 8e-9. */
		t = x / 3.75;
		t *= t;
		return x * (0.5 + t * (0.87890594 + t * (0.51498869 + t * (0.15084934
			+ t * (0.02658733 + t * (0.00301532 + t * 0.00032411))))));
	}
	/*
		otherwise: x >= 3.75
	*/
	/* Formula 9.8.4. Accuracy of the polynomial factor 2.2e-7. */
	t = 3.75 / x;   /* <= 1.0 */
	return exp (x) / sqrt (x) * (0.39894228 + t * (-0.03988024
		+ t * (-0.00362018 + t * (0.00163801 + t * (-0.01031555
		+ t * (0.02282967 + t * (-0.02895312 + t * (0.01787654
		+ t * (-0.00420059)))))))));
}

double NUMbesselI (long n, double x) {
	gsl_sf_result result;
	int status = gsl_sf_bessel_In_e (n, x, & result);
	return status == GSL_SUCCESS ? result. val : NUMundefined;
}

/* Modified Bessel function K0. Abramowicz & Stegun, p. 379. */
double NUMbessel_k0_f (double x) {
	double t;
	if (x <= 0.0) return NUMundefined;   /* Positive infinity. */
	if (x <= 2.0) {
		/* Formula 9.8.5. Accuracy 1e-8. */
		double x2 = 0.5 * x, t = x2 * x2;
		return - log (x2) * NUMbessel_i0_f (x) + (-0.57721566 + t * (0.42278420
			+ t * (0.23069756 + t * (0.03488590 + t * (0.00262698
			+ t * (0.00010750 + t * 0.00000740))))));
	}
	/*
		otherwise: 2 < x < positive infinity
	*/
	/* Formula 9.8.6. Accuracy of the polynomial factor 1.9e-7. */
	t = 2.0 / x;   /* < 1.0 */
	return exp (- x) / sqrt (x) * (1.25331414 + t * (-0.07832358
		+ t * (0.02189568 + t * (-0.01062446 + t * (0.00587872
		+ t * (-0.00251540 + t * 0.00053208))))));
}

/* Modified Bessel function K1. Abramowicz & Stegun, p. 379. */
double NUMbessel_k1_f (double x) {
	double t;
	if (x <= 0.0) return NUMundefined;   /* Positive infinity. */
	if (x <= 2.0) {
		/* Formula 9.8.7. Accuracy  of the polynomial factor 8e-9. */
		double x2 = 0.5 * x, t = x2 * x2;
		return log (x2) * NUMbessel_i1_f (x) + (1.0 / x) * (1.0 + t * (0.15443144
			+ t * (-0.67278579 + t * (-0.18156897 + t * (-0.01919402
			+ t * (-0.00110404 + t * (-0.00004686)))))));
	}
	/*
		otherwise: 2 < x < positive infinity
	*/
	/* Formula 9.8.8. Accuracy of the polynomial factor 2.2e-7. */
	t = 2.0 / x;   /* < 1.0 */
	return exp (- x) / sqrt (x) * (1.25331414 + t * (0.23498619
			 + t * (-0.03655620 + t * (0.01504268 + t * (-0.00780353
			 + t * (0.00325614 + t * (-0.00068245)))))));
}

double NUMbesselK_f (long n, double x) {
	double twoByX, besselK_min2, besselK_min1, besselK = NUMundefined;
	long i;
	Melder_assert (n >= 0 && x > 0);
	besselK_min2 = NUMbessel_k0_f (x);
	if (n == 0) return besselK_min2;
	besselK_min1 = NUMbessel_k1_f (x);
	if (n == 1) return besselK_min1;
	Melder_assert (n >= 2);
	twoByX = 2.0 / x;
	/*
		Recursion formula.
	*/
	for (i = 1; i < n; i ++) {
		besselK = besselK_min2 + twoByX * i * besselK_min1;
		besselK_min2 = besselK_min1;
		besselK_min1 = besselK;
	}
	Melder_assert (NUMdefined (besselK));
	return besselK;
}

double NUMbesselK (long n, double x) {
	gsl_sf_result result;
	int status = gsl_sf_bessel_Kn_e (n, x, & result);
	return status == GSL_SUCCESS ? result. val : NUMundefined;
}

double NUMsigmoid (double x)
	{ return x > 0.0 ? 1 / (1 + exp (- x)) : 1 - 1 / (1 + exp (x)); }

double NUMinvSigmoid (double x)
	{ return x <= 0.0 || x >= 1.0 ? NUMundefined : log (x / (1.0 - x)); }

double NUMerfcc (double x) {
	gsl_sf_result result;
	int status = gsl_sf_erfc_e (x, & result);
	return status == GSL_SUCCESS ? result. val : NUMundefined;
}

double NUMgaussP (double z) {
	return 1 - 0.5 * NUMerfcc (NUMsqrt1_2 * z);
}

double NUMgaussQ (double z) {
	return 0.5 * NUMerfcc (NUMsqrt1_2 * z);
}

double NUMincompleteGammaP (double a, double x) {
	gsl_sf_result result;
	int status = gsl_sf_gamma_inc_P_e (a, x, & result);
	return status == GSL_SUCCESS ? result. val : NUMundefined;
}

double NUMincompleteGammaQ (double a, double x) {
	gsl_sf_result result;
	int status = gsl_sf_gamma_inc_Q_e (a, x, & result);
	return status == GSL_SUCCESS ? result. val : NUMundefined;
}

double NUMchiSquareP (double chiSquare, long degreesOfFreedom) {
	if (chiSquare < 0 || degreesOfFreedom <= 0) return NUMundefined;
	return NUMincompleteGammaP (0.5 * degreesOfFreedom, 0.5 * chiSquare);
}

double NUMchiSquareQ (double chiSquare, long degreesOfFreedom) {
	if (chiSquare < 0 || degreesOfFreedom <= 0) return NUMundefined;
	return NUMincompleteGammaQ (0.5 * degreesOfFreedom, 0.5 * chiSquare);
}

double NUMcombinations (long n, long k) {
	double result = 1.0;
	long i;
	if (k > n / 2) k = n - k;
	for (i = 1; i <= k; i ++) result *= n - i + 1;
	for (i = 2; i <= k; i ++) result /= i;
	return result;
}

#define NUM_interpolate_simple_cases \
	if (nx < 1) return NUMundefined; \
	if (x > nx) return y [nx]; \
	if (x < 1) return y [1]; \
	if (x == midleft) return y [midleft]; \
	/* 1 < x < nx && x not integer: interpolate. */ \
	if (maxDepth > midright - 1) maxDepth = midright - 1; \
	if (maxDepth > nx - midleft) maxDepth = nx - midleft; \
	if (maxDepth <= NUM_VALUE_INTERPOLATE_NEAREST) return y [(long) floor (x + 0.5)]; \
	if (maxDepth == NUM_VALUE_INTERPOLATE_LINEAR) return y [midleft] + (x - midleft) * (y [midright] - y [midleft]); \
	if (maxDepth == NUM_VALUE_INTERPOLATE_CUBIC) { \
		double yl = y [midleft], yr = y [midright]; \
		double dyl = 0.5 * (yr - y [midleft - 1]), dyr = 0.5 * (y [midright + 1] - yl); \
		double fil = x - midleft, fir = midright - x; \
		return yl * fir + yr * fil - fil * fir * (0.5 * (dyr - dyl) + (fil - 0.5) * (dyl + dyr - 2 * (yr - yl))); \
	}

#if defined (sgi) || defined (sun) || defined (HPUX) || __POWERPC__
#define NUM_interpolate_sinc_MACRO(proc,type) \
double proc (type y [], long nx, double x, long maxDepth) { \
	long ix, midleft = floor (x), midright = midleft + 1, left, right; \
	double result = 0.0, a, halfsina, aa, daa, cosaa, sinaa, cosdaa, sindaa; \
	NUM_interpolate_simple_cases \
	left = midright - maxDepth, right = midleft + maxDepth; \
	a = NUMpi * (x - midleft); \
	halfsina = 0.5 * sin (a); \
	aa = a / (x - left + 1); cosaa = cos (aa); sinaa = sin (aa); \
	daa = NUMpi / (x - left + 1); cosdaa = cos (daa); sindaa = sin (daa); \
	for (ix = midleft; ix >= left; ix --) { \
		double d = halfsina / a * (1.0 + cosaa), help; \
		result += y [ix] * d; \
		a += NUMpi; \
		help = cosaa * cosdaa - sinaa * sindaa; \
		sinaa = cosaa * sindaa + sinaa * cosdaa; \
		cosaa = help; \
		halfsina = - halfsina; \
	} \
	a = NUMpi * (midright - x); \
	halfsina = 0.5 * sin (a); \
	aa = a / (right - x + 1); cosaa = cos (aa); sinaa = sin (aa); \
	daa = NUMpi / (right - x + 1); cosdaa = cos (daa); sindaa = sin (daa); \
	for (ix = midright; ix <= right; ix ++) { \
		double d = halfsina / a * (1.0 + cosaa), help; \
		result += y [ix] * d; \
		a += NUMpi;	 \
		help = cosaa * cosdaa - sinaa * sindaa; \
		sinaa = cosaa * sindaa + sinaa * cosdaa; \
		cosaa = help; \
		halfsina = - halfsina; \
	} \
	return result; \
}
#else
#define NUM_interpolate_sinc_MACRO(proc,type) \
double proc (type y [], long nx, double x, long maxDepth) { \
	long ix, midleft = floor (x), midright = midleft + 1, left, right; \
	double result = 0.0, a, halfsina, aa, daa; \
	NUM_interpolate_simple_cases \
	left = midright - maxDepth, right = midleft + maxDepth; \
	a = NUMpi * (x - midleft); \
	halfsina = 0.5 * sin (a); \
	aa = a / (x - left + 1); \
	daa = NUMpi / (x - left + 1); \
	for (ix = midleft; ix >= left; ix --) { \
		double d = halfsina / a * (1.0 + cos (aa)); \
		result += y [ix] * d; \
		a += NUMpi;	 \
		aa += daa;	 \
		halfsina = - halfsina; \
	} \
	a = NUMpi * (midright - x); \
	halfsina = 0.5 * sin (a); \
	aa = a / (right - x + 1); \
	daa = NUMpi / (right - x + 1); \
	for (ix = midright; ix <= right; ix ++) { \
		double d = halfsina / a * (1.0 + cos (aa)); \
		result += y [ix] * d; \
		a += NUMpi;	 \
		aa += daa; \
		halfsina = - halfsina; \
	} \
	return result; \
}
#endif
NUM_interpolate_sinc_MACRO (NUM_interpolate_sinc_f, float)
NUM_interpolate_sinc_MACRO (NUM_interpolate_sinc_d, double)

/********** Improving extrema **********/

struct improve_params {
	int depth;
	float *y_f;
	double *y_d;
	long ixmax;
	int isMaximum;
};

static double improve_evaluate_f (double x, void *closure) {
	struct improve_params *me = closure;
	double y = NUM_interpolate_sinc_f (my y_f, my ixmax, x, my depth);
	return my isMaximum ? - y : y;
}
static double improve_evaluate_d (double x, void *closure) {
	struct improve_params *me = closure;
	double y = NUM_interpolate_sinc_d (my y_d, my ixmax, x, my depth);
	return my isMaximum ? - y : y;
}

double NUMimproveExtremum_f (float *y, long nx, long ixmid, int interpolation, double *ixmid_real, int isMaximum) {
	struct improve_params params;
	double result;
	if (ixmid <= 1) { *ixmid_real = 1; return y [1]; }
	if (ixmid >= nx) { *ixmid_real = nx; return y [nx]; }
	if (interpolation <= NUM_PEAK_INTERPOLATE_NONE) { *ixmid_real = ixmid; return y [ixmid]; }
	if (interpolation == NUM_PEAK_INTERPOLATE_PARABOLIC) {
		double dy = 0.5 * (y [ixmid + 1] - y [ixmid - 1]);
		double d2y = 2 * y [ixmid] - y [ixmid - 1] - y [ixmid + 1];
		*ixmid_real = ixmid + dy / d2y;
		return y [ixmid] + 0.5 * dy * dy / d2y;
	}
	/* Sinc interpolation. */
	params. y_f = y;
	params. depth = interpolation == NUM_PEAK_INTERPOLATE_CUBIC ? 2 :
		interpolation == NUM_PEAK_INTERPOLATE_SINC70 ? 70 : 700;
	params. ixmax = nx;
	params. isMaximum = isMaximum;
	/*return isMaximum ?
		- NUM_minimize (ixmid - 1, ixmid, ixmid + 1, improve_evaluate_f, & params, 1e-10, 1e-7, ixmid_real) :
		NUM_minimize (ixmid - 1, ixmid, ixmid + 1, improve_evaluate_f, & params, 1e-10, 1e-7, ixmid_real);*/
	*ixmid_real = NUMminimize_brent (improve_evaluate_f, ixmid - 1, ixmid + 1, & params, 1e-8, & result);
	return isMaximum ? - result : result;
}

double NUMimproveMaximum_f (float *y, long nx, long ixmid, int interpolation, double *ixmid_real)
	{ return NUMimproveExtremum_f (y, nx, ixmid, interpolation, ixmid_real, 1); }
double NUMimproveMinimum_f (float *y, long nx, long ixmid, int interpolation, double *ixmid_real)
	{ return NUMimproveExtremum_f (y, nx, ixmid, interpolation, ixmid_real, 0); }

double NUMimproveExtremum_d (double *y, long nx, long ixmid, int interpolation, double *ixmid_real, int isMaximum) {
	struct improve_params params;
	double result;
	if (ixmid <= 1) { *ixmid_real = 1; return y [1]; }
	if (ixmid >= nx) { *ixmid_real = nx; return y [nx]; }
	if (interpolation <= NUM_PEAK_INTERPOLATE_NONE) { *ixmid_real = ixmid; return y [ixmid]; }
	if (interpolation == NUM_PEAK_INTERPOLATE_PARABOLIC) {
		double dy = 0.5 * (y [ixmid + 1] - y [ixmid - 1]);
		double d2y = 2 * y [ixmid] - y [ixmid - 1] - y [ixmid + 1];
		*ixmid_real = ixmid + dy / d2y;
		return y [ixmid] + 0.5 * dy * dy / d2y;
	}
	/* Sinc interpolation. */
	params. y_d = y;
	params. depth = interpolation == NUM_PEAK_INTERPOLATE_SINC70 ? 70 : 700;
	params. ixmax = nx;
	params. isMaximum = isMaximum;
	/*return isMaximum ?
		- NUM_minimize (ixmid - 1, ixmid, ixmid + 1, improve_evaluate_d, & params, 1e-10, 1e-11, ixmid_real) :
		NUM_minimize (ixmid - 1, ixmid, ixmid + 1, improve_evaluate_d, & params, 1e-10, 1e-11, ixmid_real);*/
	*ixmid_real = NUMminimize_brent (improve_evaluate_d, ixmid - 1, ixmid + 1, & params, 1e-10, & result);
	return isMaximum ? - result : result;
}

double NUMimproveMaximum_d (double *y, long nx, long ixmid, int interpolation, double *ixmid_real)
	{ return NUMimproveExtremum_d (y, nx, ixmid, interpolation, ixmid_real, 1); }
double NUMimproveMinimum_d (double *y, long nx, long ixmid, int interpolation, double *ixmid_real)
	{ return NUMimproveExtremum_d (y, nx, ixmid, interpolation, ixmid_real, 0); }

/********** Viterbi **********/

int NUM_viterbi (
	long numberOfFrames, long maxnCandidates,
	long (*getNumberOfCandidates) (long iframe, void *closure),
	double (*getLocalCost) (long iframe, long icand, void *closure),
	double (*getTransitionCost) (long iframe, long icand1, long icand2, void *closure),
	void (*putResult) (long iframe, long place, void *closure),
	void *closure)
{
	double **delta = NULL, maximum;
	long **psi = NULL, *numberOfCandidates = NULL, iframe, icand, icand1, icand2, place;
	if (! (delta = NUMdmatrix (1, numberOfFrames, 1, maxnCandidates))) goto end;
	if (! (psi = NUMlmatrix (1, numberOfFrames, 1, maxnCandidates))) goto end;
	if (! (numberOfCandidates = NUMlvector (1, numberOfFrames))) goto end;
	for (iframe = 1; iframe <= numberOfFrames; iframe ++) {
		numberOfCandidates [iframe] = getNumberOfCandidates (iframe, closure);
		for (icand = 1; icand <= numberOfCandidates [iframe]; icand ++)
			delta [iframe] [icand] = - getLocalCost (iframe, icand, closure);
	}
	for (iframe = 2; iframe <= numberOfFrames; iframe ++) {
		for (icand2 = 1; icand2 <= numberOfCandidates [iframe]; icand2 ++) {
			maximum = -1e300;
			place = 0;
			for (icand1 = 1; icand1 <= numberOfCandidates [iframe - 1]; icand1 ++) {
				double value = delta [iframe - 1] [icand1] + delta [iframe] [icand2]
					- getTransitionCost (iframe, icand1, icand2, closure);
				if (value > maximum) { maximum = value; place = icand1; }
			}
			delta [iframe] [icand2] = maximum;
			psi [iframe] [icand2] = place;
		}
	}
	/*
	 * Find the end of the most probable path.
	 */
	maximum = delta [numberOfFrames] [place = 1];
	for (icand = 2; icand <= numberOfCandidates [numberOfFrames]; icand ++)
		if (delta [numberOfFrames] [icand] > maximum)
			maximum = delta [numberOfFrames] [place = icand];
	/*
	 * Backtrack.
	 */
	for (iframe = numberOfFrames; iframe >= 1; iframe --) {
		putResult (iframe, place, closure);
		place = psi [iframe] [place];
	}
end:
	NUMdmatrix_free (delta, 1, 1);
	NUMlmatrix_free (psi, 1, 1);
	NUMlvector_free (numberOfCandidates, 1);
	if (Melder_hasError ()) return 0;
	return 1;
}

/******************/

struct parm2 {
	int ntrack;
	long ncomb;
	long **indices;
	double (*getLocalCost) (long iframe, long icand, int itrack, void *closure);
	double (*getTransitionCost) (long iframe, long icand1, long icand2, int itrack, void *closure);
	void (*putResult) (long iframe, long place, int itrack, void *closure);
	void *closure;
};

static long getNumberOfCandidates_n (long iframe, void *closure) {
	struct parm2 *me = closure;
	(void) iframe;
	return my ncomb;
}
static double getLocalCost_n (long iframe, long jcand, void *closure) {
	struct parm2 *me = closure;
	int itrack;
	double localCost = 0.0;
	for (itrack = 1; itrack <= my ntrack; itrack ++)
		localCost += my getLocalCost (iframe, my indices [jcand] [itrack], itrack, my closure);
	return localCost;
}
static double getTransitionCost_n (long iframe, long jcand1, long jcand2, void *closure) {
	struct parm2 *me = closure;
	int itrack;
	double transitionCost = 0.0;
	for (itrack = 1; itrack <= my ntrack; itrack ++)
		transitionCost += my getTransitionCost (iframe,
			my indices [jcand1] [itrack], my indices [jcand2] [itrack], itrack, my closure);
	return transitionCost;
}
static void putResult_n (long iframe, long jplace, void *closure) {
	struct parm2 *me = closure;
	int itrack;
	for (itrack = 1; itrack <= my ntrack; itrack ++)
		my putResult (iframe, my indices [jplace] [itrack], itrack, my closure);
}

int NUM_viterbi_multi (
	long nframe, long ncand, int ntrack,
	double (*getLocalCost) (long iframe, long icand, int itrack, void *closure),
	double (*getTransitionCost) (long iframe, long icand1, long icand2, int itrack, void *closure),
	void (*putResult) (long iframe, long place, int itrack, void *closure),
	void *closure)
{
	double ncomb;
	struct parm2 parm;
	int itrack, jtrack;
	long *icand = NULL, jcomb;
	parm.indices = NULL;

	if (ntrack > ncand) return Melder_error ("(NUMviterbin:) "
		"Number of tracks (%d) should not exceed number of candidates (%ld).", ntrack, ncand);
	ncomb = NUMcombinations (ncand, ntrack);
	if (ncomb > 10000000) return Melder_error ("(NUMviterbin:) "
		"Unrealistically high number of combinations (%ld).", ncomb);
	parm. ntrack = ntrack;
	parm. ncomb = ncomb;

	/*
	 * For ncand == 5 and ntrack == 3, parm.indices is going to contain:
	 *   1 2 3
	 *   1 2 4
	 *   1 2 5
	 *   1 3 4
	 *   1 3 5
	 *   1 4 5
	 *   2 3 4
	 *   2 3 5
	 *   2 4 5
	 *   3 4 5
	 */
	if (! (parm. indices = NUMlmatrix (1, parm. ncomb, 1, ntrack)) ||
	    ! (icand = NUMlvector (1, ntrack))) goto end;
	for (itrack = 1; itrack <= ntrack; itrack ++)
		icand [itrack] = itrack;   /* Start out with "1 2 3". */
	jcomb = 0;
	for (;;) {
		jcomb ++;
		for (itrack = 1; itrack <= ntrack; itrack ++)
			parm. indices [jcomb] [itrack] = icand [itrack];
		for (itrack = ntrack; itrack >= 1; itrack --) {
			if (++ icand [itrack] <= ncand - (ntrack - itrack)) {
				for (jtrack = itrack + 1; jtrack <= ntrack; jtrack ++)
					icand [jtrack] = icand [itrack] + jtrack - itrack;
				break;
			}
		}
		if (itrack == 0) break;
	}
	Melder_assert (jcomb == ncomb);
	parm. getLocalCost = getLocalCost;
	parm. getTransitionCost = getTransitionCost;
	parm. putResult = putResult;
	parm. closure = closure;
	if (! NUM_viterbi (nframe, ncomb, getNumberOfCandidates_n, getLocalCost_n, getTransitionCost_n,
		putResult_n, & parm)) goto end;
end:
	NUMlmatrix_free (parm. indices, 1, 1);
	NUMlvector_free (icand, 1);
	if (Melder_hasError ()) return 0;
	return 1;
}

int NUMrotationsPointInPolygon (double x0, double y0, long n, float x [], float y []) {
	long nup = 0, i;
	int upold = y [n] > y0, upnew;
	for (i = 1; i <= n; i ++) if ((upnew = y [i] > y0) != upold) {
		long j = i == 1 ? n : i - 1;
		if (x0 < x [i] + (x [j] - x [i]) * (y0 - y [i]) / (y [j] - y [i])) {
			if (upnew) nup ++; else nup --;
		}
		upold = upnew;
	}
	return nup;
}

/* End of file NUM.c */