Author: tlinnet Date: Wed Jun 25 02:14:32 2014 New Revision: 24292 URL: http://svn.gna.org/viewcvs/relax?rev=24292&view=rev Log: Inserted an extremely interesting development in NS R1rho 2site. If one do a transpose of M0, one can calculate all the matrix evolutions in the start via numpy einsum. Since M0 is in higher a dimensions, one should not do a numpy transpose, but swap/roll the outer M0 6x1 axis. 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_r1rho_2site.py Modified: branches/disp_spin_speed/lib/dispersion/ns_r1rho_2site.py URL: http://svn.gna.org/viewcvs/relax/branches/disp_spin_speed/lib/dispersion/ns_r1rho_2site.py?rev=24292&r1=24291&r2=24292&view=diff ============================================================================== --- branches/disp_spin_speed/lib/dispersion/ns_r1rho_2site.py (original) +++ branches/disp_spin_speed/lib/dispersion/ns_r1rho_2site.py Wed Jun 25 02:14:32 2014 @@ -50,7 +50,7 @@ """ # Python module imports. -from numpy import array, cos, dot, einsum, float64, log, multiply, sin, sum +from numpy import array, cos, dot, einsum, float64, log, multiply, rollaxis, sin, sum # relax module imports. from lib.float import isNaN @@ -238,31 +238,63 @@ # Magnetization evolution. Rexpo_M0_mat = einsum('...ij,...jk', Rexpo_mat, M0) - # Loop over spins. - for si in range(NS): - # Loop over the spectrometer frequencies. - for mi in range(NM): - # Loop over offsets: - for oi in range(NO): - # Extract number of points. - num_points_i = num_points[0, si, mi, oi] - - # Loop over the time points, back calculating the R2eff values. - for j in range(num_points_i): + # Transpose M0, to prepare for dot operation. Roll the last axis one back. + M0_T = rollaxis(M0, 6, 5) + + if NS != 1: + # Magnetization evolution, which include all dimensions. + MA_mat = einsum('...ij,...jk', M0_T, Rexpo_M0_mat) + + # Loop over the spectrometer frequencies. + for mi in range(NM): + # Loop over offsets: + for oi in range(NO): + # Extract number of points. + num_points_i = num_points[0, 0, mi, oi] + + # Loop over the time points, back calculating the R2eff values. + for j in range(num_points_i): + + # If the number of spins are 1, do the fastest method by dot product. Else do it by einstein summation. + if NS == 1: + # Set the spin index for one spin. + si = 0 # Extract the preformed matrix that rotate the magnetization previous to spin-lock into the weff frame. - M0_i= M0[0, si, mi, oi, j, :, 0] - - # This matrix is a propagator that will evolve the magnetization with the matrix R. - Rexpo_i = Rexpo_mat[0, si, mi, oi, j] + M0_i= M0_T[0, si, mi, oi, j] # Extract from the pre-formed Magnetization evolution matrix. - Rexpo_M0_mat_i = Rexpo_M0_mat[0, si, mi, oi, j, :, 0] + Rexpo_M0_mat_i = Rexpo_M0_mat[0, si, mi, oi, j] # Magnetization evolution. - MA = dot(M0_i, Rexpo_M0_mat_i) + MA = dot(M0_i, Rexpo_M0_mat_i)[0, 0] # The next lines calculate the R1rho using a two-point approximation, i.e. assuming that the decay is mono-exponential. if MA <= 0.0 or isNaN(MA): back_calc[0, si, mi, oi, j] = 1e99 else: back_calc[0, si, mi, oi, j]= -inv_relax_time[0, si, mi, oi, j] * log(MA) + + # If there is multiple spin a clustered analysis. + else: + # Loop over spins. + for si in range(NS): + # Extract the preformed matrix that rotate the magnetization previous to spin-lock into the weff frame. + M0_i= M0_T[0, si, mi, oi, j] + + # Extract from the pre-formed Magnetization evolution matrix. + Rexpo_M0_mat_i = Rexpo_M0_mat[0, si, mi, oi, j] + + # Magnetization evolution. + MA = dot(M0_i, Rexpo_M0_mat_i)[0, 0] + + MA_mat_i = MA_mat[0, si, mi, oi, j, 0, 0] + + # Diff + diff = MA - MA_mat_i + if diff != 0.0: + print "oh no" + + if MA_mat_i <= 0.0 or isNaN(MA_mat_i): + back_calc[0, si, mi, oi, j] = 1e99 + else: + back_calc[0, si, mi, oi, j]= -inv_relax_time[0, si, mi, oi, j] * log(MA_mat_i)