Hi, I have a strong feeling that this can be taken even further. I think that Dmitry Korzhnev's cpmg_fit software already does this, though I am not sure. It will require careful study of the formula for magentisation evolution and how to write it in a form to be used outside of the ND loop. It may require a specially constructed structure for the dispersion points and a matrix exponential. I'm not sure how to do this yet, but I think that it will be possible to eliminate all loops in the numeric models.
Regards, Edward On 23 June 2014 10:15, Troels Emtekær Linnet <[email protected]> wrote: > Hi Ed. > > I think we can kill the for loops over the spins. > > If I find a way to access the array for all spins, and then only loop over > mi, and di, to make the dot product according to the power of the tau_cpmg, > then > we have saved the costly clustering for loop. > > I have been looking into strides: > http://wiki.nmr-relax.com/Numpy_linalg#Stride_tricks > > But I haven't figured it out yet. > > Best > Troels > > > 2014-06-23 10:03 GMT+02:00 Edward d'Auvergne <[email protected]>: >> >> Interesting! If the magnetisation propagation can also be shifted out >> of the loops, then the loops will disappear and the numeric models >> will become far quicker. With a matrix exponential, maybe this is >> possible. I'll think about this for a while. >> >> Regards, >> >> Edward >> >> >> On 19 June 2014 20:52, <[email protected]> wrote: >> > Author: tlinnet >> > Date: Thu Jun 19 20:52:57 2014 >> > New Revision: 24171 >> > >> > URL: http://svn.gna.org/viewcvs/relax?rev=24171&view=rev >> > Log: >> > Moved the costly calculation of the matrix exponential out of for loops. >> > >> > It was the numpy.eig and numpy.inv which was draining power. >> > >> > This speeds up model NS R1rho 2site, by a factor 4X: >> > BEFORE: >> > Single: >> > ncalls tottime percall cumtime percall filename:lineno(function) >> > 1 0.000 0.000 32.552 32.552 <string>:1(<module>) >> > 1 0.002 0.002 32.552 32.552 >> > pf_nsr1rho2site:530(single) >> > Cluster: >> > ncalls tottime percall cumtime percall filename:lineno(function) >> > 1 0.000 0.000 33.307 33.307 <string>:1(<module>) >> > 1 0.008 0.008 33.307 33.307 >> > pf_nsr1rho2site:554(cluster) >> > >> > AFTER: >> > Single: >> > ncalls tottime percall cumtime percall filename:lineno(function) >> > 1 0.000 0.000 8.286 8.286 <string>:1(<module>) >> > 1 0.002 0.002 8.286 8.286 >> > pf_nsr1rho2site:530(single) >> > Cluster: >> > ncalls tottime percall cumtime percall filename:lineno(function) >> > 1 0.000 0.000 8.223 8.223 <string>:1(<module>) >> > 1 0.007 0.007 8.223 8.223 >> > pf_nsr1rho2site:554(cluster) >> > >> > 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=24171&r1=24170&r2=24171&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 20:52:57 2014 >> > @@ -56,7 +56,7 @@ >> > # relax module imports. >> > from lib.dispersion.ns_matrices import rr1rho_3d, rr1rho_3d_rankN >> > from lib.float import isNaN >> > -from lib.linear_algebra.matrix_exponential import matrix_exponential >> > +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): >> > @@ -106,6 +106,9 @@ >> > # The matrix that contains all the contributions to the evolution, >> > i.e. relaxation, exchange and chemical shift evolution. >> > R_mat = rr1rho_3d_rankN(R1=r1, r1rho_prime=r1rho_prime, pA=pA, >> > pB=pB, dw=dw, omega=omega, offset=offset, w1=spin_lock_fields, k_AB=k_AB, >> > k_BA=k_BA, relax_time=relax_time) >> > >> > + # This matrix is a propagator that will evolve the magnetization >> > with the matrix R. >> > + Rexpo_mat = matrix_exponential_rankN(R_mat) >> > + >> > # Loop over spins. >> > for si in range(NS): >> > # Loop over the spectrometer frequencies. >> > @@ -135,19 +138,16 @@ >> > >> > # Loop over the time points, back calculating the R2eff >> > values. >> > for j in range(num_points_i): >> > - # The matrix that contains all the contributions to >> > the evolution, i.e. relaxation, exchange and chemical shift evolution. >> > - R_mat_i = R_mat[0, si, mi, oi, j] >> > - >> > # The following lines rotate the magnetization >> > previous to spin-lock into the weff frame. >> > theta = atan2(spin_lock_fields_i[j], dA) >> > M0[0] = sin(theta) # The A state initial X >> > magnetisation. >> > M0[2] = cos(theta) # The A state initial Z >> > magnetisation. >> > >> > # This matrix is a propagator that will evolve the >> > magnetization with the matrix R. >> > - Rexpo = matrix_exponential(R_mat_i) >> > + Rexpo_i = Rexpo_mat[0, si, mi, oi, j] >> > >> > # Magnetization evolution. >> > - MA = dot(M0, dot(Rexpo, M0)) >> > + MA = dot(M0, dot(Rexpo_i, M0)) >> > >> > # 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): >> > >> > >> > _______________________________________________ >> > 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 > > _______________________________________________ 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
