pipe_control.minimise

1 ############################################################################### 2 # # 3 # Copyright (C) 2002-2005,2007-2014 Edward d'Auvergne # 4 # Copyright (C) 2006 Chris MacRaild # 5 # # 6 # This file is part of the program relax (http://www.nmr-relax.com). # 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 """Module for model minimisation/optimisation.""" 25 26 # Python module imports. 27 from numpy import float64, identity 28 import sys 29 30 # relax module imports. 31 from lib.errors import RelaxError, RelaxIntListIntError, RelaxLenError 32 from lib.float import isNaN 33 from lib.io import write_data 34 from multi import Processor_box 35 from pipe_control.mol_res_spin import return_spin, spin_loop 36 from pipe_control import pipes 37 from pipe_control.pipes import check_pipe 38 from specific_analyses.api import return_api, return_parameter_object 39 from status import Status; status = Status() 40 from user_functions.data import Uf_tables; uf_tables = Uf_tables() 41 42

43 -def assemble_scaling_matrix(scaling=True):

44 """Create and return the per-model scaling matrices. 45 46 @keyword scaling: If True, diagonal scaling is enabled during optimisation to allow the problem to be better conditioned. 47 @type scaling: bool 48 @return: The list of diagonal and square scaling matrices. 49 @rtype: list of numpy rank-2, float64 array or list of None 50 """ 51 52 # The specific analysis API object and parameter object. 53 api = return_api() 54 param_object = return_parameter_object() 55 56 # Initialise. 57 scaling_matrix = [] 58 59 # Loop over the models. 60 for model_info in api.model_loop(): 61 # No diagonal scaling. 62 if not scaling: 63 scaling_matrix.append(None) 64 continue 65 66 # Get the parameter names. 67 names = api.get_param_names(model_info) 68 69 # No parameters for this model. 70 if names == None or len(names) == 0: 71 scaling_matrix.append(None) 72 continue 73 74 # The parameter number. 75 n = len(names) 76 77 # Initialise. 78 scaling_matrix.append(identity(n, float64)) 79 i = 0 80 81 # Update the diagonal with the parameter specific scaling factor. 82 for i in range(n): 83 scaling_matrix[-1][i, i] = param_object.scaling(names[i], model_info=model_info) 84 85 # Return the matrix. 86 return scaling_matrix

87 88

89 -def calc(verbosity=1):

90 """Function for calculating the function value. 91 92 @keyword verbosity: The amount of information to print. The higher the value, the greater the verbosity. 93 @type verbosity: int 94 """ 95 96 # Test if the current data pipe exists. 97 check_pipe() 98 99 # Reset the minimisation statistics. 100 reset_min_stats(verbosity=verbosity) 101 102 # The specific analysis API object. 103 api = return_api() 104 105 # Deselect spins lacking data: 106 api.overfit_deselect() 107 108 # Create the scaling matrix. 109 scaling_matrix = assemble_scaling_matrix() 110 111 # Get the Processor box singleton (it contains the Processor instance) and alias the Processor. 112 processor_box = Processor_box() 113 processor = processor_box.processor 114 115 # Monte Carlo simulation calculation. 116 if hasattr(cdp, 'sim_state') and cdp.sim_state == 1: 117 # Loop over the simulations. 118 for i in range(cdp.sim_number): 119 # Status. 120 if status.current_analysis: 121 status.auto_analysis[status.current_analysis].mc_number = i 122 else: 123 status.mc_number = i 124 125 # Calculation. 126 api.calculate(verbosity=verbosity-1, sim_index=i, scaling_matrix=scaling_matrix) 127 128 # Print out. 129 if verbosity and not processor.is_queued(): 130 print("Simulation " + repr(i+1)) 131 132 # Unset the status. 133 if status.current_analysis: 134 status.auto_analysis[status.current_analysis].mc_number = None 135 else: 136 status.mc_number = None 137 138 # Minimisation. 139 else: 140 api.calculate(verbosity=verbosity, scaling_matrix=scaling_matrix) 141 142 # Execute any queued commands. 143 processor.run_queue()

144 145

146 -def grid_search(lower=None, upper=None, inc=None, verbosity=1, constraints=True, skip_preset=True):

