maths_fns.n_state

40 """Class containing the target function of the optimisation of the N-state model.""" 41

42 - def __init__(self, model=None, N=None, init_params=None, probs=None, full_tensors=None, red_data=None, red_errors=None, full_in_ref_frame=None, fixed_tensors=None, pcs=None, pcs_errors=None, pcs_weights=None, rdcs=None, rdc_errors=None, rdc_weights=None, xh_vect=None, temp=None, frq=None, dip_const=None, atomic_pos=None, paramag_centre=None, scaling_matrix=None, centre_fixed=True):

43 """Set up the class instance for optimisation. 44 45 The N-state models 46 ================== 47 48 All constant data required for the N-state model are initialised here. Depending on the base data used for optimisation, different data structures can be supplied. However a number of structures must be provided for the N-state model. These are: 49 50 - model, the type of N-state model. This can be '2-domain', 'population', or 'fixed'. 51 - N, the number of states (or structures). 52 - init_params, the initial parameter values. 53 - scaling_matrix, the matrix used for parameter scaling during optimisation. 54 55 56 2-domain N-state model 57 ---------------------- 58 59 If the model type is set to '2-domain', then the following data structures should be supplied: 60 61 - full_tensors, the alignment tensors in matrix form. 62 - red_data, the alignment tensors in 5D form in a rank-1 array. 63 - red_errors, the alignment tensor errors in 5D form in a rank-1 array. This data is not obligatory. 64 - full_in_ref_frame, an array of flags specifying if the tensor in the reference frame is the full or reduced tensor. 65 66 67 The population N-state model 68 ============================ 69 70 In this model, populations are optimised for each state. Additionally the alignment tensors for anisotropic data can also be optimised if they have not been supplied (through the full_tensors arg). 71 72 73 PCS base data 74 ------------- 75 76 If pseudocontact shift data is to be used for optimisation, then the following should be supplied: 77 78 - pcs, the pseudocontact shifts. 79 - pcs_errors, the optional pseudocontact shift error values. 80 - temp, the temperatures of the experiments. 81 - frq, the frequencies of the experiments. 82 83 84 PCS and PRE base data 85 --------------------- 86 87 If either pseudocontact shift or PRE data is to be used for optimisation, then the following should be supplied: 88 89 - atomic_pos, the positions of all atoms. 90 - paramag_centre, the paramagnetic centre position. 91 92 93 RDC base data 94 ------------- 95 96 If residual dipolar coupling data is to be used for optimisation, then the following should be supplied: 97 98 - rdcs, the residual dipolar couplings. 99 - rdc_errors, the optional residual dipolar coupling errors. 100 - xh_vect, the heteronucleus to proton unit vectors. 101 - dip_const, the dipolar constants. 102 103 104 @keyword model: The N-state model type. This can be one of '2-domain', 'population' or 'fixed'. 105 @type model: str 106 @keyword N: The number of states. 107 @type N: int 108 @keyword init_params: The initial parameter values. Optimisation must start at some point! 109 @type init_params: numpy float64 array 110 @keyword probs: The probabilities for each state c. The length of this list should be equal to N. 111 @type probs: list of float 112 @keyword full_tensors: An array of the {Sxx, Syy, Sxy, Sxz, Syz} values for all full tensors. The format is [Sxx1, Syy1, Sxy1, Sxz1, Syz1, Sxx2, Syy2, Sxy2, Sxz2, Syz2, ..., Sxxn, Syyn, Sxyn, Sxzn, Syzn] 113 @type full_tensors: list of rank-2, 3D numpy arrays 114 @keyword red_data: An array of the {Sxx, Syy, Sxy, Sxz, Syz} values for all reduced tensors. The format is the same as for full_tensors. 115 @type red_data: numpy float64 array 116 @keyword red_errors: An array of the {Sxx, Syy, Sxy, Sxz, Syz} errors for all reduced tensors. The array format is the same as for full_tensors. 117 @type red_errors: numpy float64 array 118 @keyword full_in_ref_frame: An array of flags specifying if the tensor in the reference frame is the full or reduced tensor. 119 @type full_in_ref_frame: numpy rank-1 array 120 @keyword fixed_tensors: An array of flags specifying if the tensor is fixed or will be optimised. 121 @type fixed_tensors: list of bool 122 @keyword pcs: The PCS lists. The first index must correspond to the different alignment media i and the second index to the spin systems j. 123 @type pcs: numpy rank-2 array 124 @keyword pcs_errors: The PCS error lists. The dimensions of this argument are the same as for 'pcs'. 125 @type pcs_errors: numpy rank-2 array 126 @keyword pcs_weights: The PCS weight lists. The dimensions of this argument are the same as for 'pcs'. 127 @type pcs_weights: numpy rank-2 array 128 @keyword rdcs: The RDC lists. The first index must correspond to the different alignment media i and the second index to the spin systems j. 129 @type rdcs: numpy rank-2 array 130 @keyword rdc_errors: The RDC error lists. The dimensions of this argument are the same as for 'rdcs'. 131 @type rdc_errors: numpy rank-2 array 132 @keyword rdc_weights: The RDC weight lists. The dimensions of this argument are the same as for 'rdcs'. 133 @type rdc_weights: numpy rank-2 array 134 @keyword xh_vect: The unit XH vector lists. The first index must correspond to the spin systems and the second index to each structure (its size being equal to the number of states). 135 @type xh_vect: numpy rank-2 array 136 @keyword temp: The temperature of each experiment, used for the PCS. 137 @type temp: numpy rank-1 array 138 @keyword frq: The frequency of each alignment, used for the PCS. 139 @type frq: numpy rank-1 array 140 @keyword dip_const: The dipolar constants for each XH vector. The indices correspond to the spin systems j. 141 @type dip_const: numpy rank-1 array 142 @keyword atomic_pos: The atomic positions of all spins. The first index is the spin systems j and the second is the structure or state c. 143 @type atomic_pos: numpy rank-3 array 144 @keyword paramag_centre: The paramagnetic centre position (or positions). 145 @type paramag_centre: numpy rank-1, 3D array or rank-2, Nx3 array 146 @keyword scaling_matrix: The square and diagonal scaling matrix. 147 @type scaling_matrix: numpy rank-2 array 148 @keyword centre_fixed: A flag which if False will cause the paramagnetic centre to be optimised. 149 @type centre_fixed: bool 150 """ 151 152 # Store the data inside the class instance namespace. 153 self.N = N 154 self.params = 1.0 * init_params # Force a copy of the data to be stored. 155 self.fixed_tensors = fixed_tensors 156 self.deltaij = pcs 157 self.Dij = rdcs 158 self.dip_vect = xh_vect 159 self.dip_const = dip_const 160 self.temp = temp 161 self.frq = frq 162 self.atomic_pos = atomic_pos 163 self.paramag_centre = paramag_centre 164 self.centre_fixed = centre_fixed 165 self.total_num_params = len(init_params) 166 167 # Initialise the function value, gradient, and Hessian. 168 self.chi2 = 0.0 169 self.dchi2 = zeros((self.total_num_params), float64) 170 self.d2chi2 = zeros((self.total_num_params, self.total_num_params), float64) 171 172 # Scaling initialisation. 173 self.scaling_matrix = scaling_matrix 174 if self.scaling_matrix != None: 175 self.scaling_flag = True 176 else: 177 self.scaling_flag = False 178 179 # The 2-domain N-state model. 180 if model == '2-domain': 181 # Some checks. 182 if red_data == None or not len(red_data): 183 raise RelaxError("The red_data argument " + repr(red_data) + " must be supplied.") 184 if red_errors == None or not len(red_errors): 185 raise RelaxError("The red_errors argument " + repr(red_errors) + " must be supplied.") 186 if full_in_ref_frame == None or not len(full_in_ref_frame): 187 raise RelaxError("The full_in_ref_frame argument " + repr(full_in_ref_frame) + " must be supplied.") 188 189 # Tensor set up. 190 self.full_tensors = array(full_tensors, float64) 191 self.num_tensors = len(self.full_tensors) / 5 192 self.red_data = red_data 193 self.red_errors = red_errors 194 self.full_in_ref_frame = full_in_ref_frame 195 196 # Alignment tensor in rank-2, 3D form. 197 self.A = zeros((self.num_tensors, 3, 3), float64) 198 for i in range(self.num_tensors): 199 to_tensor(self.A[i], self.full_tensors[5*i:5*i+5]) 200 201 # Initialise some empty numpy objects for storage of: 202 # R: the transient rotation matricies. 203 # RT: the transposes of the rotation matricies. 204 # red_bc: the back-calculated reduced alignment tensors. 205 # red_bc_vector: the back-calculated reduced alignment tensors in vector form {Sxx, Syy, Sxy, Sxz, Syz}. 206 self.R = zeros((self.N, 3, 3), float64) 207 self.RT = zeros((self.N, 3, 3), float64) 208 self.red_bc = zeros((self.num_tensors, 3, 3), float64) 209 self.red_bc_vector = zeros(self.num_tensors*5, float64) 210 211 # Set the target function. 212 self.func = self.func_2domain 213 self.dfunc = None 214 self.d2func = None 215 216 # The flexible population or equal probability N-state models. 217 elif model in ['population', 'fixed']: 218 # Set the RDC and PCS flags (indicating the presence of data). 219 self.rdc_flag = True 220 self.pcs_flag = True 221 if rdcs == None or len(rdcs) == 0: 222 self.rdc_flag = False 223 if pcs == None or len(pcs) == 0: 224 self.pcs_flag = False 225 226 # Some checks. 227 if self.rdc_flag and (xh_vect == None or not len(xh_vect)): 228 raise RelaxError("The xh_vect argument " + repr(xh_vect) + " must be supplied.") 229 if self.pcs_flag and (atomic_pos == None or not len(atomic_pos)): 230 raise RelaxError("The atomic_pos argument " + repr(atomic_pos) + " must be supplied.") 231 232 # The total number of spins. 233 self.num_spins = 0 234 if self.rdc_flag: 235 self.num_spins = len(rdcs[0]) 236 elif self.pcs_flag: 237 self.num_spins = len(pcs[0]) 238 239 # The total number of alignments. 240 self.num_align = 0 241 if self.rdc_flag: 242 self.num_align = len(rdcs) 243 elif self.pcs_flag: 244 self.num_align = len(pcs) 245 246 # Alignment tensor function and gradient matrices. 247 self.A = zeros((self.num_align, 3, 3), float64) 248 self.dA = zeros((5, 3, 3), float64) 249 250 # Fixed alignment tensors. 251 if full_tensors != None: 252 # Convert to numpy. 253 self.full_tensors = array(full_tensors, float64) 254 255 # Set up the alignment data. 256 self.num_align_params = 0 257 index = 0 258 for i in range(self.num_align): 259 # Fill the alignment tensor object with the fixed tensors. 260 if fixed_tensors[i]: 261 to_tensor(self.A[i], self.full_tensors[5*index:5*index+5]) 262 index += 1 263 264 # The number of alignment parameters. 265 if not fixed_tensors[i]: 266 self.num_align_params += 5 267 268 # Optimised alignment tensors. 269 else: 270 # The alignment tensor gradients don't change, so pre-calculate them. 271 dAi_dAxx(self.dA[0]) 272 dAi_dAyy(self.dA[1]) 273 dAi_dAxy(self.dA[2]) 274 dAi_dAxz(self.dA[3]) 275 dAi_dAyz(self.dA[4]) 276 277 # PCS errors. 278 if self.pcs_flag: 279 err = False 280 for i in xrange(len(pcs_errors)): 281 for j in xrange(len(pcs_errors[i])): 282 if not isNaN(pcs_errors[i, j]): 283 err = True 284 if err: 285 self.pcs_sigma_ij = pcs_errors 286 else: 287 # Missing errors (the values need to be small, close to ppm units, so the chi-squared value is comparable to the RDC). 288 self.pcs_sigma_ij = 0.03 * 1e-6 * ones((self.num_align, self.num_spins), float64) 289 290 # RDC errors. 291 if self.rdc_flag: 292 err = False 293 for i in xrange(len(rdc_errors)): 294 for j in xrange(len(rdc_errors[i])): 295 if not isNaN(rdc_errors[i, j]): 296 err = True 297 if err: 298 self.rdc_sigma_ij = rdc_errors 299 else: 300 # Missing errors. 301 self.rdc_sigma_ij = ones((self.num_align, self.num_spins), float64) 302 303 # Missing data matrices (RDC). 304 if self.rdc_flag: 305 self.missing_Dij = zeros((self.num_align, self.num_spins), float64) 306 307 # Missing data matrices (PCS). 308 if self.pcs_flag: 309 self.missing_deltaij = zeros((self.num_align, self.num_spins), float64) 310 311 # Clean up problematic data and put the weights into the errors.. 312 if self.rdc_flag or self.pcs_flag: 313 for i in xrange(self.num_align): 314 for j in xrange(self.num_spins): 315 if self.rdc_flag: 316 if isNaN(self.Dij[i, j]): 317 # Set the flag. 318 self.missing_Dij[i, j] = 1 319 320 # Change the NaN to zero. 321 self.Dij[i, j] = 0.0 322 323 # Change the error to one (to avoid zero division). 324 self.rdc_sigma_ij[i, j] = 1.0 325 326 # Change the weight to one. 327 rdc_weights[i, j] = 1.0 328 329 if self.pcs_flag: 330 if isNaN(self.deltaij[i, j]): 331 # Set the flag. 332 self.missing_deltaij[i, j] = 1 333 334 # Change the NaN to zero. 335 self.deltaij[i, j] = 0.0 336 337 # Change the error to one (to avoid zero division). 338 self.pcs_sigma_ij[i, j] = 1.0 339 340 # Change the weight to one. 341 pcs_weights[i, j] = 1.0 342 343 # The RDC weights. 344 if self.rdc_flag: 345 self.rdc_sigma_ij[i, j] = self.rdc_sigma_ij[i, j] / sqrt(rdc_weights[i, j]) 346 347 # The PCS weights. 348 if self.pcs_flag: 349 self.pcs_sigma_ij[i, j] = self.pcs_sigma_ij[i, j] / sqrt(pcs_weights[i, j]) 350 351 352 # The paramagnetic centre vectors and distances. 353 if self.pcs_flag: 354 # Initialise the data structures. 355 self.paramag_unit_vect = zeros(atomic_pos.shape, float64) 356 self.paramag_dist = zeros((self.num_spins, self.N), float64) 357 self.pcs_const = zeros((self.num_align, self.num_spins, self.N), float64) 358 if self.paramag_centre == None: 359 self.paramag_centre = zeros(3, float64) 360 361 # Set up the paramagnetic info. 362 self.paramag_info() 363 364 # PCS function, gradient, and Hessian matrices. 365 self.deltaij_theta = zeros((self.num_align, self.num_spins), float64) 366 self.ddeltaij_theta = zeros((self.total_num_params, self.num_align, self.num_spins), float64) 367 self.d2deltaij_theta = zeros((self.total_num_params, self.total_num_params, self.num_align, self.num_spins), float64) 368 369 # RDC function, gradient, and Hessian matrices. 370 self.Dij_theta = zeros((self.num_align, self.num_spins), float64) 371 self.dDij_theta = zeros((self.total_num_params, self.num_align, self.num_spins), float64) 372 self.d2Dij_theta = zeros((self.total_num_params, self.total_num_params, self.num_align, self.num_spins), float64) 373 374 # Set the target function, gradient, and Hessian (paramagnetic centre optimisation). 375 if not self.centre_fixed: 376 self.func = self.func_population 377 self.dfunc = None 378 self.d2func = None 379 380 # Set the standard target function, gradient, and Hessian. 381 else: 382 self.func = self.func_population 383 self.dfunc = self.dfunc_population 384 self.d2func = self.d2func_population 385 386 # Fixed probabilities. 387 if model == 'fixed': 388 # The probs are unpacked by self.func in the population model, so just override that function. 389 self.func = self.func_tensor_opt 390 self.dfunc = self.dfunc_tensor_opt 391 self.d2func = self.d2func_tensor_opt 392 393 # The zero Hessian. 394 self.zero_hessian = zeros(self.num_spins, float64) 395 396 # The probability array. 397 if probs: 398 self.probs = probs 399 400 # All structures have initial equal probability. 401 else: 402 self.probs = ones(self.N, float64) / self.N

