specific_analyses.relax

50 """Class containing functions for relaxation curve fitting.""" 51 52 # Class variable for storing the class instance (for the singleton design pattern). 53 instance = None 54

55 - def __init__(self):

56 """Initialise the class by placing API_common methods into the API.""" 57 58 # Place methods into the API. 59 self.base_data_loop = self._base_data_loop_spin 60 self.model_loop = self._model_loop_spin 61 self.return_conversion_factor = self._return_no_conversion_factor 62 self.return_value = self._return_value_general 63 self.set_error = self._set_error_spin 64 self.set_param_values = self._set_param_values_spin 65 self.set_selected_sim = self._set_selected_sim_spin 66 self.sim_init_values = self._sim_init_values_spin 67 self.sim_return_param = self._sim_return_param_spin 68 self.sim_return_selected = self._sim_return_selected_spin 69 70 # Place a copy of the parameter list object in the instance namespace. 71 self._PARAMS = Relax_fit_params()

72 73

74 - def create_mc_data(self, data_id=None):

75 """Create the Monte Carlo peak intensity data. 76 77 @keyword data_id: The spin identification string, as yielded by the base_data_loop() generator method. 78 @type data_id: str 79 @return: The Monte Carlo simulation data. 80 @rtype: list of floats 81 """ 82 83 # Initialise the MC data data structure. 84 mc_data = {} 85 86 # Get the spin container. 87 spin = return_spin(data_id) 88 89 # Skip deselected spins. 90 if not spin.select: 91 return 92 93 # Skip spins which have no data. 94 if not hasattr(spin, 'peak_intensity'): 95 return 96 97 # Test if the model is set. 98 if not hasattr(spin, 'model') or not spin.model: 99 raise RelaxNoModelError 100 101 # Loop over the spectral time points. 102 for id in list(cdp.relax_times.keys()): 103 # Back calculate the value. 104 value = back_calc(spin=spin, relax_time_id=id) 105 106 # Append the value. 107 mc_data[id] = value 108 109 # Return the MC data. 110 return mc_data

111 112

113 - def data_init(self, data_cont, sim=False):

114 """Initialise the spin specific data structures. 115 116 @param data_cont: The spin container. 117 @type data_cont: SpinContainer instance 118 @keyword sim: The Monte Carlo simulation flag, which if true will initialise the simulation data structure. 119 @type sim: bool 120 """ 121 122 # Loop over the data structure names. 123 for name in self.data_names(set='params'): 124 # Data structures which are initially empty arrays. 125 list_data = [ 'params' ] 126 if name in list_data: 127 init_data = [] 128 129 # Otherwise initialise the data structure to None. 130 else: 131 init_data = None 132 133 # If the name is not in 'data_cont', add it. 134 if not hasattr(data_cont, name): 135 setattr(data_cont, name, init_data)

136 137

138 - def grid_search(self, lower=None, upper=None, inc=None, constraints=True, verbosity=1, sim_index=None):

139 """The exponential curve fitting grid search method. 140 141 @keyword lower: The lower bounds of the grid search which must be equal to the number of parameters in the model. 142 @type lower: array of numbers 143 @keyword upper: The upper bounds of the grid search which must be equal to the number of parameters in the model. 144 @type upper: array of numbers 145 @keyword inc: The increments for each dimension of the space for the grid search. The number of elements in the array must equal to the number of parameters in the model. 146 @type inc: array of int 147 @keyword constraints: If True, constraints are applied during the grid search (eliminating parts of the grid). If False, no constraints are used. 148 @type constraints: bool 149 @keyword verbosity: A flag specifying the amount of information to print. The higher the value, the greater the verbosity. 150 @type verbosity: int 151 @keyword sim_index: The index of the simulation to apply the grid search to. If None, the normal model is optimised. 152 @type sim_index: int 153 """ 154 155 # Minimisation. 156 self.minimise(min_algor='grid', lower=lower, upper=upper, inc=inc, constraints=constraints, verbosity=verbosity, sim_index=sim_index)

157 158

159 - def minimise(self, min_algor=None, min_options=None, func_tol=None, grad_tol=None, max_iterations=None, constraints=False, scaling=True, verbosity=0, sim_index=None, lower=None, upper=None, inc=None):

