package monashts.tempo.univariate.distances;

import static monashts.tempo.univariate.distances.utils.dist;
import static monashts.tempo.univariate.distances.utils.min;
import static java.lang.Double.POSITIVE_INFINITY;
import static java.lang.Double.min;
import static java.lang.Math.abs;

public class TWE {

	static double distance(double[] lines, double[] cols, double nu, double lambda, double cutoff) {
		// Ensure that lines are longer than columns
		if (lines.length < cols.length) {
			double[] swap = lines;
			lines = cols;
			cols = swap;
		}

		// --- --- --- Declarations
		int nblines = lines.length;
		int nbcols = cols.length;

		// Setup buffers - no extra initialization required - border condition managed in the code.
		double[] buffers = new double[2 * nbcols];
		int c = 0;
		int p = nbcols;

		// Buffer holding precomputed distance between columns
		double[] distcol = new double[nbcols];

		// Line & columns indices
		int i = 0;
		int j = 0;

		// Cost accumulator in a line, also used as the "left neighbor"
		double cost = 0;

		// EAP variable: track where to start the next line, and the position of the previous pruning point.
		// Must be init to 0: index 0 is the next starting point and also the "previous pruning point"
		int next_start = 0;
		int prev_pp = 0;

		// Constants: we only consider timestamp spaced by 1, so:
		// In the "delete" case, we always have a time difference of 1, so we always have 1*nu+lambda
		final double nu_lambda = nu + lambda;
		// In the "match" case, we always have nu*(|i-j|+|(i-1)-(j-1)|) == 2*nu*|i-j|
		final double nu2 = 2 * nu;

		// --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- ---
		// Create a new tighter upper bounds using the last alignment.
		// The last alignment can only be computed if we have nbcols >= 2
		double ub = cutoff;
		if (nbcols >= 2) {
			double li = lines[nblines - 1];
			double li1 = lines[nblines - 2];
			double co = cols[nbcols - 1];
			double co1 = cols[nbcols - 2];
			double distli = dist(li1, li);
			double distco = dist(co1, co);
			double la = min(distco + nu_lambda, // "Delete_B": over the columns / Prev
					dist(li, co) + dist(li1, co1) + nu2 * (nblines - nbcols), // Match: Diag. Ok: nblines >= nbcols
					distli + nu_lambda // "Delete_A": over the lines / Top
			);
			ub = (cutoff + utils.EPSILON) - la;
		}

		// --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- ---
		// Initialisation of the first line. Deal with the line top border condition.
		{
			// Case [0,0]: special "Match case"
			cost = dist(lines[0], cols[0]);
			buffers[c + 0] = cost;
			// Distance for the first column is relative to 0 "by conventions" (from the paper, section 4.2)
			distcol[0] = dist(0, cols[0]);
			// Rest of the line: [i==0, j>=1]: "Delete_B case" (prev). We also initialize 'distcol' here.
			for (j = 1; j < nbcols; ++j) {
				double d = dist(cols[j - 1], cols[j]);
				distcol[j] = d;
				cost = cost + d + nu_lambda;
				buffers[c + j] = cost;
				if (cost <= ub) {
					prev_pp = j + 1;
				} else {
					break;
				}
			}
			// Complete the initialisation of distcol
			for (; j < nbcols; ++j) {
				double d = dist(cols[j - 1], cols[j]);
				distcol[j] = d;
			}
			// Next line.
			++i;
		}

		// --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- ---
		// Main loop, starts at the second line
		for (; i < nblines; ++i) {
			// --- --- --- Swap and variables init
			int swap = c;
			c = p;
			p = swap;
			double li = lines[i];
			double li1 = lines[i - 1];
			double distli = dist(li1, li);
			int curr_pp = next_start; // Next pruning point init at the start of the line
			j = next_start;
			// --- --- --- Stage 0: Special case for the first column. Can only look up (border on the left)
			{
				cost = buffers[p + j] + distli + nu_lambda; // "Delete_A" / Top
				buffers[c + j] = cost;
				if (cost <= ub) {
					curr_pp = j + 1;
				} else {
					++next_start;
				}
				++j;
			}
			// --- --- --- Stage 1: Up to the previous pruning point while advancing next_start: diag and top
			for (; j == next_start && j < prev_pp; ++j) {
				cost = min(buffers[p + j - 1] + dist(li, cols[j]) + dist(li1, cols[j - 1]) + nu2 * abs(i - j), // "Match": Diag
						buffers[p + j] + distli + nu_lambda // "Delete_A" / Top
				);
				buffers[c + j] = cost;
				if (cost <= ub) {
					curr_pp = j + 1;
				} else {
					++next_start;
				}
			}
			// --- --- --- Stage 2: Up to the previous pruning point without advancing next_start: left, diag and top
			for (; j < prev_pp; ++j) {
				cost = min(cost + distcol[j] + nu_lambda, // "Delete_B": over the columns / Prev
						buffers[p + j - 1] + dist(li, cols[j]) + dist(li1, cols[j - 1]) + nu2 * abs(i - j), // "Matc"h: Diag
						buffers[p + j] + distli + nu_lambda // "Delete_A": over the lines / Top
				);
				buffers[c + j] = cost;
				if (cost <= ub) {
					curr_pp = j + 1;
				}
			}
			// --- --- --- Stage 3: At the previous pruning point. Check if we are within bounds.
			if (j < nbcols) { // If so, two cases.
				if (j == next_start) { // Case 1: Advancing next start: only diag.
					cost = buffers[p + j - 1] + dist(li, cols[j]) + dist(li1, cols[j - 1]) + nu2 * abs(i - j); // "Match": Diag
					buffers[c + j] = cost;
					if (cost <= ub) {
						curr_pp = j + 1;
					} else {
						// Special case if we are on the last alignment: return the actual cost if we are <= cutoff
						if (i == nblines - 1 && j == nbcols - 1 && cost <= cutoff) {
							return cost;
						} else {
							return POSITIVE_INFINITY;
						}
					}
				} else { // Case 2: Not advancing next start: possible path in previous cells: left and diag.
					cost = min(cost + distcol[j] + nu_lambda, // "Delete_B": over the columns / Prev
							buffers[p + j - 1] + dist(li, cols[j]) + dist(li1, cols[j - 1]) + nu2 * abs(i - j) // "Match": Diag
					);
					buffers[c + j] = cost;
					if (cost <= ub) {
						curr_pp = j + 1;
					}
				}
				++j;
			} else { // Previous pruning point is out of bound: exit if we extended next start up to here.
				if (j == next_start) {
					// But only if we are above the original UB
					// Else set the next starting point to the last valid column
					if (cost > cutoff) {
						return POSITIVE_INFINITY;
					} else {
						next_start = nbcols - 1;
					}
				}
			}
			// --- --- --- Stage 4: After the previous pruning point: only prev.
			// Go on while we advance the curr_pp; if it did not advance, the rest of the line is guaranteed to be > ub.
			for (; j == curr_pp && j < nbcols; ++j) {
				cost = cost + distcol[j] + nu_lambda; // "Delete_B": over the columns / Prev
				buffers[c + j] = cost;
				if (cost <= ub) {
					++curr_pp;
				}
			}
			// --- --- ---
			prev_pp = curr_pp;
		} // End of main loop for(;i<nblines;++i)

		// --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- ---
		// Finalization
		// Check for last alignment (i==nblines implied, Stage 4 implies j<=nbcols).
		// Cost must be <= original bound.
		if (j == nbcols && cost <= cutoff) {
			return cost;
		} else {
			return POSITIVE_INFINITY;
		}

	}

}
