CIRDLES / Squid

Squid3 is being developed by the Cyber Infrastructure Research and Development Lab for the Earth Sciences (CIRDLES.org) at the College of Charleston, Charleston, SC and Geoscience Australia as a re-implementation in Java of Ken Ludwig's Squid 2.5. - please contribute your expertise!
http://cirdles.org/projects/squid/
Apache License 2.0
12 stars 24 forks source link

Non-issue: Results of comparisons for Perm3 #330

Closed sbodorkos closed 5 years ago

sbodorkos commented 5 years ago

The final tests relate to Perm3, which reflects 208Pb/232Th geochronology without any 206Pb/238U calibration. There were a total of 4 permutations to be tested:

  1. 204corrected + Thcentric
  2. 204corrected + Ucentric
  3. 207corrected + Thcentric
  4. 207corrected + Ucentric

Discrepancies arising were broadly similar to those documented for Perm1 and Perm2+Perm4, and for the same reasons (i.e. the issues are Perm-independent). In general, I will not revisit these, other than to explain "after-market" adjustments to the SQUID 2.50 workbooks, implemented universally in order to minimise disagreement with Squid3 results in a transparent way. Discrepancies specific to Perm3 will be documented in the next post.

SQUID 2.50 source files and outputs

Building blocks My comparison uses an edited version of SQUID 2.50.11.02.03 (VBA code has been modified in order to fix bugs and minimise discrepancies in the results of double-precision arithmetic between VBA and Java). If you're confident in your ability to switch versions of SQUID 2.50 in Excel 2003 and you want to try it out, it is included in this ZIP, together with the 2 relevant Tasks (spanning Thcentric/Ucentric), and my (GA Session 180050) 10-peak zircon XML test file: SQUID2.50+TwoPerm3Tasks+180050XML.zip.

Results If you'd rather just fast-forward to my SQUID 2.50 results based on that combination of files, they are here. The first set comprises 4 "live" Excel SQUID-workbooks corresponding to the 4 cases defined in the list above: 180050_SQ2.50_Perm3_notfrozen.zip. The second set comprises the same 4 SQUID-books in "frozen" format i.e. with formulae stripped, and all cells containing evaluated results: 180050_SQ2.50_Perm3_frozen.zip.

"After market" adjustments included in these Results The suffixes "BF[BE]_BQ[BP]_CI[CH]" refer to the SQUID 2.50 Excel-columns where I have made adjustments to improve agreement with Squid3. The adjustments are invariant across all 4 cases under consideration; their locations depend only on the identity of the index isotope: 204Pb[207Pb]. The location variation arisies from a quirk of SQUID 2.50 in Perm3 mode: Tasks utilising the 204-correction do tabulate 204Cor_Com206_Pct for the Reference Material, but those utilising the 207-correction do NOT tabulate 207Cor_Com206_Pct.

Column BF[BE] (magic number): In Perm3, this column contains the concentration of the "secondary parent" element (e.g. Th when the concentration RM is defined in terms of its U content, and vice versa). In Perm3 Tasks in SQUID 2.50, this involves invocation of the hard-coded 'magic number' 1.033, which turns out to be a coefficient linking present-day U abundance to present-day Th abundance by reference to present-day 238U/235U (noting that Squid3 models present-day 238U/235U as a property of the reference material, rather than a universal physical constant).

For the purpose of this comparison, in the SQUID 2.50 workbooks, I have manually replaced the magic number 1.033 with the full double-precision quotient of the analogous expression in Squid3, in exactly the same way as was implemented for the Perm1 comparison (the below graphic is copied from there):

image

Without this adjustment, secondary element (ppm) values between SQUID 2.50 and Squid3 usually diverge at the first or second decimal place; afterwards. After adjustment, the ppm difference is usually smaller than 1e-10.

Column BQ[BP] (12 significant-figure rounding): This dataset contains an interesting one-off complication that illustrates some unintended consequences of universal 12-significant-figure rounding in both SQUID 2.50 and Squid3, implemented by @bowring and I in an attempt to get "best possible" agreement between isotopic ratios and NU-switched Task equation, with the aim of minimising downstream divergence.

