Hi,
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 <[email protected]> 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 <[email protected]>:
>> 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 <[email protected]> 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 <[email protected]>:
>>>> 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 <[email protected]>:
>>>>
>>>>> 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, <[email protected]> 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
>>>>> > [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
>>>>
>>>>
_______________________________________________
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
