specific_analyses.relax_disp.estimate

1 ############################################################################### 2 # # 3 # Copyright (C) 2014 Troels E. Linnet # 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 """Target functions for relaxation exponential curve fitting with both minfx and scipy.optimize.leastsq.""" 24 25 # Python module imports. 26 from copy import deepcopy 27 from numpy import array, asarray, diag, exp, log, ones, sqrt, sum, transpose, zeros 28 from minfx.generic import generic_minimise 29 import sys 30 from warnings import warn 31 32 # relax module imports. 33 from dep_check import C_module_exp_fn, scipy_module 34 from lib.dispersion.variables import MODEL_R2EFF 35 from lib.errors import RelaxError 36 from lib.statistics import multifit_covar 37 from lib.text.sectioning import subsection 38 from lib.warnings import RelaxWarning 39 from pipe_control.mol_res_spin import generate_spin_string, spin_loop 40 from specific_analyses.relax_disp.checks import check_model_type 41 from specific_analyses.relax_disp.data import average_intensity, loop_exp_frq_offset_point, loop_time, return_param_key_from_data 42 from specific_analyses.relax_disp.parameters import disassemble_param_vector 43 from target_functions.chi2 import chi2_rankN, dchi2 44 from target_functions.relax_fit_wrapper import Relax_fit_opt 45 46 # Scipy installed. 47 if scipy_module: 48 # Import leastsq. 49 from scipy.optimize import leastsq 50 51

52 -def estimate_r2eff_err(spin_id=None, epsrel=0.0, verbosity=1):

53 """This will estimate the R2eff and i0 errors from the covariance matrix Qxx. Qxx is calculated from the Jacobian matrix and the optimised parameters. 54 55 @keyword spin_id: The spin identification string. 56 @type spin_id: str 57 @param epsrel: Any columns of R which satisfy |R_{kk}| <= epsrel |R_{11}| are considered linearly-dependent and are excluded from the covariance matrix, where the corresponding rows and columns of the covariance matrix are set to zero. 58 @type epsrel: float 59 @keyword verbosity: The amount of information to print. The higher the value, the greater the verbosity. 60 @type verbosity: int 61 """ 62 63 # Check that the C modules have been compiled. 64 if not C_module_exp_fn: 65 raise RelaxError("Relaxation curve fitting is not available. Try compiling the C modules on your platform.") 66 67 # Perform checks. 68 check_model_type(model=MODEL_R2EFF) 69 70 # Loop over the spins. 71 for cur_spin, mol_name, resi, resn, cur_spin_id in spin_loop(selection=spin_id, full_info=True, return_id=True, skip_desel=True): 72 # Generate spin string. 73 spin_string = generate_spin_string(spin=cur_spin, mol_name=mol_name, res_num=resi, res_name=resn) 74 75 # Raise Error, if not optimised. 76 if not (hasattr(cur_spin, 'r2eff') and hasattr(cur_spin, 'i0')): 77 raise RelaxError("Spin %s does not contain optimised 'r2eff' and 'i0' values. Try execute: minimise.execute(min_algor='Newton', constraints=False)"%(spin_string)) 78 79 # Raise warning, if gradient count is 0. This could point to a lack of minimisation first. 80 if hasattr(cur_spin, 'g_count'): 81 if getattr(cur_spin, 'g_count') == 0.0: 82 text = "Spin %s contains a gradient count of 0.0. Is the R2eff parameter optimised? Try execute: minimise.execute(min_algor='Newton', constraints=False)" %(spin_string) 83 warn(RelaxWarning("%s." % text)) 84 85 # Print information. 86 if verbosity >= 1: 87 # Individual spin block section. 88 top = 2 89 if verbosity >= 2: 90 top += 2 91 subsection(file=sys.stdout, text="Estimating R2eff error for spin: %s"%spin_string, prespace=top) 92 93 # Loop over each spectrometer frequency and dispersion point. 94 for exp_type, frq, offset, point, ei, mi, oi, di in loop_exp_frq_offset_point(return_indices=True): 95 # The parameter key. 96 param_key = return_param_key_from_data(exp_type=exp_type, frq=frq, offset=offset, point=point) 97 98 # Extract values. 99 r2eff = getattr(cur_spin, 'r2eff')[param_key] 100 i0 = getattr(cur_spin, 'i0')[param_key] 101 102 # Pack data 103 param_vector = [r2eff, i0] 104 105 # The peak intensities, errors and times. 106 values = [] 107 errors = [] 108 times = [] 109 for time in loop_time(exp_type=exp_type, frq=frq, offset=offset, point=point): 110 values.append(average_intensity(spin=cur_spin, exp_type=exp_type, frq=frq, offset=offset, point=point, time=time)) 111 errors.append(average_intensity(spin=cur_spin, exp_type=exp_type, frq=frq, offset=offset, point=point, time=time, error=True)) 112 times.append(time) 113 114 # Convert to numpy array. 115 values = asarray(values) 116 errors = asarray(errors) 117 times = asarray(times) 118 119 # Initialise data in C code. 120 scaling_list = [1.0, 1.0] 121 model = Relax_fit_opt(model='exp', num_params=len(param_vector), values=values, errors=errors, relax_times=times, scaling_matrix=scaling_list) 122 123 # Use the direct Jacobian from function. 124 jacobian_matrix_exp = transpose(asarray( model.jacobian(param_vector) ) ) 125 weights = 1. / errors**2 126 127 # Get the co-variance 128 pcov = multifit_covar(J=jacobian_matrix_exp, weights=weights) 129 130 # To compute one standard deviation errors on the parameters, take the square root of the diagonal covariance. 131 param_vector_error = sqrt(diag(pcov)) 132 133 # Extract values. 134 r2eff_err, i0_err = param_vector_error 135 136 # Copy r2eff dictionary, to r2eff_err dictionary. They have same keys to the dictionary, 137 if not hasattr(cur_spin, 'r2eff_err'): 138 setattr(cur_spin, 'r2eff_err', deepcopy(getattr(cur_spin, 'r2eff'))) 139 if not hasattr(cur_spin, 'i0_err'): 140 setattr(cur_spin, 'i0_err', deepcopy(getattr(cur_spin, 'i0'))) 141 142 # Set error. 143 cur_spin.r2eff_err[param_key] = r2eff_err 144 cur_spin.i0_err[param_key] = i0_err 145 146 # Get other relevant information. 147 chi2 = getattr(cur_spin, 'chi2') 148 149 # Print information. 150 print_strings = [] 151 if verbosity >= 1: 152 # Add print strings. 153 point_info = "%s at %3.1f MHz, for offset=%3.3f ppm and dispersion point %-5.1f, with %i time points." % (exp_type, frq/1E6, offset, point, len(times)) 154 print_strings.append(point_info) 155 156 par_info = "r2eff=%3.3f r2eff_err=%3.4f, i0=%6.1f, i0_err=%3.4f, chi2=%3.3f.\n" % ( r2eff, r2eff_err, i0, i0_err, chi2) 157 print_strings.append(par_info) 158 159 if verbosity >= 2: 160 time_info = ', '.join(map(str, times)) 161 print_strings.append('For time array: '+time_info+'.\n\n') 162 163 # Print info 164 if len(print_strings) > 0: 165 for print_string in print_strings: 166 print(print_string),

