/*
 * Licensed to the Apache Software Foundation (ASF) under one or more
 * contributor license agreements.  See the NOTICE file distributed with
 * this work for additional information regarding copyright ownership.
 * The ASF licenses this file to You under the Apache License, Version 2.0
 * (the "License"); you may not use this file except in compliance with
 * the License.  You may obtain a copy of the License at
 *
 *      http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
package org.apache.commons.math.transform;

import org.apache.commons.math.FunctionEvaluationException;
import org.apache.commons.math.MathRuntimeException;
import org.apache.commons.math.analysis.UnivariateRealFunction;

/**
 * Implements the
 * <a href="http://www.archive.chipcenter.com/dsp/DSP000517F1.html">Fast
 * Hadamard Transform</a> (FHT). Transformation of an input vector x to the
 * output vector y.
 * <p>
 * In addition to transformation of real vectors, the Hadamard transform can
 * transform integer vectors into integer vectors. However, this integer
 * transform cannot be inverted directly. Due to a scaling factor it may lead to
 * rational results. As an example, the inverse transform of integer vector (0,
 * 1, 0, 1) is rational vector (1/2, -1/2, 0, 0).
 * </p>
 * 
 * @version $Revision$ $Date$
 * @since 2.0
 */
public class FastHadamardTransformer implements RealTransformer {

	/** {@inheritDoc} */
	public double[] transform(double f[]) throws IllegalArgumentException {
		return fht(f);
	}

	/** {@inheritDoc} */
	public double[] transform(UnivariateRealFunction f, double min, double max, int n)
			throws FunctionEvaluationException, IllegalArgumentException {
		return fht(FastFourierTransformer.sample(f, min, max, n));
	}

	/** {@inheritDoc} */
	public double[] inversetransform(double f[]) throws IllegalArgumentException {
		return FastFourierTransformer.scaleArray(fht(f), 1.0 / f.length);
	}

	/** {@inheritDoc} */
	public double[] inversetransform(UnivariateRealFunction f, double min, double max, int n)
			throws FunctionEvaluationException, IllegalArgumentException {
		final double[] unscaled = fht(FastFourierTransformer.sample(f, min, max, n));
		return FastFourierTransformer.scaleArray(unscaled, 1.0 / n);
	}

	/**
	 * Transform the given real data set.
	 * <p>
	 * The integer transform cannot be inverted directly, due to a scaling
	 * factor it may lead to double results.
	 * </p>
	 * 
	 * @param f
	 *            the integer data array to be transformed (signal)
	 * @return the integer transformed array (spectrum)
	 * @throws IllegalArgumentException
	 *             if any parameters are invalid
	 */
	public int[] transform(int f[]) throws IllegalArgumentException {
		return fht(f);
	}

