/* PointProcess.c
 *
 * Copyright (C) 1992-2003 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 2001/03/26
 * pb 2002/03/08 speed up creation of Poisson process
 * pb 2002/03/08 GPL
 * pb 2002/03/08 getMeanPeriod, getStdevPeriod, more jitter functions
 * pb 2003/05/21 more jitter functions
 */

#include "PointProcess.h"

#include "oo_DESTROY.h"
#include "PointProcess_def.h"
#include "oo_COPY.h"
#include "PointProcess_def.h"
#include "oo_EQUAL.h"
#include "PointProcess_def.h"
#include "oo_WRITE_ASCII.h"
#include "PointProcess_def.h"
#include "oo_READ_ASCII.h"
#include "PointProcess_def.h"
#include "oo_WRITE_BINARY.h"
#include "PointProcess_def.h"
#include "oo_READ_BINARY.h"
#include "PointProcess_def.h"
#include "oo_DESCRIPTION.h"
#include "PointProcess_def.h"

static void infoPeriods (PointProcess me, double shortestPeriod, double longestPeriod, int precision) {
	long numberOfPeriods = PointProcess_getNumberOfPeriods (me, 0.0, 0.0, shortestPeriod, longestPeriod);
	double meanPeriod = PointProcess_getMeanPeriod (me, 0.0, 0.0, shortestPeriod, longestPeriod);
	double stdevPeriod = PointProcess_getStdevPeriod (me, 0.0, 0.0, shortestPeriod, longestPeriod);
	double jitter_local = PointProcess_getJitter_local (me, 0.0, 0.0, shortestPeriod, longestPeriod);
	double jitter_local_absolute = PointProcess_getJitter_local_absolute (me, 0.0, 0.0, shortestPeriod, longestPeriod);
	double jitter_rap = PointProcess_getJitter_rap (me, 0.0, 0.0, shortestPeriod, longestPeriod);
	double jitter_ppq5 = PointProcess_getJitter_ppq5 (me, 0.0, 0.0, shortestPeriod, longestPeriod);
	double jitter_ddp = PointProcess_getJitter_ddp (me, 0.0, 0.0, shortestPeriod, longestPeriod);
	MelderInfo_writeLine2 ("     Number of periods: ", Melder_integer (numberOfPeriods));
	MelderInfo_writeLine3 ("     Mean period: ", Melder_double (meanPeriod), " seconds");
	MelderInfo_writeLine3 ("     Stdev period: ", Melder_double (stdevPeriod), " seconds");
	MelderInfo_writeLine2 ("     Jitter (local): ", Melder_percent (jitter_local, precision));
	MelderInfo_writeLine3 ("     Jitter (local, absolute): ", Melder_fixedExponent (jitter_local_absolute, -6, precision), " seconds");
	MelderInfo_writeLine2 ("     Jitter (rap): ", Melder_percent (jitter_rap, precision));
	MelderInfo_writeLine2 ("     Jitter (ppq5): ", Melder_percent (jitter_ppq5, precision));
	MelderInfo_writeLine2 ("     Jitter (ddp): ", Melder_percent (jitter_ddp, precision));
	MelderInfo_close ();
}

static void info (I) {
	iam (PointProcess);
	Melder_info ("Time domain: from %.15g to %.15g seconds.", my xmin, my xmax);
	Melder_info ("Number of times: %ld", my nt);
	if (my nt) {
		Melder_info ("First time: %.15g seconds.", my t [1]);
		Melder_info ("Last time: %.15g seconds.", my t [my nt]);
	}
	Melder_info ("Periods between 0.1 ms and 20 ms (pitch between 50 and 10000 Hz):");
	infoPeriods (me, 1e-4, 20e-3, 3);
	Melder_info ("All periods:");
	infoPeriods (me, 0.0, 0.0, 6);
}

class_methods (PointProcess, Function)
	class_method_local (PointProcess, destroy)
	class_method (info)
	class_method_local (PointProcess, description)
	class_method_local (PointProcess, copy)
	class_method_local (PointProcess, equal)
	class_method_local (PointProcess, writeAscii)
	class_method_local (PointProcess, readAscii)
	class_method_local (PointProcess, writeBinary)
	class_method_local (PointProcess, readBinary)
