specific_analyses.relax_disp.estimate

1 ############################################################################### 2 # # 3 # Copyright (C) 2014 Troels E. Linnet # 4 # Copyright (C) 2014 Edward d'Auvergne # 5 # # 6 # This file is part of the program relax (http://www.nmr-relax.com). # 7 # # 8 # This program is free software: you can redistribute it and/or modify # 9 # it under the terms of the GNU General Public License as published by # 10 # the Free Software Foundation, either version 3 of the License, or # 11 # (at your option) any later version. # 12 # # 13 # This program is distributed in the hope that it will be useful, # 14 # but WITHOUT ANY WARRANTY; without even the implied warranty of # 15 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # 16 # GNU General Public License for more details. # 17 # # 18 # You should have received a copy of the GNU General Public License # 19 # along with this program. If not, see <http://www.gnu.org/licenses/>. # 20 # # 21 ############################################################################### 22 23 # Module docstring. 24 """Target functions for relaxation exponential curve fitting with both minfx and scipy.optimize.leastsq.""" 25 26 # Python module imports. 27 from copy import deepcopy 28 from numpy import array, asarray, diag, dot, exp, eye, log, multiply, ones, sqrt, sum, transpose, zeros 29 from minfx.generic import generic_minimise 30 import sys 31 from warnings import warn 32 33 # relax module imports. 34 from dep_check import C_module_exp_fn, scipy_module 35 from lib.dispersion.variables import MODEL_R2EFF 36 from lib.errors import RelaxError 37 from lib.statistics import multifit_covar 38 from lib.text.sectioning import section, subsection 39 from lib.warnings import RelaxWarning 40 from pipe_control.mol_res_spin import generate_spin_string, spin_loop 41 from specific_analyses.relax_disp.checks import check_model_type 42 from specific_analyses.relax_disp.data import average_intensity, loop_exp_frq_offset_point, loop_time, return_param_key_from_data 43 from specific_analyses.relax_disp.parameters import disassemble_param_vector 44 from target_functions.chi2 import chi2_rankN, dchi2 45 from target_functions.relax_fit_wrapper import Relax_fit_opt 46 47 # Scipy installed. 48 if scipy_module: 49 # Import leastsq. 50 from scipy.optimize import leastsq 51 52

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

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

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

172 -class Exp:

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

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

188 189

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

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

205 206

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

208 """Setup options to scipy.optimize.leastsq. 209 210 @keyword ftol: The function tolerance for the relative error desired in the sum of squares, parsed to leastsq. 211 @type ftol: float 212 @keyword xtol: The error tolerance for the relative error desired in the approximate solution, parsed to leastsq. 213 @type xtol: float 214 @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]. 215 @type maxfev: int 216 @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). 217 @type factor: float 218 """ 219 220 # Store settings. 221 self.ftol = ftol 222 self.xtol = xtol 223 self.maxfev = maxfev 224 self.factor = factor

225 226

227 - 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):

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

282 283

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

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

310 311

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

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

328 329

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

331 """Calculate the residual vector betwen measured values and the function values. 332 333 @param params: The vector of parameter values. 334 @type params: numpy rank-1 float array 335 @keyword times: The time points. 336 @type times: numpy array 337 @param values: The measured values. 338 @type values: numpy array 339 @return: The residuals. 340 @rtype: numpy array 341 """ 342 343 # 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. 344 K = values - self.func_exp(params=params, times=times) 345 346 # Return 347 return K

348 349

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

351 """Calculate the weighted residual vector betwen measured values and the function values. 352 353 @param params: The vector of parameter values. 354 @type params: numpy rank-1 float array 355 @keyword times: The time points. 356 @type times: numpy array 357 @param values: The measured values. 358 @type values: numpy array 359 @param errors: The standard deviation of the measured intensity values per time point. 360 @type errors: numpy array 361 @return: The weighted residuals. 362 @rtype: numpy array 363 """ 364 365 # 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. 366 Kw = 1. / errors * self.func_exp_residual(params=params, times=times, values=values) 367 368 # Return 369 return Kw

370 371

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

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

398 399

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

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

423 424

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

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

460 461

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

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

491 492

493 -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):

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

660 661

662 -def minimise_leastsq(E=None):

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

769 770

771 -def minimise_minfx(E=None):

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

851

Source Code for Module specific_analyses.relax_disp.estimate_r2eff