403 404

405 - def func_2domain(self, params):

406 """The target function for optimisation of the 2-domain N-state model. 407 408 This function should be passed to the optimisation algorithm. It accepts, as an array, a vector of parameter values and, using these, returns the single chi-squared value corresponding to that coordinate in the parameter space. If no tensor errors are supplied, then the SSE (the sum of squares error) value is returned instead. The chi-squared is simply the SSE normalised to unit variance (the SSE divided by the error squared). 409 410 @param params: The vector of parameter values. 411 @type params: list of float 412 @return: The chi-squared or SSE value. 413 @rtype: float 414 """ 415 416 # Scaling. 417 if self.scaling_flag: 418 params = dot(params, self.scaling_matrix) 419 420 # Reset the back-calculated the reduced tensor structure. 421 self.red_bc = self.red_bc * 0.0 422 423 # Loop over the N states. 424 for c in xrange(self.N): 425 # The rotation matrix. 426 euler_to_R_zyz(params[self.N-1+3*c], params[self.N-1+3*c+1], params[self.N-1+3*c+2], self.R[c]) 427 428 # Its transpose. 429 self.RT[c] = transpose(self.R[c]) 430 431 # The probability of state c. 432 if c < self.N-1: 433 pc = params[c] 434 435 # The probability of state N (1 minus the pc of all other states). 436 else: 437 pc = 1.0 438 for c2 in xrange(self.N-1): 439 pc = pc - params[c2] 440 441 # Back-calculate the reduced tensors for sum element c and add these to red_bc. 442 for i in xrange(self.num_tensors): 443 # Normal RT.X.R rotation. 444 if self.full_in_ref_frame[i]: 445 self.red_bc[i] = self.red_bc[i] + pc * dot(self.RT[c], dot(self.A[i], self.R[c])) 446 447 # Inverse R.X.RT rotation. 448 else: 449 self.red_bc[i] = self.red_bc[i] + pc * dot(self.R[c], dot(self.A[i], self.RT[c])) 450 451 # 5D vectorise the back-calculated tensors (create red_bc_vector from red_bc). 452 for i in xrange(self.num_tensors): 453 self.red_bc_vector[5*i] = self.red_bc[i, 0, 0] # Sxx. 454 self.red_bc_vector[5*i+1] = self.red_bc[i, 1, 1] # Syy. 455 self.red_bc_vector[5*i+2] = self.red_bc[i, 0, 1] # Sxy. 456 self.red_bc_vector[5*i+3] = self.red_bc[i, 0, 2] # Sxz. 457 self.red_bc_vector[5*i+4] = self.red_bc[i, 1, 2] # Syz. 458 459 # Return the chi-squared value. 460 return chi2(self.red_data, self.red_bc_vector, self.red_errors)