167 168 169 #### This class is only for testing. 170

171 -class Exp:

172 - def __init__(self, verbosity=1):

173 """Class for to set settings for minimisation and dispersion target functions for minimisation. 174 175 This class contains minimisation functions for both minfx and scipy.optimize.leastsq. 176 177 @keyword verbosity: The amount of information to print. The higher the value, the greater the verbosity. 178 @type verbosity: int 179 """ 180 181 # Store. 182 self.verbosity = verbosity 183 184 # Initialize standard settings. 185 self.set_settings_leastsq() 186 self.set_settings_minfx()

187 188

189 - def setup_data(self, times=None, values=None, errors=None):

190 """Setup for minimisation with minfx. 191 192 @keyword times: The time points. 193 @type times: numpy array 194 @keyword values: The measured intensity values per time point. 195 @type values: numpy array 196 @keyword errors: The standard deviation of the measured intensity values per time point. 197 @type errors: numpy array 198 """ 199 200 # Store variables. 201 self.times = times 202 self.values = values 203 self.errors = errors

204 205

206 - def set_settings_leastsq(self, ftol=1e-15, xtol=1e-15, maxfev=10000000, factor=100.0):

207 """Setup options to scipy.optimize.leastsq. 208 209 @keyword ftol: The function tolerance for the relative error desired in the sum of squares, parsed to leastsq. 210 @type ftol: float 211 @keyword xtol: The error tolerance for the relative error desired in the approximate solution, parsed to leastsq. 212 @type xtol: float 213 @keyword maxfev: The maximum number of function evaluations, parsed to leastsq. If zero, then 100*(N+1) is the maximum function calls. N is the number of elements in x0=[r2eff, i0]. 214 @type maxfev: int 215 @keyword factor: The initial step bound, parsed to leastsq. It determines the initial step bound (''factor * || diag * x||''). Should be in the interval (0.1, 100). 216 @type factor: float 217 """ 218 219 # Store settings. 220 self.ftol = ftol 221 self.xtol = xtol 222 self.maxfev = maxfev 223 self.factor = factor

