/* Vector.c
 *
 * Copyright (C) 1992-2002 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/06/04
 * pb 2002/07/16 GPL
 */

#include "Vector.h"

static double getVector (I, long icol) {
	iam (Vector);
	if (icol < 1 || icol > my nx) return 0.0;
	return my z [1] [icol];
}

static double getFunction1 (I, double x) {
	iam (Vector);
	double rcol = (x - my x1) / my dx + 1.0;
	long icol = floor (rcol);
	double dcol = rcol - icol;
	double z1 = icol < 1 || icol > my nx ? 0.0 : my z [1] [icol];
	double z2 = icol < 0 || icol >= my nx ? 0.0 : my z [1] [icol + 1];
	return (1.0 - dcol) * z1 + dcol * z2;
}

class_methods (Vector, Matrix)
	class_method (getVector)
	class_method (getFunction1)
	us -> getMatrix = NULL;
	us -> getFunction2 = NULL;
class_methods_end

/***** Get content. *****/

double Vector_getValueAtX (I, double x, int interpolation) { iam (Vector);
	double leftEdge = my x1 - 0.5 * my dx, rightEdge = leftEdge + my nx * my dx;
	if (x <  leftEdge || x > rightEdge) return NUMundefined;
	return NUM_interpolate_sinc_f (my z [1], my nx, Sampled_xToIndex (me, x),
		interpolation == 3 ? 70 : interpolation == 4 ? 700 : interpolation);
}

/***** Get shape. *****/

void Vector_getMinimumAndX (I, double xmin, double xmax, int interpolation,
	double *return_minimum, double *return_xOfMinimum)
{
	iam (Vector);
	long imin, imax, i, n = my nx;
	float *y = my z [1];
	double minimum, x;
	if (xmax <= xmin) { xmin = my xmin; xmax = my xmax; }
	if (! Sampled_getWindowSamples (me, xmin, xmax, & imin, & imax)) {
		/*
		 * No samples between xmin and xmax.
		 * Try to return the lesser of the values at these two points.
		 */
		double yleft = Vector_getValueAtX (me, xmin, interpolation > 0 ? 1 : 0);
		double yright = Vector_getValueAtX (me, xmax, interpolation > 0 ? 1 : 0);
		minimum = yleft < yright ? yleft : yright;
		x = yleft == yright ? (xmin + xmax) / 2 : yleft < yright ? xmin : xmax;
	} else {
		minimum = y [imin], x = imin;
		if (y [imax] < minimum) minimum = y [imax], x = imax;
		if (imin == 1) imin ++;
		if (imax == my nx) imax --;
		for (i = imin; i <= imax; i ++) {
			if (y [i] < y [i - 1] && y [i] <= y [i + 1]) {
				double i_real, localMinimum = NUMimproveMinimum_f (y, n, i, interpolation, & i_real);
				if (localMinimum < minimum) minimum = localMinimum, x = i_real;
			}
		}
		x = my x1 + (x - 1) * my dx;   /* Convert sample to x. */
		if (x < xmin) x = xmin; else if (x > xmax) x = xmax;
	}
	if (return_minimum) *return_minimum = minimum;
	if (return_xOfMinimum) *return_xOfMinimum = x;
}
double Vector_getMinimum (I, double xmin, double xmax, int interpolation) { iam (Vector);
	double minimum;
	Vector_getMinimumAndX (me, xmin, xmax, interpolation, & minimum, NULL);
	return minimum;
}
double Vector_getXOfMinimum (I, double xmin, double xmax, int interpolation) { iam (Vector);
	double xOfMinimum;
	Vector_getMinimumAndX (me, xmin, xmax, interpolation, NULL, & xOfMinimum);
	return xOfMinimum;
}
void Vector_getMaximumAndX (I, double xmin, double xmax, int interpolation,
	double *return_maximum, double *return_xOfMaximum)
{
	iam (Vector);
	long imin, imax, i, n = my nx;
	float *y = my z [1];
	double maximum, x;
	if (xmax <= xmin) { xmin = my xmin; xmax = my xmax; }
	if (! Sampled_getWindowSamples (me, xmin, xmax, & imin, & imax)) {
		/*
		 * No samples between xmin and xmax.
		 * Try to return the greater of the values at these two points.
		 */
		double yleft = Vector_getValueAtX (me, xmin, interpolation > 0 ? 1 : 0);
		double yright = Vector_getValueAtX (me, xmax, interpolation > 0 ? 1 : 0);
		maximum = yleft > yright ? yleft : yright;
		x = yleft == yright ? (xmin + xmax) / 2 : yleft > yright ? xmin : xmax;
	} else {
		maximum = y [imin], x = imin;
		if (y [imax] > maximum) maximum = y [imax], x = imax;
		if (imin == 1) imin ++;
		if (imax == my nx) imax --;
		for (i = imin; i <= imax; i ++) {
			if (y [i] > y [i - 1] && y [i] >= y [i + 1]) {
				double i_real, localMaximum = NUMimproveMaximum_f (y, n, i, interpolation, & i_real);
				if (localMaximum > maximum) maximum = localMaximum, x = i_real;
			}
		}
		x = my x1 + (x - 1) * my dx;   /* Convert sample to x. */
		if (x < xmin) x = xmin; else if (x > xmax) x = xmax;
	}
	if (return_maximum) *return_maximum = maximum;
	if (return_xOfMaximum) *return_xOfMaximum = x;
}
double Vector_getMaximum (I, double xmin, double xmax, int interpolation) { iam (Vector);
	double maximum;
	Vector_getMaximumAndX (me, xmin, xmax, interpolation, & maximum, NULL);
	return maximum;
}
double Vector_getXOfMaximum (I, double xmin, double xmax, int interpolation) { iam (Vector);
	double xOfMaximum;
	Vector_getMaximumAndX (me, xmin, xmax, interpolation, NULL, & xOfMaximum);
	return xOfMaximum;
}
double Vector_getAbsoluteExtremum (I, double xmin, double xmax, int interpolation) { iam (Vector);
	double minimum = fabs (Vector_getMinimum (me, xmin, xmax, interpolation));
	double maximum = fabs (Vector_getMaximum (me, xmin, xmax, interpolation));
	return minimum > maximum ? minimum : maximum;
}

