package mil.nga.giat.geowave.analytic.mapreduce.kde;

import java.util.ArrayList;
import java.util.List;

import org.apache.log4j.Logger;

public class GaussianFilter
{
	private final static Logger LOGGER = Logger.getLogger(GaussianFilter.class);
	private static final double SQRT_2_PI = Math.sqrt(2 * Math.PI);

	/**
	 * This kernel was computed with sigma = 1 for x=(-3,-2,-1,0,1,2,3)
	 */
	private static double[] majorSmoothingGaussianKernel = new double[] {
		0.006,
		0.061,
		0.242,
		0.383,
		0.242,
		0.061,
		0.006
	};

	// private static double[] minorSmoothingGaussianKernel = new double[] {
	// 0.2186801,
	// 0.531923041,
	// 0.2186801
	// };

	private static class ValueRange
	{
		private final double min;
		private final double max;

		private ValueRange(
				final double min,
				final double max ) {
			this.min = min;
			this.max = max;
		}

		public double getMin() {
			return min;
		}

		public double getMax() {
			return max;
		}
	}

	private static final ValueRange[] valueRangePerDimension = new ValueRange[] {
		new ValueRange(
				-180,
				180),
		new ValueRange(
				-90,
				90)
	};

	public static void incrementPt(
			final double lat,
			final double lon,
			final CellCounter results,
			final int numXPosts,
			final int numYPosts ) {
		incrementBBox(
				lon,
				lon,
				lat,
				lat,
				results,
				numXPosts,
				numYPosts,
				1);
	}

	public static void incrementPt(
			final double lat,
			final double lon,
			final CellCounter results,
			final int numXPosts,
			final int numYPosts,
			double contributionScaleFactor ) {
		incrementBBox(
				lon,
				lon,
				lat,
				lat,
				results,
				numXPosts,
				numYPosts,
				contributionScaleFactor);
	}

	public static void incrementPtFast(
			final double lat,
			final double lon,
			final CellCounter results,
			final int numXPosts,
			final int numYPosts ) {
		final int numDimensions = 2;
		final double[] binLocationPerDimension = new double[numDimensions];
		final int[] binsPerDimension = new int[] {
			numXPosts,
			numYPosts
		};
		final double[] valsPerDimension = new double[] {
			lon,
			lat
		};
		for (int d = 0; d < numDimensions; d++) {
			final ValueRange valueRange = valueRangePerDimension[d];
			final double span = (valueRange.getMax() - valueRange.getMin());
			binLocationPerDimension[d] = (((valsPerDimension[d] - valueRange.getMin()) / span) * binsPerDimension[d]);
		}
		final double[] gaussianKernel = getGaussianKernel(
				1,
				3);
		final int maxOffset = gaussianKernel.length / 2;
		final List<int[]> offsets = getOffsets(
				numDimensions,
				0,
				new int[numDimensions],
				gaussianKernel,
				maxOffset);
		for (final int[] offset : offsets) {
			final double blur = getBlurFromOffset(
					offset,
					gaussianKernel,
					maxOffset);
			final List<BinPositionAndContribution> positionsAndContributions = getPositionsAndContributionPt(
					numDimensions,
					0,
					binLocationPerDimension,
					blur,
					new int[numDimensions],
					binsPerDimension,
					offset);
			for (final BinPositionAndContribution positionAndContribution : positionsAndContributions) {
				results.increment(
						positionAndContribution.position,
						positionAndContribution.contribution);
			}
		}
	}