224 225

226 - def set_settings_minfx(self, scaling_matrix=None, min_algor='simplex', c_code=True, constraints=False, chi2_jacobian=False, func_tol=1e-25, grad_tol=None, max_iterations=10000000):

227 """Setup options to minfx. 228 229 @keyword scaling_matrix: The square and diagonal scaling matrix. 230 @type scaling_matrix: numpy rank-2 float array 231 @keyword min_algor: The minimisation algorithm 232 @type min_algor: string 233 @keyword c_code: If optimise with C code. 234 @type c_code: bool 235 @keyword constraints: If constraints should be used. 236 @type constraints: bool 237 @keyword chi2_jacobian: If the chi2 Jacobian should be used. 238 @type chi2_jacobian: bool 239 @keyword func_tol: The function tolerance which, when reached, terminates optimisation. Setting this to None turns of the check. 240 @type func_tol: None or float 241 @keyword grad_tol: The gradient tolerance which, when reached, terminates optimisation. Setting this to None turns of the check. 242 @type grad_tol: None or float 243 @keyword max_iterations: The maximum number of iterations for the algorithm. 244 @type max_iterations: int 245 """ 246 247 # Store variables. 248 self.scaling_matrix = scaling_matrix 249 self.c_code = c_code 250 self.chi2_jacobian = chi2_jacobian 251 252 # Scaling initialisation. 253 self.scaling_flag = False 254 if self.scaling_matrix != None: 255 self.scaling_flag = True 256 257 # Set algorithm. 258 self.min_algor = min_algor 259 260 # Define if constraints should be used. 261 self.constraints = constraints 262 263 # Settings to minfx. 264 self.func_tol = func_tol 265 self.grad_tol = grad_tol 266 self.max_iterations = max_iterations 267 268 # Options to minfx depends if contraints is set. 269 if self.constraints: 270 self.min_options = ('%s'%(self.min_algor),) 271 self.min_algor = 'Log barrier' 272 self.A = array([ [ 1., 0.], 273 [-1., 0.], 274 [ 0., 1.]] ) 275 self.b = array([ 0., -200., 0.]) 276 277 else: 278 self.min_options = () 279 self.A = None 280 self.b = None

281 282

283 - def estimate_x0_exp(self, times=None, values=None):

284 """Estimate starting parameter x0 = [r2eff_est, i0_est], by converting the exponential curve to a linear problem. Then solving by linear least squares of: ln(Intensity[j]) = ln(i0) - time[j]* r2eff. 285 286 @keyword times: The time points. 287 @type times: numpy array 288 @keyword values: The measured intensity values per time point. 289 @type values: numpy array 290 @return: The list with estimated r2eff and i0 parameter for optimisation, [r2eff_est, i0_est] 291 @rtype: list 292 """ 293 294 # Convert to linear problem. 295 w = log(values) 296 x = - 1. * times 297 n = len(times) 298 299 # Solve by linear least squares. 300 b = (sum(x*w) - 1./n * sum(x) * sum(w) ) / ( sum(x**2) - 1./n * (sum(x))**2 ) 301 a = 1./n * sum(w) - b * 1./n * sum(x) 302 303 # Convert back from linear to exp function. Best estimate for parameter. 304 r2eff_est = b 305 i0_est = exp(a) 306 307 # Return. 308 return [r2eff_est, i0_est]

309 310

311 - def func_exp(self, params=None, times=None):

312 """Calculate the function values of exponential function. 313 314 @param params: The vector of parameter values. 315 @type params: numpy rank-1 float array 316 @keyword times: The time points. 317 @type times: numpy array 318 @return: The function values. 319 @rtype: numpy array 320 """ 321 322 # Unpack 323 r2eff, i0 = params 324 325 # Calculate. 326 return i0 * exp( -times * r2eff)

327 328

329 - def func_exp_residual(self, params=None, times=None, values=None):

