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