Hi Andy, I have modified your script to produce text files of your grid, displaying the maximum R2eff difference for each grid point (https://gna.org/support/download.php?file_id=20654). This has produced a number of maps which clearly show the problem:
map_Baldwin: https://gna.org/support/download.php?file_id=20656 map_Carver: https://gna.org/support/download.php?file_id=20655 map_Meiboom: https://gna.org/support/download.php?file_id=20657 map_Nikolai: https://gna.org/support/download.php?file_id=20658 map_NikolaiFun: https://gna.org/support/download.php?file_id=20660 map_NikolaiLong: https://gna.org/support/download.php?file_id=20659 For map_NikolaiFun (and map_Nikolai), it can be seen that the problem area is located to only one region of the parameter space: pB dw kex max_delta_R2eff 0.400000000000000 1000.000000000000000 5.492802716530586 0.592752848196470 0.252382937792077 1000.000000000000000 5.492802716530586 3.845196468401902 0.159242868221399 1000.000000000000000 5.492802716530586 5.431089306584756 0.100475457260383 1000.000000000000000 5.492802716530586 2.682512609040295 0.063395727698445 1000.000000000000000 5.492802716530586 0.107557528833629 0.040000000000000 1000.000000000000000 5.492802716530586 1.862111109685991 0.025238293779208 1000.000000000000000 5.492802716530586 1.248136912461243 0.015924286822140 1000.000000000000000 5.492802716530586 1.439357156013592 0.010047545726038 1000.000000000000000 5.492802716530586 0.384699031445987 0.006339572769844 1000.000000000000000 5.492802716530586 0.165313396220368 0.004000000000000 1000.000000000000000 5.492802716530586 0.316275262141682 0.400000000000000 1000.000000000000000 2.343672911592100 11.633432971003034 0.252382937792077 1000.000000000000000 2.343672911592100 nan 0.159242868221399 1000.000000000000000 2.343672911592100 45.046754399379282 0.100475457260383 1000.000000000000000 2.343672911592100 3.357238601090618 0.063395727698445 1000.000000000000000 2.343672911592100 3.758080253732464 0.040000000000000 1000.000000000000000 2.343672911592100 2.232620263385279 0.025238293779208 1000.000000000000000 2.343672911592100 19.760618010603249 0.015924286822140 1000.000000000000000 2.343672911592100 8.862633579602273 0.010047545726038 1000.000000000000000 2.343672911592100 15.978405412625353 0.006339572769844 1000.000000000000000 2.343672911592100 8.952365343475268 0.004000000000000 1000.000000000000000 2.343672911592100 11.945565346143020 0.400000000000000 1000.000000000000000 1.000000000000000 nan 0.252382937792077 1000.000000000000000 1.000000000000000 nan 0.159242868221399 1000.000000000000000 1.000000000000000 nan 0.100475457260383 1000.000000000000000 1.000000000000000 104.735381560928317 0.063395727698445 1000.000000000000000 1.000000000000000 102.520903807959172 0.040000000000000 1000.000000000000000 1.000000000000000 65.987026895834390 0.025238293779208 1000.000000000000000 1.000000000000000 78.994037851916559 0.015924286822140 1000.000000000000000 1.000000000000000 nan 0.010047545726038 1000.000000000000000 1.000000000000000 nan 0.006339572769844 1000.000000000000000 1.000000000000000 nan 0.004000000000000 1000.000000000000000 1.000000000000000 44.553251112908285 I.e. dw of 1000 rad.s^-1 and kex being very, very low. In all other regions of the parameter space, the max_delta_R2eff values are in the order of 1e-8 or lower. Looking at the equivalent region in map_NikolaiLong: pB dw kex max_delta_R2eff 0.400000000000000 1000.000000000000000 5.492802716530586 0.000000022779017 0.252382937792077 1000.000000000000000 5.492802716530586 0.000000370744953 0.159242868221399 1000.000000000000000 5.492802716530586 0.000000004931843 0.100475457260383 1000.000000000000000 5.492802716530586 0.000000003682589 0.063395727698445 1000.000000000000000 5.492802716530586 0.000000007636758 0.040000000000000 1000.000000000000000 5.492802716530586 0.000000002312100 0.025238293779208 1000.000000000000000 5.492802716530586 0.000000003527191 0.015924286822140 1000.000000000000000 5.492802716530586 0.000000003627945 0.010047545726038 1000.000000000000000 5.492802716530586 0.000000004664540 0.006339572769844 1000.000000000000000 5.492802716530586 0.000000001267414 0.004000000000000 1000.000000000000000 5.492802716530586 0.000000006526546 0.400000000000000 1000.000000000000000 2.343672911592100 0.000000013242452 0.252382937792077 1000.000000000000000 2.343672911592100 0.000000074741680 0.159242868221399 1000.000000000000000 2.343672911592100 0.000000004875817 0.100475457260383 1000.000000000000000 2.343672911592100 0.000000123757648 0.063395727698445 1000.000000000000000 2.343672911592100 0.000000011314873 0.040000000000000 1000.000000000000000 2.343672911592100 0.000000054738139 0.025238293779208 1000.000000000000000 2.343672911592100 0.000000011630442 0.015924286822140 1000.000000000000000 2.343672911592100 0.000000063631723 0.010047545726038 1000.000000000000000 2.343672911592100 0.000000004935935 0.006339572769844 1000.000000000000000 2.343672911592100 0.000000023574356 0.004000000000000 1000.000000000000000 2.343672911592100 0.000000021224806 0.400000000000000 1000.000000000000000 1.000000000000000 0.000000577754879 0.252382937792077 1000.000000000000000 1.000000000000000 0.000000347777326 0.159242868221399 1000.000000000000000 1.000000000000000 0.000000228853165 0.100475457260383 1000.000000000000000 1.000000000000000 0.000000259100318 0.063395727698445 1000.000000000000000 1.000000000000000 0.000000179476287 0.040000000000000 1000.000000000000000 1.000000000000000 0.000000030331906 0.025238293779208 1000.000000000000000 1.000000000000000 0.000000064205007 0.015924286822140 1000.000000000000000 1.000000000000000 0.000000276358881 0.010047545726038 1000.000000000000000 1.000000000000000 0.000000110869264 0.006339572769844 1000.000000000000000 1.000000000000000 0.000000195626962 0.004000000000000 1000.000000000000000 1.000000000000000 0.000000211432798 This makes it clear where the higher numerical precision is required. However as this situation will only arise when optimisation fails and shoots off to infinity. If the optimisation reaches a dw value of 1000 rad.s^-1, I do not believe that an optimisation algorithm could ever be able to recover and return back to reasonable values. Therefore I don't think that I'll be implementing a higher precision version of Nikolai's code in relax. It will only take time and I don't see a benefit. Anyway, your model that has been labelled as 'B14' in relax (http://wiki.nmr-relax.com/B14) will soon be fully implemented and documented. As Troels mentioned, we will have two versions of your equations implemented as separate models - one with R20A == R20B for most users, and one with R20A != R20B for those who like adventure. With all of the code optimisations we can perform within relax, that Troels is currently working on, our implementation should be faster than the Python code you originally sent us (which forms the basis of the relax implementation). Troels might be able to make it up to twice as fast. So you will soon have at your disposal a top quality implementation with a very user friendly user interface! This will also give you access to a pile of powerful tools - grid search capabilities for finding a non-biased initial optimisation position, the fast and very accurate Nelder-Mead simplex optimiser, constraints implemented with the powerful Log-Barrier iterative constraint algorithm, model nesting in the auto-analysis for speed (i.e. we could take the CR72 model results as the optimisation starting point for the B14 model rather than performing a grid search), Monte Carlo simulations for error estimates, 2D Grace graph plotting abilities for parameter values, 3D optimisation space mapping capabilities (via the dx.map user function), plotting parameter values directly onto 3D structures in Molmol or PyMOL, etc. Cheers, Edward On 6 May 2014 12:21, Andrew Baldwin <[email protected]> wrote: > Okay. The proper definition of tcp for your code shows that both functions > (fun and rieko) behave similarly. Good! > > I also agree that even with the strictest tests, a long double > implementation of the algorithm comes out with errors that are much smaller > than anything we could ever hope to measure experimentally. 23 decimal > places seems to be needed to get things to machine precision. This is > overkill. The paper makes all these points in a manner that I hope is clear > (see last fig). > > Best, > > Andy. > > > >> >> The tcp in this function in relax is defined as 1/(4*nu_CPMG), but in >> your script you have "tcp=1./(2*nu_cpmg);". Apart from that, it looks >> ok. Changing to the correct tcp gives: >> >> Relative times of the various approaches: >> Baldwin 0.159632 1.0 >> Numerical 19.974315 125.127261451 >> Meiboom 0.056497 0.353920266613 >> Carver 0.080027 0.501321790117 >> Nikolai 0.453068 2.83820286659 >> NikolaiFun 0.284052 1.77941766062 >> NikolaiLong 50.020116 313.34642177 >> NikolaiLong 50.430594 315.917823494 >> Maximum error over grid (s-1): >> Meiboom: 237214256.582 >> Baldwin: 6.57168541807e-10 >> CarverRichards: 29.9287668504 >> Nikolai fun: 104.735381561 >> Nikolai dougle (9): 142.46676955 >> Nikolai long double (18): 0.063449549192819563320223 >> Nikolai long double (23): 5.7775487944482434059452e-7 >> >> That's far more reasonable for Nikolai's fun! >> >> >>> I was hoping that the modification would make it more stable. If anything >>> the reverse is true. I would have included a note to this effect in the >>> paper (and still will if you can show it's better). Note the error is >>> very >>> temperamental. You can miss it even by setting up the grid search >>> differently. Give it a narrow grid, and you won't see the error. Only >>> relatively specific, but seemingly hard to predict sets of parameters >>> (all >>> with large deltaOs) lead it to explode. >> >> It is worth looking at all deviations rather than just the single >> maximum R2eff difference. It could be that there is an edge case in >> the grid whereby a maths domain error occurs. Or an trig function >> overflow. The median deviation to avoid such outliers might also be >> of worth. By switching in your code all numpy.max() calls with >> numpy.median(), and making sure the "Nikolai long double (*)" results >> are pre-converted to a numpy float64 array, I then see an incredibly >> different picture: >> >> Relative times of the various approaches: >> Baldwin 0.152725 1.0 >> Numerical 19.596446 128.3119725 >> Meiboom 0.05604 0.366934031756 >> Carver 0.077236 0.505719430349 >> Nikolai 0.43944 2.87732853167 >> NikolaiFun 0.273977 1.7939237191 >> NikolaiLong 50.243433 328.979754461 >> NikolaiLong 50.784841 332.524740547 >> Median error over grid (s-1): >> Meiboom: 4.0064828013 >> Baldwin: 3.9630521087e-12 >> CarverRichards: 0.00661151728704 >> Nikolai fun: 6.62225829728e-12 >> Nikolai dougle (9): 6.70041799822e-12 >> Nikolai long double (18): 2.62099543458e-12 >> Nikolai long double (23): 2.27758174227e-12 >> >> ! >> >> So this is much clearer - outliers are the major difference you are >> seeing. A study of such outliers and their cause would be quite >> interesting. But in most of these cases, simple collapsed equations >> for each model should be derivable at the edges which then produce the >> expected R2eff over the entire parameter space rather than >> computational artefacts in some corners. Anyway, this result makes it >> crystal clear that the different models should not be judged by the >> single maximum outlier over the entire parameter space. >> >> Regards, >> >> Edward > > _______________________________________________ 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