	/**
	 * The FHT (Fast Hadamard Transformation) which uses only subtraction and
	 * addition. <br>
	 * Requires <b>Nlog2N = n2</b><sup>n</sup> additions. <br>
	 * <br>
	 * <b><u>Short Table of manual calculation for N=8:</u></b>
	 * <ol>
	 * <li><b>x</b> is the input vector we want to transform</li>
	 * <li><b>y</b> is the output vector which is our desired result</li>
	 * <li>a and b are just helper rows</li>
	 * </ol>
	 * 
	 * <pre>
	 * <code>
	 * +----+----------+---------+----------+
	 * | <b>x</b>  |    <b>a</b>     |    <b>b</b>    |    <b>y</b>     |
	 * +----+----------+---------+----------+
	 * | x<sub>0</sub> | a<sub>0</sub>=x<sub>0</sub>+x<sub>1</sub> | b<sub>0</sub>=a<sub>0</sub>+a<sub>1</sub> | y<sub>0</sub>=b<sub>0</sub>+b<sub>1</sub> |
	 * +----+----------+---------+----------+
	 * | x<sub>1</sub> | a<sub>1</sub>=x<sub>2</sub>+x<sub>3</sub> | b<sub>0</sub>=a<sub>2</sub>+a<sub>3</sub> | y<sub>0</sub>=b<sub>2</sub>+b<sub>3</sub> |
	 * +----+----------+---------+----------+
	 * | x<sub>2</sub> | a<sub>2</sub>=x<sub>4</sub>+x<sub>5</sub> | b<sub>0</sub>=a<sub>4</sub>+a<sub>5</sub> | y<sub>0</sub>=b<sub>4</sub>+b<sub>5</sub> |
	 * +----+----------+---------+----------+
	 * | x<sub>3</sub> | a<sub>3</sub>=x<sub>6</sub>+x<sub>7</sub> | b<sub>0</sub>=a<sub>6</sub>+a<sub>7</sub> | y<sub>0</sub>=b<sub>6</sub>+b<sub>7</sub> |
	 * +----+----------+---------+----------+
	 * | x<sub>4</sub> | a<sub>0</sub>=x<sub>0</sub>-x<sub>1</sub> | b<sub>0</sub>=a<sub>0</sub>-a<sub>1</sub> | y<sub>0</sub>=b<sub>0</sub>-b<sub>1</sub> |
	 * +----+----------+---------+----------+
	 * | x<sub>5</sub> | a<sub>1</sub>=x<sub>2</sub>-x<sub>3</sub> | b<sub>0</sub>=a<sub>2</sub>-a<sub>3</sub> | y<sub>0</sub>=b<sub>2</sub>-b<sub>3</sub> |
	 * +----+----------+---------+----------+
	 * | x<sub>6</sub> | a<sub>2</sub>=x<sub>4</sub>-x<sub>5</sub> | b<sub>0</sub>=a<sub>4</sub>-a<sub>5</sub> | y<sub>0</sub>=b<sub>4</sub>-b<sub>5</sub> |
	 * +----+----------+---------+----------+
	 * | x<sub>7</sub> | a<sub>3</sub>=x<sub>6</sub>-x<sub>7</sub> | b<sub>0</sub>=a<sub>6</sub>-a<sub>7</sub> | y<sub>0</sub>=b<sub>6</sub>-b<sub>7</sub> |
	 * +----+----------+---------+----------+
	 * </code>
	 * </pre>
	 * 
	 * <b><u>How it works</u></b>
	 * <ol>
	 * <li>Construct a matrix with N rows and n+1 columns<br>
	 * <b>hadm[n+1][N]</b> <br>
	 * <i>(If I use [x][y] it always means [row-offset][column-offset] of a
	 * Matrix with n rows and m columns. Its entries go from M[0][0] to
	 * M[n][m])</i></li>
	 * <li>Place the input vector <b>x[N]</b> in the first column of the matrix
	 * <b>hadm</b></li>
	 * <li>The entries of the submatrix D<sub>top</sub> are calculated as
	 * follows. <br>
	 * D<sub>top</sub> goes from entry [0][1] to [N/2-1][n+1]. <br>
	 * The columns of D<sub>top</sub> are the pairwise mutually exclusive sums
	 * of the previous column</li>
	 * <li>The entries of the submatrix D<sub>bottom</sub> are calculated as
	 * follows. <br>
	 * D<sub>bottom</sub> goes from entry [N/2][1] to [N][n+1]. <br>
	 * The columns of D<sub>bottom</sub> are the pairwise differences of the
	 * previous column</li>
	 * <li>How D<sub>top</sub> and D<sub>bottom</sub> you can understand best
	 * with the example for N=8 above.
	 * <li>The output vector y is now in the last column of <b>hadm</b></li>
	 * <li><i>Algorithm from:
	 * http://www.archive.chipcenter.com/dsp/DSP000517F1.html</i></li>
	 * </ol>
	 * <br>
	 * <b><u>Visually</u></b>
	 * 
	 * <pre>
	 *        +--------+---+---+---+-----+---+
	 *        |   0    | 1 | 2 | 3 | ... |n+1|
	 * +------+--------+---+---+---+-----+---+
	 * |0     | x<sub>0</sub>     |       /\            |
	 * |1     | x<sub>1</sub>     |       ||            |
	 * |2     | x<sub>2</sub>     |   <= D<sub>top</sub>  =>       |
	 * |...   | ...    |       ||            |
	 * |N/2-1 | x<sub>N/2-1</sub>  |       \/            |
	 * +------+--------+---+---+---+-----+---+
	 * |N/2   | x<sub>N/2</sub>   |       /\            |
	 * |N/2+1 | x<sub>N/2+1</sub>  |       ||            |
	 * |N/2+2 | x<sub>N/2+2</sub>  |  <= D<sub>bottom</sub>  =>      | which is in the last column of the matrix
	 * |...   | ...    |       ||            |
	 * |N     | x<sub>N/2</sub>   |        \/           |
	 * +------+--------+---+---+---+-----+---+
	 * </pre>
	 * 
	 * @param x
	 *            input vector
	 * @return y output vector
	 * @exception IllegalArgumentException
	 *                if input array is not a power of 2
	 */
	protected double[] fht(double x[]) throws IllegalArgumentException {

		// n is the row count of the input vector x
		final int n = x.length;
		final int halfN = n / 2;

		// n has to be of the form n = 2^p !!
		if (!FastFourierTransformer.isPowerOf2(n)) {
			throw MathRuntimeException.createIllegalArgumentException("{0} is not a power of 2", n);
		}

		// Instead of creating a matrix with p+1 columns and n rows
		// we will use two single dimension arrays which we will use in an
		// alternating way.
		double[] yPrevious = new double[n];
		double[] yCurrent = x.clone();

		// iterate from left to right (column)
		for (int j = 1; j < n; j <<= 1) {

			// switch columns
			final double[] yTmp = yCurrent;
			yCurrent = yPrevious;
			yPrevious = yTmp;

			// iterate from top to bottom (row)
			for (int i = 0; i < halfN; ++i) {
				// D<sub>top</sub>
				// The top part works with addition
				final int twoI = 2 * i;
				yCurrent[i] = yPrevious[twoI] + yPrevious[twoI + 1];
			}
			for (int i = halfN; i < n; ++i) {
				// D<sub>bottom</sub>
				// The bottom part works with subtraction
				final int twoI = 2 * i;
				yCurrent[i] = yPrevious[twoI - n] - yPrevious[twoI - n + 1];
			}
		}

		// return the last computed output vector y
		return yCurrent;

	}

