minfx.grid

1 ############################################################################### 2 # # 3 # Copyright (C) 2003-2014 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 U{minfx optimisation library<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, sparseness=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 sparseness: An optional argument for defining sparsity, or regions of the grid to not search over. This is a list whereby each element is a list of two parameters indices. These two parameters will then be assumed to be decoupled and the grid search space between them will be skipped. Symmetry is not observed, so if [0, 1] is sent in, then maybe [1, 0] should be as well. 64 @type sparseness: list of list of int 65 @keyword verbosity: The verbosity level. 0 corresponds to no output, 1 is standard, and higher values cause greater and greater amount of output. 66 @type verbosity: int 67 @keyword print_prefix: The text to place before the printed output. 68 @type print_prefix: str 69 @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. 70 @rtype: tuple of numpy rank-1 array, float, int, str 71 """ 72 73 # Checks. 74 if num_incs == None and incs == None: 75 raise MinfxError("Either the incs arg or the num_incs, lower, and upper args must be supplied.") 76 elif num_incs: 77 # Check that this is a list. 78 if not isinstance(num_incs, list): 79 raise MinfxError("The num_incs argument '%s' must be a list." % num_incs) 80 if not isinstance(lower, list): 81 raise MinfxError("The lower argument '%s' must be a list." % lower) 82 if not isinstance(upper, list): 83 raise MinfxError("The upper argument '%s' must be a list." % upper) 84 85 # Lengths. 86 if len(num_incs) != len(lower): 87 raise MinfxError("The '%s' num_incs and '%s' lower arguments are of different lengths" % (num_incs, lower)) 88 if len(num_incs) != len(upper): 89 raise MinfxError("The '%s' num_incs and '%s' upper arguments are of different lengths" % (num_incs, upper)) 90 91 # Catch models with zero parameters. 92 if num_incs == [] or incs == []: 93 # Print out. 94 if verbosity: 95 print("Cannot run a grid search on a model with zero parameters, directly calculating the function value.") 96 97 # Empty parameter vector. 98 x0 = zeros(0, float64) 99 100 # The function value. 101 fk = func(x0) 102 103 # The results tuple. 104 return x0, fk, 1, "No optimisation" 105 106 107 # Initialise. 108 if num_incs: 109 n = len(num_incs) 110 else: 111 n = len(incs) 112 grid_size = 0 113 total_steps = 1 114 step_num = ones(n, int) 115 min_step = ones(n, int) 116 params = zeros((n), float64) 117 min_params = zeros((n), float64) 118 sparseness_flag = False 119 if sparseness != None and len(sparseness): 120 sparseness_flag = True 121 122 # Linear grid search. 123 # The incs data structure eliminates the round-off error of summing a step size value to the parameter value. 124 if num_incs: 125 incs = [] 126 127 # Loop over the dimensions. 128 for k in range(n): 129 params[k] = lower[k] 130 min_params[k] = lower[k] 131 total_steps = total_steps * num_incs[k] 132 incs.append([]) 133 134 # Loop over the increments of dimension k. 135 for i in range(num_incs[k]): 136 # Single grid search increment in dimension k, so use the average of the lower and upper. 137 if num_incs[k] == 1: 138 incs[k].append((lower[k] + upper[k]) / 2.0) 139 140 # More than 1 increment. 141 else: 142 incs[k].append(lower[k] + i * (upper[k] - lower[k]) / (num_incs[k] - 1)) 143 144 # User supplied grid search. 145 else: 146 for k in range(n): 147 total_steps = total_steps * len(incs[k]) 148 149 # Print out. 150 if verbosity: 151 if verbosity >= 2: 152 print(print_prefix) 153 print(print_prefix) 154 print(print_prefix + "Grid search") 155 print(print_prefix + "~~~~~~~~~~~") 156 157 # Linear constraints. 158 if A != None and b != None: 159 constraint_flag = 1 160 constraint_linear = Constraint_linear(A, b) 161 c = constraint_linear.func 162 if verbosity >= 3: 163 print(print_prefix + "Linear constraint matrices.") 164 print(print_prefix + "A: " + repr(A)) 165 print(print_prefix + "b: " + repr(b)) 166 167 # Bound constraints. 168 elif l != None and u != None: 169 constraint_flag = 1 170 raise MinfxError("Bound constraints are not implemented yet.") 171 172 # General constraints. 173 elif c != None: 174 constraint_flag = 1 175 176 # No constraints. 177 else: 178 constraint_flag = 0 179 180 # Set a ridiculously large initial grid value. 181 f_min = 1e300 182 183 # Initial printout. 184 if verbosity: 185 print("\n" + print_prefix + "Searching through %s grid nodes." % total_steps) 186 187 # Test if the grid is too large. 188 if total_steps >= 1e8: 189 raise MinfxError("A grid search of size %s is too large." % total_steps) 190 191 # Search the grid. 192 k = 0 193 curr_dim = 0 194 for i in range(total_steps): 195 # The flag for skipping the grid point if outside a constraint or in a sparse region. 196 skip = False 197 198 # Check that the grid point does not violate a constraint, and if it does, skip the function call. 199 if constraint_flag: 200 ci = c(params) 201 if min(ci) < 0.0: 202 if verbosity >= 3: 203 print_iter(k, min_params, print_prefix=print_prefix) 204 print(print_prefix + "Constraint violated, skipping grid point.") 205 print(print_prefix + "ci: " + repr(ci)) 206 print("") 207 skip = True 208 209 # Check if in a sparse region to skip. 210 if sparseness_flag: 211 # Loop over the sparseness blocks. 212 for block in sparseness: 213 # The currently incremented parameter is in the block. 214 if curr_dim == block[0]: 215 if step_num[block[1]] != min_step[block[1]]: 216 skip = True 217 break 218 219 # Function call, test, and increment grid_size. 220 min_found = False 221 if not skip: 222 # Back calculate the current function value. 223 f = func(*(params,)+args) 224 225 # Test if the current function value is less than the least function value. 226 if f < f_min: 227 f_min = f 228 min_params = 1.0 * params 229 min_step = 1.0 * step_num 230 min_found = True 231 232 # Print out code. 233 if verbosity: 234 print_iter(k=k, xk=min_params, fk=f_min, print_prefix=print_prefix) 235 236 # Grid count. 237 grid_size = grid_size + 1 238 239 # Print out code. 240 if verbosity >= 2: 241 if not min_found: 242 print_iter(k=k, xk=min_params, fk=f, print_prefix=print_prefix) 243 if verbosity >= 3: 244 if min_found: 245 print(print_prefix + "Minimum found.") 246 print(print_prefix + "%-20s%-20s" % ("Increment:", step_num)) 247 print(print_prefix + "%-20s%-20s" % ("Current dimension:", curr_dim)) 248 print(print_prefix + "%-20s%-20s" % ("Params:", params)) 249 print(print_prefix + "%-20s%-20s" % ("Min params:", min_params)) 250 print(print_prefix + "%-20s%-20s" % ("Min indices:", min_step)) 251 print(print_prefix + "%-20s%-20.8g" % ("f:", f)) 252 print(print_prefix + "%-20s%-20.8g\n" % ("Min f:", f_min)) 253 254 # Increment k. 255 k = k + 1 256 257 # Increment the grid search. 258 for j in range(n): 259 if step_num[j] < len(incs[j]): 260 step_num[j] = step_num[j] + 1 261 params[j] = incs[j][step_num[j]-1] 262 if j > curr_dim: 263 curr_dim = j 264 break # Exit so that the other step numbers are not incremented. 265 else: 266 step_num[j] = 1 267 params[j] = incs[j][0] 268 269 # Return the results. 270 return min_params, f_min, grid_size, None