	public static void incrementBBox(
			final double minX,
			final double maxX,
			final double minY,
			final double maxY,
			final CellCounter results,
			final int numXPosts,
			final int numYPosts,
			double contributionScaleFactor ) {
		final int numDimensions = 2;
		final double[] minBinLocationPerDimension = new double[numDimensions];
		final double[] maxBinLocationPerDimension = new double[numDimensions];
		final int[] binsPerDimension = new int[] {
			numXPosts,
			numYPosts
		};
		final ValueRange[] valueRangePerDimension = new ValueRange[] {
			new ValueRange(
					-180,
					180),
			new ValueRange(
					-90,
					90)
		};
		final double[] minsPerDimension = new double[] {
			minX,
			minY
		};
		final double[] maxesPerDimension = new double[] {
			maxX,
			maxY
		};
		for (int d = 0; d < numDimensions; d++) {
			final ValueRange valueRange = valueRangePerDimension[d];
			final double span = (valueRange.getMax() - valueRange.getMin());
			minBinLocationPerDimension[d] = (((minsPerDimension[d] - valueRange.getMin()) / span) * binsPerDimension[d]);
			maxBinLocationPerDimension[d] = (((maxesPerDimension[d] - valueRange.getMin()) / span) * binsPerDimension[d]);
			// give it a buffer of 1 for being counted within this bounds
			// because we perform smoothing on the values anyway
			if ((maxBinLocationPerDimension[d] < -1) || (minBinLocationPerDimension[d] > binsPerDimension[d])) {
				// not in bounds
				return;
			}
			else {
				minBinLocationPerDimension[d] = Math.max(
						minBinLocationPerDimension[d],
						-1);
				maxBinLocationPerDimension[d] = Math.min(
						maxBinLocationPerDimension[d],
						binsPerDimension[d]);
			}

		}
		final double[] gaussianKernel = getGaussianKernel(
				1,
				3);
		final int maxOffset = gaussianKernel.length / 2;
		final List<int[]> offsets = getOffsets(
				numDimensions,
				0,
				new int[numDimensions],
				gaussianKernel,
				maxOffset);
		for (final int[] offset : offsets) {
			final double blur = getBlurFromOffset(
					offset,
					gaussianKernel,
					maxOffset);
			final List<BinPositionAndContribution> positionsAndContributions = getPositionsAndContribution(
					numDimensions,
					0,
					minBinLocationPerDimension,
					maxBinLocationPerDimension,
					blur,
					new int[numDimensions],
					binsPerDimension,
					offset);
			for (final BinPositionAndContribution positionAndContribution : positionsAndContributions) {
				results.increment(
						positionAndContribution.position,
						positionAndContribution.contribution * contributionScaleFactor);
			}
		}
	}

	protected static double getSigma(
			final int radius,
			final int order ) {
		return ((radius * 2.0) + 1.0) / (5.0 + (0.8 * order));
	}

	protected static double[] getGaussianKernel(
			final double sigma,
			final int radius ) {
		return majorSmoothingGaussianKernel;
		// final double[] kernel = new double[(radius * 2) + 1];
		// int index = 0;
		// for (int i = radius; i >= -radius; i--) {
		// kernel[index++] = computePDF(
		// 0,
		// sigma,
		// i);
		// }
		// return normalizeSumToOne(kernel);
	}

	protected static double computePDF(
			final double mean,
			final double sigma,
			final double sample ) {
		final double delta = sample - mean;
		return Math.exp((-delta * delta) / (2.0 * sigma * sigma)) / (sigma * SQRT_2_PI);
	}

	protected static double[] normalizeSumToOne(
			final double[] kernel ) {
		final double[] retVal = new double[kernel.length];
		double total = 0;
		for (final double element : kernel) {
			total += element;
		}
		for (int i = 0; i < kernel.length; i++) {
			retVal[i] = kernel[i] / total;
		}
		return retVal;
	}

	static private List<int[]> getOffsets(
			final int numDimensions,
			final int currentDimension,
			final int[] currentOffsetsPerDimension,
			final double[] gaussianKernel,
			final int maxOffset ) {
		final List<int[]> offsets = new ArrayList<int[]>();
		if (currentDimension == numDimensions) {
			offsets.add(currentOffsetsPerDimension.clone());
		}
		else {
			for (int i = -maxOffset; i < (gaussianKernel.length - maxOffset); i++) {
				currentOffsetsPerDimension[currentDimension] = i;
				offsets.addAll(getOffsets(
						numDimensions,
						currentDimension + 1,
						currentOffsetsPerDimension,
						gaussianKernel,
						maxOffset));
			}
		}
		return offsets;
	}

	static private double getBlurFromOffset(
			final int[] indexIntoGaussianPerDimension,
			final double[] gaussianKernel,
			final int maxOffset ) {
		double blurFactor = 1;

		for (final int index : indexIntoGaussianPerDimension) {
			blurFactor *= gaussianKernel[index + maxOffset];
		}
		return blurFactor;
	}

