///////////////////////////////////////////////////////////////////////////////
// For information as to what this class does, see the Javadoc, below.       //
// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006,       //
// 2007, 2008, 2009, 2010, 2014, 2015 by Peter Spirtes, Richard Scheines, Joseph   //
// Ramsey, and Clark Glymour.                                                //
//                                                                           //
// This program is free software; you can redistribute it and/or modify      //
// it under the terms of the GNU General Public License as published by      //
// the Free Software Foundation; either version 2 of the License, or         //
// (at your option) any later version.                                       //
//                                                                           //
// This program is distributed in the hope that it will be useful,           //
// but WITHOUT ANY WARRANTY; without even the implied warranty of            //
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the             //
// GNU General Public License for more details.                              //
//                                                                           //
// You should have received a copy of the GNU General Public License         //
// along with this program; if not, write to the Free Software               //
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA //
///////////////////////////////////////////////////////////////////////////////

package edu.cmu.tetrad.search.kernel;

import edu.cmu.tetrad.data.DataSet;
import edu.cmu.tetrad.graph.Node;

/**
 * Gaussian kernel for a given bandwidth. Default bandwidth is set using the median distance heuristic.
 *
 * @author Robert Tillman
 */
public final class KernelGaussian implements Kernel {

    /**
     * The bandwidth
     */
    private double sigma;

    /**
     * Creates a new Gaussian kernel with the given bandwidth
     *
     * @param sigma the bandwidth
     */
    public KernelGaussian(double sigma) {
        if (sigma <= 0)
            throw new IllegalArgumentException("Sigma must be > 0");
        this.sigma = sigma;
    }

    /**
     * Creates a new Gaussian kernel using the median distance between points to set the bandwidth
     *
     * @param dataset dataet containing variable used to set bandwidth
     * @param node    variable used to set bandwidth
     */
    public KernelGaussian(DataSet dataset, Node node) {
        setMedianBandwidth(dataset, node);
    }

    /**
     * @return the bandwidth
     */
    public double getBandwidth() {
        return this.sigma;
    }

    /**
     * Evaluates the kernel at two given points
     *
     * @param i first point
     * @param j second point
     */
    public double eval(double i, double j) {
        double evalKernel = Math.exp(-.5 * (Math.pow((i - j), 2) / Math.pow(sigma, 2)));
        return evalKernel;
    }

    /**
     * Default setting of bandwidth based on median distance heuristic
     */
    public void setDefaultBw(DataSet dataset, Node node) {
        setMedianBandwidth(dataset, node);
    }


    /**
     * Sets the bandwidth of the kernel to median distance between two points in the given vector
     *
     * @param dataset dataet containing variable used to set bandwidth
     * @param node    variable used to set bandwidth
     */
    public void setMedianBandwidth(DataSet dataset, Node node) {
        int col = dataset.getColumn(node);
        int m = dataset.getNumRows();

        double[] diff = new double[(int) Math.pow(m, 2) - m];
        int c = 0;
        for (int i = 0; i < (m - 1); i++) {
            for (int j = (i + 1); j < m; j++) {
                diff[c] = Math.abs(dataset.getDouble(i, col) - dataset.getDouble(j, col));
                c++;
            }
        }

        this.sigma = find(diff, 0, (m - 1));
    }

    // private method for finding median distance

    private double find(double[] a, int from, int to) {
        int low = from;
        int high = to;
        int median = (low + high) / 2;
        do {
            if (high <= low) {
                return a[median];
            }
            if (high == low + 1) {
                if (a[low] > a[high]) {
                    swap(a, low, high);
                }
                return a[median];
            }
            int middle = (low + high) / 2;
            if (a[middle] > a[high]) {
                swap(a, middle, high);
            }
            if (a[low] > a[high]) {
                swap(a, low, high);
            }
            if (a[middle] > a[low]) {
                swap(a, middle, low);
            }
            swap(a, middle, low + 1);
            int ll = low + 1;
            int hh = high;
            do {
                do {
                    ll++;
                }
                while (a[low] > a[ll]);
                do {
                    hh--;
                }
                while (a[hh] > a[low]);
                if (hh < ll) {
                    break;
                }
                swap(a, ll, hh);
            }
            while (true);
            swap(a, low, hh);
            if (hh <= median) {
                low = ll;
            }
            if (hh >= median) {
                high = hh - 1;
            }
        }
        while (true);
    }

    private void swap(double[] a, int i1, int i2) {
        double temp = a[i1];
        a[i1] = a[i2];
        a[i2] = temp;
    }


}