The checking code should have been in a different commit - or even better a different commit adding this as a unit test. Also, shift the fixed matrices out of the function for speed.
Cheers, Edward On 20 June 2014 09:18, <[email protected]> wrote: > Author: tlinnet > Date: Fri Jun 20 09:18:00 2014 > New Revision: 24189 > > URL: http://svn.gna.org/viewcvs/relax?rev=24189&view=rev > Log: > Implemented the collection of the multidimensional matrix m1 and m2 in model > ns mmq 2site. > > Inserted also a check, that the newly computed matrix is equal. > They are, to the 6 digit. > > 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_mmq_2site.py > > Modified: branches/disp_spin_speed/lib/dispersion/ns_mmq_2site.py > URL: > http://svn.gna.org/viewcvs/relax/branches/disp_spin_speed/lib/dispersion/ns_mmq_2site.py?rev=24189&r1=24188&r2=24189&view=diff > ============================================================================== > --- branches/disp_spin_speed/lib/dispersion/ns_mmq_2site.py (original) > +++ branches/disp_spin_speed/lib/dispersion/ns_mmq_2site.py Fri Jun 20 > 09:18:00 2014 > @@ -51,7 +51,7 @@ > > # Python module imports. > from math import floor > -from numpy import array, conj, dot, float64, log > +from numpy import array, conj, complex64, dot, float64, log, multiply, sum > > # relax module imports. > from lib.float import isNaN > @@ -81,6 +81,76 @@ > matrix[0, 1] = k_BA > matrix[1, 0] = k_AB > matrix[1, 1] = -k_BA + 1.j*dw - R20B > + > + > +def populate_matrix_rankN(R20A=None, R20B=None, dw=None, k_AB=None, > k_BA=None, tcp=None): > + """The Bloch-McConnell matrix for 2-site exchange, for rank > [NE][NS][NM][NO][ND][2][2]. > + > + @keyword R20A: The transverse, spin-spin relaxation rate for > state A. > + @type R20A: numpy float array of rank [NE][NS][NM][NO][ND] > + @keyword R20B: The transverse, spin-spin relaxation rate for > state B. > + @type R20B: numpy float array of rank [NE][NS][NM][NO][ND] > + @keyword dw: The combined chemical exchange difference > parameters between states A and B in rad/s. This can be any combination of > dw and dwH. > + @type dw: numpy float array of rank [NE][NS][NM][NO][ND] > + @keyword k_AB: The rate of exchange from site A to B (rad/s). > + @type k_AB: float > + @keyword k_BA: The rate of exchange from site B to A (rad/s). > + @type k_BA: float > + @keyword tcp: The tau_CPMG times (1 / 4.nu1). > + @type tcp: numpy float array of rank [NE][NS][NM][NO][ND] > + @return: The relaxation matrix. > + @rtype: numpy float array of rank > [NE][NS][NM][NO][ND][2][2] > + """ > + > + # Pre-multiply with tcp. > + r20a_tcp = R20A * tcp > + r20b_tcp = R20B * tcp > + k_AB_tcp = k_AB * tcp > + k_BA_tcp = k_BA * tcp > + # Complex dw. > + dw_tcp_C = dw * tcp * 1j > + > + # Fill in the elements. > + #matrix[0, 0] = -k_AB - R20A > + #matrix[0, 1] = k_BA > + #matrix[1, 0] = k_AB > + #matrix[1, 1] = -k_BA + 1.j*dw - R20B > + > + m_r20a = array([ > + [-1.0, 0.0], > + [0.0, 0.0],], complex64) > + > + m_r20b = array([ > + [0.0, 0.0], > + [0.0, -1.0],], complex64) > + > + m_k_AB = array([ > + [-1.0, 0.0], > + [1.0, 0.0],], complex64) > + > + m_k_BA = array([ > + [0.0, 1.0], > + [0.0, -1.0],], complex64) > + > + m_dw = array([ > + [0.0, 0.0], > + [0.0, 1.0],], complex64) > + > + # Multiply and expand. > + m_r20a_tcp = multiply.outer( r20a_tcp, m_r20a ) > + m_r20b_tcp = multiply.outer( r20b_tcp, m_r20b ) > + > + # Multiply and expand. > + m_k_AB_tcp = multiply.outer( k_AB_tcp, m_k_AB ) > + m_k_BA_tcp = multiply.outer( k_BA_tcp, m_k_BA ) > + > + # Multiply and expand. > + m_dw_tcp_C = multiply.outer( dw_tcp_C, m_dw ) > + > + # Collect matrix. > + matrix = (m_r20a_tcp + m_r20b_tcp + m_k_AB_tcp + m_k_BA_tcp + m_dw_tcp_C) > + > + return matrix > > > def r2eff_ns_mmq_2site_mq(M0=None, F_vector=array([1, 0], float64), m1=None, > m2=None, R20A=None, R20B=None, pA=None, dw=None, dwH=None, kex=None, > inv_tcpmg=None, tcp=None, back_calc=None, num_points=None, power=None): > @@ -141,6 +211,10 @@ > # Extract shape of experiment. > NS, NM, NO = num_points.shape > > + # Populate the m1 and m2 matrices (only once per function call for > speed). > + m1_mat = populate_matrix_rankN(R20A=R20A, R20B=R20B, dw=-dw - dwH, > k_AB=k_AB, k_BA=k_BA, tcp=tcp) > + m2_mat = populate_matrix_rankN(R20A=R20A, R20B=R20B, dw=dw - dwH, > k_AB=k_AB, k_BA=k_BA, tcp=tcp) > + > # Loop over spins. > for si in range(NS): > # Loop over the spectrometer frequencies. > @@ -160,6 +234,23 @@ > > # Loop over the time points, back calculating the R2eff > values. > for i in range(num_points_i): > + m1_mat_i = m1_mat[si, mi, oi, i] > + m2_mat_i = m2_mat[si, mi, oi, i] > + > + diff_m1 = abs(sum(m1*tcp[si, mi, oi, i] - m1_mat_i)) > + if diff_m1 > 1.0e-06: > + print diff_m1 > + print m1*tcp[si, mi, oi, i] > + print m1_mat_i > + print asd > + > + diff_m2 = abs(sum(m2*tcp[si, mi, oi, i] - m2_mat_i)) > + if diff_m2 > 1.0e-06: > + print diff_m2 > + print m2*tcp[si, mi, oi, i] > + print m2_mat_i > + print asd > + > # The M1 and M2 matrices. > M1 = matrix_exponential(m1*tcp[si, mi, oi, i]) # > Equivalent to D+. > M2 = matrix_exponential(m2*tcp[si, mi, oi, i]) # > Equivalent to Z-. > > > _______________________________________________ > 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