147 """The grid search function. 148 149 @keyword lower: The lower bounds of the grid search which must be equal to the number of parameters in the model. 150 @type lower: array of numbers 151 @keyword upper: The upper bounds of the grid search which must be equal to the number of parameters in the model. 152 @type upper: array of numbers 153 @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. 154 @type inc: int or list of int 155 @keyword verbosity: The amount of information to print. The higher the value, the greater the verbosity. 156 @type verbosity: int 157 @keyword constraints: If True, constraints are applied during the grid search (elinating parts of the grid). If False, no constraints are used. 158 @type constraints: bool 159 @keyword skip_preset: This argument, when True, allows any parameter which already has a value set to be skipped in the grid search. 160 @type skip_preset: bool 161 """ 162 163 # Test if the current data pipe exists. 164 check_pipe() 165 166 # The specific analysis API object. 167 api = return_api() 168 169 # Deselect models lacking data: 170 api.overfit_deselect() 171 172 # Determine the model specific grid bounds, and allow for the zooming grid search, and check the inc argument. 173 model_lower, model_upper, model_inc = grid_setup(lower, upper, inc, verbosity=verbosity, skip_preset=skip_preset) 174 175 # Create the scaling matrix. 176 scaling_matrix = assemble_scaling_matrix() 177 178 # Get the Processor box singleton (it contains the Processor instance) and alias the Processor. 179 processor_box = Processor_box() 180 processor = processor_box.processor 181 182 # Monte Carlo simulation grid search. 183 if hasattr(cdp, 'sim_state') and cdp.sim_state == 1: 184 # Loop over the simulations. 185 for i in range(cdp.sim_number): 186 # Reset the minimisation statistics. 187 reset_min_stats(sim_index=i, verbosity=verbosity) 188 189 # Status. 190 if status.current_analysis: 191 status.auto_analysis[status.current_analysis].mc_number = i 192 else: 193 status.mc_number = i 194 195 # Optimisation. 196 api.grid_search(lower=model_lower, upper=model_upper, inc=model_inc, scaling_matrix=scaling_matrix, constraints=constraints, verbosity=verbosity-1, sim_index=i) 197 198 # Print out. 199 if verbosity and not processor.is_queued(): 200 print("Simulation " + repr(i+1)) 201 202 # Unset the status. 203 if status.current_analysis: 204 status.auto_analysis[status.current_analysis].mc_number = None 205 else: 206 status.mc_number = None 207 208 # Grid search. 209 else: 210 # Reset the minimisation statistics. 211 reset_min_stats(verbosity=verbosity) 212 213 # Optimise. 214 api.grid_search(lower=model_lower, upper=model_upper, inc=model_inc, scaling_matrix=scaling_matrix, constraints=constraints, verbosity=verbosity) 215 216 # Execute any queued commands. 217 processor.run_queue()

218 219

220 -def grid_setup(lower=None, upper=None, inc=None, verbosity=1, skip_preset=True):

