// 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.optimization.fitting;

import static org.junit.Assert.assertEquals;
import static org.junit.Assert.assertTrue;

import java.util.Random;

import org.apache.commons.math.analysis.polynomials.PolynomialFunction;
import org.apache.commons.math.optimization.DifferentiableMultivariateVectorialOptimizer;
import org.apache.commons.math.optimization.OptimizationException;
import org.apache.commons.math.optimization.general.GaussNewtonOptimizer;
import org.apache.commons.math.optimization.general.LevenbergMarquardtOptimizer;
import org.junit.Test;

public class PolynomialFitterTest {

	@Test
	public void testNoError() throws OptimizationException {
		Random randomizer = new Random(64925784252l);
		for (int degree = 1; degree < 10; ++degree) {
			PolynomialFunction p = buildRandomPolynomial(degree, randomizer);

			PolynomialFitter fitter = new PolynomialFitter(degree, new LevenbergMarquardtOptimizer());
			for (int i = 0; i <= degree; ++i) {
				fitter.addObservedPoint(1.0, i, p.value(i));
			}

			PolynomialFunction fitted = fitter.fit();

			for (double x = -1.0; x < 1.0; x += 0.01) {
				double error = Math.abs(p.value(x) - fitted.value(x)) / (1.0 + Math.abs(p.value(x)));
				assertEquals(0.0, error, 1.0e-6);
			}

		}

	}

	@Test
	public void testSmallError() throws OptimizationException {
		Random randomizer = new Random(53882150042l);
		double maxError = 0;
		for (int degree = 0; degree < 10; ++degree) {
			PolynomialFunction p = buildRandomPolynomial(degree, randomizer);

			PolynomialFitter fitter = new PolynomialFitter(degree, new LevenbergMarquardtOptimizer());
			for (double x = -1.0; x < 1.0; x += 0.01) {
				fitter.addObservedPoint(1.0, x, p.value(x) + 0.1 * randomizer.nextGaussian());
			}

			PolynomialFunction fitted = fitter.fit();

			for (double x = -1.0; x < 1.0; x += 0.01) {
				double error = Math.abs(p.value(x) - fitted.value(x)) / (1.0 + Math.abs(p.value(x)));
				maxError = Math.max(maxError, error);
				assertTrue(Math.abs(error) < 0.1);
			}
		}
		assertTrue(maxError > 0.01);

	}

	@Test
	public void testRedundantSolvable() {
		// Levenberg-Marquardt should handle redundant information gracefully
		checkUnsolvableProblem(new LevenbergMarquardtOptimizer(), true);
	}

	@Test
	public void testRedundantUnsolvable() {
		// Gauss-Newton should not be able to solve redundant information
		DifferentiableMultivariateVectorialOptimizer optimizer = new GaussNewtonOptimizer(true);
		checkUnsolvableProblem(optimizer, false);
	}

	private void checkUnsolvableProblem(DifferentiableMultivariateVectorialOptimizer optimizer, boolean solvable) {
		Random randomizer = new Random(1248788532l);
		for (int degree = 0; degree < 10; ++degree) {
			PolynomialFunction p = buildRandomPolynomial(degree, randomizer);

			PolynomialFitter fitter = new PolynomialFitter(degree, optimizer);

			// reusing the same point over and over again does not bring
			// information, the problem cannot be solved in this case for
			// degrees greater than 1 (but one point is sufficient for
			// degree 0)
			for (double x = -1.0; x < 1.0; x += 0.01) {
				fitter.addObservedPoint(1.0, 0.0, p.value(0.0));
			}

			try {
				fitter.fit();
				assertTrue(solvable || (degree == 0));
			} catch (OptimizationException e) {
				assertTrue((!solvable) && (degree > 0));
			}

		}

	}

	private PolynomialFunction buildRandomPolynomial(int degree, Random randomizer) {
		final double[] coefficients = new double[degree + 1];
		for (int i = 0; i <= degree; ++i) {
			coefficients[i] = randomizer.nextGaussian();
		}
		return new PolynomialFunction(coefficients);
	}

}
