minfx.grid

1 ############################################################################### 2 # # 3 # Copyright (C) 2003-2013 Edward d'Auvergne # 4 # # 5 # This file is part of the minfx optimisation library, # 6 # https://sourceforge.net/projects/minfx # 7 # # 8 # This program is free software: you can redistribute it and/or modify # 9 # it under the terms of the GNU General Public License as published by # 10 # the Free Software Foundation, either version 3 of the License, or # 11 # (at your option) any later version. # 12 # # 13 # This program is distributed in the hope that it will be useful, # 14 # but WITHOUT ANY WARRANTY; without even the implied warranty of # 15 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # 16 # GNU General Public License for more details. # 17 # # 18 # You should have received a copy of the GNU General Public License # 19 # along with this program. If not, see <http://www.gnu.org/licenses/>. # 20 # # 21 ############################################################################### 22 23 # Module docstring. 24 """The grid search algorithm. 25 26 This file is part of the minfx optimisation library at U{https://sourceforge.net/projects/minfx}. 27 """ 28 29 # Python module imports. 30 from numpy import array, float64, ones, zeros 31 32 # Minfx module imports. 33 from minfx.base_classes import print_iter 34 from minfx.constraint_linear import Constraint_linear 35 from minfx.errors import MinfxError 36 37

38 -def grid(func, args=(), num_incs=None, lower=None, upper=None, incs=None, A=None, b=None, l=None, u=None, c=None, verbosity=0, print_prefix=""):

