///////////////////////////////////////////////////////////////////////////////
// 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, 2022 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.sem;

import edu.cmu.tetrad.data.DataUtils;
import edu.cmu.tetrad.graph.Node;
import edu.cmu.tetrad.util.Matrix;
import edu.cmu.tetrad.util.RandomUtil;
import edu.cmu.tetrad.util.TetradLogger;
import org.apache.commons.math3.util.FastMath;

import java.util.List;

/**
 * Optimizes a SEM by randomly selecting points in cubes of decreasing size about
 * a given point.
 *
 * @author josephramsey
 */
public class SemOptimizerScattershot implements SemOptimizer {
    static final long serialVersionUID = 23L;
    private int numRestarts;

    //=============================CONSTRUCTORS=========================//

    /**
     * Blank constructor.
     */
    public SemOptimizerScattershot() {
    }

    /**
     * Generates a simple exemplar of this class to test serialization.
     */
    public static SemOptimizerScattershot serializableInstance() {
        return new SemOptimizerScattershot();
    }

    //==============================PUBLIC METHODS========================//

    /**
     * Optimizes the fitting function of the given Sem using the Powell method
     * from Numerical Recipes by adjusting the freeParameters of the Sem.
     */
    public void optimize(SemIm semIm) {
        Matrix sampleCovar = semIm.getSampleCovar();

        if (sampleCovar == null) {
            throw new NullPointerException("Sample covar has not been set.");
        }

        if (DataUtils.containsMissingValue(sampleCovar)) {
            throw new IllegalArgumentException("Please remove or impute missing values.");
        }

        if (DataUtils.containsMissingValue(sampleCovar)) {
            throw new IllegalArgumentException("Please remove or impute missing values.");
        }

        if (this.numRestarts < 1) this.numRestarts = 1;

        TetradLogger.getInstance().log("info", "Trying EM...");
        TetradLogger.getInstance().log("info", "Trying scattershot...");

        double min = Double.POSITIVE_INFINITY;
        SemIm _sem = null;

        // With local search on points in the width 1 iteration, multiple iterations of the whole search
        // doesn't seem necessary.
        for (int i = 0; i < this.numRestarts + 1; i++) {
            TetradLogger.getInstance().log("details", "Trial " + (i + 1));
//            System.out.println("Trial " + (i + 1));
            SemIm _sem2 = new SemIm(semIm);
            optimize2(_sem2);
            double chisq = _sem2.getChiSquare();

            if (FastMath.abs(chisq) < min) {
                min = FastMath.abs(chisq);
                _sem = _sem2;
            }
        }

        if (_sem == null) {
            throw new NullPointerException("Minimal score SEM could not be found.");
        }

//        new SemOptimizerPalCds().optimize2(_sem);

        for (Parameter param : semIm.getFreeParameters()) {
            Node nodeA = param.getNodeA();
            Node nodeB = param.getNodeB();

            Node _nodeA = _sem.getVariableNode(nodeA.getName());
            Node _nodeB = _sem.getVariableNode(nodeB.getName());

            double value = _sem.getParamValue(_nodeA, _nodeB);
            semIm.setParamValue(param, value);
        }


//        optimize2(semIm);
    }

    @Override
    public void setNumRestarts(int numRestarts) {
        this.numRestarts = numRestarts;
    }

    @Override
    public int getNumRestarts() {
        return this.numRestarts;
    }

    public String toString() {
        return "Sem Optimizer Scattershot";
    }

    private void optimize2(SemIm semIm) {
        FittingFunction f = new SemFittingFunction(semIm);

        double[] p = semIm.getFreeParamValues();

        f.setAvoidNegativeVariances(true);
        iterateFindLowerRandom(f, p, 1.0, 1500);
        iterateFindLowerRandom(f, p, 0.5, 500);
        iterateFindLowerRandom(f, p, 0.25, 500);
        iterateFindLowerRandom(f, p, 0.1, 500);
        iterateFindLowerRandom(f, p, 0.1, 500);
        iterateFindLowerRandom(f, p, 0.05, 500);
        iterateFindLowerRandom(f, p, 0.01, 500);
        iterateFindLowerRandom(f, p, 0.005, 50);
        iterateFindLowerRandom(f, p, 0.001, 50);
        iterateFindLowerRandom(f, p, 0.0005, 50);
        iterateFindLowerRandom(f, p, 0.0001, 50);

        semIm.setFreeParamValues(p);
    }

    private void iterateFindLowerRandom(FittingFunction fcn, double[] p,
                                        double range, int iterations) {
        int t = 0;

        while (++t < 2000) {
            boolean found;

            try {
                found = findLowerRandom(fcn, p, range, iterations);
            } catch (Exception e) {
                return;
            }

            if (!found) {
                return;
            }
        }
    }

