Hi Troels, Here are some tips for more speed here - shift many calculations out of the loops, i.e. dA. It should be possible to shift theta out as well. And sin(theta) and cos(theta). Maybe even building up M0. Such changes will be of benefit in eventually eliminating all looping from this numeric model. Though the bottleneck is clearly the magnetisation evolution dot products, these changes will make the model slightly faster and it will be a great infrastructure improvement.
Then once a matrix exponential is used for magetisation evolutions, maybe: MA = (Rexpo . M0 . MA)**(di + 1) and all looping eliminated, this model will be lightning fast and much closer to the analytic models rather than the current ~300 times slower speed. Cheers, Edward On 19 June 2014 21:05, <[email protected]> wrote: > Author: tlinnet > Date: Thu Jun 19 21:05:51 2014 > New Revision: 24173 > > URL: http://svn.gna.org/viewcvs/relax?rev=24173&view=rev > Log: > Cleaned up the code of NS R1rho 2site, and removed the matrix argument to the > function. > > 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=24173&r1=24172&r2=24173&view=diff > ============================================================================== > --- branches/disp_spin_speed/lib/dispersion/ns_r1rho_2site.py (original) > +++ branches/disp_spin_speed/lib/dispersion/ns_r1rho_2site.py Thu Jun 19 > 21:05:51 2014 > @@ -59,7 +59,7 @@ > from lib.linear_algebra.matrix_exponential import matrix_exponential, > matrix_exponential_rankN > > > -def ns_r1rho_2site(M0=None, matrix=None, r1rho_prime=None, omega=None, > offset=None, r1=0.0, pA=None, dw=None, kex=None, spin_lock_fields=None, > relax_time=None, inv_relax_time=None, back_calc=None, num_points=None): > +def ns_r1rho_2site(M0=None, r1rho_prime=None, omega=None, offset=None, > r1=0.0, pA=None, dw=None, kex=None, spin_lock_fields=None, relax_time=None, > inv_relax_time=None, back_calc=None, num_points=None): > """The 2-site numerical solution to the Bloch-McConnell equation for > R1rho data. > > This function calculates and stores the R1rho values. > @@ -67,8 +67,6 @@ > > @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 > - @keyword matrix: A numpy array to be populated to create the > evolution matrix. > - @type matrix: numpy rank-2, 6D float64 array > @keyword r1rho_prime: The R1rho_prime parameter value (R1rho with > no exchange). > @type r1rho_prime: numpy float array of rank [NS][NM][NO][ND] > @keyword omega: The chemical shift for the spin in rad/s. > @@ -115,31 +113,17 @@ > for mi in range(NM): > # Loop over offsets: > for oi in range(NO): > - # Extract parameters from array. > - omega_i = omega[0, si, mi, oi, 0] > - offset_i = offset[0, si, mi, oi, 0] > - r1_i = r1[0, si, mi, oi, 0] > - dw_i = dw[0, si, mi, oi, 0] > - > - r1rho_prime_i = r1rho_prime[0, si, mi, oi] > - spin_lock_fields_i = spin_lock_fields[0, si, mi, oi] > - relax_time_i = relax_time[0, si, mi, oi] > - inv_relax_time_i = inv_relax_time[0, si, mi, oi] > - back_calc_i = back_calc[0, si, mi, oi] > + # Extract number of points. > num_points_i = num_points[0, si, mi, oi] > > # Repetitive calculations (to speed up calculations). > - Wa = omega_i # Larmor frequency [s^-1]. > - Wb = omega_i + dw_i # Larmor frequency [s^-1]. > - W = pA*Wa + pB*Wb # Population-averaged Larmor > frequency [s^-1]. > - dA = Wa - offset_i # Offset of spin-lock from A. > - dB = Wb - offset_i # Offset of spin-lock from B. > - d = W - offset_i # Offset of spin-lock from > population-average. > + # Offset of spin-lock from A. > + dA = omega[0, si, mi, oi, 0] - offset[0, si, mi, oi, 0] > > # Loop over the time points, back calculating the R2eff > values. > for j in range(num_points_i): > # The following lines rotate the magnetization previous > to spin-lock into the weff frame. > - theta = atan2(spin_lock_fields_i[j], dA) > + theta = atan2(spin_lock_fields[0, si, mi, oi, j], dA) > M0[0] = sin(theta) # The A state initial X > magnetisation. > M0[2] = cos(theta) # The A state initial Z > magnetisation. > > @@ -147,10 +131,11 @@ > Rexpo_i = Rexpo_mat[0, si, mi, oi, j] > > # Magnetization evolution. > - MA = dot(M0, dot(Rexpo_i, M0)) > + MA = dot(Rexpo_i, M0) > + MA = dot(M0, MA) > > # 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_i[j] * > log(MA) > + back_calc[0, si, mi, oi, j]= -inv_relax_time[0, si, > mi, oi, j] * log(MA) > > > _______________________________________________ > 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