	/**
	 * The FHT (Fast Hadamard Transformation) which uses only subtraction and
	 * addition.
	 * 
	 * @param x
	 *            input vector
	 * @return y output vector
	 * @exception IllegalArgumentException
	 *                if input array is not a power of 2
	 */
	protected int[] fht(int x[]) throws IllegalArgumentException {

		// n is the row count of the input vector x
		final int n = x.length;
		final int halfN = n / 2;

		// n has to be of the form n = 2^p !!
		if (!FastFourierTransformer.isPowerOf2(n)) {
			throw MathRuntimeException.createIllegalArgumentException("{0} is not a power of 2", n);
		}

		// Instead of creating a matrix with p+1 columns and n rows
		// we will use two single dimension arrays which we will use in an
		// alternating way.
		int[] yPrevious = new int[n];
		int[] yCurrent = x.clone();

		// iterate from left to right (column)
		for (int j = 1; j < n; j <<= 1) {

			// switch columns
			final int[] yTmp = yCurrent;
			yCurrent = yPrevious;
			yPrevious = yTmp;

			// iterate from top to bottom (row)
			for (int i = 0; i < halfN; ++i) {
				// D<sub>top</sub>
				// The top part works with addition
				final int twoI = 2 * i;
				yCurrent[i] = yPrevious[twoI] + yPrevious[twoI + 1];
			}
			for (int i = halfN; i < n; ++i) {
				// D<sub>bottom</sub>
				// The bottom part works with subtraction
				final int twoI = 2 * i;
				yCurrent[i] = yPrevious[twoI - n] - yPrevious[twoI - n + 1];
			}
		}

		// return the last computed output vector y
		return yCurrent;

	}

}