160 """Relaxation curve fitting minimisation method. 161 162 @keyword min_algor: The minimisation algorithm to use. 163 @type min_algor: str 164 @keyword min_options: An array of options to be used by the minimisation algorithm. 165 @type min_options: array of str 166 @keyword func_tol: The function tolerance which, when reached, terminates optimisation. Setting this to None turns of the check. 167 @type func_tol: None or float 168 @keyword grad_tol: The gradient tolerance which, when reached, terminates optimisation. Setting this to None turns of the check. 169 @type grad_tol: None or float 170 @keyword max_iterations: The maximum number of iterations for the algorithm. 171 @type max_iterations: int 172 @keyword constraints: If True, constraints are used during optimisation. 173 @type constraints: bool 174 @keyword scaling: If True, diagonal scaling is enabled during optimisation to allow the problem to be better conditioned. 175 @type scaling: bool 176 @keyword verbosity: The amount of information to print. The higher the value, the greater the verbosity. 177 @type verbosity: int 178 @keyword sim_index: The index of the simulation to optimise. This should be None if normal optimisation is desired. 179 @type sim_index: None or int 180 @keyword lower: The lower bounds of the grid search which must be equal to the number of parameters in the model. This optional argument is only used when doing a grid search. 181 @type lower: array of numbers 182 @keyword upper: The upper bounds of the grid search which must be equal to the number of parameters in the model. This optional argument is only used when doing a grid search. 183 @type upper: array of numbers 184 @keyword inc: The increments for each dimension of the space for the grid search. The number of elements in the array must equal to the number of parameters in the model. This argument is only used when doing a grid search. 185 @type inc: array of int 186 """ 187 188 # Test if sequence data is loaded. 189 if not exists_mol_res_spin_data(): 190 raise RelaxNoSequenceError 191 192 # Loop over the sequence. 193 for spin, mol_name, res_num, res_name in spin_loop(full_info=True): 194 # Skip deselected spins. 195 if not spin.select: 196 continue 197 198 # Skip spins which have no data. 199 if not hasattr(spin, 'peak_intensity'): 200 continue 201 202 # Create the initial parameter vector. 203 param_vector = assemble_param_vector(spin=spin) 204 205 # Diagonal scaling. 206 scaling_matrix = assemble_scaling_matrix(spin=spin, scaling=scaling) 207 if len(scaling_matrix): 208 param_vector = dot(inv(scaling_matrix), param_vector) 209 210 # Get the grid search minimisation options. 211 if match('^[Gg]rid', min_algor): 212 inc, lower_new, upper_new = grid_search_setup(spin=spin, param_vector=param_vector, lower=lower, upper=upper, inc=inc, scaling_matrix=scaling_matrix) 213 214 # Linear constraints. 215 if constraints: 216 A, b = linear_constraints(spin=spin, scaling_matrix=scaling_matrix) 217 else: 218 A, b = None, None 219 220 # Print out. 221 if verbosity >= 1: 222 # Get the spin id string. 223 spin_id = generate_spin_id_unique(mol_name=mol_name, res_num=res_num, res_name=res_name, spin_num=spin.num, spin_name=spin.name) 224 225 # Individual spin printout. 226 if verbosity >= 2: 227 print("\n\n") 228 229 string = "Fitting to spin " + repr(spin_id) 230 print("\n\n" + string) 231 print(len(string) * '~') 232 233 234 # Initialise the function to minimise. 235 ###################################### 236 237 # The keys. 238 keys = list(spin.peak_intensity.keys()) 239 240 # The peak intensities and times. 241 values = [] 242 errors = [] 243 times = [] 244 for key in keys: 245 # The values. 246 if sim_index == None: 247 values.append(spin.peak_intensity[key]) 248 else: 249 values.append(spin.peak_intensity_sim[sim_index][key]) 250 251 # The errors. 252 errors.append(spin.peak_intensity_err[key]) 253 254 # The relaxation times. 255 times.append(cdp.relax_times[key]) 256 257 # The scaling matrix in a diagonalised list form. 258 scaling_list = [] 259 for i in range(len(scaling_matrix)): 260 scaling_list.append(scaling_matrix[i, i]) 261 262 setup(num_params=len(spin.params), num_times=len(values), values=values, sd=errors, relax_times=times, scaling_matrix=scaling_list) 263 264 265 # Setup the minimisation algorithm when constraints are present. 266 ################################################################ 267 268 if constraints and not match('^[Gg]rid', min_algor): 269 algor = min_options[0] 270 else: 271 algor = min_algor 272 273 274 # Levenberg-Marquardt minimisation. 275 ################################### 276 277 if match('[Ll][Mm]$', algor) or match('[Ll]evenburg-[Mm]arquardt$', algor): 278 # Reconstruct the error data structure. 279 lm_error = zeros(len(spin.relax_times), float64) 280 index = 0 281 for k in range(len(spin.relax_times)): 282 lm_error[index:index+len(relax_error[k])] = relax_error[k] 283 index = index + len(relax_error[k]) 284 285 min_options = min_options + (self.relax_fit.lm_dri, lm_error) 286 287 288 # Minimisation. 289 ############### 290 291 # Grid search. 292 if search('^[Gg]rid', min_algor): 293 results = grid(func=func_wrapper, args=(), num_incs=inc, lower=lower_new, upper=upper_new, A=A, b=b, verbosity=verbosity) 294 295 # Unpack the results. 296 param_vector, chi2, iter_count, warning = results 297 f_count = iter_count 298 g_count = 0.0 299 h_count = 0.0 300 301 # Minimisation. 302 else: 303 results = generic_minimise(func=func_wrapper, dfunc=dfunc_wrapper, d2func=d2func_wrapper, args=(), x0=param_vector, min_algor=min_algor, min_options=min_options, func_tol=func_tol, grad_tol=grad_tol, maxiter=max_iterations, A=A, b=b, full_output=True, print_flag=verbosity) 304 305 # Unpack the results. 306 if results == None: 307 return 308 param_vector, chi2, iter_count, f_count, g_count, h_count, warning = results 309 310 # Scaling. 311 if scaling: 312 param_vector = dot(scaling_matrix, param_vector) 313 314 # Disassemble the parameter vector. 315 disassemble_param_vector(param_vector=param_vector, spin=spin, sim_index=sim_index) 316 317 # Monte Carlo minimisation statistics. 318 if sim_index != None: 319 # Chi-squared statistic. 320 spin.chi2_sim[sim_index] = chi2 321 322 # Iterations. 323 spin.iter_sim[sim_index] = iter_count 324 325 # Function evaluations. 326 spin.f_count_sim[sim_index] = f_count 327 328 # Gradient evaluations. 329 spin.g_count_sim[sim_index] = g_count 330 331 # Hessian evaluations. 332 spin.h_count_sim[sim_index] = h_count 333 334 # Warning. 335 spin.warning_sim[sim_index] = warning 336 337 338 # Normal statistics. 339 else: 340 # Chi-squared statistic. 341 spin.chi2 = chi2 342 343 # Iterations. 344 spin.iter = iter_count 345 346 # Function evaluations. 347 spin.f_count = f_count 348 349 # Gradient evaluations. 350 spin.g_count = g_count 351 352 # Hessian evaluations. 353 spin.h_count = h_count 354 355 # Warning. 356 spin.warning = warning