    /**
     * @return true iff a new point was found with a lower score.
     */
    private boolean findLowerRandom(FittingFunction fcn, double[] p,
                                    double width, int numPoints) {
        double fP = fcn.evaluate(p);

        if (Double.isNaN(fP)) {
            throw new IllegalArgumentException("Center point must evaluate!");
        }

        // This point will remain fixed, the center of the search.
        double[] fixedP = new double[p.length];
        System.arraycopy(p, 0, fixedP, 0, p.length);

        // This point will move around randomly. If it ever has a lower
        // score than p, it will be copied into p (and returned).
        double[] pTemp = new double[p.length];
        System.arraycopy(p, 0, pTemp, 0, p.length);

        for (int i = 0; i < numPoints; i++) {
            randomPointAboutCenter(pTemp, fixedP, width);
            double f = fcn.evaluate(pTemp);

            if (f == Double.POSITIVE_INFINITY) {
                i--;
                continue;
            }

            // Try to find the lowest reachable spot for each of the trial points at width 1.
            if (width == 1) {
                int t = 0;
                while (++t < 2000) {
                    if (!findLowerRandomLocal(fcn, pTemp)) break;
                }
            }

            if (f < fP) {
                System.arraycopy(pTemp, 0, p, 0, pTemp.length);
                TetradLogger.getInstance().log("optimization", "Cube width = " + width + " FML = " + f);
                return true;
            }
        }

        return false;
    }

    private boolean findLowerRandomLocal(FittingFunction fcn, double[] p) {
        double fP = fcn.evaluate(p);

        if (Double.isNaN(fP)) {
            throw new IllegalArgumentException("Center point must evaluate!");
        }

        // This point will remain fixed, the center of the search.
        double[] fixedP = new double[p.length];
        System.arraycopy(p, 0, fixedP, 0, p.length);

        // This point will move around randomly. If it ever has a lower
        // score than p, it will be copied into p (and returned).
        double[] pTemp = new double[p.length];
        System.arraycopy(p, 0, pTemp, 0, p.length);

        for (int i = 0; i < 10; i++) {
            randomPointAboutCenter(pTemp, fixedP, 0.2);
            double f = fcn.evaluate(pTemp);

            if (f == Double.POSITIVE_INFINITY) {
                i++;
                continue;
            }

            if (f < fP) {
                System.arraycopy(pTemp, 0, p, 0, pTemp.length);
                TetradLogger.getInstance().log("optimization", "Cube width = " + 0.2 + " FML = " + f);
                return true;
            }
        }

        return false;
    }

    private void randomPointAboutCenter(double[] pTemp, double[] fixedP, double width) {
        for (int j = 0; j < pTemp.length; j++) {
            double v = getRandom().nextDouble();
            pTemp[j] = fixedP[j] + (-width / 2.0 + width * v);
        }
    }

    /**
     * Evaluates a fitting function for an array of freeParameters.
     *
     * @author josephramsey
     */
    interface FittingFunction {

        /**
         * @return the value of the function for the given array of parameter
         * values.
         */
        double evaluate(double[] argument);

        void setAvoidNegativeVariances(boolean avoidNegativeVariances);
    }

    /**
     * Wraps a Sem for purposes of calculating its fitting function for given
     * parameter values.
     *
     * @author josephramsey
     */
    static class SemFittingFunction implements FittingFunction {

        /**
         * The wrapped Sem.
         */
        private final SemIm sem;
        private final List<Parameter> freeParameters;
        private boolean avoidNegativeVariances;

        /**
         * Constructs a new CoefFittingFunction for the given Sem.
         */
        public SemFittingFunction(SemIm sem) {
            this.sem = sem;
            this.freeParameters = sem.getFreeParameters();
        }

        /**
         * Computes the maximum likelihood function value for the given                                   G
         * freeParameters values as given by the optimizer. These values are mapped
         * to parameter values.
         */
        public double evaluate(double[] parameters) {
            this.sem.setFreeParamValues(parameters);

            for (double parameter : parameters) {
                if (Double.isNaN(parameter) || Double.isInfinite(parameter)) {
                    return Double.POSITIVE_INFINITY;
                }
            }

            double fml = this.sem.getScore();

            if (Double.isNaN(fml) || Double.isInfinite(fml)) {
                return Double.POSITIVE_INFINITY;
            }

            if (this.avoidNegativeVariances) {
                for (int i = 0; i < parameters.length; i++) {
                    if (this.freeParameters.get(i).getType() == ParamType.VAR && parameters[i] <= 0.0) {
                        return Double.POSITIVE_INFINITY;
                    }
                }
            }

            if (Double.isNaN(fml)) {
                return Double.POSITIVE_INFINITY;
            }

            if (fml < 0) {
                return Double.POSITIVE_INFINITY;
            }

            return fml;
        }

        @Override
        public void setAvoidNegativeVariances(boolean avoidNegativeVariances) {
            this.avoidNegativeVariances = avoidNegativeVariances;
        }
    }

    private RandomUtil getRandom() {
        return RandomUtil.getInstance();
    }
}



