Hi,
Such equations are common after asking for a partial derivative,
however there are a huge number of simplifications possible. The
essential key is to use the chain rule yourself before any symbolic
calculations. For example in:
R2eff = 1/2 (R20A + R20B + kex - 2nu_cpmg * acosh(Dpos * cosh(etapos)
- Dneg * cos(etaneg)))
we can say that Dpos(theta) is a function of the model parameters
theta (and Dneg(theta), etapos(theta), and etaneg(theta) as well).
Then instead of expanding this fully we can construct, for example,
the dR2eff_dkex() function which is a function of the Dpos(theta) and
dDpos_dkex(theta). Then you can handle the Dpos equation in isolation
and work out its partial derivatives for all parameters. This will
result in a series of partial derivatives for Dpos, Dneg, etapos, and
etaneg for all parameters. So it requires a lot more code in
lib.dispersion.cr72, but repetitive calculations are easier to
identify and it is much faster and it is simpler to handle. The chain
rule can also be used for eta(theta), Psi(theta) and zeta(theta).
Then the partial derivative equations will have the same R2eff, D+/-,
eta+/-, Psi, zeta equation layout as the CR72 section of the manual
(http://www.nmr-relax.com/manual/full_CR72_2_site_CPMG_model.html).
You just have one partial derivative section per parameter.
I used Maxima and Mathematica for the simplifications to the
model-free analysis gradients and Hessians. But that was after using
the chain rule to break the equations into unique parts which are now
the different sections in the "Optimisation of relaxation data -
values, gradients, and Hessians" chapter of the manual
(http://www.nmr-relax.com/manual/Optimisation_relaxation_data_values_gradients_Hessians.html).
The key there was that I identified patterns in the equations and
performed part of the simplification myself, and then let Maxima
simplify only parts of the full equation. This works much better.
Anyway, these derivations required many, many pages of my notebook to
complete. Imagine the "Optimisation of relaxation data - values,
gradients, and Hessians" chapter of the manual expanded ~20 times for
all the derivations.
Regards,
Edward
On 2 September 2014 20:06, <tlinnet@xxxxxxxxxxxxx> wrote:
Author: tlinnet
Date: Tue Sep 2 20:06:13 2014
New Revision: 25567
URL: http://svn.gna.org/viewcvs/relax?rev=25567&view=rev
Log:
Added derivation of CR72 simple, where R20A = R20B.
This is some very very very long eqations, but they will be used as a test
case.
task #7824(https://gna.org/task/index.php?7824): Model parameter ERROR
estimation from Jacobian and Co-variance matrix of dispersion models.
Added:
branches/est_par_error/test_suite/shared_data/dispersion/estimate_par_err/cr72/cr72_simple_der.log
branches/est_par_error/test_suite/shared_data/dispersion/estimate_par_err/cr72/cr72_simple_der.py
Added:
branches/est_par_error/test_suite/shared_data/dispersion/estimate_par_err/cr72/cr72_simple_der.log
URL:
http://svn.gna.org/viewcvs/relax/branches/est_par_error/test_suite/shared_data/dispersion/estimate_par_err/cr72/cr72_simple_der.log?rev=25567&view=auto
==============================================================================
---
branches/est_par_error/test_suite/shared_data/dispersion/estimate_par_err/cr72/cr72_simple_der.log
(added)
+++
branches/est_par_error/test_suite/shared_data/dispersion/estimate_par_err/cr72/cr72_simple_der.log
Tue Sep 2 20:06:13 2014
@@ -0,0 +1,12 @@
+Form the Jacobian matrix by:
+------------------------------------------------------------------------------
+
+d_f_d_r20a = 1.00000000000000
+
+d_f_d_pA = -cpmg_frqs*(8.0*dw**2*kex*(0.5*dw**2 + 0.5*kex**2)*(kex*pA -
kex*(-pA + 1.0))*cos(eta_scale*sqrt(dw**2 - kex**2 + sqrt(4.0*dw**2*(kex*pA
- kex*(-pA + 1.0))**2 + (-dw**2 +
kex**2)**2))/cpmg_frqs)/(4.0*dw**2*(kex*pA - kex*(-pA + 1.0))**2 + (-dw**2
+ kex**2)**2)**(3/2) - 8.0*dw**2*kex*(0.5*dw**2 + 0.5*kex**2)*(kex*pA -
kex*(-pA + 1.0))*cosh(eta_scale*sqrt(-dw**2 + kex**2 +
sqrt(4.0*dw**2*(kex*pA - kex*(-pA + 1.0))**2 + (-dw**2 +
kex**2)**2))/cpmg_frqs)/(4.0*dw**2*(kex*pA - kex*(-pA + 1.0))**2 + (-dw**2
+ kex**2)**2)**(3/2) + 4.0*dw**2*eta_scale*kex*(kex*pA - kex*(-pA +
1.0))*((0.5*dw**2 + 0.5*kex**2)/sqrt(4.0*dw**2*(kex*pA - kex*(-pA +
1.0))**2 + (-dw**2 + kex**2)**2) - 0.5)*sin(eta_scale*sqrt(dw**2 - kex**2 +
sqrt(4.0*dw**2*(kex*pA - kex*(-pA + 1.0))**2 + (-dw**2 +
kex**2)**2))/cpmg_frqs)/(cpmg_frqs*sqrt(4.0*dw**2*(kex*pA - kex*(-pA +
1.0))**2 + (-dw**2 + kex**2)**2)*sqrt(dw**2 - kex**2 +
sqrt(4.0*dw**2*(kex*pA - kex*(-pA + 1.0))**2 + (-dw**2 + kex**2)**2))) +
4.0*dw**2*eta_scale*kex*(kex*pA - kex*(-pA + 1.0))*((0.5*dw**2 +
0.5*kex**2)/sqrt(4.0*dw**2*(kex*pA - kex*(-pA + 1.0))**2 + (-dw**2 +
kex**2)**2) + 0.5)*sinh(eta_scale*sqrt(-dw**2 + kex**2 +
sqrt(4.0*dw**2*(kex*pA - kex*(-pA + 1.0))**2 + (-dw**2 +
kex**2)**2))/cpmg_frqs)/(cpmg_frqs*sqrt(4.0*dw**2*(kex*pA - kex*(-pA +
1.0))**2 + (-dw**2 + kex**2)**2)*sqrt(-dw**2 + kex**2 +
sqrt(4.0*dw**2*(kex*pA - kex*(-pA + 1.0))**2 + (-dw**2 +
kex**2)**2))))/sqrt((-((0.5*dw**2 + 0.5*kex**2)/sqrt(4.0*dw**2*(kex*pA -
kex*(-pA + 1.0))**2 + (-dw**2 + kex**2)**2) - 0.5)*cos(eta_scale*sqrt(dw**2
- kex**2 + sqrt(4.0*dw**2*(kex*pA - kex*(-pA + 1.0))**2 + (-dw**2 +
kex**2)**2))/cpmg_frqs) + ((0.5*dw**2 + 0.5*kex**2)/sqrt(4.0*dw**2*(kex*pA
- kex*(-pA + 1.0))**2 + (-dw**2 + kex**2)**2) +
0.5)*cosh(eta_scale*sqrt(-dw**2 + kex**2 + sqrt(4.0*dw**2*(kex*pA -
kex*(-pA + 1.0))**2 + (-dw**2 + kex**2)**2))/cpmg_frqs))**2 - 1)
+
+d_f_d_dw = -cpmg_frqs*(-(1.0*dw/sqrt(4.0*dw**2*(kex*pA - kex*(-pA +
1.0))**2 + (-dw**2 + kex**2)**2) + (0.5*dw**2 + 0.5*kex**2)*(2*dw*(-dw**2 +
kex**2) - 4.0*dw*(kex*pA - kex*(-pA + 1.0))**2)/(4.0*dw**2*(kex*pA -
kex*(-pA + 1.0))**2 + (-dw**2 +
kex**2)**2)**(3/2))*cos(eta_scale*sqrt(dw**2 - kex**2 +
sqrt(4.0*dw**2*(kex*pA - kex*(-pA + 1.0))**2 + (-dw**2 +
kex**2)**2))/cpmg_frqs) + (1.0*dw/sqrt(4.0*dw**2*(kex*pA - kex*(-pA +
1.0))**2 + (-dw**2 + kex**2)**2) + (0.5*dw**2 + 0.5*kex**2)*(2*dw*(-dw**2 +
kex**2) - 4.0*dw*(kex*pA - kex*(-pA + 1.0))**2)/(4.0*dw**2*(kex*pA -
kex*(-pA + 1.0))**2 + (-dw**2 +
kex**2)**2)**(3/2))*cosh(eta_scale*sqrt(-dw**2 + kex**2 +
sqrt(4.0*dw**2*(kex*pA - kex*(-pA + 1.0))**2 + (-dw**2 +
kex**2)**2))/cpmg_frqs) + eta_scale*(-dw + (-2*dw*(-dw**2 + kex**2) +
4.0*dw*(kex*pA - kex*(-pA + 1.0))**2)/(2*sqrt(4.0*dw**2*(kex*pA - kex*(-pA
+ 1.0))**2 + (-dw**2 + kex**2)**2)))*((0.5*dw**2 +
0.5*kex**2)/sqrt(4.0*dw**2*(kex*pA - kex*(-pA + 1.0))**2 + (-dw**2 +
kex**2)**2) + 0.5)*sinh(eta_scale*sqrt(-dw**2 + kex**2 +
sqrt(4.0*dw**2*(kex*pA - kex*(-pA + 1.0))**2 + (-dw**2 +
kex**2)**2))/cpmg_frqs)/(cpmg_frqs*sqrt(-dw**2 + kex**2 +
sqrt(4.0*dw**2*(kex*pA - kex*(-pA + 1.0))**2 + (-dw**2 + kex**2)**2))) +
eta_scale*(dw + (-2*dw*(-dw**2 + kex**2) + 4.0*dw*(kex*pA - kex*(-pA +
1.0))**2)/(2*sqrt(4.0*dw**2*(kex*pA - kex*(-pA + 1.0))**2 + (-dw**2 +
kex**2)**2)))*((0.5*dw**2 + 0.5*kex**2)/sqrt(4.0*dw**2*(kex*pA - kex*(-pA +
1.0))**2 + (-dw**2 + kex**2)**2) - 0.5)*sin(eta_scale*sqrt(dw**2 - kex**2 +
sqrt(4.0*dw**2*(kex*pA - kex*(-pA + 1.0))**2 + (-dw**2 +
kex**2)**2))/cpmg_frqs)/(cpmg_frqs*sqrt(dw**2 - kex**2 +
sqrt(4.0*dw**2*(kex*pA - kex*(-pA + 1.0))**2 + (-dw**2 +
kex**2)**2))))/sqrt((-((0.5*dw**2 + 0.5*kex**2)/sqrt(4.0*dw**2*(kex*pA -
kex*(-pA + 1.0))**2 + (-dw**2 + kex**2)**2) - 0.5)*cos(eta_scale*sqrt(dw**2
- kex**2 + sqrt(4.0*dw**2*(kex*pA - kex*(-pA + 1.0))**2 + (-dw**2 +
kex**2)**2))/cpmg_frqs) + ((0.5*dw**2 + 0.5*kex**2)/sqrt(4.0*dw**2*(kex*pA
- kex*(-pA + 1.0))**2 + (-dw**2 + kex**2)**2) +
0.5)*cosh(eta_scale*sqrt(-dw**2 + kex**2 + sqrt(4.0*dw**2*(kex*pA -
kex*(-pA + 1.0))**2 + (-dw**2 + kex**2)**2))/cpmg_frqs))**2 - 1)
+
+d_f_d_kex = -cpmg_frqs*(-(1.0*kex/sqrt(4.0*dw**2*(kex*pA - kex*(-pA +
1.0))**2 + (-dw**2 + kex**2)**2) + (0.5*dw**2 +
0.5*kex**2)*(-2.0*dw**2*(4*pA - 2.0)*(kex*pA - kex*(-pA + 1.0)) -
2*kex*(-dw**2 + kex**2))/(4.0*dw**2*(kex*pA - kex*(-pA + 1.0))**2 + (-dw**2
+ kex**2)**2)**(3/2))*cos(eta_scale*sqrt(dw**2 - kex**2 +
sqrt(4.0*dw**2*(kex*pA - kex*(-pA + 1.0))**2 + (-dw**2 +
kex**2)**2))/cpmg_frqs) + (1.0*kex/sqrt(4.0*dw**2*(kex*pA - kex*(-pA +
1.0))**2 + (-dw**2 + kex**2)**2) + (0.5*dw**2 +
0.5*kex**2)*(-2.0*dw**2*(4*pA - 2.0)*(kex*pA - kex*(-pA + 1.0)) -
2*kex*(-dw**2 + kex**2))/(4.0*dw**2*(kex*pA - kex*(-pA + 1.0))**2 + (-dw**2
+ kex**2)**2)**(3/2))*cosh(eta_scale*sqrt(-dw**2 + kex**2 +
sqrt(4.0*dw**2*(kex*pA - kex*(-pA + 1.0))**2 + (-dw**2 +
kex**2)**2))/cpmg_frqs) + eta_scale*(-kex + (2.0*dw**2*(4*pA - 2.0)*(kex*pA
- kex*(-pA + 1.0)) + 2*kex*(-dw**2 + kex**2))/(2*sqrt(4.0*dw**2*(kex*pA -
kex*(-pA + 1.0))**2 + (-dw**2 + kex**2)**2)))*((0.5*dw**2 +
0.5*kex**2)/sqrt(4.0*dw**2*(kex*pA - kex*(-pA + 1.0))**2 + (-dw**2 +
kex**2)**2) - 0.5)*sin(eta_scale*sqrt(dw**2 - kex**2 +
sqrt(4.0*dw**2*(kex*pA - kex*(-pA + 1.0))**2 + (-dw**2 +
kex**2)**2))/cpmg_frqs)/(cpmg_frqs*sqrt(dw**2 - kex**2 +
sqrt(4.0*dw**2*(kex*pA - kex*(-pA + 1.0))**2 + (-dw**2 + kex**2)**2))) +
eta_scale*(kex + (2.0*dw**2*(4*pA - 2.0)*(kex*pA - kex*(-pA + 1.0)) +
2*kex*(-dw**2 + kex**2))/(2*sqrt(4.0*dw**2*(kex*pA - kex*(-pA + 1.0))**2 +
(-dw**2 + kex**2)**2)))*((0.5*dw**2 + 0.5*kex**2)/sqrt(4.0*dw**2*(kex*pA -
kex*(-pA + 1.0))**2 + (-dw**2 + kex**2)**2) +
0.5)*sinh(eta_scale*sqrt(-dw**2 + kex**2 + sqrt(4.0*dw**2*(kex*pA -
kex*(-pA + 1.0))**2 + (-dw**2 +
kex**2)**2))/cpmg_frqs)/(cpmg_frqs*sqrt(-dw**2 + kex**2 +
sqrt(4.0*dw**2*(kex*pA - kex*(-pA + 1.0))**2 + (-dw**2 +
kex**2)**2))))/sqrt((-((0.5*dw**2 + 0.5*kex**2)/sqrt(4.0*dw**2*(kex*pA -
kex*(-pA + 1.0))**2 + (-dw**2 + kex**2)**2) - 0.5)*cos(eta_scale*sqrt(dw**2
- kex**2 + sqrt(4.0*dw**2*(kex*pA - kex*(-pA + 1.0))**2 + (-dw**2 +
kex**2)**2))/cpmg_frqs) + ((0.5*dw**2 + 0.5*kex**2)/sqrt(4.0*dw**2*(kex*pA
- kex*(-pA + 1.0))**2 + (-dw**2 + kex**2)**2) +
0.5)*cosh(eta_scale*sqrt(-dw**2 + kex**2 + sqrt(4.0*dw**2*(kex*pA -
kex*(-pA + 1.0))**2 + (-dw**2 + kex**2)**2))/cpmg_frqs))**2 - 1) + 0.5
+
+------------------------------------------------------------------------------
Added:
branches/est_par_error/test_suite/shared_data/dispersion/estimate_par_err/cr72/cr72_simple_der.py
URL:
http://svn.gna.org/viewcvs/relax/branches/est_par_error/test_suite/shared_data/dispersion/estimate_par_err/cr72/cr72_simple_der.py?rev=25567&view=auto
==============================================================================
---
branches/est_par_error/test_suite/shared_data/dispersion/estimate_par_err/cr72/cr72_simple_der.py
(added)
+++
branches/est_par_error/test_suite/shared_data/dispersion/estimate_par_err/cr72/cr72_simple_der.py
Tue Sep 2 20:06:13 2014
@@ -0,0 +1,75 @@
+from sympy import *
+
+# In contrast to other Computer Algebra Systems, in SymPy you have to
declare symbolic variables explicitly:
+r20a = Symbol('r20a')
+pA = Symbol('pA')
+dw = Symbol('dw')
+kex = Symbol('kex')
+cpmg_frqs = Symbol('cpmg_frqs')
+eta_scale = Symbol('eta_scale')
+
+r20b = r20a
+
+pB = 1.0 - pA
+dw2 = dw**2
+r20_kex = (r20a + r20b + kex) / 2.0
+k_BA = pA * kex
+k_AB = pB * kex
+# When r20b = r20a
+Psi = kex**2 - dw2
+zeta = -2.0*dw * (k_BA - k_AB)
+
+sqrt_psi2_zeta2 = sqrt(Psi**2 + zeta**2)
+D_part = (0.5*Psi + dw2) / sqrt_psi2_zeta2
+Dpos = 0.5 + D_part
+Dneg = -0.5 + D_part
+eta_fact = eta_scale / cpmg_frqs
+etapos = eta_fact * sqrt(Psi + sqrt_psi2_zeta2)
+etaneg = eta_fact * sqrt(-Psi + sqrt_psi2_zeta2)
+fact = Dpos * cosh(etapos) - Dneg * cos(etaneg)
+back_calc = cpmg_frqs * acosh(fact)
+back_calc = r20_kex - back_calc
+
+
+d_f_d_r20a = diff(back_calc, r20a)
+d_f_d_pA = diff(back_calc, pA)
+d_f_d_dw = diff(back_calc, dw)
+d_f_d_kex = diff(back_calc, kex)
+
+tfile = open("cr72_simple_der.log", 'w')
+
+text = ("""Form the Jacobian matrix by:
+------------------------------------------------------------------------------
+
+d_f_d_r20a = %s
+
+d_f_d_pA = %s
+
+d_f_d_dw = %s
+
+d_f_d_kex = %s
+
+------------------------------------------------------------------------------
+"""% (d_f_d_r20a, d_f_d_pA, d_f_d_dw, d_f_d_kex) )
+
+tfile.write(text)
+tfile.close()
+
+# Numpy chokes on this.
+#sim_d_f_d_r20a = simplify(d_f_d_r20a)
+#print sim_d_f_d_r20a
+#sim_d_f_d_r20b = simplify(d_f_d_r20b)
+#sim_d_f_d_r20b = simplify(d_f_d_pA)
+#sim_d_f_d_dw = simplify(d_f_d_dw)
+#sim_d_f_d_kex = simplify(d_f_d_kex)
+
+# The vectorial function.
+X = Matrix([back_calc])
+# What to derive for.
+Y = Matrix([r20a, r20b, pA, dw, kex])
+
+# Make the Jacobian
+Jacobian = X.jacobian(Y)
+
+#Jac2 = Jacobian.simplify()
+#simplify(Jacobian)
_______________________________________________
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
Re: r25567 - /branches/est_par_error/test_suite/shared_data/dispersion/estimate_par_err/cr72/