357 358

359 - def overfit_deselect(self, data_check=True, verbose=True):

360 """Deselect spins which have insufficient data to support minimisation. 361 362 @keyword data_check: A flag to signal if the presence of base data is to be checked for. 363 @type data_check: bool 364 @keyword verbose: A flag which if True will allow printouts. 365 @type verbose: bool 366 """ 367 368 # Print out. 369 if verbose: 370 print("\nOver-fit spin deselection:") 371 372 # Test the sequence data exists. 373 if not exists_mol_res_spin_data(): 374 raise RelaxNoSequenceError 375 376 # Loop over spin data. 377 deselect_flag = False 378 for spin, spin_id in spin_loop(return_id=True): 379 # Skip deselected spins. 380 if not spin.select: 381 continue 382 383 # Check if data exists. 384 if not hasattr(spin, 'peak_intensity'): 385 warn(RelaxDeselectWarning(spin_id, 'missing intensity data')) 386 spin.select = False 387 deselect_flag = True 388 continue 389 390 # Require 3 or more data points. 391 elif len(spin.peak_intensity) < 3: 392 warn(RelaxDeselectWarning(spin_id, 'insufficient data, 3 or more data points are required')) 393 spin.select = False 394 deselect_flag = True 395 continue 396 397 # Check that the number of relaxation times is complete. 398 if len(spin.peak_intensity) != len(cdp.relax_times): 399 raise RelaxError("The %s peak intensity points of the spin '%s' does not match the expected number of %s (the IDs %s do not match %s)." % (len(spin.peak_intensity), spin_id, len(cdp.relax_times), list(spin.peak_intensity.keys()), list(cdp.relax_times.keys()))) 400 401 # Final printout. 402 if verbose and not deselect_flag: 403 print("No spins have been deselected.")

404 405

406 - def return_data(self, data_id=None):

407 """Function for returning the peak intensity data structure. 408 409 @keyword data_id: The spin identification string, as yielded by the base_data_loop() generator method. 410 @type data_id: str 411 @return: The peak intensity data structure. 412 @rtype: list of float 413 """ 414 415 # Get the spin container. 416 spin = return_spin(data_id) 417 418 # Return the peak intensities. 419 return spin.peak_intensity

420 421

422 - def return_error(self, data_id):

423 """Return the standard deviation data structure. 424 425 @param data_id: The spin identification string, as yielded by the base_data_loop() generator 426 method. 427 @type data_id: str 428 @return: The standard deviation data structure. 429 @rtype: list of float 430 """ 431 432 # Get the spin container. 433 spin = return_spin(data_id) 434 435 # Return the error list. 436 return spin.peak_intensity_err

437 438

439 - def sim_pack_data(self, data_id, sim_data):

440 """Pack the Monte Carlo simulation data. 441 442 @param data_id: The spin identification string, as yielded by the base_data_loop() generator method. 443 @type data_id: str 444 @param sim_data: The Monte Carlo simulation data. 445 @type sim_data: list of float 446 """ 447 448 # Get the spin container. 449 spin = return_spin(data_id) 450 451 # Create the data structure. 452 spin.peak_intensity_sim = sim_data

Source Code for Module specific_analyses.relax_fit.api