Hi Edward. Thank you for the detailed description.
I currently don't have time to look into api implementation. Best Troels 2014-09-12 11:53 GMT+02:00 Edward d'Auvergne <[email protected]>: > Hi Troels, > > For the backend function, instead of duplicating code and having one > backend in the relaxation dispersion analysis and one backend in the > relaxation curve-fitting analysis - now is the time to learn the > specific analysis API! I cannot emphasis how simple this is. The > steps would be: > > 1) Create the error_analysis.covariance_matrix user function frontend > in user_functions/. This would be instead of creating a new relax_fit > user function, and would be pretty much identical anyway. > > 2a) I have renamed the pipe_control.monte_carlo module to > pipe_control.error_analysis. I have prepended the text 'monte_carlo_' > to all of the current functions. This module will be used for all > error analysis techniques: Monte Carlo simulations, covariance > matrix, Jackknife simulations, Bootstrapping (which is currently via > the Monte Carlo functions), etc. > > 2b) Create a new function in pipe_control.error_analysis called > covariance_matrix(). > > 3) Inside covariance_matrix(), call "api = return_api()", and then > call "api.covariance_matrix()". For the 'relax-fit' analysis, this > will be the covariance_matrix() method in the > specific_analyses.relax_fit.api module. > > 4) Add a stub covariance_matrix() method to > specific_analyses.api_base. This defines how the method should be > named and what arguments it takes. > > 5) Add a new covariance_matrix() method to specific_analyses.relax_fit.api. > > That's all there is too it. It is really time for this! The entire > point of the specific analysis API is to avoid what you are currently > implementing in relax. If the pattern of what you are doing is > extended to all analysis types, which might happen in the future, then > there will be huge amounts of code duplication required. This is a > bad design. You need to dive into the specific API! The steps above > are all that are required for this. Oh, a system test with a call to > the error_analysis.covariance_matrix user function would be the best > way to start this. > > Regards, > > Edward > > On 12 September 2014 11:25, <[email protected]> wrote: >> Author: tlinnet >> Date: Fri Sep 12 11:25:50 2014 >> New Revision: 25775 >> >> URL: http://svn.gna.org/viewcvs/relax?rev=25775&view=rev >> Log: >> Implemented back-end function to estimate Rx and I0 errors from Jacobian >> matrix. >> >> This is to prepare for user funcion in relax_fit, to estimate errors. >> >> Added: >> trunk/specific_analyses/relax_fit/estimate_rx_err.py >> Modified: >> trunk/specific_analyses/relax_fit/__init__.py >> >> Modified: trunk/specific_analyses/relax_fit/__init__.py >> URL: >> http://svn.gna.org/viewcvs/relax/trunk/specific_analyses/relax_fit/__init__.py?rev=25775&r1=25774&r2=25775&view=diff >> ============================================================================== >> --- trunk/specific_analyses/relax_fit/__init__.py (original) >> +++ trunk/specific_analyses/relax_fit/__init__.py Fri Sep 12 11:25:50 >> 2014 >> @@ -25,6 +25,7 @@ >> # The available modules. >> __all__ = [ >> 'api', >> + 'estimate_rx_err', >> 'optimisation', >> 'parameter_object', >> 'parameters', >> >> Added: 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=25775&view=auto >> ============================================================================== >> --- trunk/specific_analyses/relax_fit/estimate_rx_err.py (added) >> +++ trunk/specific_analyses/relax_fit/estimate_rx_err.py Fri Sep 12 >> 11:25:50 2014 >> @@ -0,0 +1,157 @@ >> +############################################################################### >> +# >> # >> +# 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), >> >> >> _______________________________________________ >> relax (http://www.nmr-relax.com) >> >> This is the relax-commits mailing list >> [email protected] >> >> To unsubscribe from this list, get a password >> reminder, or change your subscription options, >> visit the list information page at >> https://mail.gna.org/listinfo/relax-commits > > _______________________________________________ > relax (http://www.nmr-relax.com) > > This is the relax-devel mailing list > [email protected] > > To unsubscribe from this list, get a password > reminder, or change your subscription options, > visit the list information page at > https://mail.gna.org/listinfo/relax-devel _______________________________________________ relax (http://www.nmr-relax.com) This is the relax-devel mailing list [email protected] To unsubscribe from this list, get a password reminder, or change your subscription options, visit the list information page at https://mail.gna.org/listinfo/relax-devel
