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, <[email protected]> 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 > [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
