specific_analyses.n_state_model.parameters

1 ############################################################################### 2 # # 3 # Copyright (C) 2007-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 """The N-state model or structural ensemble analysis parameter handling.""" 24 25 # Python module imports. 26 from numpy import array, float64, zeros 27 from re import search 28 from warnings import warn 29 30 # relax module imports. 31 from lib.errors import RelaxNoModelError 32 from lib.warnings import RelaxWarning 33 from pipe_control import align_tensor 34 from pipe_control.pipes import check_pipe 35 from pipe_control.align_tensor import opt_uses_align_data, opt_uses_tensor 36 from specific_analyses.n_state_model.data import base_data_types 37 38

39 -def assemble_param_vector(sim_index=None):

40 """Assemble all the parameters of the model into a single array. 41 42 @param sim_index: The index of the simulation to optimise. This should be None if normal optimisation is desired. 43 @type sim_index: None or int 44 @return: The parameter vector used for optimisation. 45 @rtype: numpy array 46 """ 47 48 # Test if the model is selected. 49 if not hasattr(cdp, 'model') or not isinstance(cdp.model, str): 50 raise RelaxNoModelError 51 52 # Determine the data type. 53 data_types = base_data_types() 54 55 # Initialise the parameter vector. 56 param_vector = [] 57 58 # A RDC or PCS data type requires the alignment tensors to be at the start of the parameter vector (unless the tensors are fixed). 59 if opt_uses_align_data(): 60 for i in range(len(cdp.align_tensors)): 61 # Skip non-optimised tensors. 62 if not opt_uses_tensor(cdp.align_tensors[i]): 63 continue 64 65 # No values set. 66 if not hasattr(cdp.align_tensors[i], 'A_5D'): 67 param_vector += [None, None, None, None, None] 68 69 # Otherwise add the parameters. 70 else: 71 param_vector += list(cdp.align_tensors[i].A_5D) 72 73 # Monte Carlo simulation data structures. 74 if sim_index != None: 75 # Populations. 76 if cdp.model in ['2-domain', 'population']: 77 probs = cdp.probs_sim[sim_index] 78 79 # Euler angles. 80 if cdp.model == '2-domain': 81 alpha = cdp.alpha_sim[sim_index] 82 beta = cdp.beta_sim[sim_index] 83 gamma = cdp.gamma_sim[sim_index] 84 85 # Normal data structures. 86 else: 87 # Populations. 88 if cdp.model in ['2-domain', 'population']: 89 probs = cdp.probs 90 91 # Euler angles. 92 if cdp.model == '2-domain': 93 alpha = cdp.alpha 94 beta = cdp.beta 95 gamma = cdp.gamma 96 97 # The probabilities (exclude that of state N). 98 if cdp.model in ['2-domain', 'population']: 99 param_vector = param_vector + probs[0:-1] 100 101 # The Euler angles. 102 if cdp.model == '2-domain': 103 for i in range(cdp.N): 104 param_vector.append(alpha[i]) 105 param_vector.append(beta[i]) 106 param_vector.append(gamma[i]) 107 108 # The paramagnetic centre. 109 if hasattr(cdp, 'paramag_centre_fixed') and not cdp.paramag_centre_fixed: 110 if not hasattr(cdp, 'paramagnetic_centre'): 111 for i in range(3): 112 param_vector.append(None) 113 elif sim_index != None: 114 if cdp.paramagnetic_centre_sim[sim_index] is None: 115 for i in range(3): 116 param_vector.append(None) 117 else: 118 for i in range(3): 119 param_vector.append(cdp.paramagnetic_centre_sim[sim_index][i]) 120 else: 121 for i in range(3): 122 param_vector.append(cdp.paramagnetic_centre[i]) 123 124 # Return a numpy arrary. 125 return array(param_vector, float64)

126 127

128 -def disassemble_param_vector(param_vector=None, data_types=None, sim_index=None):

