Hi, This is an awesome improvement! Note that you should also mention the removal of the debugging code in the message, or better have it in a different commit.
Cheers, Edward On 19 June 2014 17:42, <[email protected]> wrote: > 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 > [email protected] > > 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 _______________________________________________ relax (http://www.nmr-relax.com) This is the relax-devel mailing list [email protected] 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-devel
