package monashts.tempo.univariate.distances;

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 MSM {

	static double split_merge_cost(double new_point, double xi, double yj, double c) {
		if (((xi <= new_point) && (new_point <= yj)) || ((yj <= new_point) && (new_point <= xi))) {
			return c;
		} else {
			return c + min(abs(new_point - xi), abs(new_point - yj));
		}
	}

	static double distance(double[] lines, double[] cols, double co, 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;

		// 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;

		// --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- ---
		// 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 cj = cols[nbcols - 1];
			double cj1 = cols[nbcols - 2];
			double la = min(abs(li - cj), // Diag: Move
					split_merge_cost(cj, li, cj1, co), // Previous: Split/Merge
					split_merge_cost(li, li1, cj, co) // Above: Split/Merge
			);
			ub = (cutoff + utils.EPSILON) - la;
		}

		// --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- ---
		// Initialisation: compute the first line. Required as the main loop starts at line=1, not 0.
		{
			double l0 = lines[0];
			// First cell (0,0) is a special case. Early abandon if above the cut-off point.
			{
				cost = abs(l0 - cols[0]); // Very first cell
				buffers[c + 0] = cost;
				if (cost <= ub) {
					prev_pp = 1;
				} else {
					return POSITIVE_INFINITY;
				}
			}
			// Rest of the line, a cell only depends on the previous cell. Stop when > ub, update prev_pp.
			for (j = 1; j < nbcols; ++j) {
				cost = cost + split_merge_cost(cols[j], l0, cols[j - 1], co);
				if (cost <= ub) {
					buffers[c + j] = cost;
					prev_pp = j + 1;
				} else {
					break;
				}
			}
			// Next line.
			++i;
		}

		// --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- ---
		// Main loop
		for (; i < nblines; ++i) {
			// --- --- --- Swap and variables init
			int swap = c;
			c = p;
			p = swap;
			double li = lines[i];
			double li1 = lines[i - 1];
			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] + split_merge_cost(li, li1, cols[j], co);
				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) {
				double cj = cols[j];
				cost = min(buffers[p + j - 1] + abs(li - cj), // Diag: Move
						buffers[p + j] + split_merge_cost(li, li1, cj, co) // Above: Split/Merge
				);
				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) {
				double cj = cols[j];
				cost = min(buffers[p + j - 1] + abs(li - cj), // Diag: Move
						cost + split_merge_cost(cj, li, cols[j - 1], co), // Previous: Split/Merge
						buffers[p + j] + split_merge_cost(li, li1, cj, co) // Above: Split/Merge
				);
				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.
				double cj = cols[j];
				if (j == next_start) { // Case 1: Advancing next start: only diag.
					cost = buffers[p + j - 1] + abs(li - cj); // Diag: Move
					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(buffers[p + j - 1] + abs(li - cj), // Diag: Move
							cost + split_merge_cost(cj, li, cols[j - 1], co) // Previous: Split/Merge
					);
					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 + split_merge_cost(cols[j], li, cols[j - 1], co); // Previous: Split/Merge
				buffers[c + j] = cost;
				if (cost <= ub) {
					++curr_pp;
				}
			}
			// --- --- ---
			prev_pp = curr_pp;
		} // End of main loop for(;i<nblines;++i)

		// --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- ---
		// Finalisation
		// 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;
		}

	}

}
