GopherProxy

	tdownhill_simplex.cpp - numeric - C++ library with numerical algorithms
	git clone git://src.adamsgaard.dk/numeric
	Log
	Files
	Refs
	LICENSE
	---
	tdownhill_simplex.cpp (3700B)
	---
	1 #include <iostream>
	2 #include <armadillo>
	3 #include <vector>
	4 #include <functional>
	5 #include "downhill_simplex.h"
	6
	7 using namespace std;
	8 using namespace arma;
	9
	10 // Amoeba constructor
	11 amoeba::amoeba(function<double(vec)> fun, vector<vec> simplex)
	12 : d(simplex.size()-1), f(fun), value(zeros<vec>(d)), p_ce(zeros<vec>(d…
	13 {
	14 p_ce_o = p_ce;
	15
	16 for (int i=0; i<d+1; ++i)
	17 p.push_back(simplex[i]);
	18 for (int i=0; i<d+1; ++i)
	19 value[i] = f(p[i]);
	20
	21 update();
	22
	23 }
	24
	25 // Update amoeba parameters
	26 void amoeba::update()
	27 {
	28
	29 p_ce_o = p_ce;
	30
	31 // Find highest point
	32 hi = 0;
	33 for (int i=1; i<d+1; ++i)
	34 if (value[i] > value[hi])
	35 hi = i;
	36
	37 // Find lowest point
	38 lo = 0;
	39 for (int i=1; i<d+1; ++i)
	40 if (value[i] < value[lo])
	41 lo = i;
	42
	43 // Find centroid
	44 p_ce = zeros<vec>(d);
	45 for (int i=0; i<d+1; ++i)
	46 if (i != hi)
	47 p_ce += p[i];
	48 p_ce /= d;
	49
	50 // Write amoeba position to stderr
	51 pos();
	52 }
	53
	54 // Find size of simplex
	55 double amoeba::size()
	56 {
	57 double s = 0;
	58 for (int i=0; i<d+1; ++i)
	59 if (i != lo) {
	60 double n = norm(p[i] - p[lo], 2.0f);
	61 if (n>s)
	62 s=n;
	63 }
	64 return s;
	65 }
	66
	67 void amoeba::downhill(double simplex_size_goal)
	68 {
	69 // Find operation to update position
	70 while (size() > simplex_size_goal) {
	71
	72 vec p_re = p_ce + (p_ce - p[hi]); // Try reflection
	73 double f_re = f(p_re);
	74 if (f_re < value[lo]) {
	75 vec p_ex = p_ce + 1.0f * (p_ce - p[hi]); // Try expansion
	76 double f_ex = f(p_ex);
	77 if (f_ex < f_re) {
	78 value[hi] = f_ex;
	79 p[hi] = p_ex; // Accept expansion
	80 update(); continue; // Start loop over
	81 }
	82 }
	83
	84 if (f_re < value[hi]) {
	85 value[hi] = f_re;
	86 p[hi] = p_re; // Accept reflection
	87 update(); continue; // Start loop over
	88 }
	89
	90 vec p_co = p_ce + 0.5f * (p[hi] - p_ce); // Try contraction
	91 double f_co = f(p_co);
	92 if (f_co < value[hi]) {
	93 value[hi] = f_co;
	94 p[hi] = p_co; // Accept contraction
	95 update(); continue; // Start loop over
	96 }
	97
	98 for (int i=0; i<d+1; ++i)
	99 if (i != lo) {
	100 p[i] = 0.5f * (p[i] + p[lo]); // Do reduction
	101 value[i] = f(p[i]);
	102 }
	103 update(); continue;
	104 }
	105 }
	106
	107 void amoeba::downhill_mod(double simplex_size_goal)
	108 {
	109 // Find operation to update position
	110 while (size() > simplex_size_goal) {
	111
	112 vec p_re = p_ce + (p_ce - p[hi]); // Try reflection
	113
	114 double f_re = f(p_re);
	115 if (f_re < value[lo]) {
	116 vec p_ex_n = p_ce + 1.0f * (p_ce - p[hi]); // Try expansion
	117 double f_ex_n = f(p_ex_n);
	118 vec p_ex_o = p_ce_o + 1.0f * (p_ce_o - p[hi]); // Try expansion, o…
	119 double f_ex_o = f(p_ex_o);
	120
	121 // Find out which expansion gave the lowest value
	122 vec p_ex; double f_ex;
	123 if (f_ex_n < f_ex_o) { // New val
	124 p_ex = p_ex_n;
	125 f_ex = f_ex_n;
	126 } else { // Old val
	127 p_ex = p_ex_o;
	128 f_ex = f_ex_o;
	129 }
	130
	131 if (f_ex < f_re) {
	132 value[hi] = f_ex;
	133 p[hi] = p_ex; // Accept expansion
	134 update(); continue; // Start loop over
	135 }
	136 }
	137
	138 if (f_re < value[hi]) {
	139 value[hi] = f_re;
	140 p[hi] = p_re; // Accept reflection
	141 update(); continue; // Start loop over
	142 }
	143
	144 vec p_co = p_ce + 0.5f * (p[hi] - p_ce); // Try contraction
	145 double f_co = f(p_co);
	146 if (f_co < value[hi]) {
	147 value[hi] = f_co;
	148 p[hi] = p_co; // Accept contraction
	149 update(); continue; // Start loop over
	150 }
	151
	152 for (int i=0; i<d+1; ++i)
	153 if (i != lo) {
	154 p[i] = 0.5f * (p[i] + p[lo]); // Do reduction
	155 value[i] = f(p[i]);
	156 }
	157 update(); continue;
	158 }
	159 }
	160
	161 // Write simplex points to stderr
	162 void amoeba::pos()
	163 {
	164 for (int i=0; i<d+1; ++i)
	165 std::cerr << p[i][0] << '\t'
	166 << p[i][1] << '\n';
	167 }
	168
	169 // Return lowest point of simplex
	170 vec amoeba::low()
	171 {
	172 return p[lo];
	173 }