class_methods_end

int PointProcess_init (I, double tmin, double tmax, long initialMaxnt) {
	iam (PointProcess);
	if (! Function_init (me, tmin, tmax)) return 0;
	if (initialMaxnt < 1) initialMaxnt = 1;
	my maxnt = initialMaxnt;
	my nt = 0;
	if (! (my t = NUMdvector (1, my maxnt))) return 0;
	return 1;
}

PointProcess PointProcess_create (double tmin, double tmax, long initialMaxnt) {
	PointProcess me = new (PointProcess);
	if (! me || ! PointProcess_init (me, tmin, tmax, initialMaxnt)) forget (me);
	return me;
}

PointProcess PointProcess_createPoissonProcess (double startingTime, double finishingTime, double density) {
	long nt = NUMrandomPoisson ((finishingTime - startingTime) * density), i;
	PointProcess me = PointProcess_create (startingTime, finishingTime, nt);
	if (! me) return NULL;
	my nt = nt;
	for (i = 1; i <= nt; i ++)
		my t [i] = NUMrandomUniform (startingTime, finishingTime);
	NUMsort_d (my nt, my t);
	return me;
}

long PointProcess_getLowIndex (PointProcess me, double t) {
	long left, right;
	if (my nt == 0 || t < my t [1])
		return 0;
	if (t >= my t [my nt])   /* Special case that often occurs in practice. */
		return my nt;
	Melder_assert (my nt != 1);   /* May fail if my t [1] is NaN. */
	/* Start binary search. */
	left = 1, right = my nt;
	while (left < right - 1) {
		long mid = (left + right) / 2;
		if (t >= my t [mid]) left = mid; else right = mid;
	}
	Melder_assert (right == left + 1);
	return left;
}

long PointProcess_getHighIndex (PointProcess me, double t) {
	long left, right;
	if (my nt == 0)
		return 0;
	if (t <= my t [1])
		return 1;
	if (t > my t [my nt])
		return my nt + 1;
	/* Start binary search. */
	left = 1, right = my nt;
	while (left < right - 1) {
		long mid = (left + right) / 2;
		if (t > my t [mid]) left = mid; else right = mid;
	}
	Melder_assert (right == left + 1);
	return right;
}

long PointProcess_getNearestIndex (PointProcess me, double t) {
	long left, right;
	if (my nt == 0)
		return 0;
	if (t <= my t [1])
		return 1;
	if (t >= my t [my nt])
		return my nt;
	/* Start binary search. */
	left = 1, right = my nt;
	while (left < right - 1) {
		long mid = (left + right) / 2;
		if (t >= my t [mid]) left = mid; else right = mid;
	}
	Melder_assert (right == left + 1);
	return t - my t [left] < my t [right] - t ? left : right;
}

int PointProcess_addPoint (PointProcess me, double t) {
	Melder_assert (t != NUMundefined);
	if (my nt >= my maxnt) {
		double *dum = NUMdvector (1, 2 * my maxnt);
		if (! dum) return 0;
		NUMdvector_copyElements (my t, dum, 1, my nt);
		NUMdvector_free (my t, 1);
		my t = dum;
		my maxnt *= 2;
	}
	if (my nt == 0 || t >= my t [my nt]) {   /* Special case that often occurs in practice. */
		my t [++ my nt] = t;
	} else {
		long left = PointProcess_getLowIndex (me, t), i;
		for (i = my nt; i > left; i --) my t [i + 1] = my t [i];
		my nt ++;
		my t [left + 1] = t;
	}
	return 1;
}

void PointProcess_removePoint (PointProcess me, long index) {
	long i;
	if (index < 1 || index > my nt) return;
	for (i = index; i < my nt; i ++)
		my t [i] = my t [i + 1];
	my nt --;
}

void PointProcess_removePointNear (PointProcess me, double t) {
	PointProcess_removePoint (me, PointProcess_getNearestIndex (me, t));
}

void PointProcess_removePoints (PointProcess me, long first, long last) {
	long i, distance;
	if (first < 1) first = 1;
	if (last > my nt) last = my nt;
	if ((distance = last - first + 1) <= 0) return;
	for (i = first + distance; i <= my nt; i ++)
		my t [i - distance] = my t [i];
	my nt -= distance;
}