The issue relates to the %err values associated with Uncor_208Pb232Th_CalibConst, which in this dataset encompass values both above and below 1%. After rounding, values above 1% have 11 decimal places; those below have 12. This means that when values diverge between SQUID 2.50 and Squid3 at the final digit (of the rounded data), the cell-by-cell difference is a full order of magnitude greater when the values being compared are greater than 1%, than is the case when the compared values are less than 1%. This is true even in the "1.01% vs 0.99%" case.

In the current dataset, the magnitude of SQUID 2.50 vs Squid3 "final-digit" divergences is 1e-12 for 8 of the first 23 TEMORA analyses (1, 2, 4, 11, 14, 18, 19, 21), because the compared values span the range 0.60% to 0.98%. But for the final analysis (Tem.24.1) has a value of 1.52%, so its "final-digit" divergence is 1e-11, and that analysis alone outweighs all of the other 8 in terms of its effect on the %err values for the common-Pb corrected 208Pb/232Th, which in turn influences the inverse-variance weighted mean A value, which in turn affects all the downstream "built-in" expressions. This result in widespread and contradictory (albeit minor) discrepancies when comparing the SQUID 2.50 and Squid3 datasets.

To mitigate the issue, I manually pasted the Squid3 value (1.52247762935) over the top of the SQUID 2.50 value (1.52247762934) for the %err associated with 'Uncorr Pb/Th const' in analysis Tem.24.1. This resolves the immediate problem, but it also underlines the need to abandon 12-significant-figure rounding in both SQUID 2.50 and Squid3 as soon as practicable.

Column CI[CH] (Pb76 function output): This is the same issue as outlined in Perm1 and Perm2+Perm4: the LudwigLibrary function "=Pb76(416.8)" yields a value of 0.0551127497968983 in Squid3, and 0.0551127497968982 in SQUID 2.50. In this comparison, I manually pasted the Squid3 value over the top of both SQUID 2.50 values for "Std 207Pb/206Pb" (i.e. RefRad_207Pb206Pb). This minimises (but does not completely eliminate) the systematic errors downstream.

Squid3 source files and outputs

Building Blocks For this Perm3-specific set of comparisons, I upgraded by Squid3 version to 1.2.1, in order to take advantage of the opportunity to test the bug fix to CSV reporting of the value of 204cor_208Pb232Th_Age_RM and its 1sigmaErr, as reported in issue #326 and resolved in version 1.2.1. Here are the 2 relevant Project files (spanning Thcentric/Ucentric): Sq3 121 180050 PERM3 CONCTh+U.zip.

Results If you'd rather just fast-forward to my Squid3 results based on that combination of files, they are here. This ZIP (121 180050 Sq3 Perm3 CSVs.zip) comprises 4 CSV files, corresponding to the 4 cases defined above, generated in Squid3 (after setting the index isotope in Visualisations... Reference Material... Weighted Means) via Report Tables... Reference Materials... Custom Report Tables... per Squid 2.x.

SQUID 2.50 vs Squid3 comparison: Perm3

