package com.imsl.test.example.math;

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

/**
 * <p>
 * Approximates a function with a polyharmonic spline 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 "custom" polyharmonic spline radial basis function is specified in
 * the class, {@link PolyHarmonicSpline}.
 * A radial basis approximation of the
 * noisy data is computed and 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 order of the Polyharmonic Spline, <i>k</i>=2. The
 * function is sampled at 200 random points and the error is computed at 10000
 * random points.</p>
 *
 * @see <a href="RadialBasisEx2.java">Code</a>
 * @see <a href="doc-files/RadialBasisEx2.out">Output</a>
 */
public class RadialBasisEx2 {

    // 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);
    }

    /**
     * <p>
     * RadialBasis Example 2b: Defines a polyharmonic spline radial basis
     * function.</p>
     *
     * The "custom" polyharmonic spline radial basis function is specified
     * $$\varphi(r)=\left\{
     * \begin{array}{ll}r^{k}&k=1,3,5,\ldots\\r^{k}\ln{r}&k=2,4,6,\ldots\end{array}
     * \right.$$
     */
    static public class PolyHarmonicSpline implements RadialBasis.Function {

        private int order = 3;
        private boolean isEven = false;

        public PolyHarmonicSpline(int order) {
            this.isEven = order % 2 == 0;
            this.order = order;
        }

        @Override
        public double f(double x) {
            if (this.isEven) {
                return Math.pow(x, order) * Math.log(x);
            }
            return Math.pow(x, order);
        }

        @Override
        public double g(double x) {
            if (order == 1) {
                return 1;
            }
            if (this.isEven) {
                return order * Math.pow(x, order - 1) * Math.log(x)
                        + Math.pow(x, order - 1);
            }
            return order * Math.pow(x, order - 1);
        }

        public int getOrder() {
            return order;
        }

        public boolean isEvenOrder() {
            return isEven;
        }
    }

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

        // 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(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 * 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 PolyHarmonicSpline(2));
        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 * (double) 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);
        System.out.println("Using even order equation: "
                + ((PolyHarmonicSpline) rb.getRadialFunction()).isEvenOrder());
    }
}