void PointProcess_removePointsBetween (PointProcess me, double tmin, double tmax) {
	PointProcess_removePoints (me, PointProcess_getHighIndex (me, tmin), PointProcess_getLowIndex (me, tmax));
}

void PointProcess_draw (PointProcess me, Graphics g, double tmin, double tmax, int garnish) {
	if (tmax <= tmin) { tmin = my xmin; tmax = my xmax; }
	Graphics_setWindow (g, tmin, tmax, -1, 1);
	if (my nt) {
		long imin = PointProcess_getHighIndex (me, tmin), imax = PointProcess_getLowIndex (me, tmax), i;
		int lineType = Graphics_inqLineType (g);
		Graphics_setLineType (g, Graphics_DOTTED);
		Graphics_setInner (g);
		for (i = imin; i <= imax; i ++) {
			Graphics_line (g, my t [i], -1, my t [i], 1);
		}
		Graphics_setLineType (g, lineType);
		Graphics_unsetInner (g);
	}
	if (garnish) {
		Graphics_drawInnerBox (g);
		Graphics_textBottom (g, 1, "Time (s)");
		Graphics_marksBottom (g, 2, 1, 1, 0);
	}
}

double PointProcess_getInterval (PointProcess me, double t) {
	long ileft = PointProcess_getLowIndex (me, t);
	if (ileft <= 0 || ileft >= my nt) return NUMundefined;
	return my t [ileft + 1] - my t [ileft];
}

PointProcess PointProcesses_union (PointProcess me, PointProcess thee) {
	PointProcess him = Data_copy (me);
	long i;
	if (thy xmin < my xmin) his xmin = thy xmin;
	if (thy xmax > my xmax) his xmax = thy xmax;
	for (i = 1; i <= thy nt; i ++)
		if (! PointProcess_addPoint (him, thy t [i])) { forget (him); return NULL; }
	return him;
}

long PointProcess_findPoint (PointProcess me, double t) {
	long left = 1, right = my nt;
	if (my nt == 0) return 0;
	if (t < my t [left] || t > my t [right]) return 0;
	while (left < right - 1) {
		long mid = (left + right) / 2;   /* tleft <= t <= tright */
		if (t == my t [mid]) return mid;
		if (t > my t [mid])
			left = mid;
		else
			right = mid;
	}
	if (t == my t [left]) return left;
	if (t == my t [right]) return right;
	return 0;
}

PointProcess PointProcesses_intersection (PointProcess me, PointProcess thee) {
	PointProcess him = Data_copy (me);
	long i;
	if (thy xmin > my xmin) his xmin = thy xmin;
	if (thy xmax < my xmax) his xmax = thy xmax;
	for (i = my nt; i >= 1; i --)
		if (! PointProcess_findPoint (thee, my t [i]))
			PointProcess_removePoint (him, i);
	return him;
}

PointProcess PointProcesses_difference (PointProcess me, PointProcess thee) {
	PointProcess him = Data_copy (me);
	long i;
	for (i = my nt; i >= 1; i --)
		if (PointProcess_findPoint (thee, my t [i]))
			PointProcess_removePoint (him, i);
	return him;
}

int PointProcess_fill (PointProcess me, double tmin, double tmax, double period) {
	long n, i;
	double t;
	if (tmax <= tmin) tmin = my xmin, tmax = my xmax;   /* Autowindowing. */
	n = floor ((tmax - tmin) / period);
	t = 0.5 * (tmin + tmax - n * period);
	for (i = 1, t = 0.5 * (tmin + tmax - n * period); i <= n; i ++, t += period)
		if (! PointProcess_addPoint (me, t)) return 0;
	return 1;
}

int PointProcess_voice (PointProcess me, double period, double maxT) {
	long ipointleft, ipointright;
	double beginVoiceless = my xmin, endVoiceless;
	for (ipointleft = 1; ipointleft <= my nt; ipointleft = ipointright + 1) {
		endVoiceless = my t [ipointleft];
		if (! PointProcess_fill (me, beginVoiceless, endVoiceless, period)) return 0;
		for (ipointright = ipointleft + 1; ipointright <= my nt; ipointright ++)
			if (my t [ipointright] - my t [ipointright - 1] > maxT)
				break;
		ipointright --;
		beginVoiceless = my t [ipointright] + 0.005;
	}
	endVoiceless = my xmax;
	PointProcess_fill (me, beginVoiceless, endVoiceless, period);
	return 1;
}

