mailRe: r24301 - /branches/disp_spin_speed/lib/dispersion/ns_cpmg_2site_3d.py


Others Months | Index by Date | Thread Index
>>   [Date Prev] [Date Next] [Thread Prev] [Thread Next]

Header


Content

Posted by Edward d'Auvergne on June 25, 2014 - 09:07:
This logic will really help collapse this model so the Python looping
can be eliminated completely :)  The linear algebra operations for
shifting the matrix power out of the loops might be much more
complicated though.

Cheers,

Edward


On 25 June 2014 02:14,  <tlinnet@xxxxxxxxxxxxx> wrote:
Author: tlinnet
Date: Wed Jun 25 02:14:48 2014
New Revision: 24301

URL: http://svn.gna.org/viewcvs/relax?rev=24301&view=rev
Log:
Got rid of the inner evolution of the magnetization.

If the looping over the number of CPMG elements is given by the index l, 
and the initial magnetization has
been formed, then the number of times for propagation of magnetization is l 
= power_si_mi_di-1.
If the magnetization matrix "Mint" has the index Mint_(i,k) and the 
evolution matrix has the index Evol_(k,j), i=1, k=7, j=7
then the dot product is given by: Sum_{k=1}^{k} Mint_(1,k) * Evol_(k,j) = 
D_(1, j).
The numpy einsum formula for this would be: einsum('ik,kj -> ij', Mint, 
Evol)

Following evolution will be: Sum_{k=1}^{k} D_(1, j) * Evol_(k,j) = 
Mint_(1,k) * Evol_(k,j) * Evol_(k,j).
We can then realize, that the evolution matrix can be raised to the power 
l. Evol_P = Evol**l.
It will then be: einsum('ik,kj -> ij', Mint, Evol_P)

- Get which power to raise the matrix to.
l = power_si_mi_di-1

- Raise the square evolution matrix to the power l.
evolution_matrix_T_pwer_i = matrix_power(evolution_matrix_T_i, l)

Mint_T_i = dot(Mint_T_i, evolution_matrix_T_pwer_i)
or
Mint_T_i = einsum('ik,kj -> ij', Mint_T_i, evolution_matrix_T_pwer_i)

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_cpmg_2site_3d.py

Modified: branches/disp_spin_speed/lib/dispersion/ns_cpmg_2site_3d.py
URL: 
http://svn.gna.org/viewcvs/relax/branches/disp_spin_speed/lib/dispersion/ns_cpmg_2site_3d.py?rev=24301&r1=24300&r2=24301&view=diff
==============================================================================
--- branches/disp_spin_speed/lib/dispersion/ns_cpmg_2site_3d.py (original)
+++ branches/disp_spin_speed/lib/dispersion/ns_cpmg_2site_3d.py Wed Jun 25 
02:14:48 2014
@@ -56,6 +56,7 @@

 # Python module imports.
 from numpy import array, dot, fabs, float64, einsum, isfinite, log, min, 
multiply, rollaxis, sum, tile
+from numpy.linalg import matrix_power
 from numpy.ma import fix_invalid, masked_where

 # relax module imports.
@@ -331,9 +332,23 @@
                 # This matrix is a propagator that will evolve the 
magnetization with the matrix R for a delay tcp.
                 evolution_matrix_T_i = evolution_matrix_T_mat[0, si, mi, 
0, di]

-                # Loop over the CPMG elements, propagating the 
magnetisation.
-                for j in range(power_si_mi_di - 1):
-                    Mint_T_i = dot(Mint_T_i, evolution_matrix_T_i)
+                # If the looping over the number of CPMG elements is given 
by the index l, and the initial magnetization has
+                # been formed, then the number of times for propagation of 
magnetization is l = power_si_mi_di-1.
+                # If the magnetization matrix "Mint" has the index 
Mint_(i,k) and the evolution matrix has the index Evol_(k,j), i=1, k=7, j=7
+                # then the dot product is given by: Sum_{k=1}^{k} 
Mint_(1,k) * Evol_(k,j) = D_(1, j).
+                # The numpy einsum formula for this would be: 
einsum('ik,kj -> ij', Mint, Evol)
+                #
+                # Following evolution will be: Sum_{k=1}^{k} D_(1, j) * 
Evol_(k,j) = Mint_(1,k) * Evol_(k,j) * Evol_(k,j).
+                # We can then realize, that the evolution matrix can be 
raised to the power l. Evol^l.
+
+                # Get which power to raise the matrix to.
+                l = power_si_mi_di-1
+
+                # Raise the square evolution matrix to the power l.
+                evolution_matrix_T_pwer_i = 
matrix_power(evolution_matrix_T_i, l)
+
+                #Mint_T_i = dot(Mint_T_i, evolution_matrix_T_pwer_i)
+                Mint_T_i = einsum('ik,kj -> ij', Mint_T_i, 
evolution_matrix_T_pwer_i)

                 # The next lines calculate the R2eff using a two-point 
approximation, i.e. assuming that the decay is mono-exponential.
                 Mx = Mint_T_i[0][1] / pA


_______________________________________________
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



Related Messages


Powered by MHonArc, Updated Wed Jun 25 09:20:18 2014