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/