129 """Disassemble the parameter vector and place the values into the relevant variables. 130 131 For the 2-domain N-state model, the parameters are stored in the probability and Euler angle data structures. For the population N-state model, only the probabilities are stored. If RDCs are present and alignment tensors are optimised, then these are stored as well. 132 133 @keyword data_types: The base data types used in the optimisation. This list can contain the elements 'rdc', 'pcs' or 'tensor'. 134 @type data_types: list of str 135 @keyword param_vector: The parameter vector returned from optimisation. 136 @type param_vector: numpy array 137 @keyword sim_index: The index of the simulation to optimise. This should be None if normal optimisation is desired. 138 @type sim_index: None or int 139 """ 140 141 # Test if the model is selected. 142 if not hasattr(cdp, 'model') or not isinstance(cdp.model, str): 143 raise RelaxNoModelError 144 145 # Unpack and strip off the alignment tensor parameters. 146 if ('rdc' in data_types or 'pcs' in data_types) and not align_tensor.all_tensors_fixed(): 147 # Loop over the alignments, adding the alignment tensor parameters to the tensor data container. 148 tensor_num = 0 149 for i in range(len(cdp.align_tensors)): 150 # Skip non-optimised tensors. 151 if not opt_uses_tensor(cdp.align_tensors[i]): 152 continue 153 154 # Normal tensors. 155 if sim_index == None: 156 cdp.align_tensors[i].set(param='Axx', value=param_vector[5*tensor_num]) 157 cdp.align_tensors[i].set(param='Ayy', value=param_vector[5*tensor_num+1]) 158 cdp.align_tensors[i].set(param='Axy', value=param_vector[5*tensor_num+2]) 159 cdp.align_tensors[i].set(param='Axz', value=param_vector[5*tensor_num+3]) 160 cdp.align_tensors[i].set(param='Ayz', value=param_vector[5*tensor_num+4]) 161 162 # Monte Carlo simulated tensors. 163 else: 164 cdp.align_tensors[i].set(param='Axx', value=param_vector[5*tensor_num], category='sim', sim_index=sim_index) 165 cdp.align_tensors[i].set(param='Ayy', value=param_vector[5*tensor_num+1], category='sim', sim_index=sim_index) 166 cdp.align_tensors[i].set(param='Axy', value=param_vector[5*tensor_num+2], category='sim', sim_index=sim_index) 167 cdp.align_tensors[i].set(param='Axz', value=param_vector[5*tensor_num+3], category='sim', sim_index=sim_index) 168 cdp.align_tensors[i].set(param='Ayz', value=param_vector[5*tensor_num+4], category='sim', sim_index=sim_index) 169 170 # Increase the tensor number. 171 tensor_num += 1 172 173 # Create a new parameter vector without the tensors. 174 param_vector = param_vector[5*tensor_num:] 175 176 # Alias the Monte Carlo simulation data structures. 177 if sim_index != None: 178 # Populations. 179 if cdp.model in ['2-domain', 'population']: 180 probs = cdp.probs_sim[sim_index] 181 182 # Euler angles. 183 if cdp.model == '2-domain': 184 alpha = cdp.alpha_sim[sim_index] 185 beta = cdp.beta_sim[sim_index] 186 gamma = cdp.gamma_sim[sim_index] 187 188 # Alias the normal data structures. 189 else: 190 # Populations. 191 if cdp.model in ['2-domain', 'population']: 192 probs = cdp.probs 193 194 # Euler angles. 195 if cdp.model == '2-domain': 196 alpha = cdp.alpha 197 beta = cdp.beta 198 gamma = cdp.gamma 199 200 # Set the probabilities for states 0 to N-1 in the aliased structures. 201 if cdp.model in ['2-domain', 'population']: 202 for i in range(cdp.N-1): 203 probs[i] = param_vector[i] 204 205 # The probability for state N. 206 probs[-1] = 1 - sum(probs[0:-1]) 207 208 # Set the Euler angles in the aliased structures. 209 if cdp.model == '2-domain': 210 for i in range(cdp.N): 211 alpha[i] = param_vector[cdp.N-1 + 3*i] 212 beta[i] = param_vector[cdp.N-1 + 3*i + 1] 213 gamma[i] = param_vector[cdp.N-1 + 3*i + 2] 214 215 # The paramagnetic centre. 216 if hasattr(cdp, 'paramag_centre_fixed') and not cdp.paramag_centre_fixed: 217 # Create the structure if needed. 218 if not hasattr(cdp, 'paramagnetic_centre'): 219 cdp.paramagnetic_centre = zeros(3, float64) 220 221 # The position. 222 if sim_index == None: 223 cdp.paramagnetic_centre[0] = param_vector[-3] 224 cdp.paramagnetic_centre[1] = param_vector[-2] 225 cdp.paramagnetic_centre[2] = param_vector[-1] 226 227 # Monte Carlo simulated positions. 228 else: 229 if cdp.paramagnetic_centre_sim[sim_index] is None: 230 cdp.paramagnetic_centre_sim[sim_index] = [None, None, None] 231 cdp.paramagnetic_centre_sim[sim_index][0] = param_vector[-3] 232 cdp.paramagnetic_centre_sim[sim_index][1] = param_vector[-2] 233 cdp.paramagnetic_centre_sim[sim_index][2] = param_vector[-1]

