package com.imsl.test.example.math;

import com.imsl.math.*;
import java.util.Random;

/**
 * <p>
 * RadialBasis Example 1: Approximates a function with a 
 * Hardy multiquadric radial basis function.</p>
 *
 * Data is generated from the function $$e^{-\|\vec{x}\|^2/d}$$ where <i>d</i>
 * is the dimension, is evaluated at a set of randomly chosen points. Random
 * noise is added to the values and a radial basis approximation to the noisy
 * data is computed. The radial basis fit, using the default radial basis
 * function Hardy multiquadric with \(\delta=1\), is then compared to the
 * original function at another set of randomly chosen points. Both the average
 * error and the maximum error are computed and printed.
 *
 * <p>
 * In this example, the dimension <i>d</i>=10. The function is sampled at 200
 * random points, in the \([-1,1]^d\) hyper-cube, to which what noise in the
 * range [-0.2,0.2] is added. The error is computed at 1000 random points, also
 * from the \([-1,1]^d\) hyper-cube. The compute errors are less than the added
 * noise.
 * </p>
 *
 * @see <a href="RadialBasisEx1.java">Code</a>
 * @see <a href="doc-files/RadialBasisEx1.out">Output</a>
 */
public class RadialBasisEx1 {

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

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

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

        // Compute the error at a randomly selected set of points
        int nTest = 1000;
        double maxError = 0.0;
        double aveError = 0.0;
        double x[] = new double[nDim];
        for (int k = 0; k < nTest; k++) {
            for (int i = 0; i < nDim; i++) {
                x[i] = 2.0 * rand.nextDouble() - 1.0;
            }
            double error = Math.abs(fcn(x) - rb.value(x));
            aveError += error;
            maxError = Math.max(error, maxError);
        }
        aveError /= nTest;

        System.out.println("average error is " + aveError);
        System.out.println("maximum error is " + maxError);
    }

    // The function to approximate
    static private double fcn(double x[]) {
        double sum = 0.0;
        for (int k = 0; k < x.length; k++) {
            sum += x[k] * x[k];
        }
        sum /= x.length;
        return Math.exp(-sum);
    }
}