/* * Copyright (C) 2014 Ehsan Roshani * * 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 3 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, see <http://www.gnu.org/licenses/>. */ /* * This is the optimization problem to find the best fit variables for semivariogram * */ package whitebox.stats; /** * * @author Ehsan Roshani, Ph.D. Department of Geography University of Guelph * Guelph, Ont. N1G 2W1 CANADA Phone: (519) 824-4120 x53527 Email: * eroshani@uoguelph.ca * * modified by John Lindsay */ import jmetal.core.Problem; import jmetal.core.Solution; import jmetal.encodings.solutionType.RealSolutionType; import jmetal.util.JMException; import static whitebox.stats.Kriging.SemivariogramType.EXPONENTIAL; import static whitebox.stats.Kriging.SemivariogramType.GAUSSIAN; import static whitebox.stats.Kriging.SemivariogramType.SPHERICAL; public class SemivariogramCurveFitterProblem extends Problem { // defining the lower and upper limits //Range, Sill, Nugget public static double[] LOWERLIMIT; public static double[] UPPERLIMIT; double difMin = 1000000000; Kriging.SemivariogramType SVType; double[][] Pnts; boolean Nugget; public static Kriging.Variogram var; void SetData(double[][] Points, Kriging.SemivariogramType SemiVType, boolean Nugget) { this.Pnts = Points; this.SVType = SemiVType; double Xmax = 0; double Ymax = 0; double Xmin = Double.MAX_VALUE; double Ymin = Double.MAX_VALUE; for (int i = 0; i < Points.length; i++) { if (Points[i][0] < Xmin) { Xmin = Points[i][0]; } if (Points[i][1] < Ymin) { Ymin = Points[i][1]; } if (Points[i][0] > Xmax) { Xmax = Points[i][0]; } if (Points[i][1] > Ymax) { Ymax = Points[i][1]; } } if (Nugget) { LOWERLIMIT = new double[3]; UPPERLIMIT = new double[3]; LOWERLIMIT[0] = 0; LOWERLIMIT[1] = 0; LOWERLIMIT[2] = 0; UPPERLIMIT[0] = Xmax * 1.5; UPPERLIMIT[1] = Ymax; UPPERLIMIT[2] = Ymax; } else { LOWERLIMIT = new double[2]; UPPERLIMIT = new double[2]; LOWERLIMIT[0] = 0; LOWERLIMIT[1] = 0; UPPERLIMIT[0] = Xmax * 1.5; UPPERLIMIT[1] = Ymax; } } /** * Constructor. Creates a default instance of the Water problem. * * @param solutionType The solution type must "Real" or "BinaryReal". */ public SemivariogramCurveFitterProblem(double[][] Points, Kriging.SemivariogramType SVType, boolean Nugget) { this.Nugget = Nugget; if (Nugget) { numberOfVariables_ = 3; } else { numberOfVariables_ = 2; } SetData(Points, SVType, Nugget); numberOfObjectives_ = 2; numberOfConstraints_ = 0; problemName_ = "SemiVariogramCurveFitter"; upperLimit_ = new double[numberOfVariables_]; lowerLimit_ = new double[numberOfVariables_]; upperLimit_ = new double[numberOfVariables_]; lowerLimit_ = new double[numberOfVariables_]; for (int i = 0; i < numberOfVariables_; i++) { lowerLimit_[i] = LOWERLIMIT[i]; upperLimit_[i] = UPPERLIMIT[i]; } solutionType_ = new RealSolutionType(this); } // Roshani /** * Evaluates a solution * * @param solution The solution to evaluate * @throws JMException */ @Override public void evaluate(Solution solution) throws JMException { double[] values = new double[Pnts.length]; switch (SVType) { case EXPONENTIAL: for (int i = 0; i < Pnts.length; i++) { if (Pnts[i][0] != 0) { values[i] = (Nugget ? solution.getDecisionVariables()[2].getValue() : 0) + solution.getDecisionVariables()[1].getValue() * (1 - Math.exp(-Pnts[i][0] / solution.getDecisionVariables()[0].getValue())); } else { values[i] = 0; } } break; case GAUSSIAN: for (int i = 0; i < Pnts.length; i++) { if (Pnts[i][0] != 0) { values[i] = (Nugget ? solution.getDecisionVariables()[2].getValue() : 0) + solution.getDecisionVariables()[1].getValue() * (1 - Math.exp(-(Math.pow(Pnts[i][0], 2)) / (Math.pow( solution.getDecisionVariables()[0].getValue(), 2)))); } else { values[i] = 0; } } break; case SPHERICAL: for (int i = 0; i < Pnts.length; i++) { if (Pnts[0][0] > solution.getDecisionVariables()[0].getValue()) { values[i] = (Nugget ? solution.getDecisionVariables()[2].getValue() : 0) + solution.getDecisionVariables()[1].getValue(); } else if (0 < Pnts[0][0] && Pnts[0][0] <= solution.getDecisionVariables()[0].getValue()) { values[i] = (Nugget ? solution.getDecisionVariables()[2].getValue() : 0) + solution.getDecisionVariables()[1].getValue() * (1.5 * Pnts[i][0] / solution.getDecisionVariables()[0].getValue() - 0.5 * Math.pow((Pnts[i][0] / solution.getDecisionVariables()[0].getValue()), 3)); } else { values[i] = 0; } } break; } double mse = 0; for (int i = 0; i < Pnts.length; i++) { mse += Math.pow((values[i] - Pnts[i][1]), 2); } if (mse < difMin) { Kriging k = new Kriging(); var = k.getSemivariogram(SVType, solution.getDecisionVariables()[0].getValue(), solution.getDecisionVariables()[1].getValue(), (Nugget) ? solution.getDecisionVariables()[2].getValue() : 0, false); var.mse = mse; difMin = mse; } solution.setObjective(0, mse); solution.setObjective(1, mse); } /** * NOT USED Evaluates the constraint overhead of a solution * * @param solution The solution * @throws JMException */ @Override public void evaluateConstraints(Solution solution) throws JMException { // // // double [] constraint = new double[1]; // 7 constraints // // for (int i = 0; i < k.Pairs.size(); i++) { // if (k.Pairs.get(i).Distance<= k.resolution) { // constraint[0]+=(k.resolution - k.Pairs.get(i).Distance); // } // } // // solution.setOverallConstraintViolation(constraint[0]); // solution.setNumberOfViolatedConstraint(1); // System.out.println(constraint[0]); } // evaluateConstraints }