330 """Calculate the residual vector betwen measured values and the function values. 331 332 @param params: The vector of parameter values. 333 @type params: numpy rank-1 float array 334 @keyword times: The time points. 335 @type times: numpy array 336 @param values: The measured values. 337 @type values: numpy array 338 @return: The residuals. 339 @rtype: numpy array 340 """ 341 342 # Let the vector K be the vector of the residuals. A residual is the difference between the observation and the equation calculated using the initial values. 343 K = values - self.func_exp(params=params, times=times) 344 345 # Return 346 return K

347 348

349 - def func_exp_weighted_residual(self, params=None, times=None, values=None, errors=None):

350 """Calculate the weighted residual vector betwen measured values and the function values. 351 352 @param params: The vector of parameter values. 353 @type params: numpy rank-1 float array 354 @keyword times: The time points. 355 @type times: numpy array 356 @param values: The measured values. 357 @type values: numpy array 358 @param errors: The standard deviation of the measured intensity values per time point. 359 @type errors: numpy array 360 @return: The weighted residuals. 361 @rtype: numpy array 362 """ 363 364 # Let the vector Kw be the vector of the weighted residuals. A residual is the difference between the observation and the equation calculated using the initial values. 365 Kw = 1. / errors * self.func_exp_residual(params=params, times=times, values=values) 366 367 # Return 368 return Kw

369 370

371 - def func_exp_grad(self, params=None, times=None):

372 """The gradient (Jacobian matrix) of func_exp for Co-variance calculation. 373 374 @param params: The vector of parameter values. 375 @type params: numpy rank-1 float array 376 @keyword times: The time points. 377 @type times: numpy array 378 @return: The Jacobian matrix with 'm' rows of function derivatives per 'n' columns of parameters. 379 @rtype: numpy array 380 """ 381 382 # Unpack the parameter values. 383 r2eff = params[0] 384 i0 = params[1] 385 386 # Make partial derivative, with respect to r2eff. 387 d_exp_d_r2eff = -i0 * times * exp(-r2eff * times) 388 389 # Make partial derivative, with respect to i0. 390 d_exp_d_i0 = exp(-r2eff * times) 391 392 # Define Jacobian as m rows with function derivatives and n columns of parameters. 393 jacobian_matrix_exp = transpose(array( [d_exp_d_r2eff, d_exp_d_i0] ) ) 394 395 # Return Jacobian matrix. 396 return jacobian_matrix_exp

397 398

399 - def func_exp_chi2(self, params=None, times=None, values=None, errors=None):

400 """Target function for minimising chi2 in minfx, for exponential fit. 401 402 @param params: The vector of parameter values. 403 @type params: numpy rank-1 float array 404 @keyword times: The time points. 405 @type times: numpy array 406 @param values: The measured values. 407 @type values: numpy array 408 @param errors: The standard deviation of the measured intensity values per time point. 409 @type errors: numpy array 410 @return: The chi2 value. 411 @rtype: float 412 """ 413 414 # Calculate. 415 back_calc = self.func_exp(params=params, times=times) 416 417 # Return the total chi-squared value. 418 chi2 = chi2_rankN(data=values, back_calc_vals=back_calc, errors=errors) 419 420 # Calculate and return the chi-squared value. 421 return chi2

422 423

424 - def func_exp_chi2_grad(self, params=None, times=None, values=None, errors=None):

425 """Target function for the gradient (Jacobian matrix) to minfx, for exponential fit . 426 427 @param params: The vector of parameter values. 428 @type params: numpy rank-1 float array 429 @keyword times: The time points. 430 @type times: numpy array 431 @param values: The measured values. 432 @type values: numpy array 433 @param errors: The standard deviation of the measured intensity values per time point. 434 @type errors: numpy array 435 @return: The Jacobian matrix with 'm' rows of function derivatives per 'n' columns of parameters, which have been summed together. 436 @rtype: numpy array 437 """ 438 439 # Get the back calc. 440 back_calc = self.func_exp(params=params, times=times) 441 442 # Get the Jacobian, with partial derivative, with respect to r2eff and i0. 443 exp_grad = self.func_exp_grad(params=params, times=times) 444 445 # Transpose back, to get rows. 446 exp_grad_t = transpose(exp_grad) 447 448 # n is number of fitted parameters. 449 n = len(params) 450 451 # Define array to update parameters in. 452 jacobian_chi2_minfx = zeros([n]) 453 454 # Update value elements. 455 dchi2(dchi2=jacobian_chi2_minfx, M=n, data=values, back_calc_vals=back_calc, back_calc_grad=exp_grad_t, errors=errors) 456 457 # Return Jacobian matrix. 458 return jacobian_chi2_minfx