39 """The grid search algorithm. 40 41 @param func: The target function. This should take the parameter vector as the first argument and return a single float. 42 @type func: function 43 @keyword args: A tuple of arguments to pass to the function, if needed. 44 @type args: tuple 45 @keyword num_incs: The number of linear increments to be used in the grid search. The length should be equal to the number of parameters and each element corresponds to the number of increments for the respective parameter. This is overridden if the incs argument is supplied. 46 @type num_incs: list of int 47 @keyword lower: The list of lower bounds for the linear grid search. This must be supplied if incs is not. 48 @type lower: list of float 49 @keyword upper: The list of upper bounds for the linear grid search. This must be supplied if incs is not. 50 @type upper: list of float 51 @keyword incs: The parameter increment values. This overrides the num_incs, lower, and upper arguments used in generating a linear grid. 52 @type incs: list of lists 53 @keyword A: The linear constraint matrix A, such that A.x >= b. 54 @type A: numpy rank-2 array 55 @keyword b: The linear constraint scalar vectors, such that A.x >= b. 56 @type b: numpy rank-1 array 57 @keyword l: The lower bound constraint vector, such that l <= x <= u. 58 @type l: list of float 59 @keyword u: The upper bound constraint vector, such that l <= x <= u. 60 @type u: list of float 61 @keyword c: A user supplied constraint function. 62 @type c: function 63 @keyword verbosity: The verbosity level. 0 corresponds to no output, 1 is standard, and higher values cause greater and greater amount of output. 64 @type verbosity: int 65 @keyword print_prefix: The text to place before the printed output. 66 @type print_prefix: str 67 @return: The optimisation information including the parameter vector at the best grid point, the function value at this grid point, the number of iterations (equal to the number of function calls), and a warning. 68 @rtype: tuple of numpy rank-1 array, float, int, str 69 """ 70 71 # Checks. 72 if num_incs == None and incs == None: 73 raise MinfxError("Either the incs arg or the num_incs, lower, and upper args must be supplied.") 74 elif num_incs: 75 # Check that this is a list. 76 if not isinstance(num_incs, list): 77 raise MinfxError("The num_incs argument '%s' must be a list." % num_incs) 78 if not isinstance(lower, list): 79 raise MinfxError("The lower argument '%s' must be a list." % lower) 80 if not isinstance(upper, list): 81 raise MinfxError("The upper argument '%s' must be a list." % upper) 82 83 # Lengths. 84 if len(num_incs) != len(lower): 85 raise MinfxError("The '%s' num_incs and '%s' lower arguments are of different lengths" % (num_incs, lower)) 86 if len(num_incs) != len(upper): 87 raise MinfxError("The '%s' num_incs and '%s' upper arguments are of different lengths" % (num_incs, upper)) 88 89 # Catch models with zero parameters. 90 if num_incs == [] or incs == []: 91 # Print out. 92 if verbosity: 93 print("Cannot run a grid search on a model with zero parameters, directly calculating the function value.") 94 95 # Empty parameter vector. 96 x0 = zeros(0, float64) 97 98 # The function value. 99 fk = func(x0) 100 101 # The results tuple. 102 return x0, fk, 1, "No optimisation" 103 104 105 # Initialise. 106 if num_incs: 107 n = len(num_incs) 108 else: 109 n = len(incs) 110 grid_size = 0 111 total_steps = 1 112 step_num = ones(n, int) 113 params = zeros((n), float64) 114 min_params = zeros((n), float64) 115 116 # Linear grid search. 117 # The incs data structure eliminates the round-off error of summing a step size value to the parameter value. 118 if num_incs: 119 incs = [] 120 121 # Loop over the dimensions. 122 for k in range(n): 123 params[k] = lower[k] 124 min_params[k] = lower[k] 125 total_steps = total_steps * num_incs[k] 126 incs.append([]) 127 128 # Loop over the increments of dimension k. 129 for i in range(num_incs[k]): 130 # Single grid search increment in dimension k, so use the average of the lower and upper. 131 if num_incs[k] == 1: 132 incs[k].append((lower[k] + upper[k]) / 2.0) 133 134 # More than 1 increment. 135 else: 136 incs[k].append(lower[k] + i * (upper[k] - lower[k]) / (num_incs[k] - 1)) 137 138 # User supplied grid search. 139 else: 140 for k in range(n): 141 total_steps = total_steps * len(incs[k]) 142 143 # Print out. 144 if verbosity: 145 if verbosity >= 2: 146 print(print_prefix) 147 print(print_prefix) 148 print(print_prefix + "Grid search") 149 print(print_prefix + "~~~~~~~~~~~") 150 151 # Linear constraints. 152 if A != None and b != None: 153 constraint_flag = 1 154 constraint_linear = Constraint_linear(A, b) 155 c = constraint_linear.func 156 if verbosity >= 3: 157 print(print_prefix + "Linear constraint matrices.") 158 print(print_prefix + "A: " + repr(A)) 159 print(print_prefix + "b: " + repr(b)) 160 161 # Bound constraints. 162 elif l != None and u != None: 163 constraint_flag = 1 164 raise MinfxError("Bound constraints are not implemented yet.") 165 166 # General constraints. 167 elif c != None: 168 constraint_flag = 1 169 170 # No constraints. 171 else: 172 constraint_flag = 0 173 174 # Set a ridiculously large initial grid value. 175 f_min = 1e300 176 177 # Initial printout. 178 if verbosity: 179 print("\n" + print_prefix + "Searching through %s grid nodes." % total_steps) 180 181 # Test if the grid is too large. 182 if total_steps >= 1e8: 183 raise MinfxError("A grid search of size %s is too large." % total_steps) 184 185 # Search the grid. 186 k = 0 187 for i in range(total_steps): 188 # Check that the grid point does not violate a constraint, and if it does, skip the function call. 189 skip = 0 190 if constraint_flag: 191 ci = c(params) 192 if min(ci) < 0.0: 193 if verbosity >= 3: 194 print_iter(k, min_params, print_prefix=print_prefix) 195 print(print_prefix + "Constraint violated, skipping grid point.") 196 print(print_prefix + "ci: " + repr(ci)) 197 print("") 198 skip = 1 199 200 # Function call, test, and increment grid_size. 201 if not skip: 202 # Back calculate the current function value. 203 f = func(*(params,)+args) 204 205 # Test if the current function value is less than the least function value. 206 if f < f_min: 207 f_min = f 208 min_params = 1.0 * params 209 210 # Print out code. 211 if verbosity: 212 print_iter(k=k, xk=min_params, fk=f_min, print_prefix=print_prefix) 213 214 # Grid count. 215 grid_size = grid_size + 1 216 217 # Print out code. 218 if verbosity >= 2: 219 if f != f_min: 220 print_iter(k=k, xk=min_params, fk=f, print_prefix=print_prefix) 221 if verbosity >= 3: 222 print(print_prefix + "%-20s%-20s" % ("Increment:", repr(step_num))) 223 print(print_prefix + "%-20s%-20s" % ("Params:", repr(params))) 224 print(print_prefix + "%-20s%-20s" % ("Min params:", repr(min_params))) 225 print(print_prefix + "%-20s%-20g\n" % ("f:", f)) 226 print(print_prefix + "%-20s%-20g\n" % ("Min f:", f_min)) 227 228 # Increment k. 229 k = k + 1 230 231 # Increment the grid search. 232 for j in range(n): 233 if step_num[j] < len(incs[j]): 234 step_num[j] = step_num[j] + 1 235 params[j] = incs[j][step_num[j]-1] 236 break # Exit so that the other step numbers are not incremented. 237 else: 238 step_num[j] = 1 239 params[j] = incs[j][0] 240 241 # Return the results. 242 return min_params, f_min, grid_size, None

