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 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, 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 verbosity: The verbosity level. 0 corresponds to no output, 1 is standard, and higher values cause greater and greater amount of output. 283 @type verbosity: int 284 @keyword print_prefix: The text to place before the printed output. 285 @type print_prefix: str 286 @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. 287 @rtype: tuple of numpy rank-1 array, float, int, str 288 """ 289 290 # Initialise. 291 total_steps = len(points) 292 n = len(points[0]) 293 294 # Set a ridiculously large initial function value. 295 f_min = 1e300 296 297 # Initial printout. 298 if verbosity: 299 print("\n" + print_prefix + "Searching through %s grid nodes." % total_steps) 300 301 # Test if the grid is too large. 302 if total_steps >= 1e8: 303 raise MinfxError("A grid search of size %s is too large." % total_steps) 304 305 # Search the grid. 306 for k in range(total_steps): 307 # Back calculate the current function value. 308 f = func(*(points[k],)+args) 309 310 # A better point. 311 if f < f_min: 312 # Switch to the new point. 313 f_min = f 314 min_params = 1.0 * points[k] 315 316 # Print out code. 317 if verbosity: 318 print_iter(k=k, xk=min_params, fk=f_min, print_prefix=print_prefix) 319 320 # Print out code. 321 if verbosity >= 2: 322 if f != f_min: 323 print_iter(k=k, xk=min_params, fk=f, print_prefix=print_prefix) 324 if verbosity >= 3: 325 print(print_prefix + "%-20s%-20s" % ("Params:", repr(points[k]))) 326 print(print_prefix + "%-20s%-20s" % ("Min params:", repr(min_params))) 327 print(print_prefix + "%-20s%-20g\n" % ("f:", f)) 328 print(print_prefix + "%-20s%-20g\n" % ("Min f:", f_min)) 329 330 # Return the results. 331 return min_params, f_min, total_steps, None

332 333

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

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

460

Source Code for Module minfx.grid