Hi, That's very useful to know. I didn't know that numpy could switch the axes like this, performing a localised transform. I actually implemented my own slow version for rank-4 tensors for the frame order theory in the lib.linear_algebra.kronecker_product.transpose_*() functions, but this should be significantly faster.
Regards, Edward On 25 June 2014 02:14, <[email protected]> wrote: > 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) > > > _______________________________________________ > 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