	private static List<BinPositionAndContribution> getPositionsAndContributionPt(
			final int numDimensions,
			final int currentDimension,
			final double[] locationPerDimension,
			final double currentContribution,
			final int[] finalIndexPerDimension,
			final int[] binsPerDimension,
			final int[] offset ) {
		final List<BinPositionAndContribution> positions = new ArrayList<BinPositionAndContribution>();
		if (currentDimension == numDimensions) {
			positions.add(new BinPositionAndContribution(
					getPosition(
							finalIndexPerDimension,
							binsPerDimension),
					currentContribution));
		}
		else {
			final int floorOfLocation = (int) (locationPerDimension[currentDimension]);
			final int[] floorLocation = finalIndexPerDimension;
			floorLocation[currentDimension] = floorOfLocation + offset[currentDimension];
			if ((floorLocation[currentDimension] >= 0)
					&& (floorLocation[currentDimension] < binsPerDimension[currentDimension])) {
				positions.addAll(getPositionsAndContributionPt(
						numDimensions,
						currentDimension + 1,
						locationPerDimension,
						currentContribution,
						floorLocation,
						binsPerDimension,
						offset));
			}
		}
		return positions;
	}

	private static List<BinPositionAndContribution> getPositionsAndContribution(
			final int numDimensions,
			final int currentDimension,
			final double[] minLocationPerDimension,
			final double[] maxLocationPerDimension,
			final double currentContribution,
			final int[] finalIndexPerDimension,
			final int[] binsPerDimension,
			final int[] offset ) {
		final List<BinPositionAndContribution> positions = new ArrayList<BinPositionAndContribution>();
		if (currentDimension == numDimensions) {
			positions.add(new BinPositionAndContribution(
					getPosition(
							finalIndexPerDimension,
							binsPerDimension),
					currentContribution));
		}
		else {
			final int floorOfLocation = (int) (minLocationPerDimension[currentDimension]);
			final int[] floorLocation = finalIndexPerDimension.clone();
			floorLocation[currentDimension] = floorOfLocation + offset[currentDimension];
			if ((floorLocation[currentDimension] >= 0)
					&& (floorLocation[currentDimension] < binsPerDimension[currentDimension])) {
				positions.addAll(getPositionsAndContribution(
						numDimensions,
						currentDimension + 1,
						minLocationPerDimension,
						maxLocationPerDimension,
						currentContribution,
						floorLocation,
						binsPerDimension,
						offset));
			}
			final int ceilOfLocation = (int) Math.ceil(maxLocationPerDimension[currentDimension]);
			/**
			 * the exterior cells are covered above by the floor of the min and
			 * ceil of the max, everything in between is covered below
			 */
			final int startLocation = Math.max(
					floorOfLocation + offset[currentDimension] + 1,
					0);
			final int stopLocation = Math.min(
					ceilOfLocation + offset[currentDimension],
					binsPerDimension[currentDimension]);
			if (startLocation < stopLocation) {
				for (int location = startLocation; location < stopLocation; location++) {
					final int[] middleLocation = finalIndexPerDimension.clone();
					middleLocation[currentDimension] = location;
					positions.addAll(getPositionsAndContribution(
							numDimensions,
							currentDimension + 1,
							minLocationPerDimension,
							maxLocationPerDimension,
							currentContribution,
							middleLocation,
							binsPerDimension,
							offset));
				}
			}
		}
		return positions;
	}

	private static long getPosition(
			final int[] positionPerDimension,
			final int[] binsPerDimension ) {
		long retVal = 0;
		double multiplier = 1;
		for (int d = positionPerDimension.length - 1; d >= 0; d--) {
			retVal += (positionPerDimension[d] * multiplier);
			multiplier *= binsPerDimension[d];
		}
		return retVal;
	}

	private static class BinPositionAndContribution
	{
		final private long position;
		final private double contribution;

		private BinPositionAndContribution(
				final long position,
				final double contribution ) {
			this.position = position;
			this.contribution = contribution;
		}
	}