221 """Determine the per-model grid bounds, allowing for the zooming grid search. 222 223 @keyword lower: The user supplied lower bounds of the grid search which must be equal to the number of parameters in the model. 224 @type lower: list of numbers 225 @keyword upper: The user supplied upper bounds of the grid search which must be equal to the number of parameters in the model. 226 @type upper: list of numbers 227 @keyword inc: The user supplied grid search increments. 228 @type inc: int or list of int 229 @keyword verbosity: The amount of information to print. The higher the value, the greater the verbosity. 230 @type verbosity: int 231 @keyword skip_preset: This argument, when True, allows any parameter which already has a value set to be skipped in the grid search. 232 @type skip_preset: bool 233 @return: The per-model grid upper and lower bounds. The first dimension of each structure corresponds to the model, the second the model parameters. 234 @rtype: tuple of lists of lists of float, lists of lists of float, list of lists of int 235 """ 236 237 # The specific analysis API object and parameter object. 238 api = return_api() 239 param_object = return_parameter_object() 240 241 # Initialise. 242 model_lower = [] 243 model_upper = [] 244 model_inc = [] 245 246 # Loop over the models. 247 for model_info in api.model_loop(): 248 # Get the parameter names and current values. 249 names = api.get_param_names(model_info) 250 values = api.get_param_values(model_info) 251 252 # No parameters for this model. 253 if names == None or len(names) == 0: 254 model_lower.append([]) 255 model_upper.append([]) 256 model_inc.append([]) 257 continue 258 259 # The parameter number. 260 n = len(names) 261 262 # Make sure that the length of the parameter array is > 0. 263 if n == 0: 264 raise RelaxError("Cannot run a grid search on a model with zero parameters.") 265 266 # Check that the user supplied bound lengths are ok. 267 if lower != None and len(lower) != n: 268 raise RelaxLenError('lower bounds', n) 269 if upper != None and len(upper) != n: 270 raise RelaxLenError('upper bounds', n) 271 272 # Check the user supplied increments. 273 if isinstance(inc, list) and len(inc) != n: 274 raise RelaxLenError('increment', n) 275 if isinstance(inc, list): 276 for i in range(n): 277 if not (isinstance(inc[i], int) or inc[i] == None): 278 raise RelaxIntListIntError('increment', inc) 279 elif not isinstance(inc, int): 280 raise RelaxIntListIntError('increment', inc) 281 282 # Convert to the model increment list. 283 if isinstance(inc, int): 284 model_inc.append([inc]*n) 285 else: 286 model_inc.append(inc) 287 288 # Print out the model title. 289 api.print_model_title(prefix="Grid search setup: ", model_info=model_info) 290 291 # The grid zoom level. 292 zoom = 0 293 if hasattr(cdp, 'grid_zoom_level'): 294 zoom = cdp.grid_zoom_level 295 zoom_factor = 1.0 / 2.0**zoom 296 if zoom > 0: 297 print("Zooming grid level of %s, scaling the grid size by a factor of %s.\n" % (zoom, zoom_factor)) 298 299 # Append empty lists for the bounds to be built up. 300 model_lower.append([]) 301 model_upper.append([]) 302 303 # Loop over the parameters. 304 data = [] 305 for i in range(n): 306 # A comment for user feedback. 307 comment = 'Default bounds' 308 if lower != None and upper != None: 309 comment = 'User supplied lower and upper bound' 310 elif lower != None: 311 comment = 'User supplied lower bound' 312 elif upper != None: 313 comment = 'User supplied upper bound' 314 315 # Alias the number of increments for this parameter. 316 incs = model_inc[-1][i] 317 318 # Error checking for increment values of None. 319 if incs == None and values[i] in [None, {}, []]: 320 raise RelaxError("The parameter '%s' has no preset value, therefore a grid increment of None is not valid." % names[i]) 321 322 # The lower bound for this parameter. 323 if lower != None: 324 lower_i = lower[i] 325 else: 326 lower_i = param_object.grid_lower(names[i], incs=incs, model_info=model_info) 327 328 # The upper bound for this parameter. 329 if upper != None: 330 upper_i = upper[i] 331 else: 332 upper_i = param_object.grid_upper(names[i], incs=incs, model_info=model_info) 333 334 # The skipping logic. 335 skip = False 336 if skip_preset: 337 # Override the flag if the zoom is on. 338 if zoom: 339 skip = False 340 341 # No preset value. 342 elif values[i] in [None, {}, []]: 343 skip = False 344 345 # The preset value is a NaN value due to numpy conversions of None. 346 elif isNaN(values[i]): 347 skip = False 348 349 # Ok, now the parameter can be skipped. 350 else: 351 skip = True 352 353 # Override the skip flag if the incs value is None. 354 if incs == None: 355 skip = True 356 357 # Skip preset values. 358 if skip: 359 lower_i = values[i] 360 upper_i = values[i] 361 model_inc[-1][i] = incs = 1 362 comment = 'Preset value' 363 364 # Zooming grid. 365 elif zoom: 366 # The full size and scaled size. 367 size = upper_i - lower_i 368 zoom_size = size * zoom_factor 369 half_size = zoom_size / 2.0 370 comment = 'Zoom grid width of %s %s' % (zoom_size, param_object.units(names[i])) 371 372 # The new size around the current value. 373 lower_zoom = values[i] - half_size 374 upper_zoom = values[i] + half_size 375 376 # Outside of the original lower bound, so shift the grid to fit. 377 if zoom > 0 and lower_zoom < lower_i: 378 # The amount to shift by. 379 shift = lower_i - lower_zoom 380 381 # Set the new bounds. 382 upper_i = upper_zoom + shift 383 384 # Outside of the original upper bound, so shift the grid to fit. 385 elif zoom > 0 and upper_zoom > upper_i: 386 # The amount to shift by. 387 shift = upper_i - upper_zoom 388 389 # Set the new bounds. 390 lower_i = lower_zoom + shift 391 392 # Inside the original bounds. 393 else: 394 lower_i = lower_zoom 395 upper_i = upper_zoom 396 397 # Add to the data list for printing out. 398 data.append([names[i], "%15s" % lower_i, "%15s" % upper_i, "%15s" % incs, comment]) 399 400 # Scale the bounds. 401 scaling = param_object.scaling(names[i], model_info=model_info) 402 lower_i /= scaling 403 upper_i /= scaling 404 405 # Append. 406 model_lower[-1].append(lower_i) 407 model_upper[-1].append(upper_i) 408 409 # Printout. 410 if verbosity: 411 write_data(out=sys.stdout, headings=["Parameter", "Lower bound", "Upper bound", "Increments", "Comment"], data=data) 412 sys.stdout.write('\n') 413 414 # Return the bounds. 415 return model_lower, model_upper, model_inc

