ngmf.util.cosu.NelderMead Maven / Gradle / Ivy
Go to download
Show more of this group Show more artifacts with this name
Show all versions of oms Show documentation
Show all versions of oms Show documentation
Object Modeling System (OMS) is a pure Java object-oriented framework.
OMS v3.+ is a highly interoperable and lightweight modeling framework for component-based model and simulation development on multiple platforms.
/*
* $Id: NelderMead.java 50798ee5e25c 2013-01-09 [email protected] $
*
* This file is part of the Object Modeling System (OMS),
* 2007-2012, Olaf David and others, Colorado State University.
*
* OMS is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published by
* the Free Software Foundation, version 2.1.
*
* OMS 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 Lesser General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with OMS. If not, see .
*/
package ngmf.util.cosu;
/**
*
* @author od
*/
import java.util.Arrays;
import java.util.Comparator;
import java.util.Random;
import oms3.util.Statistics;
public class NelderMead {
static class SampleComperator implements Comparator {
private int order = 1;
public SampleComperator(boolean decreasing_order) {
order = decreasing_order ? -1 : 1;
}
@Override
public int compare(Sample d1, Sample d2) {
if (d1.fx < d2.fx) {
return -1 * order;
} else if (d1.fx == d2.fx) {
return 0 * order;
} else {
return 1 * order;
}
}
}
Sample[] initialSimplex = null;
//number of parameters!!
int n;
double[] lowBound;
double[] upBound;
double x0[] = null;
int mode;
Random generator = new Random();
int currentSampleCount;
double effValue;
int maxn;
double[] parameters;
void sort(Sample[] array) {
Arrays.sort(array, new SampleComperator(false));
}
double normalizedgeometricRange(Sample x[]) {
if (x.length == 0) {
return 0;
}
double min[] = new double[n];
double max[] = new double[n];
double mean = 0;
for (int i = 0; i < n; i++) {
min[i] = Double.POSITIVE_INFINITY;
max[i] = Double.NEGATIVE_INFINITY;
for (int j = 0; j < x.length; j++) {
min[i] = Math.min(x[j].x[i], min[i]);
max[i] = Math.max(x[j].x[i], max[i]);
}
mean += Math.log(max[i] - min[i]);
}
return Math.exp(mean / n);
}
boolean feasible(double point[]) {
for (int i = 0; i < point.length; i++) {
if (point[i] < lowBound[i] || point[i] > upBound[i]) {
return false;
}
}
return true;
}
Sample getSample(double[] x) {
return new Sample(x, funct(x));
}
public double funct(double x[]) {
//RefreshDataHandles();
for (int j = 0; j < parameters.length; j++) {
try {
parameters[j] = x[j];
} catch (Exception e) {
throw new RuntimeException("Error! Parameter No. " + j + " wasn^t found" + e.toString());
}
}
// singleRun();
double value = effValue;
//sometimes its a bad idea to calculate with NaN or Infty
double bigNumber = 10000000;
effValue = Math.max(effValue, -bigNumber);
effValue = Math.min(effValue, bigNumber);
if (Double.isNaN(effValue)) {
effValue = -bigNumber;
}
currentSampleCount++;
switch (mode) {
case Statistics.MINIMIZATION:
return value;
case Statistics.MAXIMIZATION:
return -value;
case Statistics.ABSMINIMIZATION:
return Math.abs(value);
case Statistics.ABSMAXIMIZATION:
return -Math.abs(value);
default:
return 0.0;
}
}
double[] randomSample() {
double[] sample = new double[n];
for (int i = 0; i < n; i++) {
sample[i] = (lowBound[i] + generator.nextDouble() * (upBound[i] - lowBound[i]));
}
return sample;
}
public void run() {
//first draw n+1 random points
Sample[] simplex = new Sample[n + 1];
if (initialSimplex != null) {
simplex = initialSimplex;
} else {
for (int i = 0; i < n + 1; i++) {
if (i == 0 && x0 != null) {
simplex[i] = getSample(x0);
} else {
simplex[i] = getSample(randomSample());
}
}
}
int m = simplex.length;
double alpha = 1.0;
double gamma = 2.0;
double rho = 0.5;
double sigma = 0.5;
double epsilon = 0.01;
double max_restart_count = 5;
int restart_counter = 0;
int iterationcounter = 0;
while (true) {
if (iterationcounter++ > maxn) {
System.out.println("*********************************************************");
System.out.println("Maximum number of iterations reached, finished optimization");
System.out.println("Bestpoint:" + simplex[0]);
System.out.println("*********************************************************");
return;
}
if (normalizedgeometricRange(simplex) < epsilon) {
if (max_restart_count < ++restart_counter) {
System.out.println("*********************************************************");
System.out.println("Maximum number of restarts reached, finished optimization");
System.out.println("Bestpoint:" + simplex[0]);
System.out.println("*********************************************************");
return;
}
System.out.println("restart");
for (int i = 1; i < m; i++) {
simplex[i] = getSample(randomSample());
}
}
sort(simplex);
// Compute the centroid of the simplex
double centroid[] = new double[n];
for (int j = 0; j < n; j++) {
centroid[j] = 0;
for (int i = 0; i < m - 1; i++) {
centroid[j] += simplex[i].x[j] * (1.0 / (m - 1.0));
}
}
//reflect
double reflection[] = new double[n];
for (int i = 0; i < n; i++) {
reflection[i] = centroid[i] + alpha * (centroid[i] - simplex[m - 1].x[i]);
}
Sample reflection_sample = null;
if (feasible(reflection)) {
System.out.println("reflection step");
reflection_sample = getSample(reflection);
if (simplex[0].fx < reflection_sample.fx && reflection_sample.fx < simplex[m - 1].fx) {
simplex[m - 1] = reflection_sample;
continue;
}
}
//expand
if (feasible(reflection) && simplex[0].fx >= reflection_sample.fx) {
double expansion[] = new double[n];
for (int i = 0; i < n; i++) {
expansion[i] = centroid[i] + gamma * (centroid[i] - simplex[m - 1].x[i]);
}
System.out.println("expansion step");
Sample expansion_sample = getSample(expansion);
if (feasible(expansion) && expansion_sample.fx < reflection_sample.fx) {
simplex[m - 1] = expansion_sample;
} else {
simplex[m - 1] = reflection_sample;
}
continue;
}
//contraction
if (!feasible(reflection) || simplex[m - 1].fx <= reflection_sample.fx) {
double contraction[] = new double[n];
for (int i = 0; i < n; i++) {
contraction[i] = centroid[i] + rho * (centroid[i] - simplex[m - 1].x[i]);
}
System.out.println("contraction step");
//this should not happen ..
Sample contraction_sample = null;
if (!feasible(contraction)) {
System.out.println("not feasible after contraction step");
contraction_sample = getSample(randomSample());
} else {
contraction_sample = getSample(contraction);
}
if (contraction_sample.fx < simplex[m - 1].fx) {
simplex[m - 1] = contraction_sample;
continue;
}
}
//shrink
for (int i = 1; i < m; i++) {
double shrink[] = new double[n];
for (int j = 0; j < n; j++) {
shrink[j] = simplex[0].x[j] + sigma * (simplex[i].x[j] - simplex[0].x[j]);
}
System.out.println("shrink step");
simplex[i] = getSample(shrink);
}
}
}
}
© 2015 - 2025 Weber Informatics LLC | Privacy Policy