Re: r24173 - /branches/disp_spin_speed/lib/dispersion/ns_r1rho_2site.py -- June 23, 2014

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,  <tlinnet@xxxxxxxxxxxxx> 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
relax-commits@xxxxxxx

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

Re: r24173 - /branches/disp_spin_speed/lib/dispersion/ns_r1rho_2site.py

Header

Content

Related Messages