Hi Troels,
Here again you can avoid the divide by zero problems by having:
"""
from numpy import max, abs
if min(abs(denom)) == 0:
back_calc[:] = array([1e100]*num_points)
"""
Then these tests should then pass when running with --numpy-raise.
You should also not see any significant speed changes (at least I
didn't with a quick test).
Regards,
Edward
On 19 May 2014 01:19, <[email protected]> wrote:
> Author: tlinnet
> Date: Mon May 19 01:19:02 2014
> New Revision: 23221
>
> URL: http://svn.gna.org/viewcvs/relax?rev=23221&view=rev
> Log:
> Speed-up of model TP02.
>
> task #7793: (https://gna.org/task/?7793) Speed-up of dispersion models.
>
> The change for running systemtest is:
> test_curve_type_r1rho_fixed_time
> 0.057s -> 0.049s
>
> test_tp02_data_to_ns_r1rho_2site
> 10.539s -> 10.456s
>
> test_tp02_data_to_tp02
> 8.608s -> 5.727s
>
> This is won by not checking single values in the R1rho array for math domain
> errors, but calculating all steps, and in one single round check for finite
> values.
> If just one non-finite value is found, the whole array is returned with a
> large
> penalty of 1e100.
>
> This makes all calculations be the fastest numpy array way.
>
> Modified:
> branches/disp_speed/lib/dispersion/tp02.py
>
> Modified: branches/disp_speed/lib/dispersion/tp02.py
> URL:
> http://svn.gna.org/viewcvs/relax/branches/disp_speed/lib/dispersion/tp02.py?rev=23221&r1=23220&r2=23221&view=diff
> ==============================================================================
> --- branches/disp_speed/lib/dispersion/tp02.py (original)
> +++ branches/disp_speed/lib/dispersion/tp02.py Mon May 19 01:19:02 2014
> @@ -60,7 +60,7 @@
> """
>
> # Python module imports.
> -from math import atan2, sin
> +from numpy import arctan2, isfinite, sin, sum
>
>
> def r1rho_TP02(r1rho_prime=None, omega=None, offset=None, pA=None, pB=None,
> dw=None, kex=None, R1=0.0, spin_lock_fields=None, spin_lock_fields2=None,
> back_calc=None, num_points=None):
> @@ -110,34 +110,31 @@
> # The numerator.
> numer = pA * pB * dw**2 * kex
>
> - # Loop over the dispersion points, back calculating the R1rho values.
> + # We assume that A resonates at 0 [s^-1], without loss of generality.
> + waeff2 = spin_lock_fields2 + da2 # Effective field at A.
> + wbeff2 = spin_lock_fields2 + db2 # Effective field at B.
> + weff2 = spin_lock_fields2 + d2 # Effective field at pop-average.
> +
> + # The rotating frame flip angle.
> + theta = arctan2(spin_lock_fields, d)
> +
> + # Repetitive calculations (to speed up calculations).
> + sin_theta2 = sin(theta)**2
> + R1_cos_theta2 = R1 * (1.0 - sin_theta2)
> + R1rho_prime_sin_theta2 = r1rho_prime * sin_theta2
> +
> + # Denominator.
> + denom = waeff2 * wbeff2 / weff2 + kex2
> + #denom_extended = waeff2*wbeff2/weff2+kex2-2*sin_theta2*pA*pB*dw**2
> +
> + # R1rho calculation.
> + R1rho = R1_cos_theta2 + R1rho_prime_sin_theta2 + sin_theta2 * numer /
> denom
> +
> + # 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(R1rho)):
> + R1rho = array([1e100]*num_points)
> +
> + # Parse back the value to update the back_calc class object.
> for i in range(num_points):
> - # We assume that A resonates at 0 [s^-1], without loss of generality.
> - waeff2 = spin_lock_fields2[i] + da2 # Effective field at A.
> - wbeff2 = spin_lock_fields2[i] + db2 # Effective field at B.
> - weff2 = spin_lock_fields2[i] + d2 # Effective field at
> pop-average.
> -
> - # The rotating frame flip angle.
> - theta = atan2(spin_lock_fields[i], d)
> -
> - # Repetitive calculations (to speed up calculations).
> - sin_theta2 = sin(theta)**2
> - R1_cos_theta2 = R1 * (1.0 - sin_theta2)
> - R1rho_prime_sin_theta2 = r1rho_prime * sin_theta2
> -
> - # Catch zeros (to avoid pointless mathematical operations).
> - if numer == 0.0:
> - back_calc[i] = R1_cos_theta2 + R1rho_prime_sin_theta2
> - continue
> -
> - # Denominator.
> - denom = waeff2 * wbeff2 / weff2 + kex2
> - #denom_extended = waeff2*wbeff2/weff2+kex2-2*sin_theta2*pA*pB*dw**2
> -
> - # Avoid divide by zero.
> - if denom == 0.0:
> - back_calc[i] = 1e100
> - continue
> -
> - # R1rho calculation.
> - back_calc[i] = R1_cos_theta2 + R1rho_prime_sin_theta2 + sin_theta2 *
> numer / denom
> + back_calc[i] = R1rho[i]
>
>
> _______________________________________________
> 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