461 462

463 - def func_population(self, params):

464 """The target function for optimisation of the flexible population N-state model. 465 466 Description 467 =========== 468 469 This function should be passed to the optimisation algorithm. It accepts, as an array, a vector of parameter values and, using these, returns the single chi-squared value corresponding to that coordinate in the parameter space. If no RDC errors are supplied, then the SSE (the sum of squares error) value is returned instead. The chi-squared is simply the SSE normalised to unit variance (the SSE divided by the error squared). 470 471 472 Indices 473 ======= 474 475 For this calculation, five indices are looped over and used in the various data structures. These include: 476 - i, the index over alignments, 477 - j, the index over spin systems, 478 - c, the index over the N-states (or over the structures), 479 - n, the index over the first dimension of the alignment tensor n = {x, y, z}, 480 - m, the index over the second dimension of the alignment tensor m = {x, y, z}. 481 482 483 Equations 484 ========= 485 486 To calculate the function value, a chain of equations are used. This includes the chi-squared equation and the RDC equation. 487 488 489 The chi-squared equation 490 ------------------------ 491 492 The equations are:: 493 494 ___ 495 \ (Dij - Dij(theta)) ** 2 496 chi^2(theta) = > ----------------------- , 497 /__ sigma_ij ** 2 498 ij 499 500 ___ 501 \ (delta_ij - delta_ij(theta)) ** 2 502 chi^2(theta) = > --------------------------------- , 503 /__ sigma_ij ** 2 504 ij 505 506 where: 507 - theta is the parameter vector, 508 - Dij are the measured RDCs for alignment i, spin j, 509 - Dij(theta) are the back calculated RDCs for alignment i, spin j, 510 - delta_ij are the measured PCSs for alignment i, spin j, 511 - delta_ij(theta) are the back calculated PCSs for alignment i, spin j, 512 - sigma_ij are the RDC or PCS errors. 513 514 Both chi-squared values sum. 515 516 517 The RDC equation 518 ---------------- 519 520 The RDC equation is:: 521 522 _N_ 523 \ T 524 Dij(theta) = dj > pc . mu_jc . Ai . mu_jc, 525 /__ 526 c=1 527 528 where: 529 - dj is the dipolar constant for spin j, 530 - N is the total number of states or structures, 531 - pc is the weight or probability associated with state c, 532 - mu_jc is the unit vector corresponding to spin j and state c, 533 - Ai is the alignment tensor. 534 535 The dipolar constant is henceforth defined as:: 536 537 dj = 3 / (2pi) d', 538 539 where the factor of 2pi is to convert from units of rad.s^-1 to Hertz, the factor of 3 is associated with the alignment tensor and the pure dipolar constant in SI units is:: 540 541 mu0 gI.gS.h_bar 542 d' = - --- ----------- , 543 4pi r**3 544 545 where: 546 - mu0 is the permeability of free space, 547 - gI and gS are the gyromagnetic ratios of the I and S spins, 548 - h_bar is Dirac's constant which is equal to Planck's constant divided by 2pi, 549 - r is the distance between the two spins. 550 551 552 The PCS equation 553 ---------------- 554 555 The PCS equation is:: 556 557 _N_ 558 \ T 559 delta_ij(theta) = > pc . dijc . mu_jc . Ai . mu_jc, 560 /__ 561 c=1 562 563 where: 564 - djci is the PCS constant for spin j, state c and experiment or alignment i, 565 - N is the total number of states or structures, 566 - pc is the weight or probability associated with state c, 567 - mu_jc is the unit vector corresponding to spin j and state c, 568 - Ai is the alignment tensor. 569 570 The PCS constant is defined as:: 571 572 mu0 15kT 1 573 dijc = --- ----- ---- , 574 4pi Bo**2 r**3 575 576 where: 577 - mu0 is the permeability of free space, 578 - k is Boltzmann's constant, 579 - T is the absolute temperature (different for each experiment), 580 - Bo is the magnetic field strength (different for each experiment), 581 - r is the distance between the paramagnetic centre (electron spin) and the nuclear spin (different for each spin and state). 582 583 584 Stored data structures 585 ====================== 586 587 There are a number of data structures calculated by this function and stored for subsequent use in the gradient and Hessian functions. This include the back calculated RDCs and the alignment tensors. 588 589 Dij(theta) 590 ---------- 591 592 The back calculated RDCs. This is a rank-2 tensor with indices {i, j}. 593 594 delta_ij(theta) 595 --------------- 596 597 The back calculated PCS. This is a rank-2 tensor with indices {i, j}. 598 599 Ai 600 -- 601 602 The alignment tensors. This is a rank-3 tensor with indices {i, n, m}. 603 604 605 @param params: The vector of parameter values. 606 @type params: numpy rank-1 array 607 @return: The chi-squared or SSE value. 608 @rtype: float 609 """ 610 611 # Scaling. 612 if self.scaling_flag: 613 params = dot(params, self.scaling_matrix) 614 615 # Initial chi-squared (or SSE) value. 616 chi2_sum = 0.0 617 618 # Unpack the probabilities (located at the end of the parameter array). 619 if self.N > 1: 620 self.probs = params[-(self.N-1):] 621 622 # Unpack the paramagnetic centre. 623 if not self.centre_fixed: 624 # The position. 625 self.paramag_centre = params[-3:] 626 627 # Update the paramagnetic info. 628 self.paramag_info() 629 630 # Loop over each alignment. 631 index = 0 632 for i in xrange(self.num_align): 633 # Create tensor i from the parameters. 634 if not self.fixed_tensors[i]: 635 to_tensor(self.A[i], params[5*index:5*index + 5]) 636 index += 1 637 638 # Loop over the spin systems j. 639 for j in xrange(self.num_spins): 640 # The back calculated RDC. 641 if self.rdc_flag: 642 # Calculate the average RDC. 643 if not self.missing_Dij[i, j]: 644 self.Dij_theta[i, j] = ave_rdc_tensor(self.dip_const[j], self.dip_vect[j], self.N, self.A[i], weights=self.probs) 645 646 # The back calculated PCS. 647 if self.pcs_flag: 648 # Calculate the average PCS. 649 if not self.missing_deltaij[i, j]: 650 self.deltaij_theta[i, j] = ave_pcs_tensor(self.pcs_const[i, j], self.paramag_unit_vect[j], self.N, self.A[i], weights=self.probs) 651 652 # Calculate and sum the single alignment chi-squared value (for the RDC). 653 if self.rdc_flag: 654 chi2_sum = chi2_sum + chi2(self.Dij[i], self.Dij_theta[i], self.rdc_sigma_ij[i]) 655 656 # Calculate and sum the single alignment chi-squared value (for the PCS). 657 if self.pcs_flag: 658 chi2_sum = chi2_sum + chi2(self.deltaij[i], self.deltaij_theta[i], self.pcs_sigma_ij[i]) 659 660 # Return the chi-squared value. 661 return chi2_sum

