pipe_control.minimise

1 ############################################################################### 2 # # 3 # Copyright (C) 2003-2014 Edward d'Auvergne # 4 # # 5 # This file is part of the program relax (http://www.nmr-relax.com). # 6 # # 7 # This program is free software: you can redistribute it and/or modify # 8 # it under the terms of the GNU General Public License as published by # 9 # the Free Software Foundation, either version 3 of the License, or # 10 # (at your option) any later version. # 11 # # 12 # This program is distributed in the hope that it will be useful, # 13 # but WITHOUT ANY WARRANTY; without even the implied warranty of # 14 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # 15 # GNU General Public License for more details. # 16 # # 17 # You should have received a copy of the GNU General Public License # 18 # along with this program. If not, see <http://www.gnu.org/licenses/>. # 19 # # 20 ############################################################################### 21 22 # Module docstring. 23 """Module for model minimisation/optimisation.""" 24 25 # Python module imports. 26 from numpy import float64, identity 27 import sys 28 29 # relax module imports. 30 from lib.errors import RelaxError, RelaxIntListIntError, RelaxLenError 31 from lib.float import isNaN 32 from lib.io import write_data 33 from multi import Processor_box 34 from pipe_control.mol_res_spin import return_spin, spin_loop 35 from pipe_control import pipes 36 from pipe_control.pipes import check_pipe 37 from specific_analyses.api import return_api, return_parameter_object 38 from status import Status; status = Status() 39 from user_functions.data import Uf_tables; uf_tables = Uf_tables() 40 41

42 -def assemble_scaling_matrix(scaling=True):

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

86 87

88 -def calc(verbosity=1):

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

143 144

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

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

217 218

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

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

415 416

417 -def grid_zoom(level=0):

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

429 430

431 -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):

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

533 534

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

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

608 609

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

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

669

Source Code for Module pipe_control.minimise