Let me reinforce that again if you didn't quite get it yet ;) In the 'disp_speed' branch, there are a lot of cases where an array of 1e100 is returned when an array of R20 should be returned. Unit tests would instantly show these. If you put a catch at the start for all parameter values that result in no chemical exchange, then you will avoid most problems seen with the --numpy-raise flag, and the optimisation space you see with the dx.map user function might improve a lot. And these single value checks will be faster than the numpy array checks further into the functions, most of which will become redundant.
Regards, Edward On 21 May 2014 19:10, Edward d'Auvergne <[email protected]> wrote: > Hi Troels, > > You should really consider adding this check to absolutely all > lib.dispersion modules. And as close to the first line as possible > (for most it should be on line one, even if you have to replicate some > calculations just for the 'No Rex' case). It should magically make > most problems go away! Especially in the numeric dispersion models. > And then many checks will no longer be necessary > (https://gna.org/bugs/?22065 and most errors you see in the test suite > with --numpy-raise will disappear). > > Regards, > > Edward > > > On 21 May 2014 08:48, Edward d'Auvergne <[email protected]> wrote: >> Hi Troels, >> >> Actually, if you add that check to the start of all lib.dispersion modules: >> >> # Catch parameter values that will result in no exchange, returning >> flat R2eff = R20 lines. >> if dw == 0.0 or pA == 1.0 or kex == 0.0: >> return array([R20]*num_points) >> >> it would then be worth checking if the other tests are necessary. >> This one check might eliminate the need for most of the other checks >> and the less checks there are, the faster the code will be. >> >> Regards, >> >> Edward >> >> >> >> On 21 May 2014 08:37, Edward d'Auvergne <[email protected]> wrote: >>> Hi Troels, >>> >>> With some of these, you have to be very careful that the parameter >>> combination actually does produce infinite R2eff values. I think in >>> this combination, it might not be the case. What input parameter >>> value causes weff2 to be zero? To properly check, one needs to take >>> the original equation and remove the parameter that is zero, >>> collapsing the equations. What looks like a divide by zero in the >>> original equation might actually not be because of the equation >>> simplifications. For example in almost all cases of dw being zero, >>> the R2eff equations should collapse to R20. Or in this R1rho case, to >>> R1.cos^2(theta) + R1rho'.sin^2(theta). I.e. dw=0 means there is no >>> exchange. So you should never set R2eff to 1e100 in such a case, as >>> you will no longer be able to produce the flat exchange-free R2eff >>> lines. >>> >>> Note, you'll see that in some of my original checks in the >>> lib.dispersion modules, when I am looking for a numerator that is zero >>> before looking for a denominator that is zero. You can see that in >>> the original lib.dispersion.tp02 module. But check has been removed >>> from the disp_speed branch. Without such checks, the dispersion >>> models will incorrectly model no exchange as being almost infinite >>> exchange. It would be worth looking at all checks that I carefully >>> added (they are still in the trunk), because every last one was very >>> carefully thought out with the simplification it would make to the >>> equations. So every check should be in the disp_speed branch too. >>> The checks I have added are not complete though. >>> >>> Maybe there is a far better way to perform some of these checks. For >>> example, catching parameter values that will result in zero exchange: >>> >>> # Catch parameter values that will result in no exchange, returning >>> flat R2eff = R20 lines. >>> if dw == 0.0 or pA == 1.0 or kex == 0.0: >>> return array([R20]*num_points) >>> >>> This could happen at the very start of all of the lib.dispersion >>> modules. Such a check would eliminate many of the numpy warnings >>> where R2eff values end up being Inf, NaN or zero. Actually, I feel >>> rather stupid now for not having had thought of such a brilliant check >>> before :S >>> >>> Regards, >>> >>> Edward >>> >>> >>> >>> On 20 May 2014 22:29, <[email protected]> wrote: >>>> Author: tlinnet >>>> Date: Tue May 20 22:29:40 2014 >>>> New Revision: 23270 >>>> >>>> URL: http://svn.gna.org/viewcvs/relax?rev=23270&view=rev >>>> Log: >>>> Math-domain catching for model TP02. >>>> >>>> task #7793: (https://gna.org/task/?7793) Speed-up of dispersion models. >>>> >>>> This is to implement catching of math domain errors, before they occur. >>>> These can be found via the --numpy-raise function to the systemtests. >>>> >>>> To make the code look clean, the class object "back_calc" is no longer >>>> being updated per time point, but is updated in the relax_disp target >>>> function in >>>> one go. >>>> >>>> Modified: >>>> branches/disp_speed/lib/dispersion/tp02.py >>>> branches/disp_speed/target_functions/relax_disp.py >>>> >>>> Modified: branches/disp_speed/lib/dispersion/tp02.py >>>> URL: >>>> http://svn.gna.org/viewcvs/relax/branches/disp_speed/lib/dispersion/tp02.py?rev=23270&r1=23269&r2=23270&view=diff >>>> ============================================================================== >>>> --- branches/disp_speed/lib/dispersion/tp02.py (original) >>>> +++ branches/disp_speed/lib/dispersion/tp02.py Tue May 20 22:29:40 2014 >>>> @@ -60,10 +60,10 @@ >>>> """ >>>> >>>> # Python module imports. >>>> -from numpy import arctan2, array, isfinite, sin, sum >>>> +from numpy import abs, arctan2, array, isfinite, min, 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): >>>> +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, num_points=None): >>>> """Calculate the R1rho' values for the TP02 model. >>>> >>>> See the module docstring for details. This is the Trott and Palmer >>>> (2002) equation according to Korzhnev (J. Biomol. NMR (2003), 26, 39-48). >>>> @@ -89,9 +89,7 @@ >>>> @type spin_lock_fields: numpy rank-1 float array >>>> @keyword spin_lock_fields2: The R1rho spin-lock field strengths >>>> squared (in rad^2.s^-2). This is for speed. >>>> @type spin_lock_fields2: numpy rank-1 float array >>>> - @keyword back_calc: The array for holding the back calculated >>>> R1rho values. Each element corresponds to one of the spin-lock fields. >>>> - @type back_calc: numpy rank-1 float array >>>> - @keyword num_points: The number of points on the dispersion >>>> curve, equal to the length of the spin_lock_fields and back_calc arguments. >>>> + @keyword num_points: The number of points on the dispersion >>>> curve, equal to the length of the spin_lock_fields. >>>> @type num_points: int >>>> """ >>>> >>>> @@ -115,6 +113,13 @@ >>>> wbeff2 = spin_lock_fields2 + db2 # Effective field at B. >>>> weff2 = spin_lock_fields2 + d2 # Effective field at >>>> pop-average. >>>> >>>> + # Catch math domain error of dividing with 0. >>>> + # This is when weff2 = 0. >>>> + if min(abs(weff2)) == 0: >>>> + R2eff = array([1e100]*num_points) >>>> + >>>> + return R2eff >>>> + >>>> # The rotating frame flip angle. >>>> theta = arctan2(spin_lock_fields, d) >>>> >>>> @@ -127,6 +132,13 @@ >>>> denom = waeff2 * wbeff2 / weff2 + kex2 >>>> #denom_extended = waeff2*wbeff2/weff2+kex2-2*sin_theta2*pA*pB*dw**2 >>>> >>>> + # Catch math domain error of dividing with 0. >>>> + # This is when denom=0. >>>> + if min(abs(denom)) == 0: >>>> + R1rho = array([1e100]*num_points) >>>> + >>>> + return R1rho >>>> + >>>> # R1rho calculation. >>>> R1rho = R1_cos_theta2 + R1rho_prime_sin_theta2 + sin_theta2 * numer / >>>> denom >>>> >>>> @@ -135,6 +147,4 @@ >>>> 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): >>>> - back_calc[i] = R1rho[i] >>>> + return R1rho >>>> >>>> Modified: branches/disp_speed/target_functions/relax_disp.py >>>> URL: >>>> http://svn.gna.org/viewcvs/relax/branches/disp_speed/target_functions/relax_disp.py?rev=23270&r1=23269&r2=23270&view=diff >>>> ============================================================================== >>>> --- branches/disp_speed/target_functions/relax_disp.py (original) >>>> +++ branches/disp_speed/target_functions/relax_disp.py Tue May 20 >>>> 22:29:40 2014 >>>> @@ -1875,7 +1875,7 @@ >>>> # Loop over the offsets. >>>> for oi in range(self.num_offsets[0][si][mi]): >>>> # Back calculate the R1rho values. >>>> - r1rho_TP02(r1rho_prime=R20[r20_index], >>>> omega=self.chemical_shifts[0][si][mi], offset=self.offset[0][si][mi][oi], >>>> pA=pA, pB=pB, dw=dw_frq, kex=kex, R1=self.r1[si, mi], >>>> spin_lock_fields=self.spin_lock_omega1[0][mi][oi], >>>> spin_lock_fields2=self.spin_lock_omega1_squared[0][mi][oi], >>>> back_calc=self.back_calc[0][si][mi][oi], >>>> num_points=self.num_disp_points[0][si][mi][oi]) >>>> + self.back_calc[0][si][mi][oi] = >>>> r1rho_TP02(r1rho_prime=R20[r20_index], >>>> omega=self.chemical_shifts[0][si][mi], offset=self.offset[0][si][mi][oi], >>>> pA=pA, pB=pB, dw=dw_frq, kex=kex, R1=self.r1[si, mi], >>>> spin_lock_fields=self.spin_lock_omega1[0][mi][oi], >>>> spin_lock_fields2=self.spin_lock_omega1_squared[0][mi][oi], >>>> num_points=self.num_disp_points[0][si][mi][oi]) >>>> >>>> # For all missing data points, set the >>>> back-calculated value to the measured values so that it has no effect on >>>> the chi-squared value. >>>> for di in range(self.num_disp_points[0][si][mi][oi]): >>>> >>>> >>>> _______________________________________________ >>>> 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