662 663

664 - def func_tensor_opt(self, params):

665 """The target function for optimisation of the alignment tensor from RDC and/or PCS data. 666 667 Description 668 =========== 669 670 This function should be passed to the optimisation algorithm. It accepts, as an array, a vector of parameter values and, using these, returns the single chi-squared value corresponding to that coordinate in the parameter space. If no RDC or PCS errors are supplied, then the SSE (the sum of squares error) value is returned instead. The chi-squared is simply the SSE normalised to unit variance (the SSE divided by the error squared). 671 672 673 Indices 674 ======= 675 676 For this calculation, five indices are looped over and used in the various data structures. These include: 677 - i, the index over alignments, 678 - j, the index over spin systems, 679 - c, the index over the N-states (or over the structures), 680 - n, the index over the first dimension of the alignment tensor n = {x, y, z}, 681 - m, the index over the second dimension of the alignment tensor m = {x, y, z}. 682 683 684 Equations 685 ========= 686 687 To calculate the function value, a chain of equations are used. This includes the chi-squared equation and the RDC and PCS equations. 688 689 690 The chi-squared equation 691 ------------------------ 692 693 The equations are:: 694 695 ___ 696 \ (Dij - Dij(theta)) ** 2 697 chi^2(theta) = > ----------------------- , 698 /__ sigma_ij ** 2 699 ij 700 701 ___ 702 \ (delta_ij - delta_ij(theta)) ** 2 703 chi^2(theta) = > --------------------------------- , 704 /__ sigma_ij ** 2 705 ij 706 707 where: 708 - theta is the parameter vector, 709 - Dij are the measured RDCs for alignment i, spin j, 710 - Dij(theta) are the back calculated RDCs for alignment i, spin j, 711 - delta_ij are the measured PCSs for alignment i, spin j, 712 - delta_ij(theta) are the back calculated PCSs for alignment i, spin j, 713 - sigma_ij are the RDC or PCS errors. 714 715 Both chi-squared values sum. 716 717 718 The RDC equation 719 ---------------- 720 721 The RDC equation is:: 722 723 _N_ 724 dj \ T 725 Dij(theta) = -- > mu_jc . Ai . mu_jc, 726 N /__ 727 c=1 728 729 where: 730 - dj is the dipolar constant for spin j, 731 - N is the total number of states or structures, 732 - mu_jc is the unit vector corresponding to spin j and state c, 733 - Ai is the alignment tensor. 734 735 The dipolar constant is henceforth defined as:: 736 737 dj = 3 / (2pi) d', 738 739 where the factor of 2pi is to convert from units of rad.s^-1 to Hertz, the factor of 3 is associated with the alignment tensor and the pure dipolar constant in SI units is:: 740 741 mu0 gI.gS.h_bar 742 d' = - --- ----------- , 743 4pi r**3 744 745 where: 746 - mu0 is the permeability of free space, 747 - gI and gS are the gyromagnetic ratios of the I and S spins, 748 - h_bar is Dirac's constant which is equal to Planck's constant divided by 2pi, 749 - r is the distance between the two spins. 750 751 752 The PCS equation 753 ---------------- 754 755 The PCS equation is:: 756 757 _N_ 758 1 \ T 759 delta_ij(theta) = - > dijc . mu_jc . Ai . mu_jc, 760 N /__ 761 c=1 762 763 where: 764 - djci is the PCS constant for spin j, state c and experiment or alignment i, 765 - N is the total number of states or structures, 766 - mu_jc is the unit vector corresponding to spin j and state c, 767 - Ai is the alignment tensor. 768 769 The PCS constant is defined as:: 770 771 mu0 15kT 1 772 dijc = --- ----- ---- , 773 4pi Bo**2 r**3 774 775 where: 776 - mu0 is the permeability of free space, 777 - k is Boltzmann's constant, 778 - T is the absolute temperature (different for each experiment), 779 - Bo is the magnetic field strength (different for each experiment), 780 - r is the distance between the paramagnetic centre (electron spin) and the nuclear spin (different for each spin and state). 781 782 783 Stored data structures 784 ====================== 785 786 There are a number of data structures calculated by this function and stored for subsequent use in the gradient and Hessian functions. This include the back calculated RDCs and the alignment tensors. 787 788 Dij(theta) 789 ---------- 790 791 The back calculated RDCs. This is a rank-2 tensor with indices {i, j}. 792 793 delta_ij(theta) 794 --------------- 795 796 The back calculated PCS. This is a rank-2 tensor with indices {i, j}. 797 798 Ai 799 -- 800 801 The alignment tensors. This is a rank-3 tensor with indices {i, n, m}. 802 803 804 @param params: The vector of parameter values. 805 @type params: numpy rank-1 array 806 @return: The chi-squared or SSE value. 807 @rtype: float 808 """ 809 810 # Scaling. 811 if self.scaling_flag: 812 params = dot(params, self.scaling_matrix) 813 814 # Initial chi-squared (or SSE) value. 815 chi2_sum = 0.0 816 817 # Unpack the paramagnetic centre. 818 if not self.centre_fixed: 819 # The position. 820 self.paramag_centre = params[-3:] 821 822 # Update the paramagnetic info. 823 self.paramag_info() 824 825 # Loop over each alignment. 826 index = 0 827 for i in xrange(self.num_align): 828 # Create tensor i from the parameters. 829 if not self.fixed_tensors[i]: 830 to_tensor(self.A[i], params[5*index:5*index + 5]) 831 index += 1 832 833 # Loop over the spin systems j. 834 for j in xrange(self.num_spins): 835 # The back calculated RDC. 836 if self.rdc_flag: 837 # Calculate the average RDC. 838 if not self.missing_Dij[i, j]: 839 self.Dij_theta[i, j] = ave_rdc_tensor(self.dip_const[j], self.dip_vect[j], self.N, self.A[i], weights=self.probs) 840 841 # The back calculated PCS. 842 if self.pcs_flag: 843 # Calculate the average PCS. 844 if not self.missing_deltaij[i, j]: 845 self.deltaij_theta[i, j] = ave_pcs_tensor(self.pcs_const[i, j], self.paramag_unit_vect[j], self.N, self.A[i], weights=self.probs) 846 847 # Calculate and sum the single alignment chi-squared value (for the RDC). 848 if self.rdc_flag: 849 chi2_sum = chi2_sum + chi2(self.Dij[i], self.Dij_theta[i], self.rdc_sigma_ij[i]) 850 851 # Calculate and sum the single alignment chi-squared value (for the PCS). 852 if self.pcs_flag: 853 chi2_sum = chi2_sum + chi2(self.deltaij[i], self.deltaij_theta[i], self.pcs_sigma_ij[i]) 854 855 # Return the chi-squared value. 856 return chi2_sum

857 858

859 - def dfunc_population(self, params):

