Hi Ed. I think I have fixed the 1e100 to R20. That was in commit 23297 - 23304
I still have to look at: dpl94 lm63 And, I will look at the unit tests. Best Troels 2014-05-21 19:21 GMT+02:00 Edward d'Auvergne <[email protected]>: > 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 _______________________________________________ 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