234 235

236 -def elim_no_prob():

237 """Remove all structures or states which have no probability.""" 238 239 # Test if the current data pipe exists. 240 check_pipe() 241 242 # Test if the model is setup. 243 if not hasattr(cdp, 'model'): 244 raise RelaxNoModelError('N-state') 245 246 # Test if there are populations. 247 if not hasattr(cdp, 'probs'): 248 raise RelaxError("The N-state model populations do not exist.") 249 250 # Create the data structure if needed. 251 if not hasattr(cdp, 'select_models'): 252 cdp.state_select = [True] * cdp.N 253 254 # Loop over the structures. 255 for i in range(len(cdp.N)): 256 # No probability. 257 if cdp.probs[i] < 1e-5: 258 cdp.state_select[i] = False

259 260

261 -def linear_constraints(data_types=None, scaling_matrix=None):

262 """Function for setting up the linear constraint matrices A and b. 263 264 Standard notation 265 ================= 266 267 The N-state model constraints are:: 268 269 0 <= pc <= 1, 270 271 where p is the probability and c corresponds to state c. 272 273 274 Matrix notation 275 =============== 276 277 In the notation A.x >= b, where A is an matrix of coefficients, x is an array of parameter 278 values, and b is a vector of scalars, these inequality constraints are:: 279 280 | 1 0 0 | | 0 | 281 | | | | 282 |-1 0 0 | | -1 | 283 | | | | 284 | 0 1 0 | | 0 | 285 | | | p0 | | | 286 | 0 -1 0 | | | | -1 | 287 | | . | p1 | >= | | 288 | 0 0 1 | | | | 0 | 289 | | | p2 | | | 290 | 0 0 -1 | | -1 | 291 | | | | 292 |-1 -1 -1 | | -1 | 293 | | | | 294 | 1 1 1 | | 0 | 295 296 This example is for a 4-state model, the last probability pn is not included as this 297 parameter does not exist (because the sum of pc is equal to 1). The Euler angle parameters 298 have been excluded here but will be included in the returned A and b objects. These 299 parameters simply add columns of zero to the A matrix and have no effect on b. The last two 300 rows correspond to the inequality:: 301 302 0 <= pN <= 1. 303 304 As:: 305 N-1 306 \ 307 pN = 1 - > pc, 308 /__ 309 c=1 310 311 then:: 312 313 -p1 - p2 - ... - p(N-1) >= -1, 314 315 p1 + p2 + ... + p(N-1) >= 0. 316 317 318 @keyword data_types: The base data types used in the optimisation. This list can 319 contain the elements 'rdc', 'pcs' or 'tensor'. 320 @type data_types: list of str 321 @keyword scaling_matrix: The diagonal scaling matrix. 322 @type scaling_matrix: numpy rank-2 square matrix 323 @return: The matrices A and b. 324 @rtype: tuple of len 2 of a numpy rank-2, size NxM matrix and numpy 325 rank-1, size N array 326 """ 327 328 # Starting point of the populations. 329 pop_start = 0 330 if ('rdc' in data_types or 'pcs' in data_types) and not align_tensor.all_tensors_fixed(): 331 # Loop over the alignments. 332 for i in range(len(cdp.align_tensors)): 333 # Skip non-optimised tensors. 334 if not opt_uses_tensor(cdp.align_tensors[i]): 335 continue 336 337 # Add 5 parameters. 338 pop_start += 5 339 340 # Initialisation (0..j..m). 341 A = [] 342 b = [] 343 zero_array = zeros(param_num(), float64) 344 i = pop_start 345 j = 0 346 347 # Probability parameters. 348 if cdp.model in ['2-domain', 'population']: 349 # Loop over the prob parameters (N - 1, because the sum of pc is 1). 350 for k in range(cdp.N - 1): 351 # 0 <= pc <= 1. 352 A.append(zero_array * 0.0) 353 A.append(zero_array * 0.0) 354 A[j][i] = 1.0 355 A[j+1][i] = -1.0 356 b.append(0.0) 357 b.append(-1.0 / scaling_matrix[i, i]) 358 j = j + 2 359 360 # Increment i. 361 i = i + 1 362 363 # Add the inequalities for pN. 364 A.append(zero_array * 0.0) 365 A.append(zero_array * 0.0) 366 for i in range(pop_start, pop_start+cdp.N-1): 367 A[-2][i] = -1.0 368 A[-1][i] = 1.0 369 b.append(-1.0 / scaling_matrix[i, i]) 370 b.append(0.0) 371 372 # Convert to numpy data structures. 373 A = array(A, float64) 374 b = array(b, float64) 375 376 # Return the constraint objects. 377 return A, b

