Re: r23737 - /branches/disp_spin_speed/lib/dispersion/cr72.pyRe: r23737 - /branches/disp_spin_speed/lib/dispersion/cr72.pyHi,
That should work. The 'allclose(dw, zeros(dw.shape))' check is for
all dw values being 0.0. But if spin 1 has a value of 0.0 and spin 2
has a value of 0.12345, the mask is not created. I think it would be
better if there are two dw values sent into the function, the original
from the parameter array used only in the check for zero dw values,
and then the full stucture (dw_struct, dw_full, or some other name).
Then the check could be something like:
# Test if dw is zero.
mask_dw_t = False
if min(abs(dw)) == 0.0:
mask_dw_t = True
mask_dw = masked_values(dw_full, 0.0)
This should be much faster as dw will be a short rank-1 array. What
happens inside the 'if' statement doesn't matter as being slow when
one dw value is 0.0 will have almost no effect on optimisation speeds.
Regards,
Edward
On 10 June 2014 14:29, Troels Emtekær Linnet <tlinnet@xxxxxxxxxxxxx> wrote:
Hi Edward. There is a very very smart trick for that. Masked arrays! One can: mask_dw_t = False # Test if dw is zero. if allclose(dw, zeros(dw.shape)): mask_dw_t = True mask_dw = ma.masked_values(dw, 0.0) # Calculate R2eff. R2eff = r20_kex - cpmg_frqs * arccosh( fact ) # Replace data in array. if mask_dw_t: R2eff[mask_dw] = r20a back_calc[:] = R2eff 2014-06-10 13:52 GMT+02:00 Edward d'Auvergne <edward@xxxxxxxxxxxxx>:Hi, I think the 'no Rex' tests and what R2eff is set to needs to be thought out carefully. We have two cases: - The spin independent parameters pA and kex. These will affect the R2eff values for all spins of the cluster. - The spin dependent parameter dw. This will effect the R2eff values for individual spins, and catching this and dealing with it correctly might be much harder for a spin cluster. The first case is simplest, catch single pA and kex values. So numpy operations are not needed there. For the second, I'm not sure how we handle that. But more importantly we also need to work out what to do with the R2eff values in that case. How do we catch dw = 0.0 for one spin but not the others, and then set the R2eff values to R20 for just that spin? But also let the R2eff values be calculated for all other spins. Maybe it is better to not catch dw values and to deal with that value later in the lib.dispersion functions? Regards, Edward On 10 June 2014 12:05, Troels Emtekær Linnet <tlinnet@xxxxxxxxxxxxx> wrote:Here are some hints: http://stackoverflow.com/questions/10580676/comparing-two-numpy-arrays-for-equality It would be the timing of: numpy.allclose http://docs.scipy.org/doc/numpy/reference/generated/numpy.allclose.html numpy.isclose http://docs.scipy.org/doc/numpy/reference/generated/numpy.isclose.html#numpy.isclose numpy.ma http://docs.scipy.org/doc/numpy/reference/generated/numpy.ma.masked_values.html#numpy.ma.masked_values 2014-06-10 11:57 GMT+02:00 Troels Emtekær Linnet <tlinnet@xxxxxxxxxxxxx>:Hi Edward. I thought the same. One could send in a ones and zeros numpy array in function. Set them to standard to None. If they are None, then create them. One can see on the timings, that r2eff_CR72 numeric.py:2056(allclose) are the most time consuming. Maybe there is a faster way for allclose. Masking ??? 1 0.000 0.000 2.084 2.084 <string>:1(<module>) 1 0.002 0.002 2.084 2.084 profiling_cr72.py:441(cluster) 1000 0.002 0.000 2.005 0.002 profiling_cr72.py:405(calc) 1000 0.034 0.000 2.003 0.002 relax_disp.py:995(func_CR72_full) 1000 0.141 0.000 1.960 0.002 relax_disp.py:524(calc_CR72_chi2) 1300 1.100 0.001 1.676 0.001 cr72.py:101(r2eff_CR72) 4300 0.245 0.000 0.507 0.000 numeric.py:2056(allclose) 3000 0.037 0.000 0.164 0.000 shape_base.py:761(tile) 17828 0.126 0.000 0.126 0.000 {method 'reduce' of 'numpy.ufunc' objects} 4000 0.110 0.000 0.110 0.000 {method 'repeat' of 'numpy.ndarray' objects} 8609 0.011 0.000 0.086 0.000 fromnumeric.py:1762(any) 2014-06-10 11:44 GMT+02:00 Edward d'Auvergne <edward@xxxxxxxxxxxxx>:Hi Troels, For speed, maybe you should send in the single value parameters together with the array versions into the lib.dispersion modules? The check: + if np.allclose(dw, np.zeros(dw.shape)) or np.allclose(pA, np.ones(dw.shape)) or np.allclose(kex, np.zeros(dw.shape)): will be very expensive. Regards, Edward On 8 June 2014 19:48, <tlinnet@xxxxxxxxxxxxx> wrote:Author: tlinnet Date: Sun Jun 8 19:48:35 2014 New Revision: 23737 URL: http://svn.gna.org/viewcvs/relax?rev=23737&view=rev Log: Re-implemented safety checks in lib/dispersion/cr72.py. This is now implemented for both rank-1 float array and of higher dimensions. This makes the unit tests pass for multi dimensional computing. Task #7807 (https://gna.org/task/index.php?7807): Speed-up of dispersion models for Clustered analysis. Modified: branches/disp_spin_speed/lib/dispersion/cr72.py Modified: branches/disp_spin_speed/lib/dispersion/cr72.py URL: http://svn.gna.org/viewcvs/relax/branches/disp_spin_speed/lib/dispersion/cr72.py?rev=23737&r1=23736&r2=23737&view=diff ============================================================================== --- branches/disp_spin_speed/lib/dispersion/cr72.py (original) +++ branches/disp_spin_speed/lib/dispersion/cr72.py Sun Jun 8 19:48:35 2014 @@ -128,9 +128,15 @@ rank_1 = False # Catch parameter values that will result in no exchange, returning flat R2eff = R20 lines (when kex = 0.0, k_AB = 0.0). + # For rank-1 float array. if rank_1: if dw == 0.0 or pA == 1.0 or kex == 0.0: back_calc[:] = array([r20a]*num_points) + return + # For higher dimensions, return same structure. + else: + if np.allclose(dw, np.zeros(dw.shape)) or np.allclose(pA, np.ones(dw.shape)) or np.allclose(kex, np.zeros(dw.shape)): + back_calc[:] = r20a return # The B population. @@ -165,16 +171,23 @@ # Catch math domain error of cosh(val > 710). # This is when etapos > 710. - if rank_1: - if max(etapos) > 700: + if max(etapos) > 700: + if rank_1: back_calc[:] = array([r20a]*num_points) + return + # For higher dimensions, return same structure. + else: + back_calc[:] = r20a return # The arccosh argument - catch invalid values. fact = Dpos * cosh(etapos) - Dneg * cos(etaneg) - if rank_1: - if min(fact) < 1.0: + if min(fact) < 1.0: + if rank_1: back_calc[:] = array([r20_kex]*num_points) + return + else: + back_calc[:] = r20_kex return # Calculate R2eff. @@ -182,8 +195,10 @@ # Catch errors, taking a sum over array is the fastest way to check for # +/- inf (infinity) and nan (not a number). - if rank_1: - if not isfinite(sum(R2eff)): + if not isfinite(sum(R2eff)): + if rank_1: R2eff = array([1e100]*num_points) + else: + R2eff = np.ones(R2eff.shape) * 1e100 back_calc[:] = R2eff _______________________________________________ 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_______________________________________________ relax (http://www.nmr-relax.com) This is the relax-devel mailing list relax-devel@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-devel