Hi Troels,
Why not call this user function error_analysis.covariance_matrix? By
using the specific analysis API, you can easily generalise this. It
requires only a few steps, mainly just shifting your existing code
around:
1) Rename the user function.
2) Have a backend in the new module pipe_control.error_analysis (I
may then merge the pipe_control.monte_carlo module into it for a
general module for all error estimate techniques). A good example of
how to get the specific analysis API and use it is the Monte Carlo
error_analysis() function
(http://www.nmr-relax.com/api/3.2/pipe_control.monte_carlo-pysrc.html#error_analysis).
3) Add a new API method to specific_analyses.api_base which, like
most of the others, raises a RelaxImplementError (e.g.
covariance_matrix()). This allows each analysis type to implement it
as desired, and allows the user function to fail cleanly when used on
an analysis that doesn't support it.
4) Shift the covariance matrix calculation code of multifit_covar()
into the relax library, maybe in the lib.statistics module. Btw, the
docstring of this function needs a bit of work to allow the API
documentation to be compiled.
5) Create the dispersion specific covariance_matrix() method which
obtains the Jacobian and calls the relax library function. This would
just be the estimate_r2eff_err() function converted into a method.
This is how you add functionality to this API, and it would be a
useful learning exercise. If you know how to control this API, it
will be very powerful when you are adding new features to relax. In
this case, most of your code is unrelated to the dispersion analysis
and it is a perfect example of such an independent and powerful
feature. It would be a rather simple exercise, changing the new
system tests and just shifting the current functions/methods to fit
into this infrastructure.
Note that the feature of estimating the initial exponential curve
parameters by linearisation is also partly analysis independent as it
is of benefit for the relaxation curve-fitting analysis type as well.
So it could be spun out as a separate user function (maybe
exponential_curve.linear_estimate), use the specific API, and putting
the specific_analyses/relax_disp/estimate_r2eff.py
Exp.estimate_x0_exp() method into the relax library. This too would
be a simple code shifting exercise.
What do you think? The specific analysis API is one of the most
important design concepts in relax, and being able to control and
extend it will allow new ideas to be added much more quickly.
Regards,
Edward
On 27 August 2014 20:06, <tlinnet@xxxxxxxxxxxxx> wrote:
Author: tlinnet
Date: Wed Aug 27 20:06:23 2014
New Revision: 25347
URL: http://svn.gna.org/viewcvs/relax?rev=25347&view=rev
Log:
Added front-end to the new user function relax_disp.r2eff_err_estimate(),
which will estimate the R2eff errors
from a pipe and spins with optimised values of R2eff and i0.
The co-variance matrix can be calculated from the optimised parameters,
and the Jacobian.
Big care should be taken not to directly trust these results, since the
errors are quite different compared to the Monte-Carlo simulations.
This implementation, will reach the exact same error estimation as
scipy.optimize.leastsq.
But with much better control over the data, and insight into the
calculations.
task #7822(https://gna.org/task/index.php?7822): Implement user function
to estimate R2eff and associated errors for exponential curve fitting.
Modified:
trunk/user_functions/relax_disp.py
Modified: trunk/user_functions/relax_disp.py
URL:
http://svn.gna.org/viewcvs/relax/trunk/user_functions/relax_disp.py?rev=25347&r1=25346&r2=25347&view=diff
==============================================================================
--- trunk/user_functions/relax_disp.py (original)
+++ trunk/user_functions/relax_disp.py Wed Aug 27 20:06:23 2014
@@ -40,6 +40,7 @@
from pipe_control.mol_res_spin import get_spin_ids
from specific_analyses.relax_disp.catia import catia_execute, catia_input
from specific_analyses.relax_disp.cpmgfit import cpmgfit_execute,
cpmgfit_input
+from specific_analyses.relax_disp.estimate_r2eff import estimate_r2eff_err
from specific_analyses.relax_disp.data import cpmg_setup, insignificance,
plot_disp_curves, plot_exp_curves, r2eff_read, r2eff_read_spin,
relax_time, set_exp_type, r20_from_min_r2eff, spin_lock_field,
spin_lock_offset, write_disp_curves
from specific_analyses.relax_disp.data import INTERPOLATE_DISP,
INTERPOLATE_OFFSET, X_AXIS_DISP, X_AXIS_W_EFF, X_AXIS_THETA,
Y_AXIS_R2_R1RHO, Y_AXIS_R2_EFF
from specific_analyses.relax_disp.nessy import nessy_input
@@ -627,6 +628,50 @@
uf.gui_icon = "oxygen.status.object-locked"
uf.wizard_size = (800, 500)
uf.wizard_image = ANALYSIS_IMAGE_PATH + 'relax_disp_200x200.png'
+
+
+# The relax_disp.r2eff_err_estimate user function.
+uf = uf_info.add_uf('relax_disp.r2eff_err_estimate')
+uf.title = "Estimate R2eff errors by the Jacobian matrix."
+uf.title_short = "Estimate R2eff errors."
+uf.add_keyarg(
+ name = "spin_id",
+ py_type = "str",
+ arg_type = "spin ID",
+ desc_short = "spin ID to restrict value setting to",
+ desc = "The spin ID string to restrict value setting to.",
+ can_be_none = True
+)
+uf.add_keyarg(
+ name = "epsrel",
+ py_type = "float",
+ default = 0.0,
+ desc_short = "parameter to remove linear-dependent columns.",
+ desc = "The parameter to remove linear-dependent columns when J is
rank deficient.",
+ can_be_none = False
+)
+uf.add_keyarg(
+ name = "verbosity",
+ default = 1,
+ py_type = "int",
+ desc_short = "amount of information to print.",
+ desc = "The higher the value, the greater the verbosity.",
+ can_be_none = False
+)
+# Description.
+uf.desc.append(Desc_container())
+uf.desc[-1].add_paragraph("This is a new experimental feature from
version 3.3, and should only be tried out with big care.")
+uf.desc[-1].add_paragraph("This will estimate R2eff errors by using the
exponential decay Jacobian matrix 'J' to compute the covariance matrix of
the best-fit parameters.")
+uf.desc[-1].add_paragraph("This can be an huge time saving step, when
performing model fitting in R1rho. Errors of R2eff values, are normally
estimated by time-consuming Monte-Carlo simulations.")
+uf.desc[-1].add_paragraph("This method is inspired from the GNU
Scientific Library (GSL).")
+uf.desc[-1].add_paragraph("The covariance matrix is given by: covar = Qxx
= (J^T J)^{-1}.")
+uf.desc[-1].add_paragraph("Qxx is computed by QR decomposition, Qxx=QR,
Qxx^-1=R^-1 Q^T. The columns of R which satisfy: |R_{kk}| <= epsrel
|R_{11}| are considered linearly-dependent and are excluded from the
covariance matrix (the corresponding rows and columns of the covariance
matrix are set to zero).")
+uf.desc[-1].add_paragraph("The parameter 'epsrel' is used to remove
linear-dependent columns when J is rank deficient.")
+uf.backend = estimate_r2eff_err
+uf.menu_text = "&r2eff_err_estimate"
+uf.gui_icon = "relax.relax_fit"
+uf.wizard_size = (800, 800)
+uf.wizard_image = ANALYSIS_IMAGE_PATH + sep + 'blank_150x150.png'
# The relax_disp.r2eff_read user function.
_______________________________________________
relax (http://www.nmr-relax.com)
This is the relax-commits mailing list
relax-commits@xxxxxxx
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
relax-devel@xxxxxxx
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
Re: r25347 - /trunk/user_functions/relax_disp.py