271 272

273 -def grid_point_array(func, args=(), points=None, A=None, b=None, l=None, u=None, c=None, verbosity=0, print_prefix=""):

274 """The grid search algorithm. 275 276 @param func: The target function. This should take the parameter vector as the first argument and return a single float. 277 @type func: function 278 @keyword args: A tuple of arguments to pass to the function, if needed. 279 @type args: tuple 280 @keyword points: The array of grid points to search over. 281 @type points: list or array of lists of float 282 @keyword A: The linear constraint matrix A, such that A.x >= b. 283 @type A: numpy rank-2 array 284 @keyword b: The linear constraint scalar vectors, such that A.x >= b. 285 @type b: numpy rank-1 array 286 @keyword l: The lower bound constraint vector, such that l <= x <= u. 287 @type l: list of float 288 @keyword u: The upper bound constraint vector, such that l <= x <= u. 289 @type u: list of float 290 @keyword c: A user supplied constraint function. 291 @type c: function 292 @keyword verbosity: The verbosity level. 0 corresponds to no output, 1 is standard, and higher values cause greater and greater amount of output. 293 @type verbosity: int 294 @keyword print_prefix: The text to place before the printed output. 295 @type print_prefix: str 296 @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. 297 @rtype: tuple of numpy rank-1 array, float, int, str 298 """ 299 300 # Initialise. 301 total_steps = len(points) 302 303 # Skip empty grids. 304 if total_steps == 0: 305 return None, None, None, None 306 307 # The dimensionality. 308 n = len(points[0]) 309 310 # Set a ridiculously large initial function value. 311 f_min = 1e300 312 313 # Initial printout. 314 if verbosity: 315 print("\n" + print_prefix + "Searching through %s grid nodes." % total_steps) 316 317 # Test if the grid is too large. 318 if total_steps >= 1e8: 319 raise MinfxError("A grid search of size %s is too large." % total_steps) 320 321 # Linear constraints. 322 if A != None and b != None: 323 constraint_linear = Constraint_linear(A, b) 324 c = constraint_linear.func 325 if verbosity >= 3: 326 print(print_prefix + "Linear constraint matrices.") 327 print(print_prefix + "A: " + repr(A)) 328 print(print_prefix + "b: " + repr(b)) 329 330 # Bound constraints. 331 elif l != None and u != None: 332 raise MinfxError("Bound constraints are not implemented yet.") 333 334 # The constraint flag. 335 constraint_flag = False 336 if c != None: 337 constraint_flag = True 338 339 # Search the grid. 340 for k in range(total_steps): 341 # If a grid point violates a constraint, skip it. 342 if constraint_flag: 343 ci = c(params) 344 if min(ci) < 0.0: 345 if verbosity >= 3: 346 print_iter(k, min_params, print_prefix=print_prefix) 347 print(print_prefix + "Constraint violated, skipping grid point.") 348 print(print_prefix + "ci: " + repr(ci)) 349 print("") 350 continue 351 352 # Back calculate the current function value. 353 f = func(*(points[k],)+args) 354 355 # A better point. 356 if f < f_min: 357 # Switch to the new point. 358 f_min = f 359 min_params = 1.0 * points[k] 360 361 # Print out code. 362 if verbosity: 363 print_iter(k=k, xk=min_params, fk=f_min, print_prefix=print_prefix) 364 365 # Print out code. 366 if verbosity >= 2: 367 if f != f_min: 368 print_iter(k=k, xk=min_params, fk=f, print_prefix=print_prefix) 369 if verbosity >= 3: 370 print(print_prefix + "%-20s%-20s" % ("Params:", repr(points[k]))) 371 print(print_prefix + "%-20s%-20s" % ("Min params:", repr(min_params))) 372 print(print_prefix + "%-20s%-20g\n" % ("f:", f)) 373 print(print_prefix + "%-20s%-20g\n" % ("Min f:", f_min)) 374 375 # Return the results. 376 return min_params, f_min, total_steps, None

