package com.imsl.test.example.math;

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

/**
 * <p>
 * Approximates a function using a Hardy multiquadric radial basis function.</p>
 *
 * Data is generated from the function
 * $$e^\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 Hardy multiquadric radial basis function is specified
 * \(\sqrt{r^2+\delta^2}\) 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 Hardy multiquadric radial basis
 * function \(\delta = 5.5\). The function is sampled at 100 random points and
 * the error is computed at 10000 random points.</p>
 *
 * @see <a href="RadialBasisEx3.java">Code</a>
 * @see <a href="doc-files/RadialBasisEx3.out">Output</a>
 */
public class RadialBasisEx3 {

    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.HardyMultiquadric(5.5));
        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);
    }
}