Hi Troels, It's now done:
r25781: http://article.gmane.org/gmane.science.nmr.relax.scm/23532 r25782: http://article.gmane.org/gmane.science.nmr.relax.scm/23533 r25783: http://article.gmane.org/gmane.science.nmr.relax.scm/23534 r25784: http://article.gmane.org/gmane.science.nmr.relax.scm/23535 r25785: http://article.gmane.org/gmane.science.nmr.relax.scm/23536 r25786: http://article.gmane.org/gmane.science.nmr.relax.scm/23537 r25787: http://article.gmane.org/gmane.science.nmr.relax.scm/23538 r25788: http://article.gmane.org/gmane.science.nmr.relax.scm/23539 r25789: http://article.gmane.org/gmane.science.nmr.relax.scm/23540 It took a little longer as there was a bug I had to sort out in the 1 new function (pipe_control.error_analysis.covariance_matrix()). You will see that these changes mainly involve rearranging your code. I have also deleted all of your error setting code as this now uses the set_error() API method to avoid code duplication. The covariance matrix and error calculation code has shifted into the pipe_control.error_analysis.covariance_matrix() function. And your backend in the specific_analysis.relax_fit.estimate_rx_err module is now the API method covariance_matrix() in specific_analyses.relax_fit.api. If you run: $ ./relax -s Relax_fit.test_curve_fitting_height_estimate_error you will see that everything is functional, just as you had it. Regards, Edward On 12 September 2014 13:46, Edward d'Auvergne <[email protected]> wrote: > 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