377 378

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

380 """Generator method yielding arrays of grid points. 381 382 This method will loop over the grid points one-by-one, generating a list of points and yielding these for each subdivision. 383 384 385 @keyword divisions: The number of grid subdivisions. 386 @type divisions: int 387 @keyword lower: The lower bounds of the grid search which must be equal to the number of parameters in the model. 388 @type lower: array of numbers 389 @keyword upper: The upper bounds of the grid search which must be equal to the number of parameters in the model. 390 @type upper: array of numbers 391 @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. 392 @type inc: array of int 393 @keyword A: The linear constraint matrix A, such that A.x >= b. 394 @type A: numpy rank-2 array 395 @keyword b: The linear constraint scalar vectors, such that A.x >= b. 396 @type b: numpy rank-1 array 397 @keyword l: The lower bound constraint vector, such that l <= x <= u. 398 @type l: list of float 399 @keyword u: The upper bound constraint vector, such that l <= x <= u. 400 @type u: list of float 401 @keyword c: A user supplied constraint function. 402 @type c: function 403 @return: A list of grid points for each subdivision is yielded. 404 @rtype: list of list of float 405 """ 406 407 # The dimensionality. 408 n = len(inc) 409 410 # Linear constraints. 411 if A != None and b != None: 412 constraint_flag = True 413 constraint_linear = Constraint_linear(A, b) 414 c = constraint_linear.func 415 if verbosity >= 3: 416 print(print_prefix + "Linear constraint matrices.") 417 print(print_prefix + "A: " + repr(A)) 418 print(print_prefix + "b: " + repr(b)) 419 420 # Bound constraints. 421 elif l != None and u != None: 422 constraint_flag = True 423 raise MinfxError("Bound constraints are not implemented yet.") 424 425 # General constraints. 426 elif c != None: 427 constraint_flag = True 428 429 # No constraints. 430 else: 431 constraint_flag = False 432 433 # Total number of points. 434 total_pts = 1 435 for i in range(n): 436 total_pts = total_pts * inc[i] 437 438 # Init. 439 indices = zeros(n, int) 440 pts = zeros((total_pts, n), float64) 441 442 # The increment values for each dimension. 443 incs = [] 444 for k in range(n): 445 incs.append([]) 446 for i in range(inc[k]): 447 incs[k].append(lower[k] + i * (upper[k] - lower[k]) / (inc[k] - 1)) 448 449 # Construct the list of all grid points. 450 for i in range(total_pts): 451 # Loop over the dimensions. 452 for j in range(n): 453 # Add the point coordinate. 454 pts[i][j] = incs[j][indices[j]] 455 456 # Increment the step positions. 457 for j in range(n): 458 if indices[j] < inc[j]-1: 459 indices[j] += 1 460 break # Exit so that the other step numbers are not incremented. 461 else: 462 indices[j] = 0 463 464 # Eliminate points outside of constraints. 465 if constraint_flag: 466 # Loop over all points. 467 pts_trimmed = [] 468 for i in range(total_pts): 469 # The constraint values. 470 ci = c(pts[i]) 471 472 # No constraint violations, so store the point. 473 if min(ci) >= 0.0: 474 pts_trimmed.append(pts[i]) 475 476 # No point elimination. 477 else: 478 pts_trimmed = pts 479 480 # Convert to numpy. 481 pts_trimmed = array(pts_trimmed) 482 483 # The total number of points has changed. 484 total_pts = len(pts_trimmed) 485 486 # The subdivision size (round up so the last subdivision is smaller than the rest). 487 size_float = total_pts / float(divisions) 488 size = int(size_float) 489 if size_float % 1: 490 size = size + 1 491 492 # Subdivide. 493 for i in range(min(divisions, total_pts)): 494 # The start index. 495 start = i * size 496 497 # The end index. 498 if i != divisions - 1: 499 end = (i + 1) * size 500 else: 501 end = total_pts 502 503 # Yield the subdivision. 504 yield pts[start: end]

505

Source Code for Module minfx.grid