Hi, Ok, then I'll do it! I cannot support multiple covariance matrix implementations for each analysis type. The subsequent support work, which is always required especially when a code refactorisation occurs or when users find cases where the method doesn't work, is something I'll probably have to look after. Therefore I will do this now before you get too far and make a lot more work for me in the future. Multiple implementations of the same thing are difficult to support. I'm guessing it'll take me all of about 15 min to rearrange the code and create the new pipe_control.error_analysis.covariance_matrix() function, so it is still not too late. Fortunately most of the work happens in lib.statistics.multifit_covar() which does not require any changes.
Regards, Edward On 12 September 2014 12:15, Troels Emtekær Linnet <[email protected]> wrote: > 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