long PointProcess_getWindowPoints (PointProcess me, double tmin, double tmax, long *imin, long *imax) {
	*imin = PointProcess_getHighIndex (me, tmin);
	*imax = PointProcess_getLowIndex (me, tmax);
	return *imax - *imin + 1;
}

long PointProcess_getNumberOfPeriods (PointProcess me, double tmin, double tmax, double pmin, double pmax) {
	long imin, imax, numberOfPeriods, i;
	if (tmax <= tmin) tmin = my xmin, tmax = my xmax;   /* Autowindowing. */
	numberOfPeriods = PointProcess_getWindowPoints (me, tmin, tmax, & imin, & imax) - 1;
	if (numberOfPeriods < 1) return 0;
	for (i = imin; i < imax; i ++) {
		double period = my t [i + 1] - my t [i];
		if (pmin == pmax || period >= pmin && period <= pmax) {
			(void) 0;
		} else {
			numberOfPeriods --;
		}
	}
	return numberOfPeriods;
}

double PointProcess_getMeanPeriod (PointProcess me, double tmin, double tmax, double pmin, double pmax) {
	long imin, imax, numberOfPeriods, i;
	double sum = 0.0;
	if (tmax <= tmin) tmin = my xmin, tmax = my xmax;   /* Autowindowing. */
	numberOfPeriods = PointProcess_getWindowPoints (me, tmin, tmax, & imin, & imax) - 1;
	if (numberOfPeriods < 1) return NUMundefined;
	for (i = imin; i < imax; i ++) {
		double period = my t [i + 1] - my t [i];
		if (pmin == pmax || period >= pmin && period <= pmax) {
			sum += period;
		} else {
			numberOfPeriods --;
		}
	}
	return numberOfPeriods > 0 ? sum / numberOfPeriods : NUMundefined;
}

double PointProcess_getStdevPeriod (PointProcess me, double tmin, double tmax, double pmin, double pmax) {
	long imin, imax, numberOfPeriods, i;
	double sum = 0.0, sum2 = 0.0, mean;
	if (tmax <= tmin) tmin = my xmin, tmax = my xmax;   /* Autowindowing. */
	numberOfPeriods = PointProcess_getWindowPoints (me, tmin, tmax, & imin, & imax) - 1;
	if (numberOfPeriods < 2) return NUMundefined;
	for (i = imin; i < imax; i ++) {
		double period = my t [i + 1] - my t [i];
		if (pmin == pmax || period >= pmin && period <= pmax) {
			sum += period;
		} else {
			numberOfPeriods --;
		}
	}
	if (numberOfPeriods < 2) return NUMundefined;
	mean = sum / numberOfPeriods;
	for (i = imin; i < imax; i ++) {
		double period = my t [i + 1] - my t [i];
		if (pmin == pmax || period >= pmin && period <= pmax) {
			double dperiod = period - mean;
			sum2 += dperiod * dperiod;
		}
	}
	return sqrt (sum2 / (numberOfPeriods - 1));
}

double PointProcess_getJitter_local (PointProcess me, double tmin, double tmax, double pmin, double pmax) {
	long imin, imax, numberOfPeriods, i;
	double sum = 0.0;
	if (tmax <= tmin) tmin = my xmin, tmax = my xmax;   /* Autowindowing. */
	numberOfPeriods = PointProcess_getWindowPoints (me, tmin, tmax, & imin, & imax) - 1;
	if (numberOfPeriods < 2) return NUMundefined;
	for (i = imin + 1; i < imax; i ++) {
		double p1 = my t [i] - my t [i - 1], p2 = my t [i + 1] - my t [i];
		if (pmin == pmax || p1 >= pmin && p1 <= pmax && p2 >= pmin && p2 <= pmax) {
			sum += fabs (p1 - p2);
		} else {
			numberOfPeriods --;
		}
	}
	if (numberOfPeriods < 2) return NUMundefined;
	return sum / (numberOfPeriods - 1) / PointProcess_getMeanPeriod (me, tmin, tmax, pmin, pmax);
}