459 460

461 - def func_exp_chi2_grad_array(self, params=None, times=None, values=None, errors=None):

462 """Return the gradient (Jacobian matrix) of func_exp_chi2() for parameter co-variance error estimation. This needs return as array. 463 464 @param params: The vector of parameter values. 465 @type params: numpy rank-1 float array 466 @keyword times: The time points. 467 @type times: numpy array 468 @keyword values: The measured intensity values per time point. 469 @type values: numpy array 470 @keyword errors: The standard deviation of the measured intensity values per time point. 471 @type errors: numpy array 472 @return: The Jacobian matrix with 'm' rows of function derivatives per 'n' columns of parameters. 473 @rtype: numpy array 474 """ 475 476 # Unpack the parameter values. 477 r2eff = params[0] 478 i0 = params[1] 479 480 # Make partial derivative, with respect to r2eff. 481 d_chi2_d_r2eff = 2.0 * i0 * times * ( -i0 * exp( -r2eff * times) + values) * exp( -r2eff * times ) / errors**2 482 483 # Make partial derivative, with respect to i0. 484 d_chi2_d_i0 = - 2.0 * ( -i0 * exp( -r2eff * times) + values) * exp( -r2eff * times) / errors**2 485 486 # Define Jacobian as m rows with function derivatives and n columns of parameters. 487 jacobian_matrix_exp_chi2 = transpose(array( [d_chi2_d_r2eff, d_chi2_d_i0] ) ) 488 489 return jacobian_matrix_exp_chi2

490 491

492 -def estimate_r2eff(method='minfx', min_algor='simplex', c_code=True, constraints=False, chi2_jacobian=False, spin_id=None, ftol=1e-15, xtol=1e-15, maxfev=10000000, factor=100.0, verbosity=1):