/***** Get statistics. *****/

double Vector_getMean (I, double xmin, double xmax) { iam (Vector);
	long imin, imax, i, n;
	double sum = 0.0, phase;
	double leftEdge, rightEdge, range;
	if (xmax <= xmin) { xmin = my xmin; xmax = my xmax; }
	n = Sampled_getWindowSamples (me, xmin, xmax, & imin, & imax);
	if (n < 1) {
		double min = Vector_getValueAtX (me, xmin, 1);
		double max = Vector_getValueAtX (me, xmax, 1);
		return ! NUMdefined (min) || ! NUMdefined (max) ? NUMundefined : 0.5 * (min + max);
	}
	for (i = imin; i <= imax; i ++) sum += my z [1] [i];
	/*
	 * Corrections within the sample.
	 */
	leftEdge = my x1 - 0.5 * my dx, rightEdge = leftEdge + my nx * my dx;
	range = n;
	if (xmin > leftEdge) {
		phase = (my x1 + (imin - 1) * my dx - xmin) / my dx;
		sum -= 0.5 * my z [1] [imin];
		range += phase - 0.5;
		if (imin == 1) {
			sum += phase * my z [1] [imin];
		} else {
			sum += phase * (my z [1] [imin] + 0.5 * phase * (my z [1] [imin - 1] - my z [1] [imin]));
		}
	}
	if (xmax < rightEdge) {
		phase = (xmax - (my x1 + (imax - 1) * my dx)) / my dx;
		sum -= 0.5 * my z [1] [imax];
		range += phase - 0.5;
		if (imax == my nx) {
			sum += phase * my z [1] [imax];
		} else {
			sum += phase * (my z [1] [imax] + 0.5 * phase * (my z [1] [imax + 1] - my z [1] [imax]));
		}
	}
	return sum / range;
}

double Vector_getStandardDeviation (I, double xmin, double xmax) { iam (Vector);
	long imin, imax, i, n;
	double mean, sum2 = 0.0;
	if (xmax <= xmin) { xmin = my xmin; xmax = my xmax; }
	n = Sampled_getWindowSamples (me, xmin, xmax, & imin, & imax);
	if (n < 2) return NUMundefined;
	mean = Vector_getMean (me, xmin, xmax);
	for (i = imin; i <= imax; i ++) sum2 += (my z [1] [i] - mean) * (my z [1] [i] - mean);
	return sqrt (sum2 / (n - 1));
}

double Vector_getRootMeanSquare (I, double xmin, double xmax) { iam (Vector);
	long imin, imax, i, n;
	double sum2 = 0.0;
	if (xmax <= xmin) { xmin = my xmin; xmax = my xmax; }
	n = Sampled_getWindowSamples (me, xmin, xmax, & imin, & imax);
	if (n < 1) return NUMundefined;
	for (i = imin; i <= imax; i ++) sum2 += my z [1] [i] * my z [1] [i];
	return sqrt (sum2 / n);
}

double Vector_getEnergy (I, double xmin, double xmax) { iam (Vector);
	long imin, imax, i, n;
	double sum2 = 0.0;
	if (xmax <= xmin) { xmin = my xmin; xmax = my xmax; }
	n = Sampled_getWindowSamples (me, xmin, xmax, & imin, & imax);
	if (n < 1) return NUMundefined;
	for (i = imin; i <= imax; i ++) sum2 += my z [1] [i] * my z [1] [i];
	return sum2 * my dx;
}

double Vector_getPower (I, double xmin, double xmax) { iam (Vector);
	long imin, imax, i, n;
	double sum2 = 0.0;
	if (xmax <= xmin) { xmin = my xmin; xmax = my xmax; }
	n = Sampled_getWindowSamples (me, xmin, xmax, & imin, & imax);
	if (n < 1) return NUMundefined;
	for (i = imin; i <= imax; i ++) sum2 += my z [1] [i] * my z [1] [i];
	return sum2 / n;
}

/***** Modify. *****/

void Vector_addScalar (I, double scalar) { iam (Vector);
	float scalar_f = scalar, *amp = my z [1];
	long i, n = my nx;
	for (i = 1; i <= n; i ++) amp [i] += scalar_f;
}

void Vector_subtractMean (I) { iam (Vector);
	double mean = Vector_getMean (me, 0, 0);
	Vector_addScalar (me, - mean);
}

void Vector_multiplyByScalar (I, double scalar) { iam (Vector);
	float scalar_f = scalar, *amp = my z [1];
	long i, n = my nx;
	for (i = 1; i <= n; i ++) amp [i] *= scalar_f;
}

void Vector_scale (I, double scale) { iam (Vector);
	float *amp = my z [1];
	double extremum = 0.0;
	long i, n = my nx;
	for (i = 1; i <= n; i ++) if (fabs (amp [i]) > extremum) extremum = fabs (amp [i]);
	if (extremum != 0.0) {
		scale /= extremum;
		for (i = 1; i <= n; i ++) amp [i] *= scale;
	}
}

/* End of file Vector.c */