package com.imsl.test.example.math;

import com.imsl.math.*;
import com.imsl.stat.Random;

/**
 * <p>
 * Approximates a function with a Gaussian radial basis function. </p>
 *
 * Data is generated from the function
 * $$\exp(\frac{y}{2.0}\sin{x})\cos(\frac{y}{2.0})$$ where a number of
 * (<i>x,y</i>) pairs make up a set of randomly chosen points. Random noise is
 * added to the values, a Gaussian radial basis function \(e^{-ax^2}\) is
 * specified, and a radial basis approximation of the noisy data is computed.
 * The radial basis fit is then compared to the original function at another set
 * of randomly chosen points. Both the average normalized error and the maximum
 * normalized error are computed and printed.
 *
 * <p>
 * In this example, the parameter of the Gaussian radial basis function \(a =
 * 0.1\). The function is sampled at \(100\) random points and the error is
 * computed at \)10000\) random points.</p>
 *
 * @see <a href="RadialBasisEx4.java">Code</a>
 * @see <a href="doc-files/RadialBasisEx4.out">Output</a>
 */
public class RadialBasisEx4 {

    public static void main(String args[]) {
        int nDim = 2;

        // Sample, with noise, the function at 100 randomly choosen points
        int nData = 100;
        double xData[][] = new double[nData][nDim];
        double fData[] = new double[nData];
        Random rand = new Random(123457);
        rand.setMultiplier(16807);
        double noise[] = new double[nData * nDim];
        for (int k = 0; k < nData; k++) {
            for (int i = 0; i < nDim; i++) {
                noise[k * 2 + i] = 1.0d - 2.0d * (double) rand.nextDouble();
                xData[k][i] = 3 * noise[k * 2 + i];
            }
            // noisy sample
            fData[k] = fcn(xData[k]) + noise[k * 2] / 10;
        }

        // Compute the radial basis approximation using 100 centers
        int nCenters = 100;
        RadialBasis rb = new RadialBasis(nDim, nCenters);
        rb.setRadialFunction(new RadialBasis.Gaussian(.1));
        rb.update(xData, fData);

        // Compute the error at a randomly selected set of points
        int nTest = 10000;
        double maxError = 0.0;
        double aveError = 0.0;
        double maxMagnitude = 0.0;
        double x[][] = new double[nTest][nDim];
        noise = new double[nTest * nDim];

        for (int i = 0; i < nTest; i++) {
            for (int j = 0; j < nDim; j++) {
                noise[i * 2 + j] = 1.0d - 2.0d * rand.nextDouble();
                x[i][j] = 3 * noise[i * 2 + j];
            }
            double error = Math.abs(fcn(x[i]) - rb.value(x[i]));
            maxMagnitude = Math.max(Math.abs(fcn(x[i])), maxMagnitude);
            aveError += error;
            maxError = Math.max(error, maxError);
        }
        aveError /= nTest;

        System.out.println("Average normalized error is "
                + aveError / maxMagnitude);
        System.out.println("Maximum normalized error is "
                + maxError / maxMagnitude);
    }

    // The function to approximate
    static private double fcn(double x[]) {
        return Math.exp((x[1]) / 2.0) * Math.sin(x[0]) - Math.cos((x[1]) / 2.0);
    }
}