860 """The gradient function for optimisation of the flexible population N-state model. 861 862 Description 863 =========== 864 865 This function should be passed to the optimisation algorithm. It accepts, as an array, a vector of parameter values and, using these, returns the chi-squared gradient corresponding to that coordinate in the parameter space. If no RDC errors are supplied, then the SSE (the sum of squares error) gradient is returned instead. The chi-squared gradient is simply the SSE gradient normalised to unit variance (the SSE divided by the error squared). 866 867 868 Indices 869 ======= 870 871 For this calculation, six indices are looped over and used in the various data structures. These include: 872 - k, the index over all parameters, 873 - i, the index over alignments, 874 - j, the index over spin systems, 875 - c, the index over the N-states (or over the structures), 876 - m, the index over the first dimension of the alignment tensor m = {x, y, z}. 877 - n, the index over the second dimension of the alignment tensor n = {x, y, z}, 878 879 880 Equations 881 ========= 882 883 To calculate the chi-squared gradient, a chain of equations are used. This includes the chi-squared gradient, the RDC gradient and the alignment tensor gradient. 884 885 886 The chi-squared gradient 887 ------------------------ 888 889 The equation is:: 890 ___ 891 dchi^2(theta) \ / Dij - Dij(theta) dDij(theta) \ 892 ------------- = -2 > | ---------------- . ----------- | 893 dthetak /__ \ sigma_ij**2 dthetak / 894 ij 895 896 where: 897 - theta is the parameter vector, 898 - Dij are the measured RDCs or PCSs, 899 - Dij(theta) are the back calculated RDCs or PCSs, 900 - sigma_ij are the RDC or PCS errors, 901 - dDij(theta)/dthetak is the RDC or PCS gradient for parameter k. 902 903 904 The RDC gradient 905 ---------------- 906 907 This gradient is different for the various parameter types. 908 909 pc partial derivative 910 ~~~~~~~~~~~~~~~~~~~~~ 911 912 The population parameter partial derivative is:: 913 914 dDij(theta) T 915 ----------- = dj . mu_jc . Ai . mu_jc, 916 dpc 917 918 where: 919 - dj is the dipolar constant for spin j, 920 - mu_jc is the unit vector corresponding to spin j and state c, 921 - Ai is the alignment tensor. 922 923 Amn partial derivative 924 ~~~~~~~~~~~~~~~~~~~~~~ 925 926 The alignment tensor element partial derivative is:: 927 928 _N_ 929 dDij(theta) \ T dAi 930 ----------- = dj > pc . mu_jc . ---- . mu_jc, 931 dAmn /__ dAmn 932 c=1 933 934 where: 935 - dj is the dipolar constant for spin j, 936 - pc is the weight or probability associated with state c, 937 - mu_jc is the unit vector corresponding to spin j and state c, 938 - dAi/dAmn is the partial derivative of the alignment tensor with respect to element Amn. 939 940 941 The PCS gradient 942 ---------------- 943 944 This gradient is also different for the various parameter types. 945 946 pc partial derivative 947 ~~~~~~~~~~~~~~~~~~~~~ 948 949 The population parameter partial derivative is:: 950 951 ddeltaij(theta) T 952 --------------- = dijc . mu_jc . Ai . mu_jc, 953 dpc 954 955 where: 956 - djc is the pseudocontact shift constant for spin j and state c, 957 - mu_jc is the unit vector corresponding to spin j and state c, 958 - Ai is the alignment tensor. 959 960 Amn partial derivative 961 ~~~~~~~~~~~~~~~~~~~~~~ 962 963 The alignment tensor element partial derivative is:: 964 965 _N_ 966 ddelta_ij(theta) \ T dAi 967 ---------------- = > pc . djc . mu_jc . ---- . mu_jc, 968 dAmn /__ dAmn 969 c=1 970 971 where: 972 - djc is the pseudocontact shift constant for spin j and state c, 973 - pc is the weight or probability associated with state c, 974 - mu_jc is the unit vector corresponding to spin j and state c, 975 - dAi/dAmn is the partial derivative of the alignment tensor with respect to element Amn. 976 977 978 The alignment tensor gradient 979 ----------------------------- 980 981 The five unique elements of the tensor {Axx, Ayy, Axy, Axz, Ayz} give five different partial derivatives. These are:: 982 983 dAi | 1 0 0 | 984 ---- = | 0 0 0 |, 985 dAxx | 0 0 -1 | 986 987 dAi | 0 0 0 | 988 ---- = | 0 1 0 |, 989 dAyy | 0 0 -1 | 990 991 dAi | 0 1 0 | 992 ---- = | 1 0 0 |, 993 dAxy | 0 0 0 | 994 995 dAi | 0 0 1 | 996 ---- = | 0 0 0 |, 997 dAxz | 1 0 0 | 998 999 dAi | 0 0 0 | 1000 ---- = | 0 0 1 |. 1001 dAyz | 0 1 0 | 1002 1003 As these are invariant, they can be pre-calculated. 1004 1005 1006 Stored data structures 1007 ====================== 1008 1009 There are a number of data structures calculated by this function and stored for subsequent use in the Hessian function. This include the back calculated RDC and PCS gradients and the alignment tensor gradients. 1010 1011 dDij(theta)/dthetak 1012 ------------------- 1013 1014 The back calculated RDC gradient. This is a rank-3 tensor with indices {k, i, j}. 1015 1016 ddeltaij(theta)/dthetak 1017 ----------------------- 1018 1019 The back calculated PCS gradient. This is a rank-3 tensor with indices {k, i, j}. 1020 1021 dAi/dAmn 1022 -------- 1023 1024 The alignment tensor gradients. This is a rank-3 tensor with indices {5, n, m}. 1025 1026 1027 @param params: The vector of parameter values. This is unused as it is assumed that 1028 func_population() was called first. 1029 @type params: numpy rank-1 array 1030 @return: The chi-squared or SSE gradient. 1031 @rtype: numpy rank-1 array 1032 """ 1033 1034 # Scaling. 1035 if self.scaling_flag: 1036 params = dot(params, self.scaling_matrix) 1037 1038 # Initial chi-squared (or SSE) gradient. 1039 self.dchi2 = self.dchi2 * 0.0 1040 1041 # Loop over each alignment. 1042 index = 0 1043 for i in xrange(self.num_align): 1044 # Fixed tensor, so skip. 1045 if self.fixed_tensors[i]: 1046 continue 1047 1048 # Construct the Amn partial derivative components. 1049 for j in xrange(self.num_spins): 1050 # RDC. 1051 if self.fixed_tensors[i] and self.rdc_flag and not self.missing_Dij[i, j]: 1052 self.dDij_theta[i*5, i, j] = ave_rdc_tensor_dDij_dAmn(self.dip_const[j], self.dip_vect[j], self.N, self.dA[0], weights=self.probs) 1053 self.dDij_theta[i*5+1, i, j] = ave_rdc_tensor_dDij_dAmn(self.dip_const[j], self.dip_vect[j], self.N, self.dA[1], weights=self.probs) 1054 self.dDij_theta[i*5+2, i, j] = ave_rdc_tensor_dDij_dAmn(self.dip_const[j], self.dip_vect[j], self.N, self.dA[2], weights=self.probs) 1055 self.dDij_theta[i*5+3, i, j] = ave_rdc_tensor_dDij_dAmn(self.dip_const[j], self.dip_vect[j], self.N, self.dA[3], weights=self.probs) 1056 self.dDij_theta[i*5+4, i, j] = ave_rdc_tensor_dDij_dAmn(self.dip_const[j], self.dip_vect[j], self.N, self.dA[4], weights=self.probs) 1057 1058 # PCS. 1059 if self.fixed_tensors[i] and self.pcs_flag and not self.missing_deltaij[i, j]: 1060 self.ddeltaij_theta[i*5, i, j] = ave_pcs_tensor_ddeltaij_dAmn(self.pcs_const[i, j], self.paramag_unit_vect[j], self.N, self.dA[0], weights=self.probs) 1061 self.ddeltaij_theta[i*5+1, i, j] = ave_pcs_tensor_ddeltaij_dAmn(self.pcs_const[i, j], self.paramag_unit_vect[j], self.N, self.dA[1], weights=self.probs) 1062 self.ddeltaij_theta[i*5+2, i, j] = ave_pcs_tensor_ddeltaij_dAmn(self.pcs_const[i, j], self.paramag_unit_vect[j], self.N, self.dA[2], weights=self.probs) 1063 self.ddeltaij_theta[i*5+3, i, j] = ave_pcs_tensor_ddeltaij_dAmn(self.pcs_const[i, j], self.paramag_unit_vect[j], self.N, self.dA[3], weights=self.probs) 1064 self.ddeltaij_theta[i*5+4, i, j] = ave_pcs_tensor_ddeltaij_dAmn(self.pcs_const[i, j], self.paramag_unit_vect[j], self.N, self.dA[4], weights=self.probs) 1065 1066 # Construct the pc partial derivative gradient components, looping over each state. 1067 for c in xrange(self.N - 1): 1068 # Index in the parameter array. 1069 param_index = self.num_align_params + c 1070 1071 # Loop over the spins. 1072 for j in xrange(self.num_spins): 1073 # Calculate the RDC for state c (this is the pc partial derivative). 1074 if self.rdc_flag and not self.missing_Dij[i, j]: 1075 self.dDij_theta[param_index, i, j] = rdc_tensor(self.dip_const[j], self.dip_vect[j, c], self.A[i]) 1076 1077 # Calculate the PCS for state c (this is the pc partial derivative). 1078 if self.pcs_flag and not self.missing_deltaij[i, j]: 1079 self.ddeltaij_theta[param_index, i, j] = pcs_tensor(self.pcs_const[i, j, c], self.paramag_unit_vect[j, c], self.A[i]) 1080 1081 # Construct the chi-squared gradient element for parameter k, alignment i. 1082 for k in xrange(self.total_num_params): 1083 # RDC part of the chi-squared gradient. 1084 if self.rdc_flag: 1085 self.dchi2[k] = self.dchi2[k] + dchi2_element(self.Dij[i], self.Dij_theta[i], self.dDij_theta[k, i], self.rdc_sigma_ij[i]) 1086 1087 # PCS part of the chi-squared gradient. 1088 if self.pcs_flag: 1089 self.dchi2[k] = self.dchi2[k] + dchi2_element(self.deltaij[i], self.deltaij_theta[i], self.ddeltaij_theta[k, i], self.pcs_sigma_ij[i]) 1090 1091 # Increment the index. 1092 index += 1 1093 1094 # Diagonal scaling. 1095 if self.scaling_flag: 1096 self.dchi2 = dot(self.dchi2, self.scaling_matrix) 1097 1098 # The gradient. 1099 return self.dchi2