378 379

380 -def number_of_states(N=None):

381 """Set the number of states in the N-state model. 382 383 @param N: The number of states. 384 @type N: int 385 """ 386 387 # Test if the current data pipe exists. 388 check_pipe() 389 390 # Set the value of N. 391 cdp.N = N 392 393 # Update the model, if set. 394 if hasattr(cdp, 'model'): 395 update_model()

396 397

398 -def param_model_index(param):

399 """Return the N-state model index for the given parameter. 400 401 @param param: The N-state model parameter. 402 @type param: str 403 @return: The N-state model index. 404 @rtype: str 405 """ 406 407 # Probability. 408 if search('^p[0-9]*$', param): 409 return int(param[1:]) 410 411 # Alpha Euler angle. 412 if search('^alpha', param): 413 return int(param[5:]) 414 415 # Beta Euler angle. 416 if search('^beta', param): 417 return int(param[4:]) 418 419 # Gamma Euler angle. 420 if search('^gamma', param): 421 return int(param[5:]) 422 423 # Model independent parameter. 424 return None

425 426

427 -def param_num():

428 """Determine the number of parameters in the model. 429 430 @return: The number of model parameters. 431 @rtype: int 432 """ 433 434 # Determine the data type. 435 data_types = base_data_types() 436 437 # Init. 438 num = 0 439 440 # Alignment tensor params. 441 if ('rdc' in data_types or 'pcs' in data_types) and not align_tensor.all_tensors_fixed(): 442 # Loop over the alignments. 443 for i in range(len(cdp.align_tensors)): 444 # Skip non-optimised tensors. 445 if not opt_uses_tensor(cdp.align_tensors[i]): 446 continue 447 448 # Add 5 tensor parameters. 449 num += 5 450 451 # Populations. 452 if cdp.model in ['2-domain', 'population']: 453 num = num + (cdp.N - 1) 454 455 # Euler angles. 456 if cdp.model == '2-domain': 457 num = num + 3*cdp.N 458 459 # The paramagnetic centre. 460 if hasattr(cdp, 'paramag_centre_fixed') and not cdp.paramag_centre_fixed: 461 num = num + 3 462 463 # Return the param number. 464 return num

465 466

467 -def ref_domain(ref=None):