double PointProcess_getJitter_local_absolute (PointProcess me, double tmin, double tmax, double pmin, double pmax) {
	long imin, imax, numberOfPeriods, i;
	double sum = 0.0;
	if (tmax <= tmin) tmin = my xmin, tmax = my xmax;   /* Autowindowing. */
	numberOfPeriods = PointProcess_getWindowPoints (me, tmin, tmax, & imin, & imax) - 1;
	if (numberOfPeriods < 2) return NUMundefined;
	for (i = imin + 1; i < imax; i ++) {
		double p1 = my t [i] - my t [i - 1], p2 = my t [i + 1] - my t [i];
		if (pmin == pmax || p1 >= pmin && p1 <= pmax && p2 >= pmin && p2 <= pmax) {
			sum += fabs (p1 - p2);
		} else {
			numberOfPeriods --;
		}
	}
	if (numberOfPeriods < 2) return NUMundefined;
	return sum / (numberOfPeriods - 1);
}

double PointProcess_getJitter_rap (PointProcess me, double tmin, double tmax, double pmin, double pmax) {
	long imin, imax, numberOfPeriods, i;
	double sum = 0.0;
	if (tmax <= tmin) tmin = my xmin, tmax = my xmax;   /* Autowindowing. */
	numberOfPeriods = PointProcess_getWindowPoints (me, tmin, tmax, & imin, & imax) - 1;
	if (numberOfPeriods < 3) return NUMundefined;
	for (i = imin + 2; i < imax; i ++) {
		double p1 = my t [i - 1] - my t [i - 2], p2 = my t [i] - my t [i - 1], p3 = my t [i + 1] - my t [i];
		if (pmin == pmax || p1 >= pmin && p1 <= pmax && p2 >= pmin && p2 <= pmax && p3 >= pmin && p3 <= pmax) {
			sum += fabs (p2 - (p1 + p2 + p3) / 3.0);
		} else {
			numberOfPeriods --;
		}
	}
	if (numberOfPeriods < 3) return NUMundefined;
	return sum / (numberOfPeriods - 2) / PointProcess_getMeanPeriod (me, tmin, tmax, pmin, pmax);
}

double PointProcess_getJitter_ppq5 (PointProcess me, double tmin, double tmax, double pmin, double pmax) {
	long imin, imax, numberOfPeriods, i;
	double sum = 0.0;
	if (tmax <= tmin) tmin = my xmin, tmax = my xmax;   /* Autowindowing. */
	numberOfPeriods = PointProcess_getWindowPoints (me, tmin, tmax, & imin, & imax) - 1;
	if (numberOfPeriods < 5) return NUMundefined;
	for (i = imin + 5; i <= imax; i ++) {
		double
			p1 = my t [i - 4] - my t [i - 5],
			p2 = my t [i - 3] - my t [i - 4],
			p3 = my t [i - 2] - my t [i - 3],
			p4 = my t [i - 1] - my t [i - 2],
			p5 = my t [i] - my t [i - 1];
		if (pmin == pmax || p1 >= pmin && p1 <= pmax && p2 >= pmin && p2 <= pmax && p3 >= pmin && p3 <= pmax &&
			p4 >= pmin && p4 <= pmax && p5 >= pmin && p5 <= pmax)
		{
			sum += fabs (p3 - (p1 + p2 + p3 + p4 + p5) / 5.0);
		} else {
			numberOfPeriods --;
		}
	}
	if (numberOfPeriods < 3) return NUMundefined;
	return sum / (numberOfPeriods - 2) / PointProcess_getMeanPeriod (me, tmin, tmax, pmin, pmax);
}

double PointProcess_getJitter_ddp (PointProcess me, double tmin, double tmax, double pmin, double pmax) {
	double rap = PointProcess_getJitter_rap (me, tmin, tmax, pmin, pmax);
	return NUMdefined (rap) ? 3.0 * rap : NUMundefined;
}

/* End of file PointProcess.c */