Author: bugman
Date: Tue Oct 15 13:27:37 2013
New Revision: 21113
URL: http://svn.gna.org/viewcvs/relax?rev=21113&view=rev
Log:
Replaced all usage of scipy.linalg.expm() with
lib.linear_algebra.matrix_exponential.matrix_exponential().
This is for the functions of the lib.dispersion package used for the
relaxation dispersion numeric
solution models. The change eliminates a bug in the scipy function which
uses the Pade
approximation which fails horribly for the complex part of the matrix. The
real part looks good,
but the complex part looks to have nasty truncation artefacts which is
propagated and amplified
through the Bloch-McConnell equations.
Modified:
branches/relax_disp/lib/dispersion/mq_ns_cpmg_2site.py
branches/relax_disp/lib/dispersion/ns_cpmg_2site_3d.py
branches/relax_disp/lib/dispersion/ns_cpmg_2site_expanded.py
branches/relax_disp/lib/dispersion/ns_cpmg_2site_star.py
branches/relax_disp/lib/dispersion/ns_r1rho_2site.py
Modified: branches/relax_disp/lib/dispersion/mq_ns_cpmg_2site.py
URL:
http://svn.gna.org/viewcvs/relax/branches/relax_disp/lib/dispersion/mq_ns_cpmg_2site.py?rev=21113&r1=21112&r2=21113&view=diff
==============================================================================
--- branches/relax_disp/lib/dispersion/mq_ns_cpmg_2site.py (original)
+++ branches/relax_disp/lib/dispersion/mq_ns_cpmg_2site.py Tue Oct 15
13:27:37 2013
@@ -36,12 +36,11 @@
from math import fabs, log
from numpy import array, conj, dot, float64
from numpy.linalg import matrix_power
-if dep_check.scipy_module:
- from scipy.linalg import expm
# relax module imports.
from lib.dispersion.ns_matrices import rcpmg_3d
from lib.float import isNaN
+from lib.linear_algebra.matrix_exponential import matrix_exponential
def populate_matrix(matrix=None, r20=None, dw=None, dwH=None, k_AB=None,
k_BA=None):
@@ -119,8 +118,9 @@
# Loop over the time points, back calculating the R2eff values.
for i in range(num_points):
# The M1 and M2 matrices.
- M1 = expm(m1*tcp[i])
- M2 = expm(m2*tcp[i])
+ M1 = matrix_exponential(m1*tcp[i])
+ M2 = matrix_exponential(m2*tcp[i])
+
# The complex conjugates M1* and M2*
M1_star = conj(M1)
Modified: branches/relax_disp/lib/dispersion/ns_cpmg_2site_3d.py
URL:
http://svn.gna.org/viewcvs/relax/branches/relax_disp/lib/dispersion/ns_cpmg_2site_3d.py?rev=21113&r1=21112&r2=21113&view=diff
==============================================================================
--- branches/relax_disp/lib/dispersion/ns_cpmg_2site_3d.py (original)
+++ branches/relax_disp/lib/dispersion/ns_cpmg_2site_3d.py Tue Oct 15
13:27:37 2013
@@ -36,12 +36,11 @@
# Python module imports.
from math import fabs, log
from numpy import dot
-if dep_check.scipy_module:
- from scipy.linalg import expm
# relax module imports.
from lib.dispersion.ns_matrices import rcpmg_3d
from lib.float import isNaN
+from lib.linear_algebra.matrix_exponential import matrix_exponential
def r2eff_ns_cpmg_2site_3D(r180x=None, M0=None, r10a=0.0, r10b=0.0,
r20a=None, r20b=None, pA=None, pB=None, dw=None, k_AB=None, k_BA=None,
inv_tcpmg=None, tcp=None, back_calc=None, num_points=None, power=None):
@@ -93,7 +92,7 @@
Mint = M0
# This matrix is a propagator that will evolve the magnetization
with the matrix R for a delay tcp.
- Rexpo = expm(R*tcp[i])
+ Rexpo = matrix_exponential(R*tcp[i])
# Loop over the CPMG elements, propagating the magnetisation.
for j in range(2*power[i]):
Modified: branches/relax_disp/lib/dispersion/ns_cpmg_2site_expanded.py
URL:
http://svn.gna.org/viewcvs/relax/branches/relax_disp/lib/dispersion/ns_cpmg_2site_expanded.py?rev=21113&r1=21112&r2=21113&view=diff
==============================================================================
--- branches/relax_disp/lib/dispersion/ns_cpmg_2site_expanded.py (original)
+++ branches/relax_disp/lib/dispersion/ns_cpmg_2site_expanded.py Tue Oct 15
13:27:37 2013
@@ -47,8 +47,6 @@
from math import log
from numpy import add, complex, conj, dot, exp, power, real, sqrt
from numpy.linalg import matrix_power
-if dep_check.scipy_module:
- from scipy.linalg import expm
# relax module imports.
from lib.float import isNaN
Modified: branches/relax_disp/lib/dispersion/ns_cpmg_2site_star.py
URL:
http://svn.gna.org/viewcvs/relax/branches/relax_disp/lib/dispersion/ns_cpmg_2site_star.py?rev=21113&r1=21112&r2=21113&view=diff
==============================================================================
--- branches/relax_disp/lib/dispersion/ns_cpmg_2site_star.py (original)
+++ branches/relax_disp/lib/dispersion/ns_cpmg_2site_star.py Tue Oct 15
13:27:37 2013
@@ -37,11 +37,10 @@
from math import log
from numpy import add, complex, conj, dot
from numpy.linalg import matrix_power
-if dep_check.scipy_module:
- from scipy.linalg import expm
# relax module imports.
from lib.float import isNaN
+from lib.linear_algebra.matrix_exponential import matrix_exponential
def r2eff_ns_cpmg_2site_star(Rr=None, Rex=None, RCS=None, R=None, M0=None,
r20a=None, r20b=None, dw=None, inv_tcpmg=None, tcp=None, back_calc=None,
num_points=None, power=None):
@@ -95,10 +94,10 @@
# Loop over the time points, back calculating the R2eff values.
for i in range(num_points):
# This matrix is a propagator that will evolve the magnetization
with the matrix R for a delay tcp.
- expm_R_tcp = expm(R*tcp[i])
+ eR_tcp = matrix_exponential(R*tcp[i])
# This is the propagator for an element of [delay tcp; 180 deg
pulse; 2 times delay tcp; 180 deg pulse; delay tau], i.e. for 2 times
tau-180-tau.
- prop_2 = dot(dot(expm_R_tcp, expm(cR2*tcp[i])), expm_R_tcp)
+ prop_2 = dot(dot(eR_tcp, matrix_exponential(cR2*tcp[i])), eR_tcp)
# Now create the total propagator that will evolve the magnetization
under the CPMG train, i.e. it applies the above tau-180-tau-tau-180-tau so
many times as required for the CPMG frequency under consideration.
prop_total = matrix_power(prop_2, power[i])
Modified: branches/relax_disp/lib/dispersion/ns_r1rho_2site.py
URL:
http://svn.gna.org/viewcvs/relax/branches/relax_disp/lib/dispersion/ns_r1rho_2site.py?rev=21113&r1=21112&r2=21113&view=diff
==============================================================================
--- branches/relax_disp/lib/dispersion/ns_r1rho_2site.py (original)
+++ branches/relax_disp/lib/dispersion/ns_r1rho_2site.py Tue Oct 15 13:27:37
2013
@@ -33,12 +33,11 @@
# Python module imports.
from math import atan, cos, log, pi, sin, sqrt
from numpy import dot
-if dep_check.scipy_module:
- from scipy.linalg import expm
# relax module imports.
from lib.dispersion.ns_matrices import rr1rho_3d
from lib.float import isNaN
+from lib.linear_algebra.matrix_exponential import matrix_exponential
def ns_r1rho_2site(M0=None, r1rho_prime=None, omega=None, offset=None,
r1=0.0, pA=None, pB=None, dw=None, k_AB=None, k_BA=None,
spin_lock_fields=None, relax_time=None, inv_relax_time=None, back_calc=None,
num_points=None):
@@ -102,7 +101,7 @@
M0[5] = cos(thetaB) * pB
# This matrix is a propagator that will evolve the magnetization
with the matrix R.
- Rexpo = expm(R*relax_time)
+ Rexpo = matrix_exponential(R*relax_time)
# Magnetization evolution.
Moft = dot(Rexpo, M0)
r21113 - /branches/relax_disp/lib/dispersion/