minimise.simplex

44 - def __init__(self, func, args, x0, func_tol, maxiter, full_output, print_flag, print_prefix):

45 """Class for downhill simplex minimisation specific functions. 46 47 Unless you know what you are doing, you should call the function 'simplex' rather than using 48 this class. 49 """ 50 51 # Function arguments. 52 self.func = func 53 self.args = args 54 self.xk = x0 55 self.func_tol = func_tol 56 self.maxiter = maxiter 57 self.full_output = full_output 58 self.print_flag = print_flag 59 self.print_prefix = print_prefix 60 61 # Initialise the function, gradient, and Hessian evaluation counters. 62 self.f_count = 0 63 self.g_count = 0 64 self.h_count = 0 65 66 # Initialise the warning string. 67 self.warning = None 68 69 # Initialise some constants. 70 self.n = len(self.xk) 71 self.m = self.n + 1 72 73 # Create the simplex 74 self.simplex = zeros((self.m, self.n), Float64) 75 self.simplex_vals = zeros(self.m, Float64) 76 77 self.simplex[0] = self.xk * 1.0 78 self.simplex_vals[0], self.f_count = apply(self.func, (self.xk,)+self.args), self.f_count + 1 79 80 for i in xrange(self.n): 81 j = i + 1 82 self.simplex[j] = self.xk 83 if self.xk[i] == 0.0: 84 self.simplex[j, i] = 2.5 * 1e-4 85 else: 86 self.simplex[j, i] = 1.05 * self.simplex[j, i] 87 self.simplex_vals[j], self.f_count = apply(self.func, (self.simplex[j],)+self.args), self.f_count + 1 88 89 # Order the simplex. 90 self.order_simplex() 91 92 # Set xk and fk as the vertex of the simplex with the lowest function value. 93 self.xk = self.simplex[0] * 1.0 94 self.fk = self.simplex_vals[0] 95 96 # Find the center of the simplex. 97 self.center = average(self.simplex)

98 99

100 - def new_param_func(self):

101 """The new parameter function. 102 103 Simplex movement. 104 """ 105 106 self.reflect_flag = 1 107 self.shrink_flag = 0 108 109 self.pivot_point = average(self.simplex[:-1]) 110 111 self.reflect() 112 if self.reflect_val <= self.simplex_vals[0]: 113 self.extend() 114 elif self.reflect_val >= self.simplex_vals[-2]: 115 self.reflect_flag = 0 116 if self.reflect_val < self.simplex_vals[-1]: 117 self.contract() 118 else: 119 self.contract_orig() 120 if self.reflect_flag: 121 self.simplex[-1], self.simplex_vals[-1] = self.reflect_vector, self.reflect_val 122 if self.shrink_flag: 123 self.shrink() 124 125 self.order_simplex() 126 127 # Update values. 128 self.xk_new = self.simplex[0] 129 self.fk_new = self.simplex_vals[0] 130 self.dfk_new = None 131 132 # Find the center of the simplex and calculate the distance moved. 133 self.center_new = average(self.simplex) 134 self.dist = sum(abs(self.center_new - self.center))

135 136

137 - def contract(self):

138 """Contraction step.""" 139 140 self.contract_vector = 1.5 * self.pivot_point - 0.5 * self.simplex[-1] 141 self.contract_val, self.f_count = apply(self.func, (self.contract_vector,)+self.args), self.f_count + 1 142 if self.contract_val < self.reflect_val: 143 self.simplex[-1], self.simplex_vals[-1] = self.contract_vector, self.contract_val 144 else: 145 self.shrink_flag = 1

146 147

148 - def contract_orig(self):

149 """Contraction of the original simplex.""" 150 151 self.contract_orig_vector = 0.5 * (self.pivot_point + self.simplex[-1]) 152 self.contract_orig_val, self.f_count = apply(self.func, (self.contract_orig_vector,)+self.args), self.f_count + 1 153 if self.contract_orig_val < self.simplex_vals[-1]: 154 self.simplex[-1], self.simplex_vals[-1] = self.contract_orig_vector, self.contract_orig_val 155 else: 156 self.shrink_flag = 1

157 158

159 - def extend(self):

160 """Extension step.""" 161 162 self.extend_vector = 3.0 * self.pivot_point - 2.0 * self.simplex[-1] 163 self.extend_val, self.f_count = apply(self.func, (self.extend_vector,)+self.args), self.f_count + 1 164 if self.extend_val < self.reflect_val: 165 self.simplex[-1], self.simplex_vals[-1] = self.extend_vector, self.extend_val 166 self.reflect_flag = 0

167 168

169 - def order_simplex(self):

170 """Order the vertecies of the simplex according to accending function values.""" 171 sorted = argsort(self.simplex_vals) 172 self.simplex = take(self.simplex, sorted) 173 self.simplex_vals = take(self.simplex_vals, sorted)

174 175

176 - def reflect(self):

177 """Reflection step.""" 178 179 self.reflect_vector = 2.0 * self.pivot_point - self.simplex[-1] 180 self.reflect_val, self.f_count = apply(self.func, (self.reflect_vector,)+self.args), self.f_count + 1

181 182

183 - def shrink(self):

184 """Shrinking step.""" 185 186 for i in xrange(self.n): 187 j = i + 1 188 self.simplex[j] = 0.5 * (self.simplex[0] + self.simplex[j]) 189 self.simplex_vals[j], self.f_count = apply(self.func, (self.simplex[j],)+self.args), self.f_count + 1

190 191

192 - def conv_test(self, *args):

193 """Convergence test. 194 195 Finish minimising when the function difference between the highest and lowest simplex 196 vertecies is insignificant or if the simplex doesn't move. 197 """ 198 199 if self.print_flag >= 2: 200 print self.print_prefix + "diff = " + `self.simplex_vals[-1] - self.simplex_vals[0]` 201 print self.print_prefix + "|diff| = " + `abs(self.simplex_vals[-1] - self.simplex_vals[0])` 202 print self.print_prefix + "f_tol = " + `self.func_tol` 203 print self.print_prefix + "center = " + `self.pivot_point` 204 try: 205 print self.print_prefix + "old center = " + `self.old_pivot` 206 print self.print_prefix + "center diff = " + `self.pivot_point - self.old_pivot` 207 except AttributeError: 208 pass 209 self.old_pivot = 1.0 * self.pivot_point 210 if abs(self.simplex_vals[-1] - self.simplex_vals[0]) <= self.func_tol: 211 if self.print_flag >= 2: 212 print "\n" + self.print_prefix + "???Function tolerance reached." 213 print self.print_prefix + "simplex_vals[-1]: " + `self.simplex_vals[-1]` 214 print self.print_prefix + "simplex_vals[0]: " + `self.simplex_vals[0]` 215 print self.print_prefix + "|diff|: " + `abs(self.simplex_vals[-1] - self.simplex_vals[0])` 216 print self.print_prefix + "tol: " + `self.func_tol` 217 self.xk_new = self.simplex[0] 218 self.fk_new = self.simplex_vals[0] 219 return 1 220 221 # Test if simplex has not moved. 222 if self.dist == 0.0: 223 self.warning = "Simplex has not moved." 224 self.xk_new = self.simplex[0] 225 self.fk_new = self.simplex_vals[0] 226 return 1

227 228

229 - def update(self):

230 """Update function.""" 231 232 self.xk = self.xk_new * 1.0 233 self.fk = self.fk_new 234 self.center = self.center_new * 1.0

Source Code for Module minimise.simplex