Finally, the results of the actual comparison are here: RM_Comparers_Perm3_Th+Ucentric.zip. There are 4 XLS files, one for each case defined above, and each file comprises five worksheets:

  1. "250_StdData": This is where the SQUID 2.50 dataset is Pasted and described. Cell A3 is the permutation number ("3"), Cell A4 indicates the index isotope ("4" = 204Pb, "7" = 207Pb), Cell A5 indicates the parent element by which the CONCENTRATION reference material is defined ("U" or "Th"). These values need to be entered manually, because SQUID 2.50 output varies depending on these settings. At cell A7, the rectangular array of SQUID 2.50 StandardData values (including the single row of column-headers) is Pasted (via sequential Paste as Values, then Paste Formats).
  2. "250_Comprehensive": This sheet has no user input. It maps the data columns from "250_StdData" to the superset of Unique Expression-Names defined for Reference Materials by @NicoleRayner and I (see column L of ExpressionNames.AcrossPerms_v5_NR+SB.xlsx). In general, yellow cells in row 7 represent entries from that 'unique name list' that contain 'row-by-row' data (as opposed to 'summary of column' data). Green cells in row 7 indicate uncertainties associated with the yellow cells (where those exist: note that the XLSX above does not explicitly consider the uncertainties associated with many (but not all) of the uniquely named expressions), with uncertainties labelled with magnitude and whether percentage or absolute. Orange cells represent Unique Expression names that pertain to 'summary of column' data (which I have postponed comparing, because of course those results will depend on the veracity or otherwise of the row-by-row data, which should therefore be considered first).
  3. "Sq3_RMReport": This is where the Squid3 dataset is Pasted and described. As for "250_StdData", cell A3 is the permutation number ("3"), Cell A4 indicates the index isotope ("4" = 204Pb, "7" = 207Pb), Cell A5 indicates the parent element by which the CONCENTRATION reference material is defined ("U" or "Th"). The rectangular array of Squid3 data is Pasted in TWO PARTS (via sequential Paste as Values, then Paste Formats): (a) At cell A7, the single-row header and data for column A, and (b) at cell B3, the multi-row headers (starting with 'Spot Fundamentals' and rectangular array of Squid3 ReferenceMaterial data. This two-step process should yield a rectangular array of numeric data.
  4. "Sq3_Comprehensive": This sheet has no user input. As for "250_Comprehensive", it maps the data columns from "Sq3_RMReport" to the superset of Unique Expression-Names defined for Reference Materials.
  5. "CellByCell": This is where the magic happens... a cell-by-cell comparison of the two "...Comprehensive" sheets. Where the two counterpart cells both contain numbers, the SQUID 2.50 value is subtracted from the corresponding Squid3 number.

Discrepancies specific to Perm3 will be documented in the next post.

sbodorkos commented 5 years ago

New issues specific to Perm3 stem fundamentally from the SQUID 2.50 StandardData column labelled "206*/238", which in Perm1, Perm2 and Perm4 is equivalent to Squid3 "unique name" Xcor_206Pb238U_RM, where X = 4 or 7, as the index isotope identifier for the common Pb correction. In each of those three Perms, the value is derived from a starting Uncor_206Pb238U_CalibConst (i.e. "Uncorr Pb/U const" in SQUID 2.50).

The difficulty in SQUID 2.50 Perm3 is that no "Uncorr Pb/U const" exists, so there is no consistent basis for a "206*/238" calculation. Instead, in this StandardData column in Perm3, SQUID 2.50 executes:

Xcorr_208Pb232Th_CalibConst / WtdAv_Xcor_208Pb232Th_CalibConst[0] * RefRad_208Pb232Th

which would be functionally equivalent to the theoretical (but in SQUID 2.50 and current Squid3 practice, non-existent) expressions "Xcor_208Pb232Th_RM" (Squid3), and "208*/232" in SQUID 2.50.

Arithmetically, this is fine; the practical problem is that SQUID 2.50 uses this calculation as the basis for the further calculation of all the expressions and parameters needed to construct a conventional (Wetherill) concordia diagram for StandardData, and commencing with "208*/232" renders that process meaningless.

This issue was recognised during coding of the built-in expressions for Squid3, and an alternative approach was devised, which has no SQUID 2.50 equivalent in StandardData, but which is closely based on SQUID 2.50 calculations for SampleData. In Squid3, these calculations are implemented across all Perms for completeness, but only in Perm3 do they come in to "front line" use.

Essentially, Squid3 calculates values for "Xcor_Total_20YPb23ZParent_RM" corresponding to every "Uncor_20YPb23ZParent_CalibConst" expression present in the Perm, where X identifies the index isotope for the common Pb correction (4, 7, or [in the case of Perm1 only] 8), Y identifies the daughter radiogenic Pb isotope (6 or 8) , and ZParent defines the parent actinide ("8U" or "2Th"). In Perm2 and Perm4, these calculations go no further. In Perm1, these calculations could theoretically be extended in order to calculate Xcor_208Pb232Th_RM, but in practice, neither SQUID 2.50 nor Squid3 (yet!) does this.

But in Perm3, the extension of the analogous calculations yields a plausible (albeit indirectly-calibrated) value for Xcor_206Pb238U_RM, which is important as a starting point for the generation of a Wetherill concordia array. For this reason, the Perm3-specific calculations are pursued in the Squid3 implementation.

(The scientific question of whether we should have bothered is a separate issue! It could be argued that Perm3 users, having specified a Task embodying primary 208Pb/232Th calibration and no direct calculation of a 206Pb/238U calibration, would not be interested in the representation of their Reference Material data on a conventional Wetherill diagram. But in practice, I was concerned about Squid3 delivering a Reference Material dataset that was discernibly "less" than that delivered by SQUID 2.50, notwithstanding the analogous calculations in SQUID 2.50 being isotopically incoherent as demonstrated above.

_Plus I wanted to see whether the calculation "worked", because if so, it also defines a path to the "all-Perms" delivery of calibrated Xcor_208Pb232Th_RM data, with particular respect to Perm1, should we ever decide that Squid3 should do that. In my opinion, its current omission is incongruous given the effort expended on Xcor_206Pb238U_RM, Xcor_207Pb235U_RM, and Xcor_207Pb206PbRM.)

So for Perm3 specifically, Xcor_Total_208Pb232Th_RM is calculated as:

Uncor_208Pb232Th_CalibConst / WtdAv_Xcor_208Pb232Th_CalibConst[0] * RefRad_208Pb232Th

which is similar to the "Xcor_208Pb232Th..." expression above, except "Xcor_Total_208Pb232Th..." involves "Uncor_208Pb232Th_CalibConst" rather than "Xcor_208Pb232Th_CalibConst".

From there, Xcor_Total_206Pb238U_RM is calculated as:

Xcor_Total_208Pb232Th_RM / ["208/206"] * 232Th238U_RM

where ["208/206"] is the spot-by-spot measured (uncorrected) 208Pb/206Pb, and 232Th238U_RM is the output of the (index isotope-independent) Special U-Th-Pb expression governing 232Th/238U.

At this point, the calculations branch depending on the identity of the index isotope, analogous to SQUID 2.50 procedure in SampleData calculations.

If the index isotope is 204Pb, 4cor_206Pb238U_RM is calculated directly:

4cor_Total_206Pb238U_RM * ( 1 - DefCom_206Pb204Pb * ["204/206"] )

where DefCom_206Pb204Pb is the assumed (default) common 206Pb/204Pb, and ["204/206"] is the spot-by-spot measured 204Pb/206Pb.

If the index isotope is 207Pb, the calculation involves intermediate evaluation of 7cor_206Pb238U_Age_RM, via the (Squid3 adaptation of the) LudwigLibrary function Age7corr:

Age7corr( 7cor_Total_206Pb238U_RM )

after which the corresponding ratio 7cor_206Pb238U_RM is derived via the age equation:

exp( Lambda238 * 7cor_206Pb238U_Age_RM ) - 1

where Lambda238 is the decay constant for 238U.

So, bearing in mind this is "new" code in Squid3 (i.e. beyond direct replication of SQUID 2.50), I think the arithmetical steps are logical enough, especially when compared directly with SQUID 2.50 implementation of the analogous calculations in SampleData. The actual correctness of the values obtained will depend heavily on the (very rarely evaluated!) applicability of the Special U-Th-Pb expression (controlling the index isotope-independent value of 232Th/238U) to the system being measured.

With respect to SQUID 2.50, the other "columns" affected by the Perm3-specific isotopic mismatch in calculating the column labelled 206*/238 (and their Squid3 correlatives) are:

  1. 206*/238 %err (i.e. Xcor_206Pb238U_RM)
  2. 207/235 and 207235 %err (i.e. 4cor_207Pb235U_RM)
  3. err corr (i.e. 4cor_ErrorCorrel)

Note also that in the case of 207Pb-corrected data, the new SQUID 2.50-related problems with 207/235, 207/235 %err, and err corr are separate and additional to the infidelity noted in other Perms (and summarised in the "universal" issues in the previous post in this thread), arising from arithmetic that mixes common Pb corrections in order to calculate a "207*/235" value despite the use of 207Pb as the index isotope.

The good news is that I can't spot any Perm3-specific discrepancies not covered by the above description, so TL; DR = no further changes yet needed in Squid3.

I can't think of anything else, so @bowring you can close this when ready. I guess the next step is to start regathering the threads regarding documentation of SampleData-specific 204Pb-overcount correction, common-Pb correction, and the Grouping algorithm, which will involve revisiting the WtdAv routine to document the branch ("third round of outlier rejection") I deliberately skipped about this time last year, when we were under time pressure to be ready for the SHRIMP Workshop in Korea.

sbodorkos commented 5 years ago

@bowring I must have lost you before you reached the end of my epistle in the last post! Close this "non-issue" whenever you're ready.

bowring commented 5 years ago

@sbodorkos - Great work! We have reached an important milestone.