416 417

418 -def grid_zoom(level=0):

419 """Store the grid zoom level. 420 421 @keyword level: The zoom level. 422 @type level: int 423 """ 424 425 # Test if the current data pipe exists. 426 check_pipe() 427 428 # Store the values. 429 cdp.grid_zoom_level = level

430 431

432 -def minimise(min_algor=None, line_search=None, hessian_mod=None, hessian_type=None, func_tol=None, grad_tol=None, max_iter=None, constraints=True, scaling=True, verbosity=1, sim_index=None):

433 """Minimisation function. 434 435 @keyword min_algor: The minimisation algorithm to use. 436 @type min_algor: str 437 @keyword line_search: The line search algorithm which will only be used in combination with the line search and conjugate gradient methods. This will default to the More and Thuente line search. 438 @type line_search: str or None 439 @keyword hessian_mod: The Hessian modification. This will only be used in the algorithms which use the Hessian, and defaults to Gill, Murray, and Wright modified Cholesky algorithm. 440 @type hessian_mod: str or None 441 @keyword hessian_type: The Hessian type. This will only be used in a few trust region algorithms, and defaults to BFGS. 442 @type hessian_type: str or None 443 @keyword func_tol: The function tolerance which, when reached, terminates optimisation. Setting this to None turns of the check. 444 @type func_tol: None or float 445 @keyword grad_tol: The gradient tolerance which, when reached, terminates optimisation. Setting this to None turns of the check. 446 @type grad_tol: None or float 447 @keyword max_iter: The maximum number of iterations for the algorithm. 448 @type max_iter: int 449 @keyword constraints: If True, constraints are used during optimisation. 450 @type constraints: bool 451 @keyword scaling: If True, diagonal scaling is enabled during optimisation to allow the problem to be better conditioned. 452 @type scaling: bool 453 @keyword verbosity: The amount of information to print. The higher the value, the greater the verbosity. 454 @type verbosity: int 455 @keyword sim_index: The index of the simulation to optimise. This should be None if normal optimisation is desired. 456 @type sim_index: None or int 457 """ 458 459 # Test if the current data pipe exists. 460 check_pipe() 461 462 # The specific analysis API object. 463 api = return_api() 464 465 # Re-package the minimisation algorithm, options, and constraints for the generic_minimise() calls within the specific code. 466 if constraints: 467 min_options = [min_algor] 468 469 # Determine the constraint algorithm to use. 470 min_algor = api.constraint_algorithm() 471 else: 472 min_options = [] 473 if line_search != None: 474 min_options.append(line_search) 475 if hessian_mod != None: 476 min_options.append(hessian_mod) 477 if hessian_type != None: 478 min_options.append(hessian_type) 479 min_options = tuple(min_options) 480 481 # Deselect spins lacking data: 482 api.overfit_deselect() 483 484 # Create the scaling matrix. 485 scaling_matrix = assemble_scaling_matrix(scaling) 486 487 # Get the Processor box singleton (it contains the Processor instance) and alias the Processor. 488 processor_box = Processor_box() 489 processor = processor_box.processor 490 491 # Single Monte Carlo simulation. 492 if sim_index != None: 493 # Reset the minimisation statistics. 494 reset_min_stats(sim_index=sim_index, verbosity=verbosity) 495 496 # Optimise. 497 api.minimise(min_algor=min_algor, min_options=min_options, func_tol=func_tol, grad_tol=grad_tol, max_iterations=max_iter, constraints=constraints, scaling_matrix=scaling_matrix, verbosity=verbosity, sim_index=sim_index) 498 499 # Monte Carlo simulation minimisation. 500 elif hasattr(cdp, 'sim_state') and cdp.sim_state == 1: 501 for i in range(cdp.sim_number): 502 # Reset the minimisation statistics. 503 reset_min_stats(sim_index=i, verbosity=verbosity) 504 505 # Status. 506 if status.current_analysis: 507 status.auto_analysis[status.current_analysis].mc_number = i 508 else: 509 status.mc_number = i 510 511 # Optimisation. 512 api.minimise(min_algor=min_algor, min_options=min_options, func_tol=func_tol, grad_tol=grad_tol, max_iterations=max_iter, constraints=constraints, scaling_matrix=scaling_matrix, verbosity=verbosity-1, sim_index=i) 513 514 # Print out. 515 if verbosity and not processor.is_queued(): 516 print("Simulation " + repr(i+1)) 517 518 # Unset the status. 519 if status.current_analysis: 520 status.auto_analysis[status.current_analysis].mc_number = None 521 else: 522 status.mc_number = None 523 524 # Standard minimisation. 525 else: 526 # Reset the minimisation statistics. 527 reset_min_stats(verbosity=verbosity) 528 529 # Optimise. 530 api.minimise(min_algor=min_algor, min_options=min_options, func_tol=func_tol, grad_tol=grad_tol, max_iterations=max_iter, constraints=constraints, scaling_matrix=scaling_matrix, verbosity=verbosity) 531 532 # Execute any queued commands. 533 processor.run_queue()