493 """Estimate r2eff and errors by exponential curve fitting with scipy.optimize.leastsq or minfx. 494 495 THIS IS ONLY FOR TESTING. 496 497 scipy.optimize.leastsq is a wrapper around MINPACK's lmdif and lmder algorithms. 498 499 MINPACK is a FORTRAN90 library which solves systems of nonlinear equations, or carries out the least squares minimization of the residual of a set of linear or nonlinear equations. 500 501 Errors are calculated by taking the square root of the reported co-variance. 502 503 This can be an huge time saving step, when performing model fitting in R1rho. 504 Errors of R2eff values, are normally estimated by time-consuming Monte-Carlo simulations. 505 506 Initial guess for the starting parameter x0 = [r2eff_est, i0_est], is by converting the exponential curve to a linear problem. 507 Then solving initial guess by linear least squares of: ln(Intensity[j]) = ln(i0) - time[j]* r2eff. 508 509 510 @keyword method: The method to minimise and estimate errors. Options are: 'minfx' or 'scipy.optimize.leastsq'. 511 @type method: string 512 @keyword min_algor: The minimisation algorithm 513 @type min_algor: string 514 @keyword c_code: If optimise with C code. 515 @type c_code: bool 516 @keyword constraints: If constraints should be used. 517 @type constraints: bool 518 @keyword chi2_jacobian: If the chi2 Jacobian should be used. 519 @type chi2_jacobian: bool 520 @keyword spin_id: The spin identification string. 521 @type spin_id: str 522 @keyword ftol: The function tolerance for the relative error desired in the sum of squares, parsed to leastsq. 523 @type ftol: float 524 @keyword xtol: The error tolerance for the relative error desired in the approximate solution, parsed to leastsq. 525 @type xtol: float 526 @keyword maxfev: The maximum number of function evaluations, parsed to leastsq. If zero, then 100*(N+1) is the maximum function calls. N is the number of elements in x0=[r2eff, i0]. 527 @type maxfev: int 528 @keyword factor: The initial step bound, parsed to leastsq. It determines the initial step bound (''factor * || diag * x||''). Should be in the interval (0.1, 100). 529 @type factor: float 530 @keyword verbosity: The amount of information to print. The higher the value, the greater the verbosity. 531 @type verbosity: int 532 """ 533 534 # Perform checks. 535 check_model_type(model=MODEL_R2EFF) 536 537 # Check that the C modules have been compiled. 538 if not C_module_exp_fn and method == 'minfx': 539 raise RelaxError("Relaxation curve fitting is not available. Try compiling the C modules on your platform.") 540 541 # Set class scipy setting. 542 E = Exp(verbosity=verbosity) 543 E.set_settings_leastsq(ftol=ftol, xtol=xtol, maxfev=maxfev, factor=factor) 544 545 # Check if intensity errors have already been calculated by the user. 546 precalc = True 547 for cur_spin, mol_name, resi, resn, cur_spin_id in spin_loop(selection=spin_id, full_info=True, return_id=True, skip_desel=True): 548 # No structure. 549 if not hasattr(cur_spin, 'peak_intensity_err'): 550 precalc = False 551 break 552 553 # Determine if a spectrum ID is missing from the list. 554 for id in cdp.spectrum_ids: 555 if id not in cur_spin.peak_intensity_err: 556 precalc = False 557 break 558 559 # Loop over the spins. 560 for cur_spin, mol_name, resi, resn, cur_spin_id in spin_loop(selection=spin_id, full_info=True, return_id=True, skip_desel=True): 561 # Generate spin string. 562 spin_string = generate_spin_string(spin=cur_spin, mol_name=mol_name, res_num=resi, res_name=resn) 563 564 # Print information. 565 if E.verbosity >= 1: 566 # Individual spin block section. 567 top = 2 568 if E.verbosity >= 2: 569 top += 2 570 subsection(file=sys.stdout, text="Fitting with %s to: %s"%(method, spin_string), prespace=top) 571 if method == 'minfx': 572 subsection(file=sys.stdout, text="min_algor='%s', c_code=%s, constraints=%s, chi2_jacobian?=%s"%(min_algor, c_code, constraints, chi2_jacobian), prespace=0) 573 574 # Loop over each spectrometer frequency and dispersion point. 575 for exp_type, frq, offset, point, ei, mi, oi, di in loop_exp_frq_offset_point(return_indices=True): 576 # The parameter key. 577 param_key = return_param_key_from_data(exp_type=exp_type, frq=frq, offset=offset, point=point) 578 579 # The peak intensities, errors and times. 580 values = [] 581 errors = [] 582 times = [] 583 for time in loop_time(exp_type=exp_type, frq=frq, offset=offset, point=point): 584 values.append(average_intensity(spin=cur_spin, exp_type=exp_type, frq=frq, offset=offset, point=point, time=time)) 585 errors.append(average_intensity(spin=cur_spin, exp_type=exp_type, frq=frq, offset=offset, point=point, time=time, error=True)) 586 times.append(time) 587 588 # Convert to numpy array. 589 values = asarray(values) 590 errors = asarray(errors) 591 times = asarray(times) 592 593 # Initialise data. 594 E.setup_data(values=values, errors=errors, times=times) 595 596 # Get the result based on method. 597 if method == 'scipy.optimize.leastsq': 598 # Acquire results. 599 results = minimise_leastsq(E=E) 600 601 elif method == 'minfx': 602 # Set settings. 603 E.set_settings_minfx(min_algor=min_algor, c_code=c_code, chi2_jacobian=chi2_jacobian, constraints=constraints) 604 605 # Acquire results. 606 results = minimise_minfx(E=E) 607 else: 608 raise RelaxError("Method for minimisation not known. Try setting: method='scipy.optimize.leastsq'.") 609 610 # Unpack results 611 param_vector, param_vector_error, chi2, iter_count, f_count, g_count, h_count, warning = results 612 613 # Extract values. 614 r2eff = param_vector[0] 615 i0 = param_vector[1] 616 r2eff_err = param_vector_error[0] 617 i0_err = param_vector_error[1] 618 619 # Disassemble the parameter vector. 620 disassemble_param_vector(param_vector=param_vector, spins=[cur_spin], key=param_key) 621 622 # Errors. 623 if not hasattr(cur_spin, 'r2eff_err'): 624 setattr(cur_spin, 'r2eff_err', deepcopy(getattr(cur_spin, 'r2eff'))) 625 if not hasattr(cur_spin, 'i0_err'): 626 setattr(cur_spin, 'i0_err', deepcopy(getattr(cur_spin, 'i0'))) 627 628 # Set error. 629 cur_spin.r2eff_err[param_key] = r2eff_err 630 cur_spin.i0_err[param_key] = i0_err 631 632 # Chi-squared statistic. 633 cur_spin.chi2 = chi2 634 635 # Iterations. 636 cur_spin.f_count = f_count 637 638 # Warning. 639 cur_spin.warning = warning 640 641 # Print information. 642 print_strings = [] 643 if E.verbosity >= 1: 644 # Add print strings. 645 point_info = "%s at %3.1f MHz, for offset=%3.3f ppm and dispersion point %-5.1f, with %i time points." % (exp_type, frq/1E6, offset, point, len(times)) 646 print_strings.append(point_info) 647 648 par_info = "r2eff=%3.3f r2eff_err=%3.4f, i0=%6.1f, i0_err=%3.4f, chi2=%3.3f.\n" % ( r2eff, r2eff_err, i0, i0_err, chi2) 649 print_strings.append(par_info) 650 651 if E.verbosity >= 2: 652 time_info = ', '.join(map(str, times)) 653 print_strings.append('For time array: '+time_info+'.\n\n') 654 655 # Print info 656 if len(print_strings) > 0: 657 for print_string in print_strings: 658 print(print_string),

