code: package dotnet54.measures.multivariate;
  
  import dotnet54.measures.Measure;
  import dotnet54.tscore.data.Dataset;
  
  import java.util.Random;
  
  /**
   * Some classes in this package may contain borrowed code from the timeseriesweka project (Bagnall, 2017), 
   * we might have modified (fixes, and improvements for efficiency) the original classes.
   * 
   */
  
  public class TWE extends MultivariateMeasure {
  	public static final String stiffnessParam = "nu";
  	public static final String costParam = "lambda";
  	protected double nu;
  	protected double lambda;
  
  	public TWE(Measure measureName, boolean dependentDimensions, int[] dimensionsToUse, double nu, double lambda) {
  		super(measureName, dependentDimensions, dimensionsToUse);
  		this.nu = nu;
  		this.lambda = lambda;
  	}
  
  	public double distance(double[][] series1, double[][] series2, double bsf, double nu, double lambda) throws Exception {
  		this.setNu(nu);
  		this.setLambda(lambda);
  		return distance(series1, series2, bsf);
  	}
  
  	@Override
  	public double distanceIndep(double[][] series1, double[][] series2, int dimension, double bsf) {
  		double[] vector1 = series1[dimension];
  		double[] vector2 = series2[dimension];
  		int m = vector1.length;
  		int n = vector2.length;
  		int maxLength = Math.max(m, n);
  		double dist, disti1, distj1;
  
  		int r = vector1.length; // this is just m?!
  		int c = vector2.length; // so is this, but surely it should actually
  				   // be n anyway
  		
  		int i, j;
  		/*
  		 * allocations in c double **D = (double **)calloc(r+1,
  		 * sizeof(double*)); double *Di1 = (double *)calloc(r+1,
  		 * sizeof(double)); double *Dj1 = (double *)calloc(c+1,
  		 * sizeof(double)); for(i=0; i<=r; i++) { D[i]=(double
  		 * *)calloc(c+1, sizeof(double)); }
  		 */
  
  		double[][]D = MemorySpaceProvider.getInstance(maxLength+1).getDoubleMatrix();
  		double[]Di1 = MemorySpaceProvider.getInstance(maxLength+1).getDoubleArray();
  		double[]Dj1 = MemorySpaceProvider.getInstance(maxLength+1).getDoubleArray();
  //		double[][] D = MemoryManager.getInstance().getDoubleMatrix(0);
  //		double[] Di1 = MemoryManager.getInstance().getDoubleArray(0);
  //		double[] Dj1 = MemoryManager.getInstance().getDoubleArray(1);
  
  		// FPH adding initialisation given that using matrices as fields
  		Di1[0] = 0.0;
  		Dj1[0] = 0.0;
  		// local costs initializations
  		for (j = 1; j <= c; j++) {
  			distj1 = 0;
  			if (j > 1) {
  				// CHANGE AJB 8/1/16: Only use power of
  				// 2 for speed
  				distj1 += (vector2[j - 2] - vector2[j - 1]) * (vector2[j - 2] - vector2[j - 1]);
  				// OLD VERSION
  				// distj1+=Math.pow(Math.abs(series2[j-2][k]-series2[j-1][k]),degree);
  				// in c:
  				// distj1+=pow(fabs(series2[j-2][k]-series2[j-1][k]),degree);
  			} else {
  				distj1 += vector2[j - 1] * vector2[j - 1];
  			}
  			// OLD distj1+=Math.pow(Math.abs(series2[j-1][k]),degree);
  			Dj1[j] = (distj1);
  		}
  
  		for (i = 1; i <= r; i++) {
  			disti1 = 0;
  			if (i > 1) {
  				disti1 += (vector1[i - 2] - vector1[i - 1]) * (vector1[i - 2] - vector1[i - 1]);
  			} // OLD
  			  // disti1+=Math.pow(Math.abs(series1[i-2][k]-series1[i-1][k]),degree);
  			else {
  				disti1 += (vector1[i - 1]) * (vector1[i - 1]);
  			}
  			// OLD disti1+=Math.pow(Math.abs(series1[i-1][k]),degree);
  
  			Di1[i] = (disti1);
  
  			for (j = 1; j <= c; j++) {
  				dist = 0;
  				dist += (vector1[i - 1] - vector2[j - 1]) * (vector1[i - 1] - vector2[j - 1]);
  				// dist+=Math.pow(Math.abs(series1[i-1][k]-series2[j-1][k]),degree);
  				if (i > 1 && j > 1) {
  					dist += (vector1[i - 2] - vector2[j - 2]) * (vector1[i - 2] - vector2[j - 2]);
  				}
  				// dist+=Math.pow(Math.abs(series1[i-2][k]-series2[j-2][k]),degree);
  				D[i][j] = (dist);
  			}
  		} // for i
  
  		// border of the cost matrix initialization
  		D[0][0] = 0;
  		for (i = 1; i <= r; i++) {
  			D[i][0] = Double.POSITIVE_INFINITY;
  		}
  		for (j = 1; j <= c; j++) {
  			D[0][j] = Double.POSITIVE_INFINITY;
  		}
  
  		double dmin, htrans, dist0;
  
  		for (i = 1; i <= r; i++) {
  			for (j = 1; j <= c; j++) {
  				htrans = Math.abs(i- j);
  				if (j > 1 && i > 1) {
  					htrans += Math.abs((i-1) - (j-1));
  				}
  				dist0 = D[i - 1][j - 1] + nu * htrans + D[i][j];
  				dmin = dist0;
  				if (i > 1) {
  					htrans = 1;
  				} else {
  					htrans = i;
  				}
  				dist = Di1[i] + D[i - 1][j] + lambda + nu * htrans;
  				if (dmin > dist) {
  					dmin = dist;
  				}
  				if (j > 1) {
  					htrans = 1;
  				} else {
  					htrans = j;
  				}
  				dist = Dj1[j] + D[i][j - 1] + lambda + nu * htrans;
  				if (dmin > dist) {
  					dmin = dist;
  				}
  				D[i][j] = dmin;
  			}
  		}
  
  		dist = D[r][c];
  		MemorySpaceProvider.getInstance().returnDoubleMatrix(D);
  		MemorySpaceProvider.getInstance().returnDoubleArray(Di1);
  		MemorySpaceProvider.getInstance().returnDoubleArray(Dj1);
  		return dist;
  	}
  
  	@Override
  	public double distanceDep(double[][] series1, double[][] series2, double bsf){
  		int m = series1[0].length;
  		int n = series2[0].length;
  		int maxLength = Math.max(m, n);
  		double dist, disti1, distj1;
  
  		int r = series1[0].length; // this is just m?!
  		int c = series2[0].length; // so is this, but surely it should actually
  		// be n anyway
  
  		int i, j;
  		/*
  		 * allocations in c double **D = (double **)calloc(r+1,
  		 * sizeof(double*)); double *Di1 = (double *)calloc(r+1,
  		 * sizeof(double)); double *Dj1 = (double *)calloc(c+1,
  		 * sizeof(double)); for(i=0; i<=r; i++) { D[i]=(double
  		 * *)calloc(c+1, sizeof(double)); }
  		 */
  
  		double[][]D = MemorySpaceProvider.getInstance(maxLength+1).getDoubleMatrix();
  		double[]Di1 = MemorySpaceProvider.getInstance(maxLength+1).getDoubleArray();
  		double[]Dj1 = MemorySpaceProvider.getInstance(maxLength+1).getDoubleArray();
  //		double[][] D = MemoryManager.getInstance().getDoubleMatrix(0);
  //		double[] Di1 = MemoryManager.getInstance().getDoubleArray(0);
  //		double[] Dj1 = MemoryManager.getInstance().getDoubleArray(1);
  
  		// FPH adding initialisation given that using matrices as fields
  		Di1[0] = 0.0;
  		Dj1[0] = 0.0;
  		// local costs initializations
  		for (j = 1; j <= c; j++) {
  			distj1 = 0;
  			if (j > 1) {
  				// CHANGE AJB 8/1/16: Only use power of
  				// 2 for speed
  				distj1 += adjustedSquaredEuclideanDistanceFromPointToPoint(series2, j-2, series2, j-1);
  //				distj1 += (series2[j - 2] - series2[j - 1]) * (series2[j - 2] - series2[j - 1]);
  				// OLD VERSION
  				// distj1+=Math.pow(Math.abs(series2[j-2][k]-series2[j-1][k]),degree);
  				// in c:
  				// distj1+=pow(fabs(series2[j-2][k]-series2[j-1][k]),degree);
  			} else {
  				distj1 += adjustedSquaredEuclideanNormOfPoint(series2, j-1);
  //				distj1 += series2[j - 1] * series2[j - 1];
  			}
  			// OLD distj1+=Math.pow(Math.abs(series2[j-1][k]),degree);
  			Dj1[j] = (distj1);
  		}
  
  		for (i = 1; i <= r; i++) {
  			disti1 = 0;
  			if (i > 1) {
  				disti1 += adjustedSquaredEuclideanDistanceFromPointToPoint(series1, i-2, series1, i-1);
  //				disti1 += (series1[i - 2] - series1[i - 1]) * (series1[i - 2] - series1[i - 1]);
  			} // OLD
  			// disti1+=Math.pow(Math.abs(series1[i-2][k]-series1[i-1][k]),degree);
  			else {
  				disti1 += adjustedSquaredEuclideanNormOfPoint(series1, i-1);
  //				disti1 += (series1[i - 1]) * (series1[i - 1]);
  			}
  			// OLD disti1+=Math.pow(Math.abs(series1[i-1][k]),degree);
  
  			Di1[i] = (disti1);
  
  			for (j = 1; j <= c; j++) {
  				dist = 0;
  				dist += adjustedSquaredEuclideanDistanceFromPointToPoint(series1, i-1, series2, j-1);
  //				dist += (series1[i - 1] - series2[j - 1]) * (series1[i - 1] - series2[j - 1]);
  				// dist+=Math.pow(Math.abs(series1[i-1][k]-series2[j-1][k]),degree);
  				if (i > 1 && j > 1) {
  					dist +=  adjustedSquaredEuclideanDistanceFromPointToPoint(series1, i-2, series2, j-2);
  //					dist += (series1[i - 2] - series2[j - 2]) * (series1[i - 2] - series2[j - 2]);
  				}
  				// dist+=Math.pow(Math.abs(series1[i-2][k]-series2[j-2][k]),degree);
  				D[i][j] = (dist);
  			}
  		} // for i
  
  //		// border of the cost matrix initialization
  //		D[0][0] = 0;
  //		for (i = 1; i <= r; i++) {
  //			D[i][0] = D[i - 1][0] + Di1[i];
  //		}
  //		for (j = 1; j <= c; j++) {
  //			D[0][j] = D[0][j - 1] + Dj1[j];
  //		}
  		//note fix by Mat
  //		https://github.com/uea-machine-learning/tsml/commit/bd0b2af1bbfdc3bd9ff078b468882ea493270185#diff-c7dc635e73d5e2013a87938ba4eb5a22
  		D[0][0] = 0;
  		for (i = 1; i <= r; i++) {
  			D[i][0] = Double.POSITIVE_INFINITY;
  		}
  		for (j = 1; j <= c; j++) {
  			D[0][j] = Double.POSITIVE_INFINITY;
  		}
  
  		double dmin, htrans, dist0;
  
  		for (i = 1; i <= r; i++) {
  			for (j = 1; j <= c; j++) {
  				htrans = Math.abs(i- j);
  				if (j > 1 && i > 1) {
  					htrans += Math.abs((i-1) - (j-1));
  				}
  				dist0 = D[i - 1][j - 1] + nu * htrans + D[i][j];
  				dmin = dist0;
  				if (i > 1) {
  					htrans = 1;
  				} else {
  					htrans = i;
  				}
  				dist = Di1[i] + D[i - 1][j] + lambda + nu * htrans;
  				if (dmin > dist) {
  					dmin = dist;
  				}
  				if (j > 1) {
  					htrans = 1;
  				} else {
  					htrans = j;
  				}
  				dist = Dj1[j] + D[i][j - 1] + lambda + nu * htrans;
  				if (dmin > dist) {
  					dmin = dist;
  				}
  				D[i][j] = dmin;
  			}
  		}
  
  		dist = D[r][c];
  		MemorySpaceProvider.getInstance().returnDoubleMatrix(D);
  		MemorySpaceProvider.getInstance().returnDoubleArray(Di1);
  		MemorySpaceProvider.getInstance().returnDoubleArray(Dj1);
  		return dist;
  	}
  
  	// Unused - TWE_D now uses the adjusted method
  	protected final double squaredEuclideanNormOfPoint(double[][] series1, int i) {
  		double total = 0;
  		double element;
  		for (int dimension : dimensionsToUse) {
  			element = series1[dimension][i];
  			total += element * element;
  		}
  		return total;
  	}
  
  	// Used by TWE_D
  	protected final double adjustedSquaredEuclideanNormOfPoint(double[][] series1, int i) {
  		double total = 0;
  		double element;
  		for (int dimension : dimensionsToUse) {
  			element = series1[dimension][i];
  			total += element * element;
  		}
  		return total / (2 * dimensionsToUse.length);
  	}
  
  
  	// Used by TWE_D
  	protected final double adjustedSquaredEuclideanDistanceFromPointToPoint(double[][] series1, int i,
  																			double[][] series2, int j) {
  		double total = 0;
  		double diff;
  		for (int dimension : dimensionsToUse) {
  			diff = series1[dimension][i] - series2[dimension][j];
  			total += diff * diff;
  		}
  		return total / (2 * dimensionsToUse.length);
  	}
  
  	@Override
  	public void setRandomParams(Dataset trainData, Random rand){
  		setNu(getRandomNu(trainData, rand));
  		setLambda(getRandomLambda(trainData, rand));
  	}
  
  	@Override
  	public void initParamsByID(Dataset trainData, Random rand) {
  		//pass
  	}
  
  	@Override
  	public void setParamsByID(int paramID, int seriesLength, Random rand) {
  		setNu(twe_nuParams[paramID / 10]);
  		setLambda(twe_lamdaParams[paramID % 10]);
  	}
  
  	@Override
  	public void setRandomParamByID(Dataset trainData, Random rand) {
  		int paramID = rand.nextInt(100 );
  		setParamsByID(paramID, trainData.length(), rand);
  	}
  
  	public static double getRandomNu(Dataset dataset, Random rand) {
  		return twe_nuParams[rand.nextInt(twe_nuParams.length)];
  	} 	
  	
  	public static double getRandomLambda(Dataset dataset, Random rand) {
  		return twe_lamdaParams[rand.nextInt(twe_lamdaParams.length)];
  	}
  
  	public double getNu() {
  		return nu;
  	}
  
  	public void setNu(double nu) {
  		this.nu = nu;
  	}
  
  	public double getLambda() {
  		return lambda;
  	}
  
  	public void setLambda(double lambda) {
  		this.lambda = lambda;
  	}
  
  	@Override
  	public String toString(){
  		return "TWE[,n="+nu+",l="+lambda+","
  				+super.toString()+",]";
  	}
  
  	public static final double[] twe_nuParams = { 0.00001, 0.0001, 0.0005, 0.001, 0.005, 0.01, 0.05, 0.1, 0.5, 1 };
  
  	public static final double[] twe_lamdaParams = { 0, 0.011111111, 0.022222222, 0.033333333, 0.044444444, 0.055555556, 0.066666667,
  			0.077777778, 0.088888889, 0.1 };
  
  
  //	shifaz old code backup 18/3/2020
  //	public synchronized double distance_indep(double[] ta, double[] tb, double bsf, double nu, double lambda) {
  //
  //		int m = ta.length;
  //		int n = tb.length;
  //		int maxLength = Math.max(m, n);
  //
  //		double dist, disti1, distj1;
  //
  //		int r = ta.length; // this is just m?!
  //		int c = tb.length; // so is this, but surely it should actually
  //		// be n anyway
  //
  //		int i, j;
  //		/*
  //		 * allocations in c double **D = (double **)calloc(r+1,
  //		 * sizeof(double*)); double *Di1 = (double *)calloc(r+1,
  //		 * sizeof(double)); double *Dj1 = (double *)calloc(c+1,
  //		 * sizeof(double)); for(i=0; i<=r; i++) { D[i]=(double
  //		 * *)calloc(c+1, sizeof(double)); }
  //		 */
  //
  //		double[][]D = MemorySpaceProvider.getInstance(maxLength+1).getDoubleMatrix();
  //		double[]Di1 = MemorySpaceProvider.getInstance(maxLength+1).getDoubleArray();
  //		double[]Dj1 = MemorySpaceProvider.getInstance(maxLength+1).getDoubleArray();
  ////		double[][] D = MemoryManager.getInstance().getDoubleMatrix(0);
  ////		double[] Di1 = MemoryManager.getInstance().getDoubleArray(0);
  ////		double[] Dj1 = MemoryManager.getInstance().getDoubleArray(1);
  //
  //		// FPH adding initialisation given that using matrices as fields
  //		Di1[0] = 0.0;
  //		Dj1[0] = 0.0;
  //		// local costs initializations
  //		for (j = 1; j <= c; j++) {
  //			distj1 = 0;
  //			if (j > 1) {
  //				// CHANGE AJB 8/1/16: Only use power of
  //				// 2 for speed
  //				distj1 += (tb[j - 2] - tb[j - 1]) * (tb[j - 2] - tb[j - 1]);
  //				// OLD VERSION
  //				// distj1+=Math.pow(Math.abs(tb[j-2][k]-tb[j-1][k]),degree);
  //				// in c:
  //				// distj1+=pow(fabs(tb[j-2][k]-tb[j-1][k]),degree);
  //			} else {
  //				distj1 += tb[j - 1] * tb[j - 1];
  //			}
  //			// OLD distj1+=Math.pow(Math.abs(tb[j-1][k]),degree);
  //			Dj1[j] = (distj1);
  //		}
  //
  //		for (i = 1; i <= r; i++) {
  //			disti1 = 0;
  //			if (i > 1) {
  //				disti1 += (ta[i - 2] - ta[i - 1]) * (ta[i - 2] - ta[i - 1]);
  //			} // OLD
  //			// disti1+=Math.pow(Math.abs(ta[i-2][k]-ta[i-1][k]),degree);
  //			else {
  //				disti1 += (ta[i - 1]) * (ta[i - 1]);
  //			}
  //			// OLD disti1+=Math.pow(Math.abs(ta[i-1][k]),degree);
  //
  //			Di1[i] = (disti1);
  //
  //			for (j = 1; j <= c; j++) {
  //				dist = 0;
  //				dist += (ta[i - 1] - tb[j - 1]) * (ta[i - 1] - tb[j - 1]);
  //				// dist+=Math.pow(Math.abs(ta[i-1][k]-tb[j-1][k]),degree);
  //				if (i > 1 && j > 1) {
  //					dist += (ta[i - 2] - tb[j - 2]) * (ta[i - 2] - tb[j - 2]);
  //				}
  //				// dist+=Math.pow(Math.abs(ta[i-2][k]-tb[j-2][k]),degree);
  //				D[i][j] = (dist);
  //			}
  //		} // for i
  //
  //		// border of the cost matrix initialization
  //		D[0][0] = 0;
  //		for (i = 1; i <= r; i++) {
  //			D[i][0] = D[i - 1][0] + Di1[i];
  //		}
  //		for (j = 1; j <= c; j++) {
  //			D[0][j] = D[0][j - 1] + Dj1[j];
  //		}
  //
  //		double dmin, htrans, dist0;
  //
  //		for (i = 1; i <= r; i++) {
  //			for (j = 1; j <= c; j++) {
  //				htrans = Math.abs(i- j);
  //				if (j > 1 && i > 1) {
  //					htrans += Math.abs((i-1) - (j-1));
  //				}
  //				dist0 = D[i - 1][j - 1] + nu * htrans + D[i][j];
  //				dmin = dist0;
  //				if (i > 1) {
  //					htrans = 1;
  //				} else {
  //					htrans = i;
  //				}
  //				dist = Di1[i] + D[i - 1][j] + lambda + nu * htrans;
  //				if (dmin > dist) {
  //					dmin = dist;
  //				}
  //				if (j > 1) {
  //					htrans = 1;
  //				} else {
  //					htrans = j;
  //				}
  //				dist = Dj1[j] + D[i][j - 1] + lambda + nu * htrans;
  //				if (dmin > dist) {
  //					dmin = dist;
  //				}
  //				D[i][j] = dmin;
  //			}
  //		}
  //
  //		dist = D[r][c];
  //		MemorySpaceProvider.getInstance().returnDoubleMatrix(D);
  //		MemorySpaceProvider.getInstance().returnDoubleArray(Di1);
  //		MemorySpaceProvider.getInstance().returnDoubleArray(Dj1);
  //		return dist;
  //	}
  
  	
  }
  