1100 1101

1102 - def dfunc_tensor_opt(self, params):

1103 """The gradient function for optimisation of the alignment tensor from RDC and/or PCS data. 1104 1105 Description 1106 =========== 1107 1108 This function should be passed to the optimisation algorithm. It accepts, as an array, a vector of parameter values and, using these, returns the chi-squared gradient corresponding to that coordinate in the parameter space. If no RDC errors are supplied, then the SSE (the sum of squares error) gradient is returned instead. The chi-squared gradient is simply the SSE gradient normalised to unit variance (the SSE divided by the error squared). 1109 1110 Indices 1111 ======= 1112 1113 For this calculation, six indices are looped over and used in the various data structures. These include: 1114 - k, the index over all parameters, 1115 - i, the index over alignments, 1116 - j, the index over spin systems, 1117 - c, the index over the N-states (or over the structures), 1118 - m, the index over the first dimension of the alignment tensor m = {x, y, z}. 1119 - n, the index over the second dimension of the alignment tensor n = {x, y, z}, 1120 1121 1122 Equations 1123 ========= 1124 1125 To calculate the chi-squared gradient, a chain of equations are used. This includes the chi-squared gradient, the RDC gradient and the alignment tensor gradient. 1126 1127 1128 The chi-squared gradient 1129 ------------------------ 1130 1131 The equation is:: 1132 ___ 1133 dchi^2(theta) \ / Dij - Dij(theta) dDij(theta) \ 1134 ------------- = -2 > | ---------------- . ----------- | 1135 dthetak /__ \ sigma_ij**2 dthetak / 1136 ij 1137 1138 where: 1139 - theta is the parameter vector, 1140 - Dij are the measured RDCs or PCSs, 1141 - Dij(theta) are the back calculated RDCs or PCSs, 1142 - sigma_ij are the RDC or PCS errors, 1143 - dDij(theta)/dthetak is the RDC or PCS gradient for parameter k. 1144 1145 1146 The RDC gradient 1147 ---------------- 1148 1149 The only parameters are the tensor components. 1150 1151 Amn partial derivative 1152 ~~~~~~~~~~~~~~~~~~~~~~ 1153 1154 The alignment tensor element partial derivative is:: 1155 1156 _N_ 1157 dDij(theta) dj \ T dAi 1158 ----------- = -- > mu_jc . ---- . mu_jc, 1159 dAmn N /__ dAmn 1160 c=1 1161 1162 where: 1163 - dj is the dipolar constant for spin j, 1164 - N is the total number of states or structures, 1165 - mu_jc is the unit vector corresponding to spin j and state c, 1166 - dAi/dAmn is the partial derivative of the alignment tensor with respect to element Amn. 1167 1168 1169 The PCS gradient 1170 ---------------- 1171 1172 Amn partial derivative 1173 ~~~~~~~~~~~~~~~~~~~~~~ 1174 1175 The alignment tensor element partial derivative is:: 1176 1177 _N_ 1178 ddelta_ij(theta) 1 \ T dAi 1179 ---------------- = - > djc . mu_jc . ---- . mu_jc, 1180 dAmn N /__ dAmn 1181 c=1 1182 1183 where: 1184 - djc is the pseudocontact shift constant for spin j and state c, 1185 - N is the total number of states or structures, 1186 - mu_jc is the unit vector corresponding to spin j and state c, 1187 - dAi/dAmn is the partial derivative of the alignment tensor with respect to element Amn. 1188 1189 1190 The alignment tensor gradient 1191 ----------------------------- 1192 1193 The five unique elements of the tensor {Axx, Ayy, Axy, Axz, Ayz} give five different partial derivatives. These are:: 1194 1195 dAi | 1 0 0 | 1196 ---- = | 0 0 0 |, 1197 dAxx | 0 0 -1 | 1198 1199 dAi | 0 0 0 | 1200 ---- = | 0 1 0 |, 1201 dAyy | 0 0 -1 | 1202 1203 dAi | 0 1 0 | 1204 ---- = | 1 0 0 |, 1205 dAxy | 0 0 0 | 1206 1207 dAi | 0 0 1 | 1208 ---- = | 0 0 0 |, 1209 dAxz | 1 0 0 | 1210 1211 dAi | 0 0 0 | 1212 ---- = | 0 0 1 |. 1213 dAyz | 0 1 0 | 1214 1215 As these are invariant, they can be pre-calculated. 1216 1217 1218 Stored data structures 1219 ====================== 1220 1221 There are a number of data structures calculated by this function and stored for subsequent use in the Hessian function. This include the back calculated RDC and PCS gradients and the alignment tensor gradients. 1222 1223 dDij(theta)/dthetak 1224 ------------------- 1225 1226 The back calculated RDC gradient. This is a rank-3 tensor with indices {k, i, j}. 1227 1228 ddeltaij(theta)/dthetak 1229 ----------------------- 1230 1231 The back calculated PCS gradient. This is a rank-3 tensor with indices {k, i, j}. 1232 1233 dAi/dAmn 1234 -------- 1235 1236 The alignment tensor gradients. This is a rank-3 tensor with indices {5, n, m}. 1237 1238 1239 @param params: The vector of parameter values. This is unused as it is assumed that func_population() was called first. 1240 @type params: numpy rank-1 array 1241 @return: The chi-squared or SSE gradient. 1242 @rtype: numpy rank-1 array 1243 """ 1244 1245 # Scaling. 1246 if self.scaling_flag: 1247 params = dot(params, self.scaling_matrix) 1248 1249 # Initial chi-squared (or SSE) gradient. 1250 self.dchi2 = self.dchi2 * 0.0 1251 1252 # Loop over each alignment. 1253 index = 0 1254 for i in xrange(self.num_align): 1255 # Fixed tensor, so skip. 1256 if self.fixed_tensors[i]: 1257 continue 1258 1259 # Construct the Amn partial derivative components. 1260 for j in xrange(self.num_spins): 1261 # RDC. 1262 if self.rdc_flag and not self.missing_Dij[i, j]: 1263 self.dDij_theta[index*5, index, j] = ave_rdc_tensor_dDij_dAmn(self.dip_const[j], self.dip_vect[j], self.N, self.dA[0], weights=self.probs) 1264 self.dDij_theta[index*5+1, index, j] = ave_rdc_tensor_dDij_dAmn(self.dip_const[j], self.dip_vect[j], self.N, self.dA[1], weights=self.probs) 1265 self.dDij_theta[index*5+2, index, j] = ave_rdc_tensor_dDij_dAmn(self.dip_const[j], self.dip_vect[j], self.N, self.dA[2], weights=self.probs) 1266 self.dDij_theta[index*5+3, index, j] = ave_rdc_tensor_dDij_dAmn(self.dip_const[j], self.dip_vect[j], self.N, self.dA[3], weights=self.probs) 1267 self.dDij_theta[index*5+4, index, j] = ave_rdc_tensor_dDij_dAmn(self.dip_const[j], self.dip_vect[j], self.N, self.dA[4], weights=self.probs) 1268 1269 # PCS. 1270 if self.pcs_flag and not self.missing_deltaij[i, j]: 1271 self.ddeltaij_theta[index*5, index, j] = ave_pcs_tensor_ddeltaij_dAmn(self.pcs_const[index, j], self.paramag_unit_vect[j], self.N, self.dA[0], weights=self.probs) 1272 self.ddeltaij_theta[index*5+1, index, j] = ave_pcs_tensor_ddeltaij_dAmn(self.pcs_const[index, j], self.paramag_unit_vect[j], self.N, self.dA[1], weights=self.probs) 1273 self.ddeltaij_theta[index*5+2, index, j] = ave_pcs_tensor_ddeltaij_dAmn(self.pcs_const[index, j], self.paramag_unit_vect[j], self.N, self.dA[2], weights=self.probs) 1274 self.ddeltaij_theta[index*5+3, index, j] = ave_pcs_tensor_ddeltaij_dAmn(self.pcs_const[index, j], self.paramag_unit_vect[j], self.N, self.dA[3], weights=self.probs) 1275 self.ddeltaij_theta[index*5+4, index, j] = ave_pcs_tensor_ddeltaij_dAmn(self.pcs_const[index, j], self.paramag_unit_vect[j], self.N, self.dA[4], weights=self.probs) 1276 1277 # Construct the chi-squared gradient element for parameter k, alignment i. 1278 for k in xrange(self.total_num_params): 1279 # RDC part of the chi-squared gradient. 1280 if self.rdc_flag: 1281 self.dchi2[k] = self.dchi2[k] + dchi2_element(self.Dij[index], self.Dij_theta[index], self.dDij_theta[k, index], self.rdc_sigma_ij[index]) 1282 1283 # PCS part of the chi-squared gradient. 1284 if self.pcs_flag: 1285 self.dchi2[k] = self.dchi2[k] + dchi2_element(self.deltaij[index], self.deltaij_theta[index], self.ddeltaij_theta[k, index], self.pcs_sigma_ij[index]) 1286 1287 # Increment the index. 1288 index += 1 1289 1290 # Diagonal scaling. 1291 if self.scaling_flag: 1292 self.dchi2 = dot(self.dchi2, self.scaling_matrix) 1293 1294 # The gradient. 1295 return self.dchi2