659 660

661 -def minimise_leastsq(E=None):

662 """Estimate r2eff and errors by exponential curve fitting with scipy.optimize.leastsq. 663 664 @keyword E: The Exponential function class, which contain data and functions. 665 @type E: class 666 @return: Packed list with optimised parameter, estimated parameter error, chi2, iter_count, f_count, g_count, h_count, warning 667 @rtype: list 668 """ 669 670 # Check that scipy.optimize.leastsq is available. 671 if not scipy_module: 672 raise RelaxError("scipy module is not available.") 673 674 # Initial guess for minimisation. Solved by linear least squares. 675 x0 = E.estimate_x0_exp(times=E.times, values=E.values) 676 677 # Define function to minimise. Use errors as weights in the minimisation. 678 use_weights = True 679 680 if use_weights: 681 func = E.func_exp_weighted_residual 682 # If 'sigma'/erros describes one standard deviation errors of the input data points. The estimated covariance in 'pcov' is based on these values. 683 absolute_sigma = True 684 else: 685 func = E.func_exp_residual 686 absolute_sigma = False 687 688 # All args to function. Params are packed out through function, then other parameters. 689 args=(E.times, E.values, E.errors) 690 691 # Call scipy.optimize.leastsq. 692 popt, pcov, infodict, errmsg, ier = leastsq(func=func, x0=x0, args=args, full_output=True, ftol=E.ftol, xtol=E.xtol, maxfev=E.maxfev, factor=E.factor) 693 694 # Catch errors: 695 if ier in [1, 2, 3, 4]: 696 warning = None 697 elif ier in [9]: 698 warn(RelaxWarning("%s." % errmsg)) 699 warning = errmsg 700 elif ier in [0, 5, 6, 7, 8]: 701 raise RelaxError("scipy.optimize.leastsq raises following error: '%s'." % errmsg) 702 703 # 'nfev' number of function calls. 704 f_count = infodict['nfev'] 705 706 # 'fvec': function evaluated at the output. 707 fvec = infodict['fvec'] 708 #fvec_test = func(popt, times, values) 709 710 # 'fjac': A permutation of the R matrix of a QR factorization of the final approximate Jacobian matrix, stored column wise. Together with ipvt, the covariance of the estimate can be approximated. 711 # fjac = infodict['fjac'] 712 713 # 'qtf': The vector (transpose(q) * fvec). 714 # qtf = infodict['qtf'] 715 716 # 'ipvt' An integer array of length N which defines a permutation matrix, p, such that fjac*p = q*r, where r is upper triangular 717 # with diagonal elements of nonincreasing magnitude. Column j of p is column ipvt(j) of the identity matrix. 718 719 # 'pcov': The estimated covariance matrix of popt. 720 # The diagonals provide the variance of the parameter estimate. 721 722 # The reduced chi square: Take each "difference element, which could have been weighted" (O - E) and put to order 2. Sum them, and divide by number of degrees of freedom. 723 # Calculated the (weighted) chi2 value. 724 chi2 = sum( fvec**2 ) 725 726 # N is number of observations. 727 N = len(E.values) 728 # p is number of fitted parameters. 729 p = len(x0) 730 # n is number of degrees of freedom 731 #n = N - p - 1 732 n = N - p 733 734 # The reduced chi square. 735 chi2_red = chi2 / n 736 737 # chi2_red >> 1 : indicates a poor model fit. 738 # chi2_red > 1 : indicates that the fit has not fully captured the data (or that the error variance has been underestimated) 739 # chi2_red = 1 : indicates that the extent of the match between observations and estimates is in accord with the error variance. 740 # chi2_red < 1 : indicates that the model is 'over-fitting' the data: either the model is improperly fitting noise, or the error variance has been overestimated. 741 742 if pcov is None: 743 # indeterminate covariance 744 pcov = zeros((len(popt), len(popt)), dtype=float) 745 pcov.fill(inf) 746 elif not absolute_sigma: 747 if N > p: 748 pcov = pcov * chi2_red 749 else: 750 pcov.fill(inf) 751 752 # To compute one standard deviation errors on the parameters, take the square root of the diagonal covariance. 753 perr = sqrt(diag(pcov)) 754 755 # Return as standard from minfx. 756 param_vector = popt 757 param_vector_error = perr 758 759 iter_count = 0 760 g_count = 0 761 h_count = 0 762 763 # Pack to list. 764 results = [param_vector, param_vector_error, chi2, iter_count, f_count, g_count, h_count, warning] 765 766 # Return, including errors. 767 return results