534 535

536 -def reset_min_stats(data_pipe=None, sim_index=None, verbosity=1):

537 """Function for resetting all minimisation statistics. 538 539 @keyword data_pipe: The name of the data pipe to reset the minimisation statistics of. This defaults to the current data pipe. 540 @type data_pipe: str 541 @keyword sim_index: The optional Monte Carlo simulation index. 542 @type sim_index: int 543 @keyword verbosity: The amount of information to print. The higher the value, the greater the verbosity. 544 @type verbosity: int 545 """ 546 547 # The data pipe. 548 if data_pipe == None: 549 data_pipe = pipes.cdp_name() 550 551 # Get the data pipe. 552 dp = pipes.get_pipe(data_pipe) 553 554 # The objects to reset to None. 555 names = [ 556 'chi2', 557 'iter', 558 'f_count', 559 'g_count', 560 'h_count', 561 'warning' 562 ] 563 564 # Loop over the objects. 565 flag = False 566 for name in names: 567 # The simulation name. 568 sim_name = name + '_sim' 569 570 # Global minimisation statistics. 571 if sim_index == None: 572 # Reset the object to None if it exists. 573 if hasattr(dp, name): 574 setattr(dp, name, None) 575 flag = False 576 577 # Global MC minimisation statistics. 578 else: 579 if hasattr(dp, sim_name): 580 # The MC simulation object. 581 sim_obj = getattr(dp, sim_name) 582 583 # Reset the object to None if possible. 584 if sim_index < len(sim_obj): 585 sim_obj[sim_index] = None 586 587 # Loop over all spins. 588 for spin in spin_loop(skip_desel=False): 589 # The minimisation statistics. 590 if sim_index == None: 591 # Reset the object to None if it exists. 592 if hasattr(spin, name): 593 setattr(spin, name, None) 594 flag = True 595 596 # The MC minimisation statistics. 597 else: 598 if hasattr(spin, sim_name): 599 # The MC simulation object. 600 sim_obj = getattr(spin, sim_name) 601 602 # Reset the object to None if possible. 603 if sim_index < len(sim_obj): 604 sim_obj[sim_index] = None 605 606 # Printout. 607 if verbosity and flag and sim_index == None: 608 print("Resetting the minimisation statistics.")

609 610

611 -def set(val=None, error=None, param=None, scaling=None, spin_id=None):

612 """Set global or spin specific minimisation parameters. 613 614 @keyword val: The parameter values. 615 @type val: number 616 @keyword param: The parameter names. 617 @type param: str 618 @keyword scaling: Unused. 619 @type scaling: float 620 @keyword spin_id: The spin identification string. 621 @type spin_id: str 622 """ 623 624 # Global minimisation stats. 625 if spin_id == None: 626 # Chi-squared. 627 if param == 'chi2': 628 cdp.chi2 = val 629 630 # Iteration count. 631 elif param == 'iter': 632 cdp.iter = val 633 634 # Function call count. 635 elif param == 'f_count': 636 cdp.f_count = val 637 638 # Gradient call count. 639 elif param == 'g_count': 640 cdp.g_count = val 641 642 # Hessian call count. 643 elif param == 'h_count': 644 cdp.h_count = val 645 646 # Residue specific minimisation. 647 else: 648 # Get the spin. 649 spin = return_spin(spin_id=spin_id) 650 651 # Chi-squared. 652 if param == 'chi2': 653 spin.chi2 = val 654 655 # Iteration count. 656 elif param == 'iter': 657 spin.iter = val 658 659 # Function call count. 660 elif param == 'f_count': 661 spin.f_count = val 662 663 # Gradient call count. 664 elif param == 'g_count': 665 spin.g_count = val 666 667 # Hessian call count. 668 elif param == 'h_count': 669 spin.h_count = val

670

Source Code for Module pipe_control.minimise