243 244

245 -def grid_point_array(func, args=(), points=None, verbosity=0, print_prefix=""):

246 """The grid search algorithm. 247 248 @param func: The target function. This should take the parameter vector as the first argument and return a single float. 249 @type func: function 250 @keyword args: A tuple of arguments to pass to the function, if needed. 251 @type args: tuple 252 @keyword points: The array of grid points to search over. 253 @type points: list or array of lists of float 254 @keyword verbosity: The verbosity level. 0 corresponds to no output, 1 is standard, and higher values cause greater and greater amount of output. 255 @type verbosity: int 256 @keyword print_prefix: The text to place before the printed output. 257 @type print_prefix: str 258 @return: The optimisation information including the parameter vector at the best grid point, the function value at this grid point, the number of iterations (equal to the number of function calls), and a warning. 259 @rtype: tuple of numpy rank-1 array, float, int, str 260 """ 261 262 # Initialise. 263 total_steps = len(points) 264 n = len(points[0]) 265 266 # Set a ridiculously large initial function value. 267 f_min = 1e300 268 269 # Initial printout. 270 if verbosity: 271 print("\n" + print_prefix + "Searching through %s grid nodes." % total_steps) 272 273 # Test if the grid is too large. 274 if total_steps >= 1e8: 275 raise MinfxError("A grid search of size %s is too large." % total_steps) 276 277 # Search the grid. 278 for k in range(total_steps): 279 # Back calculate the current function value. 280 f = func(*(points[k],)+args) 281 282 # A better point. 283 if f < f_min: 284 # Switch to the new point. 285 f_min = f 286 min_params = 1.0 * points[k] 287 288 # Print out code. 289 if verbosity: 290 print_iter(k=k, xk=min_params, fk=f_min, print_prefix=print_prefix) 291 292 # Print out code. 293 if verbosity >= 2: 294 if f != f_min: 295 print_iter(k=k, xk=min_params, fk=f, print_prefix=print_prefix) 296 if verbosity >= 3: 297 print(print_prefix + "%-20s%-20s" % ("Params:", repr(points[k]))) 298 print(print_prefix + "%-20s%-20s" % ("Min params:", repr(min_params))) 299 print(print_prefix + "%-20s%-20g\n" % ("f:", f)) 300 print(print_prefix + "%-20s%-20g\n" % ("Min f:", f_min)) 301 302 # Return the results. 303 return min_params, f_min, total_steps, None

304 305

306 -def grid_split(divisions=None, lower=None, upper=None, inc=None, A=None, b=None, l=None, u=None, c=None):