1296 1297

1298 - def d2func_population(self, params):

1299 """The Hessian function for optimisation of the flexible population N-state model. 1300 1301 Description 1302 =========== 1303 1304 This function should be passed to the optimisation algorithm. It accepts, as an array, a vector of parameter values and, using these, returns the chi-squared Hessian corresponding to that coordinate in the parameter space. If no RDC/PCS errors are supplied, then the SSE (the sum of squares error) Hessian is returned instead. The chi-squared Hessian is simply the SSE Hessian normalised to unit variance (the SSE divided by the error squared). 1305 1306 1307 Indices 1308 ======= 1309 1310 For this calculation, six indices are looped over and used in the various data structures. These include: 1311 - k, the index over all parameters, 1312 - i, the index over alignments, 1313 - j, the index over spin systems, 1314 - c, the index over the N-states (or over the structures), 1315 - m, the index over the first dimension of the alignment tensor m = {x, y, z}. 1316 - n, the index over the second dimension of the alignment tensor n = {x, y, z}, 1317 1318 1319 Equations 1320 ========= 1321 1322 To calculate the chi-squared gradient, a chain of equations are used. This includes the chi-squared gradient, the RDC gradient and the alignment tensor gradient. 1323 1324 1325 The chi-squared Hessian 1326 ----------------------- 1327 1328 The equation is:: 1329 ___ 1330 d2chi^2(theta) \ 1 / dDij(theta) dDij(theta) d2Dij(theta) \ 1331 --------------- = 2 > ---------- | ----------- . ----------- - (Dij-Dij(theta)) . --------------- |. 1332 dthetaj.dthetak /__ sigma_i**2 \ dthetaj dthetak dthetaj.dthetak / 1333 ij 1334 1335 where: 1336 - theta is the parameter vector, 1337 - Dij are the measured RDCs or PCSs, 1338 - Dij(theta) are the back calculated RDCs or PCSs, 1339 - sigma_ij are the RDC or PCS errors, 1340 - dDij(theta)/dthetak is the RDC or PCS gradient for parameter k. 1341 - d2Dij(theta)/dthetaj.dthetak is the RDC or PCS Hessian for parameters j and k. 1342 1343 1344 The RDC Hessian 1345 --------------- 1346 1347 pc-pd second partial derivatives 1348 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 1349 1350 The probability parameter second partial derivative is:: 1351 1352 d2Dij(theta) 1353 ------------ = 0. 1354 dpc.dpd 1355 1356 1357 pc-Anm second partial derivatives 1358 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 1359 1360 The probability parameter-tensor element second partial derivative is:: 1361 1362 d2Dij(theta) T dAi 1363 ------------ = dj . mu_jc . ---- . mu_jc. 1364 dpc.dAmn dAmn 1365 1366 1367 Amn-Aop second partial derivatives 1368 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 1369 1370 The alignment tensor element second partial derivative is:: 1371 1372 d2Dij(theta) 1373 ------------ = 0. 1374 dAmn.dAop 1375 1376 1377 The PCS Hessian 1378 --------------- 1379 1380 pc-pd second partial derivatives 1381 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 1382 1383 The probability parameter second partial derivative is:: 1384 1385 d2delta_ij(theta) 1386 ----------------- = 0. 1387 dpc.dpd 1388 1389 1390 pc-Anm second partial derivatives 1391 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 1392 1393 The probability parameter-tensor element second partial derivative is:: 1394 1395 d2delta_ij(theta) T dAi 1396 ----------------- = djc . mu_jc . ---- . mu_jc. 1397 dpc.dAmn dAmn 1398 1399 1400 Amn-Aop second partial derivatives 1401 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 1402 1403 The alignment tensor element second partial derivative is:: 1404 1405 d2delta_ij(theta) 1406 ----------------- = 0 1407 dAmn.dAop 1408 1409 1410 The alignment tensor Hessian 1411 ---------------------------- 1412 1413 The five unique elements of the tensor {Axx, Ayy, Axy, Axz, Ayz} all have the same second partial derivative of:: 1414 1415 d2Ai | 0 0 0 | 1416 --------- = | 0 0 0 |. 1417 dAmn.dAop | 0 0 0 | 1418 1419 1420 @param params: The vector of parameter values. This is unused as it is assumed that func_population() was called first. 1421 @type params: numpy rank-1 array 1422 @return: The chi-squared or SSE Hessian. 1423 @rtype: numpy rank-2 array 1424 """ 1425 1426 # Scaling. 1427 if self.scaling_flag: 1428 params = dot(params, self.scaling_matrix) 1429 1430 # Initial chi-squared (or SSE) Hessian. 1431 self.d2chi2 = self.d2chi2 * 0.0 1432 1433 # Loop over each alignment. 1434 index = 0 1435 for i in xrange(self.num_align): 1436 # Fixed tensor, so skip. 1437 if self.fixed_tensors[i]: 1438 continue 1439 1440 # Construct the pc-Amn second partial derivative Hessian components. 1441 for c in xrange(self.N - 1): 1442 # Index in the parameter array. 1443 pc_index = self.num_align_params + c 1444 1445 # Loop over the spins. 1446 for j in xrange(self.num_spins): 1447 # Calculate the RDC Hessian component. 1448 if self.fixed_tensors[i] and self.rdc_flag and not self.missing_Dij[i, j]: 1449 self.d2Dij_theta[pc_index, i*5+0, i, j] = self.d2Dij_theta[i*5+0, pc_index, i, j] = rdc_tensor(self.dip_const[j], self.dip_vect[j, c], self.dA[0]) 1450 self.d2Dij_theta[pc_index, i*5+1, i, j] = self.d2Dij_theta[i*5+1, pc_index, i, j] = rdc_tensor(self.dip_const[j], self.dip_vect[j, c], self.dA[1]) 1451 self.d2Dij_theta[pc_index, i*5+2, i, j] = self.d2Dij_theta[i*5+2, pc_index, i, j] = rdc_tensor(self.dip_const[j], self.dip_vect[j, c], self.dA[2]) 1452 self.d2Dij_theta[pc_index, i*5+3, i, j] = self.d2Dij_theta[i*5+3, pc_index, i, j] = rdc_tensor(self.dip_const[j], self.dip_vect[j, c], self.dA[3]) 1453 self.d2Dij_theta[pc_index, i*5+4, i, j] = self.d2Dij_theta[i*5+4, pc_index, i, j] = rdc_tensor(self.dip_const[j], self.dip_vect[j, c], self.dA[4]) 1454 1455 # Calculate the PCS Hessian component. 1456 if self.fixed_tensors[i] and self.pcs_flag and not self.missing_deltaij[i, j]: 1457 self.d2deltaij_theta[pc_index, i*5+0, i, j] = self.d2deltaij_theta[i*5+0, pc_index, i, j] = pcs_tensor(self.pcs_const[i, j, c], self.paramag_unit_vect[j, c], self.dA[0]) 1458 self.d2deltaij_theta[pc_index, i*5+1, i, j] = self.d2deltaij_theta[i*5+1, pc_index, i, j] = pcs_tensor(self.pcs_const[i, j, c], self.paramag_unit_vect[j, c], self.dA[1]) 1459 self.d2deltaij_theta[pc_index, i*5+2, i, j] = self.d2deltaij_theta[i*5+2, pc_index, i, j] = pcs_tensor(self.pcs_const[i, j, c], self.paramag_unit_vect[j, c], self.dA[2]) 1460 self.d2deltaij_theta[pc_index, i*5+3, i, j] = self.d2deltaij_theta[i*5+3, pc_index, i, j] = pcs_tensor(self.pcs_const[i, j, c], self.paramag_unit_vect[j, c], self.dA[3]) 1461 self.d2deltaij_theta[pc_index, i*5+4, i, j] = self.d2deltaij_theta[i*5+4, pc_index, i, j] = pcs_tensor(self.pcs_const[i, j, c], self.paramag_unit_vect[j, c], self.dA[4]) 1462 1463 # Increment the index. 1464 index += 1 1465 1466 # Loop over each alignment. 1467 for i in xrange(self.num_align): 1468 # Construct the chi-squared Hessian element for parameters j and k, alignment i. 1469 for j in xrange(self.total_num_params): 1470 for k in xrange(self.total_num_params): 1471 # RDC part of the chi-squared gradient. 1472 if self.fixed_tensors[i] and self.rdc_flag: 1473 self.d2chi2[j, k] = self.d2chi2[j, k] + d2chi2_element(self.Dij[i], self.Dij_theta[i], self.dDij_theta[j, i], self.dDij_theta[k, i], self.d2Dij_theta[j, k, i], self.rdc_sigma_ij[i]) 1474 1475 # PCS part of the chi-squared gradient. 1476 if self.fixed_tensors[i] and self.pcs_flag: 1477 self.d2chi2[j, k] = self.d2chi2[j, k] + d2chi2_element(self.deltaij[i], self.deltaij_theta[i], self.ddeltaij_theta[j, i], self.ddeltaij_theta[k, i], self.d2deltaij_theta[j, k, i], self.pcs_sigma_ij[i]) 1478 1479 # Diagonal scaling. 1480 if self.scaling_flag: 1481 self.d2chi2 = dot(self.d2chi2, self.scaling_matrix) 1482 1483 # The gradient. 1484 return self.d2chi2

