tdownhill_simplex.cpp - numeric - C++ library with numerical algorithms
HTML git clone git://src.adamsgaard.dk/numeric
DIR Log
DIR Files
DIR Refs
DIR 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, old val
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 }