Author: bugman
Date: Fri Sep 12 14:14:11 2014
New Revision: 25787
URL: http://svn.gna.org/viewcvs/relax?rev=25787&view=rev
Log:
Shifted the contents of the specific_analysis.relax_fit.estimate_rx_err
module into the API.
The estimate_rx_err() function is now the covariance_matrix() method of the
specific API. The code
for calculating the covariance matrix and errors are now in the function
pipe_control.error_analysis.covariance_matrix(), so this has been removed.
And the error setting is
performed by the set_errors() API method, so that code has been deleted as
well.
Removed:
trunk/specific_analyses/relax_fit/estimate_rx_err.py
Modified:
trunk/specific_analyses/relax_fit/api.py
Modified: trunk/specific_analyses/relax_fit/api.py
URL:
http://svn.gna.org/viewcvs/relax/trunk/specific_analyses/relax_fit/api.py?rev=25787&r1=25786&r2=25787&view=diff
==============================================================================
--- trunk/specific_analyses/relax_fit/api.py (original)
+++ trunk/specific_analyses/relax_fit/api.py Fri Sep 12 14:14:11 2014
@@ -1,6 +1,7 @@
###############################################################################
#
#
# Copyright (C) 2004-2014 Edward d'Auvergne
#
+# Copyright (C) 2014 Troels E. Linnet
#
#
#
# This file is part of the program relax (http://www.nmr-relax.com).
#
#
#
@@ -25,14 +26,16 @@
# Python module imports.
from minfx.generic import generic_minimise
from minfx.grid import grid
-from numpy import dot, float64, zeros
+from numpy import asarray, dot, float64, transpose, zeros
from numpy.linalg import inv
from re import match, search
+import sys
from warnings import warn
# relax module imports.
from dep_check import C_module_exp_fn
from lib.errors import RelaxError, RelaxNoModelError, RelaxNoSequenceError
+from lib.text.sectioning import subsection
from lib.warnings import RelaxDeselectWarning
from pipe_control.mol_res_spin import exists_mol_res_spin_data, return_spin,
spin_loop
from specific_analyses.api_base import API_base
@@ -43,7 +46,7 @@
# C modules.
if C_module_exp_fn:
- from target_functions.relax_fit import setup
+ from target_functions.relax_fit import jacobian, setup
class Relax_fit(API_base, API_common):
@@ -70,6 +73,82 @@
# Place a copy of the parameter list object in the instance
namespace.
self._PARAMS = Relax_fit_params()
+
+
+ def covariance_matrix(self, model_info=None, verbosity=1):
+ """Return the Jacobian and weights required for parameter errors via
the covariance matrix.
+
+ @keyword model_info: The spin container and the spin ID string
from the _model_loop_spin() method.
+ @type model_info: SpinContainer instance, str
+ @keyword verbosity: The amount of information to print. The
higher the value, the greater the verbosity.
+ @type verbosity: int
+ @return: The Jacobian and weight matrices for the
given model.
+ @rtype: numpy rank-2 array, numpy rank-2 array
+ """
+
+ # Unpack the data.
+ spin, spin_id = model_info
+
+ # Check that the C modules have been compiled.
+ if not C_module_exp_fn:
+ raise RelaxError("Relaxation curve fitting is not available.
Try compiling the C modules on your platform.")
+
+ # Perform checks.
+ if not cdp.curve_type == 'exp':
+ raise RelaxError("Only curve type of 'exp' is allowed for error
estimation. Set by: relax_fit.select_model('exp').")
+
+ # Raise Error, if not optimised.
+ if not (hasattr(spin, 'rx') and hasattr(spin, 'i0')):
+ raise RelaxError("Spin '%s' does not contain optimised 'rx' and
'i0' values. Try execute: minimise.execute(min_algor='Newton',
constraints=False)"%(spin_id))
+
+ # Raise warning, if gradient count is 0. This could point to a lack
of minimisation first.
+ if hasattr(spin, 'g_count'):
+ if getattr(spin, 'g_count') == 0.0:
+ text = "Spin %s contains a gradient count of 0.0. Is the rx
parameter optimised? Try execute: minimise.execute(min_algor='Newton',
constraints=False)" %(spin_id)
+ warn(RelaxWarning("%s." % text))
+
+ # Print information.
+ if verbosity >= 1:
+ # Individual spin block section.
+ top = 2
+ if verbosity >= 2:
+ top += 2
+ subsection(file=sys.stdout, text="Estimating rx error for spin:
%s"%spin_id, prespace=top)
+
+ # The keys.
+ keys = list(spin.peak_intensity.keys())
+
+ # The peak intensities and times.
+ values = []
+ errors = []
+ times = []
+ for key in keys:
+ values.append(spin.peak_intensity[key])
+ errors.append(spin.peak_intensity_err[key])
+ times.append(cdp.relax_times[key])
+
+ # Convert to numpy array.
+ values = asarray(values)
+ errors = asarray(errors)
+ times = asarray(times)
+
+ # Extract values.
+ rx = getattr(spin, 'rx')
+ i0 = getattr(spin, 'i0')
+
+ # Pack data
+ param_vector = [rx, i0]
+
+ # Initialise data in C code.
+ scaling_list = [1.0, 1.0]
+ setup(num_params=len(param_vector), num_times=len(times),
values=values, sd=errors, relax_times=times, scaling_matrix=scaling_list)
+
+ # Use the direct Jacobian from function.
+ jacobian_matrix_exp = transpose(asarray( jacobian(param_vector) ) )
+ weights = 1. / errors**2
+
+ # Return the matrices.
+ return jacobian_matrix_exp, weights
def create_mc_data(self, data_id=None):
Removed: trunk/specific_analyses/relax_fit/estimate_rx_err.py
URL:
http://svn.gna.org/viewcvs/relax/trunk/specific_analyses/relax_fit/estimate_rx_err.py?rev=25786&view=auto
==============================================================================
--- trunk/specific_analyses/relax_fit/estimate_rx_err.py (original)
+++ trunk/specific_analyses/relax_fit/estimate_rx_err.py (removed)
@@ -1,157 +0,0 @@
-###############################################################################
-#
#
-# Copyright (C) 2014 Troels E. Linnet
#
-#
#
-# This file is part of the program relax (http://www.nmr-relax.com).
#
-#
#
-# This program is free software: you can redistribute it and/or modify
#
-# it under the terms of the GNU General Public License as published by
#
-# the Free Software Foundation, either version 3 of the License, or
#
-# (at your option) any later version.
#
-#
#
-# This program is distributed in the hope that it will be useful,
#
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
#
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
#
-# GNU General Public License for more details.
#
-#
#
-# You should have received a copy of the GNU General Public License
#
-# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
-#
#
-###############################################################################
-
-# Module docstring.
-"""Estimation of rx error, from Co-variance matrix."""
-
-# Python module imports.
-from copy import deepcopy
-from numpy import asarray, diag, sqrt, transpose
-import sys
-from warnings import warn
-
-# relax module imports.
-from dep_check import C_module_exp_fn
-from lib.errors import RelaxError
-from lib.statistics import multifit_covar
-from lib.text.sectioning import subsection
-from lib.warnings import RelaxWarning
-from pipe_control.mol_res_spin import generate_spin_string, spin_loop
-
-# C modules.
-if C_module_exp_fn:
- from target_functions.relax_fit import jacobian, jacobian_chi2, setup
-
-
-def estimate_rx_err(spin_id=None, epsrel=0.0, verbosity=2):
- """This will estimate the rx and i0 errors from the covariance matrix
Qxx. Qxx is calculated from the Jacobian matrix and the optimised parameters.
-
- @keyword spin_id: The spin identification string.
- @type spin_id: str
- @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.
- @type epsrel: float
- @keyword verbosity: The amount of information to print. The higher
the value, the greater the verbosity.
- @type verbosity: int
- """
-
- # Check that the C modules have been compiled.
- if not C_module_exp_fn:
- raise RelaxError("Relaxation curve fitting is not available. Try
compiling the C modules on your platform.")
-
- # Perform checks.
- if not cdp.curve_type == 'exp':
- raise RelaxError("Only curve type of 'exp' is allowed for error
estimation. Set by: relax_fit.select_model('exp').")
-
- # Loop over the spins.
- 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):
- # Generate spin string.
- spin_string = generate_spin_string(spin=cur_spin, mol_name=mol_name,
res_num=resi, res_name=resn)
-
- # Raise Error, if not optimised.
- if not (hasattr(cur_spin, 'rx') and hasattr(cur_spin, 'i0')):
- raise RelaxError("Spin '%s' does not contain optimised 'rx' and
'i0' values. Try execute: minimise.execute(min_algor='Newton',
constraints=False)"%(spin_string))
-
- # Raise warning, if gradient count is 0. This could point to a lack
of minimisation first.
- if hasattr(cur_spin, 'g_count'):
- if getattr(cur_spin, 'g_count') == 0.0:
- text = "Spin %s contains a gradient count of 0.0. Is the rx
parameter optimised? Try execute: minimise.execute(min_algor='Newton',
constraints=False)" %(spin_string)
- warn(RelaxWarning("%s." % text))
-
- # Print information.
- if verbosity >= 1:
- # Individual spin block section.
- top = 2
- if verbosity >= 2:
- top += 2
- subsection(file=sys.stdout, text="Estimating rx error for spin:
%s"%spin_string, prespace=top)
-
- # The keys.
- keys = list(cur_spin.peak_intensity.keys())
-
- # The peak intensities and times.
- values = []
- errors = []
- times = []
- for key in keys:
- values.append(cur_spin.peak_intensity[key])
- errors.append(cur_spin.peak_intensity_err[key])
- times.append(cdp.relax_times[key])
-
- # Convert to numpy array.
- values = asarray(values)
- errors = asarray(errors)
- times = asarray(times)
-
- # Extract values.
- rx = getattr(cur_spin, 'rx')
- i0 = getattr(cur_spin, 'i0')
-
- # Pack data
- param_vector = [rx, i0]
-
- # Initialise data in C code.
- scaling_list = [1.0, 1.0]
- setup(num_params=len(param_vector), num_times=len(times),
values=values, sd=errors, relax_times=times, scaling_matrix=scaling_list)
-
- # Use the direct Jacobian from function.
- jacobian_matrix_exp = transpose(asarray( jacobian(param_vector) ) )
- weights = 1. / errors**2
-
- # Get the co-variance
- pcov = multifit_covar(J=jacobian_matrix_exp, weights=weights)
-
- # To compute one standard deviation errors on the parameters, take
the square root of the diagonal covariance.
- param_vector_error = sqrt(diag(pcov))
-
- # Extract values.
- rx_err, i0_err = param_vector_error
-
- # Copy rx, to rx_err, if not exists.
- if not hasattr(cur_spin, 'rx_err'):
- setattr(cur_spin, 'rx_err', deepcopy(getattr(cur_spin, 'rx')))
- if not hasattr(cur_spin, 'i0_err'):
- setattr(cur_spin, 'i0_err', deepcopy(getattr(cur_spin, 'i0')))
-
- # Set error.
- cur_spin.rx_err = rx_err
- cur_spin.i0_err = i0_err
-
- # Get other relevant information.
- chi2 = getattr(cur_spin, 'chi2')
-
- # Print information.
- print_strings = []
- if verbosity >= 1:
- # Add print strings.
- point_info = "Spin: '%s', with %i time points." % (spin_string,
len(times))
- print_strings.append(point_info)
-
- par_info = "rx=%3.3f rx_err=%3.4f, i0=%6.1f, i0_err=%3.4f,
chi2=%3.3f.\n" % ( rx, rx_err, i0, i0_err, chi2)
- print_strings.append(par_info)
-
- if verbosity >= 2:
- time_info = ', '.join(map(str, times))
- print_strings.append('For time array: '+time_info+'.\n\n')
-
- # Print info
- if len(print_strings) > 0:
- for print_string in print_strings:
- print(print_string),
r25787 - in /trunk/specific_analyses/relax_fit: api.py estimate_rx_err.py