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 )