[task #7807] Speed-up of dispersion models for Clustered analysis -- June 19, 2014 - 19:18

Follow-up Comment #325, task #7807 (project relax):

Storage of idea: NS CPMG 2site 3d

Skipping of the spin lloop

    # Loop over the spins
    #for si in range(NS):
    #    # Loop over the spectrometer frequencies.
    for mi in range(NM):
        # Extract number of points.
        num_points_si_mi = int(num_points[0, 0, mi, 0])

        # Loop over the time points, back calculating the R2eff values.
        for di in range(num_points_si_mi):
            # Extract the values from the higher dimensional arrays.
            inv_tcpmg_si_mi_di = inv_tcpmg[0, 0, mi, 0, di]
            power_si_mi_di = int(power[0, 0, mi, 0, di])
            #r20a_si_mi_di = r20a[0, 0, mi, 0, di]

            # Initial magnetisation.
            Mint_i = Mint_mat[:, :, mi, Ellipsis]

            # This matrix is a propagator that will evolve the magnetization
with the matrix R for a delay tcp.
            evolution_matrix_i = evolution_matrix_mat[:, :, mi, Ellipsis]

            # Loop over the CPMG elements, propagating the magnetisation.
            for j in range(power_si_mi_di):
                Mint_i = einsum('...ij,...jk', evolution_matrix_i, Mint_i)

            # Store the result.
            Mint_mat[:, :, mi, Ellipsis] =  Mint_i


    print Mint_mat[:, :, :, :, :, 1, 0]
    print asd

    # The next lines calculate the R2eff using a two-point approximation, i.e.
assuming that the decay is mono-exponential.
    #Mx = Mint_i[1][0] / pA
    # Extract data from ei, si, mi, oi, di, the [1][0]
    Mx = Mint_mat[:, :, :, :, :, 1, 0] / pA

    mask_inv = masked_invalid(Mx)
    if any(mask_inv.mask):
        Mx[mask_inv.mask] = r20a[mask_inv.mask]

    back_calc[:] = Mx


    #if Mx <= 0.0 or isNaN(Mx):
    #    back_calc[0, si, mi, 0, di] = r20a_si_mi_di
    #else:
    #    back_calc[0, si, mi, 0, di] = - inv_tcpmg_si_mi_di * log(Mx)

    # Replace data in array.
    # If dw is zero.
    if t_dw_zero:
        back_calc[mask_dw_zero.mask] = r20a[mask_dw_zero.mask]

    # Catch errors, taking a sum over array is the fastest way to check for
    # +/- inf (infinity) and nan (not a number).
    if not isfinite(sum(back_calc)):
        # Replaces nan, inf, etc. with fill value.
        fix_invalid(back_calc, copy=False, fill_value=1e100)

    _______________________________________________________

Reply to this item at:

  <http://gna.org/task/?7807>

_______________________________________________
  Message sent via/by Gna!
  http://gna.org/