code: package monashts.tempo.univariate.distances;
  
  import static monashts.tempo.univariate.distances.utils.dist;
  import static monashts.tempo.univariate.distances.utils.min;
  
  import java.util.Arrays;
  
  import static java.lang.Double.POSITIVE_INFINITY;
  import static java.lang.Double.min;
  import static java.lang.Integer.max;
  import static java.lang.Integer.min;
  
  public class ERP {
  	
  	static double distance(double[] lines, double[] cols, double gValue, int w, double cutoff) {
  		
  		// Ensure that lines are longer than columns
  		if(lines.length < cols.length) {
  			double [] swap = lines;
  			lines = cols;
  			cols = swap;
  		}
  		
          // Cap the windows and check that, given the constraint, an alignment is possible
          if (w > lines.length) { w = lines.length; }
          if (lines.length - cols.length > w) { return POSITIVE_INFINITY; }
  		
  		// --- --- --- Declarations
  		int nblines = lines.length;
  		int nbcols = cols.length;
  				
  		// Setup buffers - get an extra cell for border condition. Init to +INF.
  		double[] buffers = new double[(1+nbcols)*2];
  		Arrays.fill(buffers, POSITIVE_INFINITY);
  		int c = 0+1; 		// Start of current line in buffer - account for the extra cell
  		int p = nbcols+2;   // Start of previous line in buffer - account for two extra cells
  		
  		// 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 (most commonly used in the code).
          // First, take the "next float" after "cutoff" to deal with numerical instability.
          // Then, subtract the cost of the last alignment.
  		double la = min(
                      dist(gValue, cols[nbcols - 1]),             // Previous
                      dist(lines[nblines - 1], cols[nbcols - 1]), // Diagonal
                      dist(lines[nblines - 1], gValue)            // Above
              );
  		double ub = cutoff + utils.EPSILON - la;
  
  
          // --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- ---
          // Initialisation of the top border
          {   // Matrix Border - Top diagonal
              buffers[c - 1] = 0;
              // Matrix Border - First line
              int jStop =  min(i+w+1, nbcols);
              for (j = 0; buffers[c + j - 1] <= ub && j < jStop; ++j) {
                  buffers[c + j] = buffers[c + j - 1] + dist(gValue, cols[j]);
              }
              // Pruning point set to first +INF value (or out of bound)
              prev_pp = j;
          }
  
          // --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- ---
          // Part 1: Loop with computed left border.
          {   // The left border has a computed value while it's within the window and its value bv <= ub
              // No "front pruning" (next_start) and no early abandoning can occur while in this loop.
              int iStop =  min(i+w+1, nblines);
              for (; i < iStop; ++i) {
                  // --- --- --- Variables init
                  double li = lines[i];
                  int jStart = 0;
                  int jStop = min(i+w+1, nbcols);
                  j = jStart;
                  int curr_pp = jStart; // Next pruning point init at the start of the line
                  // --- --- --- Stage 0: Initialise the left border
                  {
                      // We haven't swap yet, so the 'top' cell is still indexed by 'c-1'.
                      cost = buffers[c - 1] + dist(li, gValue);
                      if (cost > ub) { break; }
                      else {
                      	int swap = c;
                      	c = p;
                      	p = swap;
                          buffers[c - 1] = cost;
                      }
                  }
                  // --- --- --- Stage 1: Up to the previous pruning point while advancing next_start: diag and top
                  // No stage 1 here.
                  // --- --- --- Stage 2: Up to the previous pruning point without advancing next_start: left, diag and top
                  for (; j < prev_pp; ++j) {
                      cost = min(
                              cost + dist(gValue, cols[j]),               // Previous
                              buffers[p + j - 1] + dist(li, cols[j]),     // Diagonal
                              buffers[p + j] + dist(li, gValue)           // Above
                      );
                      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 < jStop) { // Possible path in previous cells: left and diag.
                      cost = min(
                              cost + dist(gValue, cols[j]),               // Previous
                              buffers[p + j - 1] + dist(li, cols[j])      // Diagonal
                      );
                      buffers[c + j] = cost;
                      if (cost <= ub) { curr_pp = j + 1; }
                      ++j;
                  }
                  // --- --- --- 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 < jStop; ++j) {
                      cost = cost + dist(gValue, cols[j]);  // Previous
                      buffers[c + j] = cost;
                      if (cost <= ub) { ++curr_pp; }
                  }
                  // --- --- ---
                  prev_pp = curr_pp;
              }
          }
  
          // --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- ---
          // Part 2: Loop with +INF left border
          {
              for (; i < nblines; ++i) {
                  // --- --- --- Swap and variables init
              	int swap = c;
              	c = p;
              	p = swap;
                  double li = lines[i];
                  int jStart = max(i-w, next_start);
                  int jStop = min(i+w+1, nbcols); 
                  j = jStart;
                  next_start = jStart;
                  int curr_pp = jStart; // Next pruning point init at the start of the line
                  // --- --- --- Stage 0: Initialise the left border
                  {
                      cost = POSITIVE_INFINITY;
                      buffers[c + jStart - 1] = cost;
                  }
                  // --- --- --- 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]),     // Diagonal
                              buffers[p + j] + dist(li, gValue)           // Above
                      );
                      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 + dist(gValue, cols[j]),               // Previous
                              buffers[p + j - 1] + dist(li, cols[j]),     // Diagonal
                              buffers[p + j] + dist(li, gValue)           // Above
                      );
                      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 < jStop) { // If so, two cases.
                      if (j == next_start) { // Case 1: Advancing next start: only diag.
                          cost = buffers[p + j - 1] + dist(li, cols[j]);     // Diagonal
                          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 + dist(gValue, cols[j]),               // Previous
                                  buffers[p + j - 1] + dist(li, cols[j])      // Diagonal
                          );
                          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 < jStop; ++j) {
                      cost = cost + dist(gValue, cols[j]);
                      buffers[c + j] = cost;
                      if (cost <= ub) { ++curr_pp; }
                  }
                  // --- --- ---
                  prev_pp = curr_pp;
              }
          }
  
          // --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- ---
          // 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; }
      }
  
  }
  