minfx.grid

1 ############################################################################### 2 # # 3 # Copyright (C) 2003-2015 Edward d'Auvergne # 4 # # 5 # This file is part of the minfx optimisation library, # 6 # https://gna.org/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://gna.org/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 is not None and b is not 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=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 is not None and b is not 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=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, verbosity=0, print_prefix=""):

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 @keyword verbosity: The verbosity level. 0 corresponds to no output, 1 is standard, and higher values cause greater and greater amount of output. 404 @type verbosity: int 405 @keyword print_prefix: The text to place before the printed output. 406 @type print_prefix: str 407 @return: A list of grid points for each subdivision is yielded. 408 @rtype: list of list of float 409 """ 410 411 # Printout. 412 if verbosity and divisions > 1: 413 print("%sSplitting up the grid points." % print_prefix) 414 415 # The dimensionality. 416 n = len(inc) 417 418 # Linear constraints. 419 if A is not None and b is not None: 420 constraint_flag = True 421 constraint_linear = Constraint_linear(A, b) 422 c = constraint_linear.func 423 if verbosity >= 3: 424 print(print_prefix + "Linear constraint matrices.") 425 print(print_prefix + "A: " + repr(A)) 426 print(print_prefix + "b: " + repr(b)) 427 428 # Bound constraints. 429 elif l != None and u != None: 430 constraint_flag = True 431 raise MinfxError("Bound constraints are not implemented yet.") 432 433 # General constraints. 434 elif c != None: 435 constraint_flag = True 436 437 # No constraints. 438 else: 439 constraint_flag = False 440 441 # Total number of points. 442 total_pts = 1 443 for i in range(n): 444 total_pts = total_pts * inc[i] 445 446 # Init. 447 indices = zeros(n, int) 448 pts = zeros((total_pts, n), float64) 449 450 # The increment values for each dimension. 451 incs = [] 452 for k in range(n): 453 incs.append([]) 454 if inc[k] == 1: 455 incs[k].append((lower[k] + upper[k])/2.0) 456 else: 457 for i in range(inc[k]): 458 incs[k].append(lower[k] + i * (upper[k] - lower[k]) / (inc[k] - 1)) 459 460 # Construct the list of all grid points. 461 for i in range(total_pts): 462 # Loop over the dimensions. 463 for j in range(n): 464 # Add the point coordinate. 465 pts[i][j] = incs[j][indices[j]] 466 467 # Increment the step positions. 468 for j in range(n): 469 if indices[j] < inc[j]-1: 470 indices[j] += 1 471 break # Exit so that the other step numbers are not incremented. 472 else: 473 indices[j] = 0 474 475 # Eliminate points outside of constraints. 476 if constraint_flag: 477 # Loop over all points. 478 pts_trimmed = [] 479 for i in range(total_pts): 480 # The constraint values. 481 ci = c(pts[i]) 482 483 # No constraint violations, so store the point. 484 if min(ci) >= 0.0: 485 pts_trimmed.append(pts[i]) 486 487 # No point elimination. 488 else: 489 pts_trimmed = pts 490 491 # Convert to numpy. 492 pts_trimmed = array(pts_trimmed) 493 494 # The total number of points has changed. 495 total_pts = len(pts_trimmed) 496 497 # The subdivision size (round up so the last subdivision is smaller than the rest). 498 size_float = total_pts / float(divisions) 499 size = int(size_float) 500 if size_float % 1: 501 size = size + 1 502 503 # Printout. 504 if verbosity and divisions > 1: 505 print("%s Total number of grid points: %20i" % (print_prefix, total_pts)) 506 print("%s Number of divisions: %20i" % (print_prefix, divisions)) 507 print("%s Subdivision size: %20i" % (print_prefix, size)) 508 509 # Subdivide. 510 for i in range(min(divisions, total_pts)): 511 # The start index. 512 start = i * size 513 514 # The end index. 515 if i != divisions - 1: 516 end = (i + 1) * size 517 else: 518 end = total_pts 519 520 # Yield the subdivision. 521 yield pts_trimmed[start: end]

522 523

524 -def grid_split_array(divisions=None, points=None, A=None, b=None, l=None, u=None, c=None, verbosity=0, print_prefix=""):

525 """Generator method yielding arrays of grid points. 526 527 This method will loop over the grid points one-by-one, generating a list of points and yielding these for each subdivision. 528 529 530 @keyword divisions: The number of grid subdivisions. 531 @type divisions: int 532 @keyword points: The array of grid points to split up. 533 @type points: list of lists or rank-2 numpy array 534 @keyword A: The linear constraint matrix A, such that A.x >= b. 535 @type A: numpy rank-2 array 536 @keyword b: The linear constraint scalar vectors, such that A.x >= b. 537 @type b: numpy rank-1 array 538 @keyword l: The lower bound constraint vector, such that l <= x <= u. 539 @type l: list of float 540 @keyword u: The upper bound constraint vector, such that l <= x <= u. 541 @type u: list of float 542 @keyword c: A user supplied constraint function. 543 @type c: function 544 @keyword verbosity: The verbosity level. 0 corresponds to no output, 1 is standard, and higher values cause greater and greater amount of output. 545 @type verbosity: int 546 @keyword print_prefix: The text to place before the printed output. 547 @type print_prefix: str 548 @return: A list of grid points for each subdivision is yielded. 549 @rtype: list of list of float 550 """ 551 552 # Printout. 553 if verbosity and divisions > 1: 554 print("%sSplitting up the array of grid points." % print_prefix) 555 556 # Linear constraints. 557 if A is not None and b is not None: 558 constraint_flag = True 559 constraint_linear = Constraint_linear(A, b) 560 c = constraint_linear.func 561 if verbosity >= 3: 562 print(print_prefix + "Linear constraint matrices.") 563 print(print_prefix + "A: " + repr(A)) 564 print(print_prefix + "b: " + repr(b)) 565 566 # Bound constraints. 567 elif l != None and u != None: 568 constraint_flag = True 569 raise MinfxError("Bound constraints are not implemented yet.") 570 571 # General constraints. 572 elif c != None: 573 constraint_flag = True 574 575 # No constraints. 576 else: 577 constraint_flag = False 578 579 # Eliminate points outside of constraints. 580 if constraint_flag: 581 # Loop over all points. 582 pts_trimmed = [] 583 for i in range(len(points)): 584 # The constraint values. 585 ci = c(points[i]) 586 587 # No constraint violations, so store the point. 588 if min(ci) >= 0.0: 589 pts_trimmed.append(points[i]) 590 591 # No point elimination. 592 else: 593 pts_trimmed = points 594 595 # Convert to numpy. 596 pts_trimmed = array(pts_trimmed) 597 598 # Total number of points remaining. 599 total_pts = len(pts_trimmed) 600 601 # The subdivision size (round up so the last subdivision is smaller than the rest). 602 size_float = total_pts / float(divisions) 603 size = int(size_float) 604 if size_float % 1: 605 size = size + 1 606 607 # Printout. 608 if verbosity and divisions > 1: 609 print("%s Total number of grid points: %20i" % (print_prefix, total_pts)) 610 print("%s Number of divisions: %20i" % (print_prefix, divisions)) 611 print("%s Subdivision size: %20i" % (print_prefix, size)) 612 613 # Subdivide. 614 for i in range(min(divisions, total_pts)): 615 # The start index. 616 start = i * size 617 618 # The end index. 619 if i != divisions - 1: 620 end = (i + 1) * size 621 else: 622 end = total_pts 623 624 # Yield the subdivision. 625 yield pts_trimmed[start: end]

626

Source Code for Module minfx.grid