Fix inexact slope solver pH from TA #9
Conversation
|
Gentle bump! :) |
|
Thanks for this, @briochemc! And sorry I just recently noticed it. I'm not a frequent or (more importantly) adept GitHub user, and I don't seem to get email notifications for new pull requests. Doing my best to change all those things. This is great and particularly beneficial I'm sure for those in the modeling community who really don't want CO2SYS to ever hang when doing a lot of automated computations. I've updated my main branch since this pull request was opened. Will merging this PR revert those commits? Do I need to push those new commits to your master branch before merging? Again, not a particularly adept user :) so let me know what to do! And sorry again for missing this. |
|
Hi @jonathansharp, No worries! I don't think merging would have reverted anything on your side (unless there is a conflict) but I just updated my branch (I pulled your changes into my branch), so you can now safely merge this if you are happy with it :) Cheers! |
|
Perfect, thanks! I'm out on vacation this week but I'll review the changes when I return. |
|
Here is a gentle bump after quite a while! |
|
Thanks for the reminder! Obviously this fell off my radar. I remember checking this out after returning from that vacation (two and a half years ago 😑), and I wasn't able to replicate the table you shared. I'm getting differences outside the realm of machine precision, but I haven't been able to figure out why. I was getting an error at the end of 'compare_against_master.m', so I made a few changes to make it run and to specify which values are compared and then used to calculate the difference metrics. I'll paste my table below. Let me know if you're getting this same result or something different! >> compare_against_master
Running CO2SYS v3.2.0
Note: switching to 'v3.2.0'.
You are in 'detached HEAD' state. You can look around, make experimental
changes and commit them, and you can discard any commits you make in this
state without impacting any branches by switching back to a branch.
If you want to create a new branch to retain commits you create, you may
do so (now or later) by using -c with the switch command. Example:
git switch -c <new-branch-name>
Or undo this operation with:
git switch -
Turn off this advice by setting config variable advice.detachedHead to false
HEAD is now at 40e49f9 Merge pull request #8 from jonathansharp/Development
Elapsed time is 0.241042 seconds.
Running CO2SYS current commit/head
Previous HEAD position was 40e49f9 Merge pull request #8 from jonathansharp/Development
Switched to branch 'pr/9'
Your branch is up to date with 'github-desktop-briochemc/master'.
Elapsed time is 0.251568 seconds.
Relative change vs v3.2.0:
Variable Mean rel. change Min rel. change Max rel. change # of samples
TAlk 0 0 0 4320
TCO2 0 0 0 4320
pHin 8e-06 0 0.0011 4288
pCO2in 0.00015 0 0.022 4288
fCO2in 0.00015 0 0.022 4288
HCO3in 6.5e-06 0 0.00092 4288
CO3in 0.00014 0 0.019 4288
CO2in 0.00015 0 0.022 4288
BAlkin 0.00013 0 0.017 4032
OHin 0.00011 0 0.0055 4000
PAlkin 1.1e-05 0 0.00048 3744
SiAlkin 0.00012 0 0.0054 3744
AmmAlkin NaN 0
HSAlkin NaN 0
Hfreein 0.00015 0 0.021 4288
RFin 0.0052 1.2e-10 0.031 4032
OmegaCAin 0.00015 0 0.019 4032
OmegaARin 0.00015 0 0.019 4032
xCO2in 0.00015 0 0.022 4288
SIRin 0.00014 0 0.019 4288
pHout 0.00011 1.1e-12 0.00098 4032
pCO2out 0.0023 2.1e-11 0.019 4032
fCO2out 0.0023 2.1e-11 0.019 4032
HCO3out 0.00018 8.7e-13 0.0011 4032
CO3out 0.0019 1.9e-11 0.017 4032
CO2out 0.0023 2.1e-11 0.019 4032
BAlkout 0.0016 1.7e-11 0.015 4032
OHout 0.0021 2e-11 0.012 3744
PAlkout 0.0002 1.5e-12 0.0011 3744
SiAlkout 0.002 2e-11 0.011 3744
AmmAlkout NaN 0
HSAlkout NaN 0
Hfreeout 0.0021 2e-11 0.018 4032
RFout 0.0036 5.4e-10 0.024 4032
OmegaCAout 0.0019 1.9e-11 0.017 4032
OmegaARout 0.0019 1.9e-11 0.017 4032
xCO2out 0.0023 2.1e-11 0.019 4032
SIRout 0.0019 1.9e-11 0.017 4032
pHinTOTAL 8.1e-06 0 0.0011 4288
pHinSWS 8.1e-06 0 0.0011 4288
pHinFREE 8e-06 0 0.0011 4288
pHinNBS 7.9e-06 0 0.0011 4288
pHoutTOTAL 0.00011 1.1e-12 0.00099 4032
pHoutSWS 0.00011 1.1e-12 0.00099 4032
pHoutFREE 0.00011 1.1e-12 0.00098 4032
pHoutNBS 0.00011 1.1e-12 0.00098 4032
TEMPIN 0 0 0 4320
TEMPOUT 0 0 0 4320
PRESIN 0 0 0 4320
PRESOUT 0 0 0 4320
PAR1TYPE 0 0 0 4320
PAR2TYPE 0 0 0 4320
K1K2CONSTANTS 0 0 0 4320
KSO4CONSTANT 0 0 0 4320
KFCONSTANT 0 0 0 4320
BORON 0 0 0 4320
pHSCALEIN 0 0 0 4320
SAL 0 0 0 4032
PO4 0 0 0 4320
SI 0 0 0 4320
NH4 NaN 0
H2S NaN 0
K0input 0 0 0 4320
K1input 0 0 0 4320
K2input 0 0 0 4320
pK1input 0 0 0 4320
pK2input 0 0 0 4320
KWinput 0 0 0 4032
KBinput 0 0 0 4032
KFinput 0 0 0 4320
KSinput 0 0 0 4320
KP1input 0 0 0 3744
KP2input 0 0 0 3744
KP3input 0 0 0 3744
KSiinput 0 0 0 3744
KNH4input 0 0 0 3456
KH2Sinput 0 0 0 3456
K0output 0 0 0 4320
K1output 0 0 0 4320
K2output 0 0 0 4320
pK1output 0 0 0 4320
pK2output 0 0 0 4320
KWoutput 0 0 0 4032
KBoutput 0 0 0 4032
KFoutput 0 0 0 4320
KSoutput 0 0 0 4320
KP1output 0 0 0 3744
KP2output 0 0 0 3744
KP3output 0 0 0 3744
KSioutput 0 0 0 3744
KNH4output 0 0 0 3456
KH2Soutput 0 0 0 3456
TB 0 0 0 4032
TF 0 0 0 4032
TS 0 0 0 4032
TP 0 0 0 3744
TSi 0 0 0 3744
TNH4 NaN 0
TH2S NaN 0 |
|
No worries, I perfectly understand! I'm not at a computer right now but the thing I would do first is make sure you are comparing to v3.2.1 since it's the new reference now! I.e., change that version in the script that checks out the 2 branches. |
|
Does that change anything? |
|
Results were still the same. But I believe the issue was related to the derivative of Hfree (and the constituents that are calculated from Hfree) with respect to H for cases where the input pH scale is the free scale. With the scaling factor equal to zero in those cases, the derivative was calculated as infinity. I added a line that makes the derivative zero when the scaling factor is zero, and here is what I'm seeing now. Let me know what you think! >> compare_against_master
Running CO2SYS v3.2.1
Note: switching to 'v3.2.1'.
You are in 'detached HEAD' state. You can look around, make experimental
changes and commit them, and you can discard any commits you make in this
state without impacting any branches by switching back to a branch.
If you want to create a new branch to retain commits you create, you may
do so (now or later) by using -c with the switch command. Example:
git switch -c <new-branch-name>
Or undo this operation with:
git switch -
Turn off this advice by setting config variable advice.detachedHead to false
HEAD is now at 6596a8f Update README.md
Elapsed time is 0.242479 seconds.
Running CO2SYS current commit/head
Previous HEAD position was 6596a8f Update README.md
Switched to branch 'pr/9'
Your branch is ahead of 'github-desktop-briochemc/master' by 1 commit.
(use "git push" to publish your local commits)
Elapsed time is 0.263956 seconds.
Relative change vs v3.2.1:
Variable Mean rel. change Min rel. change Max rel. change # of samples
TAlk 0 0 0 4320
TCO2 0 0 0 4320
pHin 2.1e-10 0 1.4e-08 4320
pCO2in 4.1e-09 0 2.7e-07 4320
fCO2in 4.1e-09 0 2.7e-07 4320
HCO3in 1.8e-10 0 1.2e-08 4320
CO3in 3.7e-09 0 2.5e-07 4320
CO2in 4.1e-09 0 2.7e-07 4320
BAlkin 3.6e-09 0 2.3e-07 4032
OHin 3.3e-09 0 2.6e-07 4032
PAlkin 3.1e-10 0 2.3e-08 3744
SiAlkin 3.5e-09 0 2.6e-07 3744
AmmAlkin NaN 0
HSAlkin NaN 0
Hfreein 3.9e-09 0 2.6e-07 4320
RFin 3.4e-07 8.2e-14 1.3e-06 4320
OmegaCAin 4e-09 0 2.5e-07 4032
OmegaARin 4e-09 0 2.5e-07 4032
xCO2in 4.1e-09 0 2.7e-07 4320
SIRin 3.7e-09 0 2.5e-07 4320
pHout 9.4e-10 7.6e-14 9.9e-09 4320
pCO2out 1.9e-08 1.7e-12 1.9e-07 4320
fCO2out 1.9e-08 1.7e-12 1.9e-07 4320
HCO3out 1.5e-09 3.2e-14 7.8e-09 4320
CO3out 1.6e-08 1.5e-12 1.7e-07 4320
CO2out 1.9e-08 1.7e-12 1.9e-07 4320
BAlkout 1.4e-08 1.7e-11 1.5e-07 4032
OHout 1.8e-08 1.6e-12 1.1e-07 4032
PAlkout 1.8e-09 1.5e-12 9.8e-09 3744
SiAlkout 1.8e-08 2e-11 1e-07 3744
AmmAlkout NaN 0
HSAlkout NaN 0
Hfreeout 1.8e-08 1.6e-12 1.8e-07 4320
RFout 8.4e-08 1.1e-11 6.6e-07 4320
OmegaCAout 1.7e-08 1.9e-11 1.7e-07 4032
OmegaARout 1.7e-08 1.9e-11 1.7e-07 4032
xCO2out 1.9e-08 1.7e-12 1.9e-07 4320
SIRout 1.6e-08 1.5e-12 1.7e-07 4320
pHinTOTAL 2.2e-10 0 1.5e-08 4320
pHinSWS 2.2e-10 0 1.5e-08 4320
pHinFREE 2.1e-10 0 1.4e-08 4320
pHinNBS 2.1e-10 0 1.4e-08 4320
pHoutTOTAL 9.5e-10 7.6e-14 9.9e-09 4320
pHoutSWS 9.5e-10 7.6e-14 9.9e-09 4320
pHoutFREE 9.4e-10 7.6e-14 9.9e-09 4320
pHoutNBS 9.3e-10 7.6e-14 9.8e-09 4320
TEMPIN 0 0 0 4320
TEMPOUT 0 0 0 4320
PRESIN 0 0 0 4320
PRESOUT 0 0 0 4320
PAR1TYPE 0 0 0 4320
PAR2TYPE 0 0 0 4320
K1K2CONSTANTS 0 0 0 4320
KSO4CONSTANT 0 0 0 4320
KFCONSTANT 0 0 0 4320
BORON 0 0 0 4320
pHSCALEIN 0 0 0 4320
SAL 0 0 0 4032
PO4 0 0 0 4320
SI 0 0 0 4320
NH4 NaN 0
H2S NaN 0
K0input 0 0 0 4320
K1input 0 0 0 4320
K2input 0 0 0 4320
pK1input 0 0 0 4320
pK2input 0 0 0 4320
KWinput 0 0 0 4032
KBinput 0 0 0 4032
KFinput 0 0 0 4320
KSinput 0 0 0 4320
KP1input 0 0 0 3744
KP2input 0 0 0 3744
KP3input 0 0 0 3744
KSiinput 0 0 0 3744
KNH4input 0 0 0 3456
KH2Sinput 0 0 0 3456
K0output 0 0 0 4320
K1output 0 0 0 4320
K2output 0 0 0 4320
pK1output 0 0 0 4320
pK2output 0 0 0 4320
KWoutput 0 0 0 4032
KBoutput 0 0 0 4032
KFoutput 0 0 0 4320
KSoutput 0 0 0 4320
KP1output 0 0 0 3744
KP2output 0 0 0 3744
KP3output 0 0 0 3744
KSioutput 0 0 0 3744
KNH4output 0 0 0 3456
KH2Soutput 0 0 0 3456
TB 0 0 0 4032
TF 0 0 0 4032
TS 0 0 0 4032
TP 0 0 0 3744
TSi 0 0 0 3744
TNH4 NaN 0
TH2S NaN 0 |
|
That looks good to me? (I'll admit I don't remember what the scaling does.) But I would revert your double sign change for the slope and Newton step! Note that I had fixed these: Since you compute the Newton step from: H0 + slope * deltapH == 0There is a minus in front of deltapH = -slope \ H0 |
|
Okay, I'll revert that. I think the effect is the same, since the alternative with the double sign change is -H0 - slope * deltapH == 0which can be treated in the same way. But your change does bring this function more into alignment with the other iterative functions, so happy to make the change! |
2 participants
This PR is a follow-up from CO2SYS-MATLAB PR #7. It computes the exact slope used in the Newton solver for determining pH from TA.
This PR also includes some extra script called
compare_against_masterfor comparing between versions of CO2-System-Extd. It is similar tocompare_versionsexcept it only runs the CO2-System-Extd version, with the difference that it uses git to switch between the v3.2.0 tag and the current commit. (Figuring out how to run git inside a script and have it actually take effect took me a while to figure out, hence the numerous commits, which can all be squashed into a single merge commit of course.) Here's the output I get (max relative change since v3.2.0 seems to be of order 1e-7):