468 """Set the reference domain for the '2-domain' N-state model. 469 470 @param ref: The reference domain. 471 @type ref: str 472 """ 473 474 # Test if the current data pipe exists. 475 check_pipe() 476 477 # Test if the model is setup. 478 if not hasattr(cdp, 'model'): 479 raise RelaxNoModelError('N-state') 480 481 # Test that the correct model is set. 482 if cdp.model != '2-domain': 483 raise RelaxError("Setting the reference domain is only possible for the '2-domain' N-state model.") 484 485 # Test if the reference domain exists. 486 exists = False 487 for tensor_cont in cdp.align_tensors: 488 if tensor_cont.domain == ref: 489 exists = True 490 if not exists: 491 raise RelaxError("The reference domain cannot be found within any of the loaded tensors.") 492 493 # Set the reference domain. 494 cdp.ref_domain = ref 495 496 # Update the model. 497 update_model()

498 499

500 -def select_model(model=None):

501 """Select the N-state model type. 502 503 @param model: The N-state model type. Can be one of '2-domain', 'population', or 'fixed'. 504 @type model: str 505 """ 506 507 # Test if the current data pipe exists. 508 check_pipe() 509 510 # Test if the model name exists. 511 if not model in ['2-domain', 'population', 'fixed']: 512 raise RelaxError("The model name " + repr(model) + " is invalid.") 513 514 # Test if the model is setup. 515 if hasattr(cdp, 'model'): 516 warn(RelaxWarning("The N-state model has already been set up. Switching from model '%s' to '%s'." % (cdp.model, model))) 517 518 # Set the model 519 cdp.model = model 520 521 # Initialise the list of model parameters. 522 cdp.params = [] 523 524 # Update the model. 525 update_model()

526 527

528 -def update_model():

529 """Update the model parameters as necessary.""" 530 531 # Initialise the list of model parameters. 532 if not hasattr(cdp, 'params'): 533 cdp.params = [] 534 535 # Determine the number of states (loaded as structural models), if not already set. 536 if not hasattr(cdp, 'N'): 537 # Set the number. 538 if hasattr(cdp, 'structure'): 539 cdp.N = cdp.structure.num_models() 540 541 # Otherwise return as the rest cannot be updated without N. 542 else: 543 return 544 545 # Set up the parameter arrays. 546 if not cdp.params: 547 # Add the probability or population weight parameters. 548 if cdp.model in ['2-domain', 'population']: 549 for i in range(cdp.N-1): 550 cdp.params.append('p' + repr(i)) 551 552 # Add the Euler angle parameters. 553 if cdp.model == '2-domain': 554 for i in range(cdp.N): 555 cdp.params.append('alpha' + repr(i)) 556 cdp.params.append('beta' + repr(i)) 557 cdp.params.append('gamma' + repr(i)) 558 559 # Initialise the probability and Euler angle arrays. 560 if cdp.model in ['2-domain', 'population']: 561 if not hasattr(cdp, 'probs'): 562 cdp.probs = [None] * cdp.N 563 if cdp.model == '2-domain': 564 if not hasattr(cdp, 'alpha'): 565 cdp.alpha = [None] * cdp.N 566 if not hasattr(cdp, 'beta'): 567 cdp.beta = [None] * cdp.N 568 if not hasattr(cdp, 'gamma'): 569 cdp.gamma = [None] * cdp.N 570 571 # Determine the data type. 572 data_types = base_data_types() 573 574 # Set up tensors for each alignment. 575 if hasattr(cdp, 'align_ids'): 576 for id in cdp.align_ids: 577 # No tensors initialised. 578 if not hasattr(cdp, 'align_tensors'): 579 align_tensor.init(tensor=id, align_id=id) 580 581 # Find if the tensor corresponding to the id exists. 582 exists = False 583 for tensor in cdp.align_tensors: 584 if id == tensor.align_id: 585 exists = True 586 587 # Initialise the tensor. 588 if not exists: 589 align_tensor.init(tensor=id, align_id=id)

590

Source Code for Module specific_analyses.n_state_model.parameters