307 """Generator method yielding arrays of grid points. 308 309 This method will loop over the grid points one-by-one, generating a list of points and yielding these for each subdivision. 310 311 312 @keyword divisions: The number of grid subdivisions. 313 @type divisions: int 314 @keyword lower: The lower bounds of the grid search which must be equal to the number of parameters in the model. 315 @type lower: array of numbers 316 @keyword upper: The upper bounds of the grid search which must be equal to the number of parameters in the model. 317 @type upper: array of numbers 318 @keyword inc: The increments for each dimension of the space for the grid search. The number of elements in the array must equal to the number of parameters in the model. 319 @type inc: array of int 320 @keyword A: The linear constraint matrix A, such that A.x >= b. 321 @type A: numpy rank-2 array 322 @keyword b: The linear constraint scalar vectors, such that A.x >= b. 323 @type b: numpy rank-1 array 324 @keyword l: The lower bound constraint vector, such that l <= x <= u. 325 @type l: list of float 326 @keyword u: The upper bound constraint vector, such that l <= x <= u. 327 @type u: list of float 328 @keyword c: A user supplied constraint function. 329 @type c: function 330 @return: A list of grid points for each subdivision is yielded. 331 @rtype: list of list of float 332 """ 333 334 # The dimensionality. 335 n = len(inc) 336 337 # Linear constraints. 338 if A != None and b != None: 339 constraint_flag = True 340 constraint_linear = Constraint_linear(A, b) 341 c = constraint_linear.func 342 if verbosity >= 3: 343 print(print_prefix + "Linear constraint matrices.") 344 print(print_prefix + "A: " + repr(A)) 345 print(print_prefix + "b: " + repr(b)) 346 347 # Bound constraints. 348 elif l != None and u != None: 349 constraint_flag = True 350 raise MinfxError("Bound constraints are not implemented yet.") 351 352 # General constraints. 353 elif c != None: 354 constraint_flag = True 355 356 # No constraints. 357 else: 358 constraint_flag = False 359 360 # Total number of points. 361 total_pts = 1 362 for i in range(n): 363 total_pts = total_pts * inc[i] 364 365 # Init. 366 indices = zeros(n, int) 367 pts = zeros((total_pts, n), float64) 368 369 # The increment values for each dimension. 370 incs = [] 371 for k in range(n): 372 incs.append([]) 373 for i in range(inc[k]): 374 incs[k].append(lower[k] + i * (upper[k] - lower[k]) / (inc[k] - 1)) 375 376 # Construct the list of all grid points. 377 for i in range(total_pts): 378 # Loop over the dimensions. 379 for j in range(n): 380 # Add the point coordinate. 381 pts[i][j] = incs[j][indices[j]] 382 383 # Increment the step positions. 384 for j in range(n): 385 if indices[j] < inc[j]-1: 386 indices[j] += 1 387 break # Exit so that the other step numbers are not incremented. 388 else: 389 indices[j] = 0 390 391 # Eliminate points outside of constraints. 392 if constraint_flag: 393 # Loop over all points. 394 pts_trimmed = [] 395 for i in range(total_pts): 396 # The constraint values. 397 ci = c(pts[i]) 398 399 # No constraint violations, so store the point. 400 if min(ci) >= 0.0: 401 pts_trimmed.append(pts[i]) 402 403 # No point elimination. 404 else: 405 pts_trimmed = pts 406 407 # Convert to numpy. 408 pts_trimmed = array(pts_trimmed) 409 410 # The total number of points has changed. 411 total_pts = len(pts_trimmed) 412 413 # The subdivision size (round up so the last subdivision is smaller than the rest). 414 size_float = total_pts / float(divisions) 415 size = int(size_float) 416 if size_float % 1: 417 size = size + 1 418 419 # Subdivide. 420 for i in range(min(divisions, total_pts)): 421 # The start index. 422 start = i * size 423 424 # The end index. 425 if i != divisions - 1: 426 end = (i + 1) * size 427 else: 428 end = total_pts 429 430 # Yield the subdivision. 431 yield pts[start: end]

432

Source Code for Module minfx.grid