768 769

770 -def minimise_minfx(E=None):

771 """Estimate r2eff and errors by minimising with minfx. 772 773 @keyword E: The Exponential function class, which contain data and functions. 774 @type E: class 775 @return: Packed list with optimised parameter, parameter error set to 'inf', chi2, iter_count, f_count, g_count, h_count, warning 776 @rtype: list 777 """ 778 779 # Check that the C modules have been compiled. 780 if not C_module_exp_fn: 781 raise RelaxError("Relaxation curve fitting is not available. Try compiling the C modules on your platform.") 782 783 # Initial guess for minimisation. Solved by linear least squares. 784 x0 = asarray( E.estimate_x0_exp(times=E.times, values=E.values) ) 785 786 if E.c_code: 787 # Minimise with C code. 788 789 # Initialise the function to minimise. 790 scaling_list = [1.0, 1.0] 791 model = Relax_fit_opt(model='exp', num_params=len(x0), values=E.values, errors=E.errors, relax_times=E.times, scaling_matrix=scaling_list) 792 793 # Define function to minimise for minfx. 794 t_func = model.func 795 t_dfunc = model.dfunc 796 t_d2func = model.d2func 797 args=() 798 799 else: 800 # Minimise with minfx. 801 # Define function to minimise for minfx. 802 t_func = E.func_exp_chi2 803 t_dfunc = E.func_exp_chi2_grad 804 t_d2func = None 805 # All args to function. Params are packed out through function, then other parameters. 806 args=(E.times, E.values, E.errors) 807 808 # Minimise. 809 results_minfx = generic_minimise(func=t_func, dfunc=t_dfunc, d2func=t_d2func, args=args, x0=x0, min_algor=E.min_algor, min_options=E.min_options, func_tol=E.func_tol, grad_tol=E.grad_tol, maxiter=E.max_iterations, A=E.A, b=E.b, full_output=True, print_flag=0) 810 811 # Unpack 812 param_vector, chi2, iter_count, f_count, g_count, h_count, warning = results_minfx 813 814 # Extract. 815 r2eff, i0 = param_vector 816 817 # Get the Jacobian. 818 if E.c_code == True: 819 if E.chi2_jacobian: 820 # Use the chi2 Jacobian from C. 821 jacobian_matrix_exp = transpose(asarray( model.jacobian_chi2(param_vector) ) ) 822 weights = ones(E.errors.shape) 823 824 else: 825 # Use the direct Jacobian from C. 826 jacobian_matrix_exp = transpose(asarray( model.jacobian(param_vector) ) ) 827 weights = 1. / E.errors**2 828 829 elif E.c_code == False: 830 if E.chi2_jacobian: 831 # Use the chi2 Jacobian from python. 832 jacobian_matrix_exp = E.func_exp_chi2_grad_array(params=param_vector, times=E.times, values=E.values, errors=E.errors) 833 weights = ones(E.errors.shape) 834 835 else: 836 # Use the direct Jacobian from python. 837 jacobian_matrix_exp = E.func_exp_grad(params=param_vector, times=E.times) 838 weights = 1. / E.errors**2 839 840 pcov = multifit_covar(J=jacobian_matrix_exp, weights=weights) 841 842 # To compute one standard deviation errors on the parameters, take the square root of the diagonal covariance. 843 param_vector_error = sqrt(diag(pcov)) 844 845 # Pack to list. 846 results = [param_vector, param_vector_error, chi2, iter_count, f_count, g_count, h_count, warning] 847 848 # Return, including errors. 849 return results

850

Source Code for Module specific_analyses.relax_disp.estimate_r2eff