Package specific_analyses :: Package n_state_model :: Module optimisation
[hide private]
[frames] | no frames]

Source Code for Module specific_analyses.n_state_model.optimisation

  1  ############################################################################### 
  2  #                                                                             # 
  3  # Copyright (C) 2007-2015 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 optimisation functions.""" 
 24   
 25  # Python module imports. 
 26  from numpy import array, dot, float64, ones, zeros 
 27  from numpy.linalg import inv 
 28  from numpy.ma import fix_invalid 
 29   
 30  # relax module imports. 
 31  from lib.errors import RelaxError, RelaxNoModelError 
 32  from pipe_control import align_tensor 
 33  from pipe_control.align_tensor import opt_uses_align_data, opt_uses_tensor 
 34  from pipe_control.interatomic import interatomic_loop 
 35  from pipe_control.mol_res_spin import return_spin, spin_loop 
 36  from pipe_control.pcs import return_pcs_data 
 37  from pipe_control.rdc import check_rdcs, return_rdc_data 
 38  from specific_analyses.n_state_model.data import base_data_types, tensor_loop 
 39  from specific_analyses.n_state_model.parameters import assemble_param_vector, update_model 
 40  from target_functions.n_state_model import N_state_opt 
 41   
 42   
43 -def minimise_bc_data(model, sim_index=None):
44 """Extract and unpack the back calculated data. 45 46 @param model: The instantiated class containing the target function. 47 @type model: class instance 48 @keyword sim_index: The optional Monte Carlo simulation index. 49 @type sim_index: None or int 50 """ 51 52 # No alignment tensors, so nothing to do. 53 if not hasattr(cdp, 'align_tensors'): 54 return 55 56 # Simulation flag. 57 sim_flag = (sim_index != None) 58 59 # Loop over each alignment. 60 align_index = 0 61 for i in range(len(cdp.align_ids)): 62 # Skip non-optimised tensors. 63 if not opt_uses_tensor(cdp.align_tensors[i]): 64 continue 65 66 # The alignment ID. 67 align_id = cdp.align_ids[i] 68 69 # Data flags 70 rdc_flag = False 71 if hasattr(cdp, 'rdc_ids') and align_id in cdp.rdc_ids: 72 rdc_flag = True 73 pcs_flag = False 74 if hasattr(cdp, 'pcs_ids') and align_id in cdp.pcs_ids: 75 pcs_flag = True 76 77 # Spin loop. 78 pcs_index = 0 79 for spin in spin_loop(): 80 # Skip deselected spins. 81 if not spin.select: 82 continue 83 84 # Spins with PCS data. 85 if pcs_flag and hasattr(spin, 'pcs'): 86 # Initialise the data structure if necessary. 87 if sim_flag: 88 if not hasattr(spin, 'pcs_sim_bc'): 89 spin.pcs_sim_bc = {} 90 if align_id not in spin.pcs_sim_bc: 91 spin.pcs_sim_bc[align_id] = [None] * cdp.sim_number 92 else: 93 if not hasattr(spin, 'pcs_bc'): 94 spin.pcs_bc = {} 95 96 # Add the back calculated PCS (in ppm). 97 if sim_flag: 98 spin.pcs_sim_bc[align_id][sim_index] = model.deltaij_theta[align_index, pcs_index] * 1e6 99 else: 100 spin.pcs_bc[align_id] = model.deltaij_theta[align_index, pcs_index] * 1e6 101 102 # Increment the data index if the spin container has data. 103 pcs_index = pcs_index + 1 104 105 # Interatomic data container loop. 106 rdc_index = 0 107 for interatom in interatomic_loop(): 108 # Get the spins. 109 spin1 = return_spin(spin_hash=interatom._spin_hash1) 110 spin2 = return_spin(spin_hash=interatom._spin_hash2) 111 112 # RDC checks. 113 if not check_rdcs(interatom): 114 continue 115 116 # Containers with RDC data. 117 if rdc_flag and hasattr(interatom, 'rdc'): 118 # Initialise the data structure if necessary. 119 if sim_flag: 120 if not hasattr(interatom, 'rdc_sim_bc'): 121 interatom.rdc_sim_bc = {} 122 if align_id not in interatom.rdc_sim_bc: 123 interatom.rdc_sim_bc[align_id] = [0.0] * cdp.sim_number 124 else: 125 if not hasattr(interatom, 'rdc_bc'): 126 interatom.rdc_bc = {} 127 128 # Append the back calculated PCS. 129 if sim_flag: 130 interatom.rdc_sim_bc[align_id][sim_index] = model.rdc_theta[align_index, rdc_index] 131 else: 132 interatom.rdc_bc[align_id] = model.rdc_theta[align_index, rdc_index] 133 134 # Increment the data index if the interatom container has data. 135 rdc_index = rdc_index + 1 136 137 # Increment the alignment index (for the optimised tensors). 138 align_index += 1
139 140
141 -def minimise_setup_atomic_pos(sim_index=None):
142 """Set up the atomic position data structures for optimisation using PCSs and PREs as base data sets. 143 144 @keyword sim_index: The index of the simulation to optimise. This should be None if normal optimisation is desired. 145 @type sim_index: None or int 146 @return: The atomic positions (the first index is the spins, the second is the structures, and the third is the atomic coordinates) and the paramagnetic centre. 147 @rtype: numpy rank-3 array, numpy rank-1 array. 148 """ 149 150 # Initialise. 151 atomic_pos = [] 152 153 # Store the atomic positions. 154 for spin in spin_loop(): 155 # Skip deselected spins. 156 if not spin.select: 157 continue 158 159 # Only use spins with alignment/paramagnetic data. 160 if not hasattr(spin, 'pcs') and not hasattr(spin, 'pre'): 161 continue 162 163 # The position list. 164 if type(spin.pos[0]) in [float, float64]: 165 atomic_pos.append([spin.pos]) 166 else: 167 atomic_pos.append(spin.pos) 168 169 # Convert to numpy objects. 170 atomic_pos = array(atomic_pos, float64) 171 172 # The paramagnetic centre. 173 if not hasattr(cdp, 'paramagnetic_centre'): 174 paramag_centre = zeros(3, float64) 175 elif sim_index != None and not cdp.paramag_centre_fixed: 176 if not hasattr(cdp, 'paramagnetic_centre_sim') or cdp.paramagnetic_centre_sim[sim_index] is None: 177 paramag_centre = zeros(3, float64) 178 else: 179 paramag_centre = array(cdp.paramagnetic_centre_sim[sim_index]) 180 else: 181 paramag_centre = array(cdp.paramagnetic_centre) 182 183 # Return the data structures. 184 return atomic_pos, paramag_centre
185 186
187 -def minimise_setup_tensors(sim_index=None):
188 """Set up the data structures for optimisation using alignment tensors as base data sets. 189 190 @keyword sim_index: The index of the simulation to optimise. This should be None if 191 normal optimisation is desired. 192 @type sim_index: None or int 193 @return: The assembled data structures for using alignment tensors as the base 194 data for optimisation. These include: 195 - full_tensors, the data of the full alignment tensors. 196 - red_tensor_elem, the tensors as concatenated rank-1 5D arrays. 197 - red_tensor_err, the tensor errors as concatenated rank-1 5D 198 arrays. 199 - full_in_ref_frame, flags specifying if the tensor in the reference 200 frame is the full or reduced tensor. 201 @rtype: tuple of (list, numpy rank-1 array, numpy rank-1 array, numpy rank-1 202 array) 203 """ 204 205 # Initialise. 206 n = len(cdp.align_tensors.reduction) 207 full_tensors = zeros(n*5, float64) 208 red_tensors = zeros(n*5, float64) 209 red_err = ones(n*5, float64) * 1e-5 210 full_in_ref_frame = zeros(n, float64) 211 212 # Loop over the full tensors. 213 for i, tensor in tensor_loop(red=False): 214 # The full tensor. 215 full_tensors[5*i + 0] = tensor.Axx 216 full_tensors[5*i + 1] = tensor.Ayy 217 full_tensors[5*i + 2] = tensor.Axy 218 full_tensors[5*i + 3] = tensor.Axz 219 full_tensors[5*i + 4] = tensor.Ayz 220 221 # The full tensor corresponds to the frame of reference. 222 if cdp.ref_domain == tensor.domain: 223 full_in_ref_frame[i] = 1 224 225 # Loop over the reduced tensors. 226 for i, tensor in tensor_loop(red=True): 227 # The reduced tensor (simulation data). 228 if sim_index != None: 229 red_tensors[5*i + 0] = tensor.Axx_sim[sim_index] 230 red_tensors[5*i + 1] = tensor.Ayy_sim[sim_index] 231 red_tensors[5*i + 2] = tensor.Axy_sim[sim_index] 232 red_tensors[5*i + 3] = tensor.Axz_sim[sim_index] 233 red_tensors[5*i + 4] = tensor.Ayz_sim[sim_index] 234 235 # The reduced tensor. 236 else: 237 red_tensors[5*i + 0] = tensor.Axx 238 red_tensors[5*i + 1] = tensor.Ayy 239 red_tensors[5*i + 2] = tensor.Axy 240 red_tensors[5*i + 3] = tensor.Axz 241 red_tensors[5*i + 4] = tensor.Ayz 242 243 # The reduced tensor errors. 244 if hasattr(tensor, 'Axx_err'): 245 red_err[5*i + 0] = tensor.Axx_err 246 red_err[5*i + 1] = tensor.Ayy_err 247 red_err[5*i + 2] = tensor.Axy_err 248 red_err[5*i + 3] = tensor.Axz_err 249 red_err[5*i + 4] = tensor.Ayz_err 250 251 # Return the data structures. 252 return full_tensors, red_tensors, red_err, full_in_ref_frame
253 254
255 -def minimise_setup_fixed_tensors():
256 """Set up the data structures for the fixed alignment tensors. 257 258 @return: The assembled data structures for the fixed alignment tensors. 259 @rtype: numpy rank-1 array. 260 """ 261 262 # Initialise. 263 n = align_tensor.num_tensors(skip_fixed=False) - align_tensor.num_tensors(skip_fixed=True) 264 tensors = zeros(n*5, float64) 265 266 # Nothing to do. 267 if n == 0: 268 return None 269 270 # Loop over the tensors. 271 index = 0 272 for i in range(len(cdp.align_tensors)): 273 # Skip non-optimised data. 274 if not opt_uses_align_data(cdp.align_tensors[i].name): 275 continue 276 277 # No parameters have been set. 278 if not hasattr(cdp.align_tensors[i], 'Axx'): 279 continue 280 281 # The real tensors. 282 tensors[5*index + 0] = cdp.align_tensors[i].Axx 283 tensors[5*index + 1] = cdp.align_tensors[i].Ayy 284 tensors[5*index + 2] = cdp.align_tensors[i].Axy 285 tensors[5*index + 3] = cdp.align_tensors[i].Axz 286 tensors[5*index + 4] = cdp.align_tensors[i].Ayz 287 288 # Increment the index. 289 index += 1 290 291 # Return the data structure. 292 return tensors
293 294
295 -def target_fn_setup(sim_index=None, scaling_matrix=None, verbosity=0):
296 """Initialise the target function for optimisation or direct calculation. 297 298 @keyword sim_index: The index of the simulation to optimise. This should be None if normal optimisation is desired. 299 @type sim_index: None or int 300 @keyword scaling_matrix: The diagonal and square scaling matrix. 301 @type scaling_matrix: numpy rank-2, float64 array or None 302 @keyword verbosity: A flag specifying the amount of information to print. The higher the value, the greater the verbosity. 303 @type verbosity: int 304 """ 305 306 # Test if the N-state model has been set up. 307 if not hasattr(cdp, 'model'): 308 raise RelaxNoModelError('N-state') 309 310 # '2-domain' model setup tests. 311 if cdp.model == '2-domain': 312 # The number of states. 313 if not hasattr(cdp, 'N'): 314 raise RelaxError("The number of states has not been set.") 315 316 # The reference domain. 317 if not hasattr(cdp, 'ref_domain'): 318 raise RelaxError("The reference domain has not been set.") 319 320 # Update the model parameters if necessary. 321 update_model() 322 323 # Create the initial parameter vector. 324 param_vector = assemble_param_vector(sim_index=sim_index) 325 326 # Replace all NaNs with 0.0. 327 fix_invalid(param_vector, copy=False, fill_value=0.0) 328 329 # Determine if alignment tensors or RDCs are to be used. 330 data_types = base_data_types() 331 332 # The probabilities. 333 probs = None 334 if hasattr(cdp, 'probs') and len(cdp.probs) and cdp.probs[0] != None: 335 probs = cdp.probs 336 337 # Diagonal scaling. 338 if len(param_vector) and scaling_matrix is not None: 339 param_vector = dot(inv(scaling_matrix), param_vector) 340 341 # Get the data structures for optimisation using the tensors as base data sets. 342 full_tensors, red_tensor_elem, red_tensor_err, full_in_ref_frame = None, None, None, None 343 if 'tensor' in data_types: 344 full_tensors, red_tensor_elem, red_tensor_err, full_in_ref_frame = minimise_setup_tensors(sim_index=sim_index) 345 346 # Get the data structures for optimisation using PCSs as base data sets. 347 pcs, pcs_err, pcs_weight, temp, frq, pcs_pseudo_flags = None, None, None, None, None, None 348 if 'pcs' in data_types: 349 pcs, pcs_err, pcs_weight, temp, frq, pcs_pseudo_flags = return_pcs_data(sim_index=sim_index, verbosity=verbosity) 350 351 # Get the data structures for optimisation using RDCs as base data sets. 352 rdcs, rdc_err, rdc_weight, rdc_vector, rdc_dj, absolute_rdc, T_flags, j_couplings, rdc_pseudo_flags = None, None, None, None, None, None, None, None, None 353 if 'rdc' in data_types: 354 # The data. 355 rdcs, rdc_err, rdc_weight, rdc_vector, rdc_dj, absolute_rdc, T_flags, j_couplings, rdc_pseudo_flags = return_rdc_data(sim_index=sim_index, verbosity=verbosity) 356 357 # Get the fixed tensors. 358 fixed_tensors = None 359 if 'rdc' in data_types or 'pcs' in data_types: 360 full_tensors = minimise_setup_fixed_tensors() 361 362 # The flag list. 363 fixed_tensors = [] 364 for i in range(len(cdp.align_tensors)): 365 # Skip non-optimised data. 366 if not opt_uses_align_data(cdp.align_tensors[i].name): 367 continue 368 369 if cdp.align_tensors[i].fixed: 370 fixed_tensors.append(True) 371 else: 372 fixed_tensors.append(False) 373 374 # Get the atomic_positions. 375 atomic_pos, paramag_centre, centre_fixed = None, None, True 376 if 'pcs' in data_types or 'pre' in data_types: 377 atomic_pos, paramag_centre = minimise_setup_atomic_pos(sim_index=sim_index) 378 379 # Optimisation of the centre. 380 if hasattr(cdp, 'paramag_centre_fixed'): 381 centre_fixed = cdp.paramag_centre_fixed 382 383 # Set up the class instance containing the target function. 384 model = N_state_opt(model=cdp.model, N=cdp.N, init_params=param_vector, probs=probs, full_tensors=full_tensors, red_data=red_tensor_elem, red_errors=red_tensor_err, full_in_ref_frame=full_in_ref_frame, fixed_tensors=fixed_tensors, pcs=pcs, rdcs=rdcs, pcs_errors=pcs_err, rdc_errors=rdc_err, T_flags=T_flags, j_couplings=j_couplings, rdc_pseudo_flags=rdc_pseudo_flags, pcs_pseudo_flags=pcs_pseudo_flags, pcs_weights=pcs_weight, rdc_weights=rdc_weight, rdc_vect=rdc_vector, temp=temp, frq=frq, dip_const=rdc_dj, absolute_rdc=absolute_rdc, atomic_pos=atomic_pos, paramag_centre=paramag_centre, scaling_matrix=scaling_matrix, centre_fixed=centre_fixed) 385 386 # Return the data. 387 return model, param_vector, data_types
388