Re: r25513 - /branches/est_par_error/lib/dispersion/tsmfk01.py -- September 02, 2014

Hi Edward.

I currently do not aim to implement minimisation with BFGS or similar.
This is only for estimating the error.

But I think I have it.

dw is returned in rad/s, and I need to scale i back to ppm.

Best
Troels

2014-09-02 9:24 GMT+02:00 Edward d'Auvergne <edward@xxxxxxxxxxxxx>:

Hi Troels,

If you are interested in speed for these functions - something which
may not be worth your effort right now - you should minimise
repetitive calculations.  For example, you have:

  dw*tcp = 5 times,
  k_AB / dw = 2 times,
  sin(dw*tcp) = 2 times.

If you pre-calculate these in the calling function and then add them
as new arguments, the code will not look as clean but it will be much
faster.  This is done in the model-free code, but instead of arguments
I have a special data container holding all the structures which is
then passed into each function, gradient, and Hessian function.

Regards,

Edward



On 1 September 2014 20:51,  <tlinnet@xxxxxxxxxxxxx> wrote:

Author: tlinnet
Date: Mon Sep  1 20:51:45 2014
New Revision: 25513

URL: http://svn.gna.org/viewcvs/relax?rev=25513&view=rev
Log:
Added partial derivatives for model TSMFK01 and the Jacobian function.

task #7824(https://gna.org/task/index.php?7824): Model parameter ERROR 
estimation from Jacobian and Co-variance matrix of dispersion models.

Modified:
    branches/est_par_error/lib/dispersion/tsmfk01.py

Modified: branches/est_par_error/lib/dispersion/tsmfk01.py
URL: 
http://svn.gna.org/viewcvs/relax/branches/est_par_error/lib/dispersion/tsmfk01.py?rev=25513&r1=25512&r2=25513&view=diff
==============================================================================
--- branches/est_par_error/lib/dispersion/tsmfk01.py    (original)
+++ branches/est_par_error/lib/dispersion/tsmfk01.py    Mon Sep  1 
20:51:45 2014
@@ -66,7 +66,7 @@
 """

 # Python module imports.
-from numpy import fabs, min, sin, isfinite, sum
+from numpy import array, cos, fabs, min, sin, isfinite, ones, sum, 
transpose
 from numpy.ma import fix_invalid, masked_where


@@ -129,3 +129,72 @@
     if not isfinite(sum(back_calc)):
         # Replaces nan, inf, etc. with fill value.
         fix_invalid(back_calc, copy=False, fill_value=1e100)
+
+
+def r2eff_TSMFK01_jacobian(r20a=None, dw=None, k_AB=None, tcp=None):
+    """The Jacobian matrix of TSMFK01.
+
+    @keyword r20a:          The R20 parameter value of state A (R2 with 
no exchange).
+    @type r20a:             numpy float array of rank [NE][NS][NM][NO][ND]
+    @keyword dw:            The chemical exchange difference between 
states A and B in rad/s.
+    @type dw:               numpy float array of rank [NE][NS][NM][NO][ND]
+    @keyword k_AB:          The k_AB parameter value (the forward 
exchange rate in rad/s).
+    @type k_AB:             float
+    @keyword tcp:           The tau_CPMG times (1 / 4.nu1).
+    @type tcp:              numpy float array of rank [NE][NS][NM][NO][ND]
+    """
+
+    # Get the partial derivatives.
+    get_d_f_d_r20a = d_f_d_r20a(r20a=r20a, dw=dw, k_AB=k_AB, tcp=tcp)
+    get_d_f_d_dw = d_f_d_dw(r20a=r20a, dw=dw, k_AB=k_AB, tcp=tcp)
+    get_d_f_d_k_AB = d_f_d_k_AB(r20a=r20a, dw=dw, k_AB=k_AB, tcp=tcp)
+
+    return transpose(array( [get_d_f_d_r20a , get_d_f_d_dw, 
get_d_f_d_k_AB] ) )
+
+
+def d_f_d_r20a(r20a=None, dw=None, k_AB=None, tcp=None):
+    """Partial derivative with respect to r20a.
+
+    @keyword r20a:          The R20 parameter value of state A (R2 with 
no exchange).
+    @type r20a:             numpy float array of rank [NE][NS][NM][NO][ND]
+    @keyword dw:            The chemical exchange difference between 
states A and B in rad/s.
+    @type dw:               numpy float array of rank [NE][NS][NM][NO][ND]
+    @keyword k_AB:          The k_AB parameter value (the forward 
exchange rate in rad/s).
+    @type k_AB:             float
+    @keyword tcp:           The tau_CPMG times (1 / 4.nu1).
+    @type tcp:              numpy float array of rank [NE][NS][NM][NO][ND]
+    """
+
+    return ones(dw.shape)
+
+
+def d_f_d_dw(r20a=None, dw=None, k_AB=None, tcp=None):
+    """Partial derivative with respect to dw.
+
+    @keyword r20a:          The R20 parameter value of state A (R2 with 
no exchange).
+    @type r20a:             numpy float array of rank [NE][NS][NM][NO][ND]
+    @keyword dw:            The chemical exchange difference between 
states A and B in rad/s.
+    @type dw:               numpy float array of rank [NE][NS][NM][NO][ND]
+    @keyword k_AB:          The k_AB parameter value (the forward 
exchange rate in rad/s).
+    @type k_AB:             float
+    @keyword tcp:           The tau_CPMG times (1 / 4.nu1).
+    @type tcp:              numpy float array of rank [NE][NS][NM][NO][ND]
+    """
+
+    return -k_AB * cos( dw * tcp) / dw + k_AB * sin(dw * tcp) / ( dw**2 * 
tcp)
+
+
+def d_f_d_k_AB(r20a=None, dw=None, k_AB=None, tcp=None):
+    """Partial derivative with respect to k_AB.
+
+    @keyword r20a:          The R20 parameter value of state A (R2 with 
no exchange).
+    @type r20a:             numpy float array of rank [NE][NS][NM][NO][ND]
+    @keyword dw:            The chemical exchange difference between 
states A and B in rad/s.
+    @type dw:               numpy float array of rank [NE][NS][NM][NO][ND]
+    @keyword k_AB:          The k_AB parameter value (the forward 
exchange rate in rad/s).
+    @type k_AB:             float
+    @keyword tcp:           The tau_CPMG times (1 / 4.nu1).
+    @type tcp:              numpy float array of rank [NE][NS][NM][NO][ND]
+    """
+
+    return 1. - sin( dw * tcp) / (dw * tcp)


_______________________________________________
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: r25513 - /branches/est_par_error/lib/dispersion/tsmfk01.py

Header

Content

Related Messages