1485 1486

1487 - def d2func_tensor_opt(self, params):

1488 """The Hessian function for optimisation of the alignment tensor from RDC and/or PCS data. 1489 1490 Description 1491 =========== 1492 1493 This function should be passed to the optimisation algorithm. It accepts, as an array, a vector of parameter values and, using these, returns the chi-squared Hessian corresponding to that coordinate in the parameter space. If no RDC/PCS errors are supplied, then the SSE (the sum of squares error) Hessian is returned instead. The chi-squared Hessian is simply the SSE Hessian normalised to unit variance (the SSE divided by the error squared). 1494 1495 Indices 1496 ======= 1497 1498 For this calculation, six indices are looped over and used in the various data structures. These include: 1499 - k, the index over all parameters, 1500 - i, the index over alignments, 1501 - j, the index over spin systems, 1502 - c, the index over the N-states (or over the structures), 1503 - m, the index over the first dimension of the alignment tensor m = {x, y, z}. 1504 - n, the index over the second dimension of the alignment tensor n = {x, y, z}, 1505 1506 1507 Equations 1508 ========= 1509 1510 To calculate the chi-squared gradient, a chain of equations are used. This includes the chi-squared gradient, the RDC gradient and the alignment tensor gradient. 1511 1512 1513 The chi-squared Hessian 1514 ----------------------- 1515 1516 The equation is:: 1517 ___ 1518 d2chi^2(theta) \ 1 / dDij(theta) dDij(theta) d2Dij(theta) \ 1519 --------------- = 2 > ---------- | ----------- . ----------- - (Dij-Dij(theta)) . --------------- |. 1520 dthetaj.dthetak /__ sigma_i**2 \ dthetaj dthetak dthetaj.dthetak / 1521 ij 1522 1523 where: 1524 - theta is the parameter vector, 1525 - Dij are the measured RDCs or PCSs, 1526 - Dij(theta) are the back calculated RDCs or PCSs, 1527 - sigma_ij are the RDC or PCS errors, 1528 - dDij(theta)/dthetak is the RDC or PCS gradient for parameter k. 1529 - d2Dij(theta)/dthetaj.dthetak is the RDC or PCS Hessian for parameters j and k. 1530 1531 1532 The RDC Hessian 1533 --------------- 1534 1535 The only parameters are the tensor components. 1536 1537 1538 Amn-Aop second partial derivatives 1539 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 1540 1541 The alignment tensor element second partial derivative is:: 1542 1543 d2Dij(theta) 1544 ------------ = 0. 1545 dAmn.dAop 1546 1547 1548 The PCS Hessian 1549 --------------- 1550 1551 Amn-Aop second partial derivatives 1552 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 1553 1554 The alignment tensor element second partial derivative is:: 1555 1556 d2delta_ij(theta) 1557 ----------------- = 0 1558 dAmn.dAop 1559 1560 1561 The alignment tensor Hessian 1562 ---------------------------- 1563 1564 The five unique elements of the tensor {Axx, Ayy, Axy, Axz, Ayz} all have the same second partial derivative of:: 1565 1566 d2Ai | 0 0 0 | 1567 --------- = | 0 0 0 |. 1568 dAmn.dAop | 0 0 0 | 1569 1570 1571 @param params: The vector of parameter values. This is unused as it is assumed that func_population() was called first. 1572 @type params: numpy rank-1 array 1573 @return: The chi-squared or SSE Hessian. 1574 @rtype: numpy rank-2 array 1575 """ 1576 1577 # Scaling. 1578 if self.scaling_flag: 1579 params = dot(params, self.scaling_matrix) 1580 1581 # Initial chi-squared (or SSE) Hessian. 1582 self.d2chi2 = self.d2chi2 * 0.0 1583 1584 # Loop over each alignment. 1585 index = 0 1586 for i in xrange(self.num_align): 1587 # Fixed tensor, so skip. 1588 if self.fixed_tensors[i]: 1589 continue 1590 1591 # Construct the chi-squared Hessian element for parameters j and k, alignment i. 1592 for j in xrange(self.total_num_params): 1593 for k in xrange(self.total_num_params): 1594 # RDC part of the chi-squared gradient. 1595 if self.rdc_flag: 1596 self.d2chi2[j, k] = self.d2chi2[j, k] + d2chi2_element(self.Dij[i], self.Dij_theta[i], self.dDij_theta[j, i], self.dDij_theta[k, i], self.zero_hessian, self.rdc_sigma_ij[i]) 1597 1598 # PCS part of the chi-squared gradient. 1599 if self.pcs_flag: 1600 self.d2chi2[j, k] = self.d2chi2[j, k] + d2chi2_element(self.deltaij[i], self.deltaij_theta[i], self.ddeltaij_theta[j, i], self.ddeltaij_theta[k, i], self.zero_hessian, self.pcs_sigma_ij[i]) 1601 1602 # Increment the index. 1603 index += 1 1604 1605 # Diagonal scaling. 1606 if self.scaling_flag: 1607 self.d2chi2 = dot(self.d2chi2, self.scaling_matrix) 1608 1609 # The gradient. 1610 return self.d2chi2

1611 1612

1613 - def paramag_info(self):

1614 """Calculate the paramagnetic centre to spin vectors, distances and constants.""" 1615 1616 # Get the vectors and distances. 1617 if rank(self.paramag_centre) == 1: 1618 vectors_single_centre(self.atomic_pos, self.paramag_centre, self.paramag_unit_vect, self.paramag_dist) 1619 else: 1620 vectors_centre_per_state(self.atomic_pos, self.paramag_centre, self.paramag_unit_vect, self.paramag_dist) 1621 1622 # The PCS constants. 1623 for i in range(self.num_align): 1624 for j in range(self.num_spins): 1625 for c in range(self.N): 1626 self.pcs_const[i, j, c] = pcs_constant(self.temp[i], self.frq[i], self.paramag_dist[j, c])

Source Code for Module maths_fns.n_state_model