Author: tlinnet
Date: Thu Jun 19 17:42:05 2014
New Revision: 24158
URL: http://svn.gna.org/viewcvs/relax?rev=24158&view=rev
Log:
Implemented double speed of model NS CPMG 2site 3D:
This is done by moving the costly calculation of the matrix exponential out
of the for loops.
The trick was to find a method to do dot product of higher dimensions.
Thiw was done with numpy.einsum:
Example at:
http://wiki.nmr-relax.com/Numpy_linalg#Ellipsis_broadcasting_in_numpy.einsum
Example:
dot_V_W = einsum('...ij,...jk', V, W_exp_diag)
Where V, and W_exp_diag has shape: [NE][NS][NM][NO][ND][7][7]
The profiling script shows a 2X speed up.
----BEFORE:
SINGLE
ncalls tottime percall cumtime percall filename:lineno(function)
1 0.000 0.000 18.811 18.811 <string>:1(<module>)
1 0.002 0.002 18.811 18.811 pf_3d:407(single)
CLUSTER
ncalls tottime percall cumtime percall filename:lineno(function)
1 0.000 0.000 18.315 18.315 <string>:1(<module>)
1 0.001 0.001 18.315 18.315 pf_3d:431(cluster)
-----AFTER:
SINGLE
ncalls tottime percall cumtime percall filename:lineno(function)
1 0.000 0.000 8.818 8.818 <string>:1(<module>)
1 0.002 0.002 8.818 8.818 pf_3d:407(single)
CLUSTER
ncalls tottime percall cumtime percall filename:lineno(function)
1 0.000 0.000 9.082 9.082 <string>:1(<module>)
1 0.001 0.001 9.082 9.082 pf_3d:431(cluster)
Task #7807 (https://gna.org/task/index.php?7807): Speed-up of dispersion
models for Clustered analysis.
Modified:
branches/disp_spin_speed/lib/dispersion/ns_cpmg_2site_3d.py
Modified: branches/disp_spin_speed/lib/dispersion/ns_cpmg_2site_3d.py
URL:
http://svn.gna.org/viewcvs/relax/branches/disp_spin_speed/lib/dispersion/ns_cpmg_2site_3d.py?rev=24158&r1=24157&r2=24158&view=diff
==============================================================================
--- branches/disp_spin_speed/lib/dispersion/ns_cpmg_2site_3d.py (original)
+++ branches/disp_spin_speed/lib/dispersion/ns_cpmg_2site_3d.py Thu Jun 19
17:42:05 2014
@@ -55,7 +55,7 @@
"""
# Python module imports.
-from numpy import asarray, dot, fabs, isfinite, log, min, sum, tile
+from numpy import asarray, dot, fabs, isfinite, log, min, newaxis, sum,
tile
from numpy.ma import fix_invalid, masked_where
@@ -134,32 +134,10 @@
# This matrix is a propagator that will evolve the magnetization with
the matrix R for a delay tcp.
Rexpo_mat = matrix_exponential_rankN(R_mat)
- # Loop over the spins
- for si in range(NS):
- # Loop over the spectrometer frequencies.
- for mi in range(NM):
- # Extract number of points.
- num_points_si_mi = int(num_points[0, si, mi, 0])
-
- # Loop over the time points, back calculating the R2eff values.
- for di in range(num_points_si_mi):
-
- # This matrix is a propagator that will evolve the
magnetization with the matrix R for a delay tcp.
- R_mat_i = R_mat[0, si, mi, 0, di]
-
- # Store the Rexpo. Is it not possible to find the matrix
exponential of higher dimensional data?
- Rexpo = matrix_exponential(R_mat_i)
- Rexpo_t = Rexpo_mat[0, si, mi, 0, di]
-
- diff = Rexpo - Rexpo_t
- if abs(sum(diff)) > 1e-15:
- print abs(sum(diff))
- import sys
- sys.exit()
-
# Initial magnetisation.
- Mint_mat = tile(M0[None, None, None, None, None, :, None], (NE, NS,
NM, NO, ND, 1, 1) )
- r180x_mat = tile(r180x[None, None, None, None, None, :, :], (NE, NS,
NM, NO, ND, 1, 1) )
+ # Expand axis, and tile up to dimensions.
+ Mint_mat = tile(M0[newaxis, newaxis, newaxis, newaxis, newaxis, :,
newaxis], (NE, NS, NM, NO, ND, 1, 1) )
+ r180x_mat = tile(r180x[newaxis, newaxis, newaxis, newaxis, newaxis,
...], (NE, NS, NM, NO, ND, 1, 1) )
# Loop over the spins
for si in range(NS):
@@ -180,13 +158,13 @@
Mint = Mint_mat[0, si, mi, 0, di]
# This matrix is a propagator that will evolve the
magnetization with the matrix R for a delay tcp.
- Rexpo = Rexpo_mat[0, si, mi, 0, di]
+ Rexpo_i = Rexpo_mat[0, si, mi, 0, di]
r180x_i = r180x_mat[0, si, mi, 0, di]
# The essential evolution matrix.
# This is the first round.
- evolution_matrix = dot(Rexpo, r180x_i)
- evolution_matrix = dot(evolution_matrix, Rexpo)
+ evolution_matrix = dot(Rexpo_i, r180x_i)
+ evolution_matrix = dot(evolution_matrix, Rexpo_i)
# The second round.
evolution_matrix = dot(evolution_matrix, evolution_matrix )
_______________________________________________
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: r24158 - /branches/disp_spin_speed/lib/dispersion/ns_cpmg_2site_3d.py