Hi, These changes also significantly affect the Relax_disp system test speed. In the trunk, these tests take 221.226s on one system. In the branch, it is 175.890s. But there are 6 new system tests in the branch for a total of 14.82s. So the system tests are now ~1.4 times faster :)
Cheers, Edward On 25 June 2014 02:14, <[email protected]> wrote: > Author: tlinnet > Date: Wed Jun 25 02:14:38 2014 > New Revision: 24295 > > URL: http://svn.gna.org/viewcvs/relax?rev=24295&view=rev > Log: > Speeded up the code of NS r1rho 3site. > > This was essential done to numpy einsum, and doing the dot operations in > multiple dimensions. > It was though necessary to realize, that to do the proper dot product > operations, the outer two > axis if M0 should be swapped, by rolling the outer axis one back. > > 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_3site.py > branches/disp_spin_speed/target_functions/relax_disp.py > > Modified: branches/disp_spin_speed/lib/dispersion/ns_r1rho_3site.py > URL: > http://svn.gna.org/viewcvs/relax/branches/disp_spin_speed/lib/dispersion/ns_r1rho_3site.py?rev=24295&r1=24294&r2=24295&view=diff > ============================================================================== > --- branches/disp_spin_speed/lib/dispersion/ns_r1rho_3site.py (original) > +++ branches/disp_spin_speed/lib/dispersion/ns_r1rho_3site.py Wed Jun 25 > 02:14:38 2014 > @@ -56,11 +56,10 @@ > """ > > # Python module imports. > -from math import atan2 > -from numpy import array, cos, dot, einsum, float64, log, multiply, sin, sum > +from numpy import array, einsum, float64, isfinite, log, min, multiply, sin, > sum > +from numpy.ma import fix_invalid, masked_less > > # relax module imports. > -from lib.float import isNaN > from lib.dispersion.matrix_exponential import > matrix_exponential_rank_NE_NS_NM_NO_ND_x_x > > # Repetitive calculations (to speed up calculations). > @@ -332,16 +331,17 @@ > return matrix > > > -def ns_r1rho_3site(M0=None, r1rho_prime=None, omega=None, offset=None, > r1=0.0, pA=None, pB=None, dw_AB=None, dw_BC=None, kex_AB=None, kex_BC=None, > kex_AC=None, spin_lock_fields=None, relax_time=None, inv_relax_time=None, > back_calc=None, num_points=None): > +def ns_r1rho_3site(M0=None, M0_T=None, r1rho_prime=None, omega=None, > offset=None, r1=0.0, pA=None, pB=None, dw_AB=None, dw_BC=None, kex_AB=None, > kex_BC=None, kex_AC=None, spin_lock_fields=None, relax_time=None, > inv_relax_time=None, back_calc=None, num_points=None): > """The 3-site numerical solution to the Bloch-McConnell equation for > R1rho data. > > This function calculates and stores the R1rho values. > > > @keyword M0: This is a vector that contains the initial > magnetizations corresponding to the A and B state transverse magnetizations. > - @type M0: numpy float64, rank-1, 7D array > + @type M0: numpy float array of rank > [NE][NS][NM][NO][ND][9][1] > + @keyword M0_T: This is a vector that contains the initial > magnetizations corresponding to the A and B state transverse magnetizations, > where the outer two axis has been swapped for efficient dot operations. > @keyword r1rho_prime: The R1rho_prime parameter value (R1rho with > no exchange). > - @type r1rho_prime: numpy float array of rank [NS][NM][NO][ND] > + @type r1rho_prime: numpy float array of rank > [NE][NS][NM][NO][ND][1][9] > @keyword omega: The chemical shift for the spin in rad/s. > @type omega: numpy float array of rank [NS][NM][NO][ND] > @keyword offset: The spin-lock offsets for the data. > @@ -399,31 +399,20 @@ > # 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 parameters from array. > - 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): > - # 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] > - > - # Extract from the pre-formed Magnetization evolution > matrix. > - Rexpo_M0_mat_i = Rexpo_M0_mat[0, si, mi, oi, j, :, 0] > - > - # Magnetization evolution. > - MA = dot(M0_i, Rexpo_M0_mat_i) > - > - # 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) > + # Magnetization evolution, which include all dimensions. > + MA_mat = einsum('...ij,...jk', M0_T, Rexpo_M0_mat)[:, :, :, :, :, 0, 0] > + > + # Insert safe checks. > + if min(MA_mat) < 0.0: > + mask_min_MA_mat = masked_less(MA_mat, 0.0) > + # Fill with high values. > + MA_mat[mask_min_MA_mat.mask] = 1e100 > + > + # Do back calculation. > + back_calc[:] = -inv_relax_time * log(MA_mat) > + > + # Catch errors, taking a sum over array is the fastest way to check for > + # +/- inf (infinity) and nan (not a number). > + if not isfinite(sum(back_calc)): > + # Replaces nan, inf, etc. with fill value. > + fix_invalid(back_calc, copy=False, fill_value=1e100) > > Modified: branches/disp_spin_speed/target_functions/relax_disp.py > URL: > http://svn.gna.org/viewcvs/relax/branches/disp_spin_speed/target_functions/relax_disp.py?rev=24295&r1=24294&r2=24295&view=diff > ============================================================================== > --- branches/disp_spin_speed/target_functions/relax_disp.py (original) > +++ branches/disp_spin_speed/target_functions/relax_disp.py Wed Jun 25 > 02:14:38 2014 > @@ -432,6 +432,8 @@ > # The A state initial Z magnetisation. > M0_cos = multiply.outer( cos(theta_mat), M0_2 ) > self.M0 = M0_sin + M0_cos > + # Transpose M0, to prepare for dot operation. Roll the last axis > one back, corresponds to a transpose for the outer two axis. > + self.M0_T = rollaxis(self.M0, 6, 5) > > # Set up the model. > if model == MODEL_NOREX: > @@ -772,7 +774,7 @@ > self.r20_struct[:] = multiply.outer( r1rho_prime.reshape(self.NE, > self.NS, self.NM), self.no_nd_ones ) > > # Back calculate the R2eff values for each experiment type. > - ns_r1rho_3site(M0=self.M0, r1rho_prime=self.r20_struct, > omega=self.chemical_shifts, offset=self.offset, r1=self.r1, pA=pA, pB=pB, > dw_AB=self.dw_AB_struct, dw_BC=self.dw_BC_struct, kex_AB=kex_AB, > kex_BC=kex_BC, kex_AC=kex_AC, spin_lock_fields=self.spin_lock_omega1, > relax_time=self.relax_times, inv_relax_time=self.inv_relax_times, > back_calc=self.back_calc, num_points=self.num_disp_points) > + ns_r1rho_3site(M0=self.M0, M0_T=self.M0_T, > r1rho_prime=self.r20_struct, omega=self.chemical_shifts, offset=self.offset, > r1=self.r1, pA=pA, pB=pB, dw_AB=self.dw_AB_struct, dw_BC=self.dw_BC_struct, > kex_AB=kex_AB, kex_BC=kex_BC, kex_AC=kex_AC, > spin_lock_fields=self.spin_lock_omega1, relax_time=self.relax_times, > inv_relax_time=self.inv_relax_times, back_calc=self.back_calc, > num_points=self.num_disp_points) > > # Clean the data for all values, which is left over at the end of > arrays. > self.back_calc = self.back_calc*self.disp_struct > > > _______________________________________________ > 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
