Re: r24004 - /branches/disp_spin_speed/lib/dispersion/ns_mmq_2site.py -- June 17, 2014 - 10:24

Hi Troels,

Now that you have converted almost all data structures to numpy
arrays, you should change how these are accessed.  For example a[1, 1]
is cleaner than a[1][1] and is the standard way to index a numpy
array.  It is also faster:

$ python -m timeit -n 1000 -s "import numpy" "a = numpy.eye(5)" "for i
in xrange(10000): a[1][1] = 5.0"
1000 loops, best of 3: 2.2 msec per loop
$ python -m timeit -n 1000 -s "import numpy" "a = numpy.eye(5)" "for i
in xrange(10000): a[1, 1] = 5.0"
1000 loops, best of 3: 838 usec per loop
$

Though this will probably have almost zero effect on the speed of the
target function, it will be more correct.

Regards,

Edward





On 16 June 2014 22:11,  <tlinnet@xxxxxxxxxxxxx> wrote:
> Author: tlinnet
> Date: Mon Jun 16 22:11:40 2014
> New Revision: 24004
>
> URL: http://svn.gna.org/viewcvs/relax?rev=24004&view=rev
> Log:
> Fixed the use of higher dimensional data in mmq 2site sq dq zq.
>
> 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=24004&r1=24003&r2=24004&view=diff
> ==============================================================================
> --- branches/disp_spin_speed/lib/dispersion/ns_mmq_2site.py     (original)
> +++ branches/disp_spin_speed/lib/dispersion/ns_mmq_2site.py     Mon Jun 16 
> 22:11:40 2014
> @@ -281,26 +281,42 @@
>      @type power:            numpy int16, rank-1 array
>      """
>
> -    # Populate the m1 and m2 matrices (only once per function call for 
> speed).
> -    populate_matrix(matrix=m1, R20A=R20A, R20B=R20B, dw=dw, k_AB=k_AB, 
> k_BA=k_BA)
> -    populate_matrix(matrix=m2, R20A=R20A, R20B=R20B, dw=-dw, k_AB=k_AB, 
> k_BA=k_BA)
> -
> -    # Loop over the time points, back calculating the R2eff values.
> -    for i in range(num_points):
> -        # The A+/- matrices.
> -        A_pos = matrix_exponential(m1*tcp[i])
> -        A_neg = matrix_exponential(m2*tcp[i])
> -
> -        # The evolution for one n.
> -        evol_block = dot(A_pos, dot(A_neg, dot(A_neg, A_pos)))
> -
> -        # The full evolution.
> -        evol = square_matrix_power(evol_block, power[i])
> -
> -        # The next lines calculate the R2eff using a two-point 
> approximation, i.e. assuming that the decay is mono-exponential.
> -        Mx = dot(F_vector, dot(evol, M0))
> -        Mx = Mx.real
> -        if Mx <= 0.0 or isNaN(Mx):
> -            back_calc[i] = 1e99
> -        else:
> -            back_calc[i] = -inv_tcpmg[i] * log(Mx / pA)
> +
> +    # Extract shape of experiment.
> +    NS, NM, NO = num_points.shape
> +
> +    # Loop over spins.
> +    for si in range(NS):
> +        # Loop over the spectrometer frequencies.
> +        for mi in range(NM):
> +            # Loop over offsets:
> +            for oi in range(NO):
> +
> +                r20a_si_mi_oi = R20A[si][mi][oi][0]
> +                r20b_si_mi_oi = R20B[si][mi][oi][0]
> +                dw_si_mi_oi = dw[si][mi][oi][0]
> +                num_points_si_mi_oi = num_points[si][mi][oi]
> +
> +                # Populate the m1 and m2 matrices (only once per function 
> call for speed).
> +                populate_matrix(matrix=m1, R20A=r20a_si_mi_oi , 
> R20B=r20b_si_mi_oi, dw=dw_si_mi_oi, k_AB=k_AB, k_BA=k_BA)
> +                populate_matrix(matrix=m2, R20A=r20a_si_mi_oi , 
> R20B=r20b_si_mi_oi, dw=-dw_si_mi_oi, k_AB=k_AB, k_BA=k_BA)
> +
> +                # Loop over the time points, back calculating the R2eff 
> values.
> +                for i in range(num_points_si_mi_oi):
> +                    # The A+/- matrices.
> +                    A_pos = matrix_exponential(m1*tcp[si][mi][oi][i])
> +                    A_neg = matrix_exponential(m2*tcp[si][mi][oi][i])
> +
> +                    # The evolution for one n.
> +                    evol_block = dot(A_pos, dot(A_neg, dot(A_neg, A_pos)))
> +
> +                    # The full evolution.
> +                    evol = square_matrix_power(evol_block, 
> power[si][mi][oi][i])
> +
> +                    # The next lines calculate the R2eff using a two-point 
> approximation, i.e. assuming that the decay is mono-exponential.
> +                    Mx = dot(F_vector, dot(evol, M0))
> +                    Mx = Mx.real
> +                    if Mx <= 0.0 or isNaN(Mx):
> +                        back_calc[si][mi][oi][i] = 1e99
> +                    else:
> +                        back_calc[si][mi][oi][i] = 
> -inv_tcpmg[si][mi][oi][i] * log(Mx / 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