	/*
	 * protected void incrementCount( final double minx, final double maxx,
	 * final double miny, final double maxy, final int count ) { final double[]
	 * minsPerDimension = new double[]{minx, miny}; final double[]
	 * maxesPerDimension = new double[]{maxx,maxy};
	 * 
	 * for (final BoundsAndCounts counts : statistics.boundsWithCounts) {
	 * boolean inBounds = true; final double[] minBinLocationPerDimension = new
	 * double[2]; final double[] maxBinLocationPerDimension = new double[2]; for
	 * (int d = 0; d < 2; d++) { final ValueRange valueRange =
	 * counts.valueRangePerDimension[d]; final double span =
	 * (valueRange.getMax() - valueRange.getMin());
	 * minBinLocationPerDimension[d] = (((minsPerDimension[d] -
	 * valueRange.getMin()) / span) * counts.binsPerDimension[d]);
	 * maxBinLocationPerDimension[d] = (((maxesPerDimension[d] -
	 * valueRange.getMin()) / span) * counts.binsPerDimension[d]); // give it a
	 * buffer of 1 for being counted within this bounds // because we perform
	 * smoothing on the values anyway if ((maxBinLocationPerDimension[d] < -1)
	 * || (minBinLocationPerDimension[d] > counts.binsPerDimension[d])) {
	 * inBounds = false; break; } else { minBinLocationPerDimension[d] =
	 * Math.max( minBinLocationPerDimension[d], -1);
	 * maxBinLocationPerDimension[d] = Math.min( maxBinLocationPerDimension[d],
	 * counts.binsPerDimension[d]); }
	 * 
	 * } if (inBounds) { final double[] gaussianKernel
	 * =majorSmoothingGaussianKernel; final int maxOffset =
	 * gaussianKernel.length / 2; final List<int[]> offsets = getOffsets( 2, 0,
	 * new int[2], gaussianKernel, maxOffset); for (final int[] offset :
	 * offsets) { final double blur = getBlurFromOffset( offset, gaussianKernel,
	 * maxOffset); final List<BinPositionAndContribution>
	 * positionsAndContributions = getPositionsAndContribution( 2, 0,
	 * minBinLocationPerDimension, maxBinLocationPerDimension, blur, new int[2],
	 * counts.binsPerDimension, offset); for (final BinPositionAndContribution
	 * positionAndContribution : positionsAndContributions) {
	 * counts.incrementCount( positionAndContribution.position,
	 * positionAndContribution.contribution * count); } } } } }
	 * 
	 * static private List<int[]> getOffsets( final int numDimensions, final int
	 * currentDimension, final int[] currentOffsetsPerDimension, final double[]
	 * gaussianKernel, final int maxOffset ) { final List<int[]> offsets = new
	 * ArrayList<int[]>(); if (currentDimension == numDimensions) {
	 * offsets.add(currentOffsetsPerDimension.clone()); } else { for (int i =
	 * -maxOffset; i < (gaussianKernel.length - maxOffset); i++) {
	 * currentOffsetsPerDimension[currentDimension] = i;
	 * offsets.addAll(getOffsets( numDimensions, currentDimension + 1,
	 * currentOffsetsPerDimension, gaussianKernel, maxOffset)); } } return
	 * offsets; }
	 * 
	 * static private double getBlurFromOffset( final int[]
	 * indexIntoGaussianPerDimension, final double[] gaussianKernel, final int
	 * maxOffset ) { double blurFactor = 1;
	 * 
	 * for (final int index : indexIntoGaussianPerDimension) { blurFactor *=
	 * gaussianKernel[index + maxOffset]; } return blurFactor; }
	 * 
	 * private List<BinPositionAndContribution> getPositionsAndContribution(
	 * final int numDimensions, final int currentDimension, final double[]
	 * minLocationPerDimension, final double[] maxLocationPerDimension, final
	 * double currentContribution, final int[] finalIndexPerDimension, final
	 * int[] binsPerDimension, final int[] offset ) { final
	 * List<BinPositionAndContribution> positions = new
	 * ArrayList<BinPositionAndContribution>(); if (currentDimension ==
	 * numDimensions) { positions.add(new BinPositionAndContribution(
	 * getPosition( finalIndexPerDimension, binsPerDimension),
	 * currentContribution)); } else { final int floorOfLocation = (int)
	 * (minLocationPerDimension[currentDimension]); final int[] floorLocation =
	 * finalIndexPerDimension.clone(); floorLocation[currentDimension] =
	 * floorOfLocation + offset[currentDimension]; if
	 * ((floorLocation[currentDimension] >= 0) &&
	 * (floorLocation[currentDimension] < binsPerDimension[currentDimension])) {
	 * positions.addAll(getPositionsAndContribution( numDimensions,
	 * currentDimension + 1, minLocationPerDimension, maxLocationPerDimension,
	 * currentContribution, floorLocation, binsPerDimension, offset)); } final
	 * int ceilOfLocation = (int)
	 * Math.ceil(maxLocationPerDimension[currentDimension]);
	 *//**
	 * the exterior cells are covered above by the floor of the min and ceil
	 * of the max, everything in between is covered below
	 */
	/*
	 * final int startLocation = Math.max( floorOfLocation +
	 * offset[currentDimension] + 1, 0); final int stopLocation = Math.min(
	 * ceilOfLocation + offset[currentDimension],
	 * binsPerDimension[currentDimension]); if (startLocation < stopLocation) {
	 * for (int location = startLocation; location < stopLocation; location++) {
	 * final int[] middleLocation = finalIndexPerDimension.clone();
	 * middleLocation[currentDimension] = location;
	 * positions.addAll(getPositionsAndContribution( numDimensions,
	 * currentDimension + 1, minLocationPerDimension, maxLocationPerDimension,
	 * currentContribution, middleLocation, binsPerDimension, offset)); } } }
	 * return positions; }
	 * 
	 * private static int getPosition( final int[] positionPerDimension, final
	 * int[] binsPerDimension ) { int retVal = 0; double multiplier = 1; for
	 * (int d = 0; d < positionPerDimension.length; d++) { retVal +=
	 * (positionPerDimension[d] * multiplier); multiplier *=
	 * binsPerDimension[d]; } return retVal; }
	 * 
	 * protected static int[] getPositionPerDimension( final int position, final
	 * int[] binsPerDimension ) { int multiplier = 1;
	 * 
	 * final int[] positionPerDimension = new int[binsPerDimension.length]; for
	 * (int d = 0; d < positionPerDimension.length; d++) {
	 * positionPerDimension[d] = (position / multiplier) % binsPerDimension[d];
	 * multiplier *= binsPerDimension[d]; } return positionPerDimension; }
	 * 
	 * private static class BinPositionAndContribution { final private int
	 * position; final private double contribution;
	 * 
	 * private BinPositionAndContribution( final int position, final double
	 * contribution ) { this.position = position; this.contribution =
	 * contribution; } }
	 * 
	 * protected static class BoundsAndCounts { public final double minx; public
	 * final double maxx; public final double miny; public final double maxy;
	 * public final Double[] counts; public final int[] binsPerDimension;
	 * 
	 * public BoundsAndCounts( final ValueRange[] valueRangePerDimension, final
	 * Double[] counts, final int[] binsPerDimension ) {
	 * this.valueRangePerDimension = valueRangePerDimension; this.counts =
	 * counts; this.binsPerDimension = binsPerDimension; }
	 * 
	 * private void incrementCount( final int position, final double increment )
	 * { if (counts.length > position) { synchronized (counts) { if
	 * (counts[position] == null) { counts[position] = new Double( 0); } }
	 * counts[position] += increment; } else {
	 * logger.warn("position of count summary outside of bounds"); } } }
	 * 
	 * protected static class SummaryStatistics { final public
	 * List<BoundsAndCounts> boundsWithCounts;
	 * 
	 * public SummaryStatistics() { boundsWithCounts =
	 * Collections.synchronizedList(new ArrayList<BoundsAndCounts>()); }
	 * 
	 * public SummaryStatistics( final List<BoundsAndCounts> boundsWithCounts) {
	 * this.boundsWithCounts = boundsWithCounts; }
	 * 
	 * public ValueRange getCountMinMax() { double min = Double.MAX_VALUE;
	 * double max = -Double.MAX_VALUE; for (final BoundsAndCounts<RowImplType>
	 * boundsAndCount : boundsWithCounts) { for (final Double count :
	 * boundsAndCount.counts) { if (count != null) { min = Math.min( min,
	 * count); max = Math.max( max, count); } } } return new ValueRange( min,
	 * max); } }
	 */
}
