Author: tlinnet Date: Tue Jun 24 15:36:01 2014 New Revision: 24283 URL: http://svn.gna.org/viewcvs/relax?rev=24283&view=rev Log: Speeded up ns mmq 3site, by moving the forming of evolution matrix out of the for loops, and preform it. 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_mmq_3site.py Modified: branches/disp_spin_speed/lib/dispersion/ns_mmq_3site.py URL: http://svn.gna.org/viewcvs/relax/branches/disp_spin_speed/lib/dispersion/ns_mmq_3site.py?rev=24283&r1=24282&r2=24283&view=diff ============================================================================== --- branches/disp_spin_speed/lib/dispersion/ns_mmq_3site.py (original) +++ branches/disp_spin_speed/lib/dispersion/ns_mmq_3site.py Tue Jun 24 15:36:01 2014 @@ -57,7 +57,7 @@ # Python module imports. from math import floor -from numpy import array, conj, dot, float64, log, multiply, sum +from numpy import array, conj, dot, einsum, float64, log, multiply, sum # relax module imports. from lib.float import isNaN @@ -297,8 +297,20 @@ M2_mat = matrix_exponential_rank_NS_NM_NO_ND_x_x(m2_mat) # The complex conjugates M1* and M2* + # Equivalent to D+*. M1_star_mat = conj(M1_mat) + # Equivalent to Z-*. M2_star_mat = conj(M2_mat) + + # Repetitive dot products (minimised for speed). + M1_M2_mat = einsum('...ij,...jk', M1_mat, M2_mat) + M2_M1_mat = einsum('...ij,...jk', M2_mat, M1_mat) + M1_M2_M2_M1_mat = einsum('...ij,...jk', M1_M2_mat, M2_M1_mat) + M2_M1_M1_M2_mat = einsum('...ij,...jk', M2_M1_mat, M1_M2_mat) + M1_M2_star_mat = einsum('...ij,...jk', M1_star_mat, M2_star_mat) + M2_M1_star_mat = einsum('...ij,...jk', M2_star_mat, M1_star_mat) + M1_M2_M2_M1_star_mat = einsum('...ij,...jk', M1_M2_star_mat, M2_M1_star_mat) + M2_M1_M1_M2_star_mat = einsum('...ij,...jk', M2_M1_star_mat, M1_M2_star_mat) # Loop over spins. for si in range(NS): @@ -311,79 +323,68 @@ # Loop over the time points, back calculating the R2eff values. for i in range(num_points_i): - # The M1 and M2 matrices. - # Equivalent to D+. - M1_i = M1_mat[si, mi, oi, i] - # Equivalent to Z-. - M2_i = M2_mat[si, mi, oi, i] - - # The complex conjugates M1* and M2* - # Equivalent to D+*. - M1_star_i = M1_star_mat[si, mi, oi, i] - # Equivalent to Z-*. - M2_star_i = M2_star_mat[si, mi, oi, i] - - # Repetitive dot products (minimised for speed). - M1_M2 = dot(M1_i, M2_i) - M2_M1 = dot(M2_i, M1_i) - M1_M2_M2_M1 = dot(M1_M2, M2_M1) - M2_M1_M1_M2 = dot(M2_M1, M1_M2) - M1_M2_star = dot(M1_star_i, M2_star_i) - M2_M1_star = dot(M2_star_i, M1_star_i) - M1_M2_M2_M1_star = dot(M1_M2_star, M2_M1_star) - M2_M1_M1_M2_star = dot(M2_M1_star, M1_M2_star) + # Extract data from array. + power_i = power[si, mi, oi, i] + M1_M2_i = M1_M2_mat[si, mi, oi, i] + M1_M2_star_i = M1_M2_star_mat[si, mi, oi, i] + M2_M1_i = M2_M1_mat[si, mi, oi, i] + M2_M1_star_i = M2_M1_star_mat[si, mi, oi, i] + M1_M2_M2_M1_i = M1_M2_M2_M1_mat[si, mi, oi, i] + M2_M1_M1_M2_star_i = M2_M1_M1_M2_star_mat[si, mi, oi, i] + M2_M1_M1_M2_i = M2_M1_M1_M2_mat[si, mi, oi, i] + M1_M2_M2_M1_star_i = M1_M2_M2_M1_star_mat[si, mi, oi, i] # Special case of 1 CPMG block - the power is zero. - if power[si, mi, oi, i] == 1: + if power_i == 1: # M1.M2. - A = M1_M2 + A = M1_M2_i # M1*.M2*. - B = M1_M2_star + B = M1_M2_star_i # M2.M1. - C = M2_M1 + C = M2_M1_i # M2*.M1*. - D = M2_M1_star + D = M2_M1_star_i # Matrices for even number of CPMG blocks. - elif power[si, mi, oi, i] % 2 == 0: + elif power_i % 2 == 0: # The power factor (only calculate once). - fact = int(floor(power[si, mi, oi, i] / 2)) + fact = int(floor(power_i / 2)) # (M1.M2.M2.M1)^(n/2). - A = square_matrix_power(M1_M2_M2_M1, fact) + A = square_matrix_power(M1_M2_M2_M1_i, fact) # (M2*.M1*.M1*.M2*)^(n/2). - B = square_matrix_power(M2_M1_M1_M2_star, fact) + B = square_matrix_power(M2_M1_M1_M2_star_i, fact) # (M2.M1.M1.M2)^(n/2). - C = square_matrix_power(M2_M1_M1_M2, fact) + C = square_matrix_power(M2_M1_M1_M2_i, fact) # (M1*.M2*.M2*.M1*)^(n/2). - D = square_matrix_power(M1_M2_M2_M1_star, fact) + D = square_matrix_power(M1_M2_M2_M1_star_i, fact) # Matrices for odd number of CPMG blocks. else: # The power factor (only calculate once). - fact = int(floor((power[si, mi, oi, i] - 1) / 2)) + fact = int(floor((power_i - 1) / 2)) # (M1.M2.M2.M1)^((n-1)/2).M1.M2. - A = square_matrix_power(M1_M2_M2_M1, fact) - A = dot(A, M1_M2) + A = square_matrix_power(M1_M2_M2_M1_i, fact) + A = dot(A, M1_M2_i) # (M1*.M2*.M2*.M1*)^((n-1)/2).M1*.M2*. - B = square_matrix_power(M1_M2_M2_M1_star, fact) - B = dot(B, M1_M2_star) + B = square_matrix_power(M1_M2_M2_M1_star_i, fact) + B = dot(B, M1_M2_star_i) # (M2.M1.M1.M2)^((n-1)/2).M2.M1. - C = square_matrix_power(M2_M1_M1_M2, fact) - C = dot(C, M2_M1) + C = square_matrix_power(M2_M1_M1_M2_i, fact) + C = dot(C, M2_M1_i) # (M2*.M1*.M1*.M2*)^((n-1)/2).M2*.M1*. - D = square_matrix_power(M2_M1_M1_M2_star, fact) - D = dot(D, M2_M1_star) + D = square_matrix_power(M2_M1_M1_M2_star_i, fact) + D = dot(D, M2_M1_star_i) # The next lines calculate the R2eff using a two-point approximation, i.e. assuming that the decay is mono-exponential. A_B = dot(A, B) @@ -476,6 +477,11 @@ A_pos_mat = matrix_exponential_rank_NS_NM_NO_ND_x_x(m1_mat) A_neg_mat = matrix_exponential_rank_NS_NM_NO_ND_x_x(m2_mat) + # The evolution for one n. + evol_block_mat = einsum('...ij,...jk', A_neg_mat, A_pos_mat) + evol_block_mat = einsum('...ij,...jk', A_neg_mat, evol_block_mat) + evol_block_mat = einsum('...ij,...jk', A_pos_mat, evol_block_mat) + # Loop over spins. for si in range(NS): # Loop over the spectrometer frequencies. @@ -487,15 +493,12 @@ # Loop over the time points, back calculating the R2eff values. for i in range(num_points_i): - # The A+/- matrices. - A_pos_i = A_pos_mat[si, mi, oi, i] - A_neg_i = A_neg_mat[si, mi, oi, i] - - # The evolution for one n. - evol_block = dot(A_pos_i, dot(A_neg_i, dot(A_neg_i, A_pos_i))) + # Extract data from array. + power_i = power[si, mi, oi, i] + evol_block_i = evol_block_mat[si, mi, oi, i] # The full evolution. - evol = square_matrix_power(evol_block, power[si, mi, oi, i]) + evol = square_matrix_power(evol_block_i, power_i) # The next lines calculate the R2eff using a two-point approximation, i.e. assuming that the decay is mono-exponential. Mx = dot(F_vector, dot(evol, M0))