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)
r25567 - /branches/est_par_error/test_suite/shared_data/dispersion/estimate_par_err/cr72/