search

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 27 forks source link

Non-issue: Results of Perm1 comparisons #320

Closed sbodorkos closed 5 years ago

sbodorkos commented 5 years ago

Not an issue: just a stub to remind myself that the results of Perm1 U-centric comparisons between SQUID 2.50 (11.02.03) and Squid3 (1.1.8) look good across all index isotopes. I am offline for most of the next 48 hours, but will upload the comparison files to this thread on (my) Wednesday.

There are a couple of interesting features/systematic offsets that I should chase, but overall the comparison is excellent, which is good news because Perm1 U-centric Tasks are the bread and butter of U-Pb zircon work (and most U-Pb geochronology, I would think)...

sbodorkos commented 5 years ago

For anyone tempted to try replicating my SQUID 2.50 vs Squid3 comparisons with their own test dataset, you should be aware that I have expended little programmatic effort: instead, I have performed cell-by-cell subtraction based comparisons of CSV-style output, following a process to ensure that the mapped cells are comparing apples with apples. Unfortunately, this approach means that the template developed for comparison purposes is quite prescriptive:

  1. Zircon data with 10 (ideal) or more (manual editing needed) mass-stations (although mass-station identities are not prescribed), and SQUID 2.50 CPS columns switched ON for exactly 10 mass-stations. (If you have more than 10 mass-stations, you can still proceed, but you'll still need exactly 10 CPS columns ticked. When generating the corresponding Squid3 output, you'll need to open the ReferenceMaterialReport CSV and manually delete the column(s) corresponding to the excess CPS columns unticked in SQUID 2.50).
  2. Exactly 10 ratios of interest, defined alphabetically, which ought to include every ratio you might use for data processing. Their identities are not prescribed (although in SQUID 2.50, Bkrd should not be used), but they need to be defined ALPHABETICALLY in SQUID 2.50, primarily by numerator mass and secondarily by denominator mass (i.e. the first ratio should comprise the lowest-mass numerator of interest and its lowest-mass of its denominators of interest, and so on, until the 10th ratio comprises the highest-mass numerator of interest and its highest-mass denominator of interest). This measure is to "auto-match" output ordering implemented for the Squid3 CSV.
  3. Exactly 2 user-defined expressions for which ST = YES and NU = YES. Their fundamental purpose is to be column placeholders; there is no need for your data-processing to actually use them, but they must not be omitted. It's an opportunity to verify that the numerical evaluation performed in Squid3 matches that performed by SQUID 2.50. My expressions here are just ln(["206/238"]) and ln(["254/238"]).
  4. NO other user-defined expressions of any kind. The aim is to test Squid3 replication of "under-the-hood" SQUID 2.50 functionality. Comprehensive implementation and testing of user-defined functions in Squid3 will come later.

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 relevant Task, and my (GA Session 180050) 10-peak zircon XML test file: SQUID2.50+Perm1ConcUTask+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 three "live" Excel SQUID-workbooks corresponding to SQUID 2.50 data-reductions using common Pb index isotope = 204Pb, 207Pb and 208Pb respectively (180050_SQ2.50_Perm1_U_notfrozen.zip); the second set comprises the same three SQUID-books in "frozen" format i.e. with formulae stripped, and all cells containing evaluated results: 180050_SQ2.50_Perm1_U_frozen.zip).

Squid3 source files and outputs

Building Blocks The version of Squid3 I used is 1.1.8, but I am not sure how best to supply the relevant files. So in addition to the Project file itself (118 180050 PERM1 CONCU), I have also supplied copies of the four ParameterModel XML files I constructed for comparison purposes:

  1. Isotopic Reference Material = "TEMORA2 dates only v1.0"
  2. Concentration Reference Material = "M127 values only v1.0"
  3. Common Pb = "Stacey-Kramers Common Pb at 416.8 Ma v1.0"
  4. Physical Constants = "Squid 2.5 Default Physical Constants Model v1.0"

The first two are hopefully self-explanatory; the third defines the "default" common Pb isotope compositions used in Squid3's analogues of the StandardData and SampleData sheets. The fourth is important because it forces Squid3 to use present-day 238U/235U = 137.88 (after Steiger & Jager 1977) to match SQUID 2.50, in preference to Squid3's other Physical Constants models, which instead use 137.818 after Hiess et al. (2012).

All five files are included in this ZIP (Squid3_ProjectFile+4ParameterModels.zip), because the extent to which my Project file incorporates them is not clear to me. This is particularly true of the parameter model files in persistent use by the Project file: Squid3 clearly knows where the files are, but I have no idea where three of the four are stored on my PC. It was necessary for me to re-export those three to known locations on my computer, in order to bundle them in this ZIP.

Results If you'd rather just fast-forward to my Squid3 results based on that combination of files, they are here. This ZIP (180050_Squid3_Perm1_ConcU_RM.zip) comprises three XLS files, derived from 'ReferenceMaterialReportTable' CSV files 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. There is one XLS file per index isotope.

SQUID 2.50 vs Squid3 comparison: Perm1 ConcU

Finally, the results of the actual comparison are here: RM_Comparers_Perm1_Ucentric.zip. There are three XLS files, one for each index isotope, 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 ("Perm1"), Cell A4 indicates the index isotope ("4" = 204Pb, "7" = 207Pb, "8" = 208Pb), 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 ("Perm1"), Cell A4 indicates the index isotope ("4" = 204Pb, "7" = 207Pb, "8" = 208Pb), 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.

There are a number of reasons why the Excel-based comparison of SQUID 2.50 and Squid3 values is not exact, which I will not go into here. In general I am looking for very small discrepancies between the two values (ideally more than 10 orders of magnitude smaller than the values themselves). Within individual columns, I am also looking for both positive and negative values, to try and reassure myself that the discrepancies are arising from numerical "noise" beyond my control, rather than very small bugs in "the code" (either SQUID 2.50 or Squid3).

There are some interesting results in the three comparisons thus far... but I'll write about those separately.

sbodorkos commented 5 years ago

A brief description of the comparisons, with initial reference to sheet CellByCell in XLS RM_Comparer_Perm1-204corr_Ucentric:

Columns A–I: Only column C involves a calculation and it is trivial.

Columns J–AQ: In general, we get perfect agreement on the fundamental calculations (mass-station CPS and calculation of the isotopic ratios) because the individual steps of the underlying arithmetic are simple (although very numerous). In this case, it helps that the ratios are being evaluated in 'spot average' mode, because the arithmetic associated with 'linear regression to burn mid-time' is not quite perfectly robust (with respect to intra-Excel 2003 solutions between VBA code and an Excel spreadsheet tasked with performing the identical calculations).

Columns AR–AU, BB–BC, BK–BL, EX: Divergence usually begins to appear in the evaluation of NU-switched Task expressions (AR–AU = user-defined general expressions; others are part of the "Special U-Th-Pb Task equations"), because a recurring step in the VBA evaluation involves converting full double-precision binary numbers into a decimal representation (15 significant figures, scientific notation) that is converted to a text-string as part of the evaluation. This process has proven impossible to reliably replicate, even inside Excel 2003, because the rounding rules for VBA and spreadsheet-based arithmetic are different.

Expressions involving a single ratio (columns AR–AU) are more easily replicated than those involving multiple ratios (BB–BC, BK–BL, and EX), so while agreement is perfect for both value and uncertainty in the former, the latter pairs of columns (where the corresponding expression contains two isotopic ratios in each case) exhibits the beginnings of occasional divergence in the calculated uncertainties.

Column AW onwards: Agreement is overall very good. Disagreements are very small, usually at least 10 orders of magnitude smaller than the values being compared. Furthermore, in any given column, there are both positive and negative values, consistent with discrepancies arising from largely random "numerical noise", perhaps reflecting downstream propagation of divergence that began appearing in column BC especially, as those are quite far-reaching. But there are some exceptions:

Column BH (204Pb206Pb_OvCtCorFrom207Pb): Although all the numbers here are very small, there is very little dispersion and the numbers are all negative. This suggests a diagnosable difference between row-by-row SQUID 2.50 and Squid3 calculations, and the consistency of the offset suggests it invo;lves a constant, at some level. My suspicion is that some part of the Common Pb values have a different number of significant figures in SQUID 2.50 as opposed to my "Stacey-Kramers at 416.8" Squid3 parameter model, but I have not yet investigated in detail.

Column ES (4cor_Th_Concen_RM): In Perm1, calculation of Th (ppm) values in SQUID 2.50 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).

In SQUID 2.50, using the rounded value 1.033, Th (ppm) in Tem.1.1 is 366.279... ppm (see sheet 250_StdData, cell BF8). Whereas in Squid3, where the coefficient is evaluated at full double precision per the below screenshot ("Quotient" is the result), Th (ppm) in Tem.1.1 is 366.135... ppm (see sheet Sq3_RMReport, cell AY8).

image

A large part of the discrepancy is presumably attributable to the value of the coefficient used, but it is not clear whether that is the entire problem. So I performed an "after-market" repair to the SQUID 2.50 data, by Copying the Quotient value from the Squid3 "peek", and Pasting it in the RM_Comparer... XLS (sheet 250_StdData, red cell BF1). Then I modified the arithmetic in sheet 250_Comprehensive so that the SQUID 2.50 data are divided by [Quotient/1.033] in order to obtain a result that should be arithmetically identical to that produced by Squid3, if the value of the coefficient is the sole source of the discrepancy.

And it looks like it is: for Tem.1.1, the cell-by-cell comparison was -1.445e-01 before the adjustment, and +1.137e-13 afterwards. The values being compared are 3.66e+02, so the latter discrepancy is right at the limit of Excel-based double-precision arithmetic. So I am marking that one "solved".

sbodorkos commented 5 years ago

Additional notes specific to the 207-based comparisons, with initial reference to sheet CellByCell in XLS RM_Comparer_Perm1-207corr_Ucentric:

Column BH (204Pb206Pb_OvCtCorFrom207Pb): Note that this problem persists from the 204corr comparison (and see column DK below).

Columns BZ, CA and CJ (4cor_207Pb235U_RM, its %err, and 4cor_ErrCorrel_RM): These comparisons are poor because SQUID 2.50 compromises the calculation of 4cor_207Pb235U_RM when the index isotope is not 204Pb. Essentially, SQUID 2.50 performs the following (for column BZ):

4cor_207Pb235U_RM = (4cor_207Pb206Pb_RM) (Xcor_206Pb238U_RM) (Ref_238U235U)

which yields a compromised ("hybrid") result when X =7 or 8. In Squid3, we do not permit 4cor_207Pb235U_RM to be populated with hybrid expressions of this type; instead, its evaluation is contingent upon valid evaluation of 4cor_206Pb238U_RM. Columns CA and CJ represent the downstream effects of the column-BZ calculation, and none of these discrepancies need addressing (remembering that all the comparisons were fine in the previous RM_Comparer_Perm1-204corr_Ucentric comparison, where the index isotope is 204Pb and the analogous calculations are "clean".

Column DK (7cor_Com206Pb_Pct_RM): Extremely similar pattern to column BH; in fact, the undispersed negative discrepancies in column DK are even proportional to those in BH. It seems clear that the underlying cause is the same; time to investigate!

sbodorkos commented 5 years ago

Columns BH and DK above arise from "edge effects" in the double-precision arithmetic, relating to the evaluation of RefRad_207Pb206Pb i.e. reference 207Pb/206Pb of TEMORA2, via the LudwigLibrary function of the form "=Pb76(416.8)".

0.0551127497968983 is the Squid3 solution (via Peek in Expression Manager) 0.0551127497968982 is the SQUID 2.50 solution in Excel 2003

I suspect this is the "banker's rounding" effect in Microsoft VBA i.e. when rounding of decimal representations of double-precision numbers is demanded procedurally, the final digit is rounded to the nearest even number (rather than to the nearest integer, as per "symmetrical arithmetic" rounding as implemented in Excel spreadsheets, and most people's brains). But it is very encouraging that the "downstream" effects of this type of discrepancy manifest as systematic offsets; makes diagnosis a lot easier!

208-corrected 206Pb/238U data show a similar, apparently inconsistent, whereby discrepancies in the 208-corrected calibration constants are distributed either side of zero (as we would expect), as are discrepancies in the 208-corrected spot ages derived directly from them. However, spot-by-spot 208-corrected 206Pb/238U is systematically offset between SQUID 2.50 and Squid3 (every spot value in Squid3 is very slightly higher than its SQUID 2.50 equivalent).

The cause is a similar issue with the WtdMeanA value. In this case, it is conceivable that the discrepancy is fuelled by minor inconsistencies in the "upstream" spot-by-spot uncertainties of the input calibration constants (rather than being purely a rounding effect in the final value), but the difference is still very small:

0.00447238288341184 in the Squid3 Expressions peek (although I was interested to see that the corresponding value in the Squid3 Visualisation... Reference Material... Weighted Mean contains an extra digit (another 4) on the end. 0.00447238288341185 in SQUID 2.50 in Excel 2003

Again, it is encouraging to see such an insignificant offset generating uniform discrepancies in downstream calculations - this confirms my bias that a mix of positive and negative discrepancies is a relatively robust indicator of correctness in our arithmetic.

sbodorkos commented 5 years ago

Perm1 Th-centric

These permutations did not add anything new. Discrepancies relate to:

  1. The "magic number" (L1033 in Squid3) involved in generating the secondary element (ppm) from the primary element (ppm) and the 232Th/238U. SQUID 2.50 uses a "magic number" value rounded to 4 significant digits; Squid3 uses the full double-precision value; substituting the latter into the former resolves the issue.
  2. The output of Pb76(416.8) between SQUID 2.50 and Squid, as per previous post.

SQUID 2.50 Task: Same XML and Preferences as Perm1 ConcU, in conjunction with this Task: SquidTask_Zrn 10pk Perm1 ConcTh.GA.zip

SQUID 2.50 results are ZIPs each comprising sets of three Excel workbooks (one for each index isotope). One set 'unfrozen', so SQUID 2.50 users can see the detail of the calculations, and one lot 'frozen' so non-SQUID 2.50 users have full access to the calculated values:

180050_SQ2.50_Perm1_Th_notfrozen.zip 180050_SQ2.50_Perm1_Th_frozen.zip

Squid3 building blocks: Here is the Squid3 Project file: 118 180050 PERM1 CONCTH.zip. Note that this includes both the Task based on "Zrn 10pk Perm1 ConcTh", and an updated value for "Default Common 208Pb206Pb" (2.0971, not 2.0941 as it was in the Perm1 ConcU Project file).

[As luck would have it, the Perm1 ConcU results are unaffected by the typo, because the common-Pb arithmetic only utilises 206Pb204Pb and 208Pb204Pb, and I had typed both of those values correctly.]

Squid3 results: Here a ZIP containing the three RM ReportTable CSVs (one for each index isotope): 180050_Sq3_Perm1_ConcTh_RM.zip

SQUID 2.50 vs Squid3 comparisons: A ZIP containing three Excel workbooks, each comprising five worksheets, as per their equivalent files in Perm1 ConcU: RM_Comparers_Perm1_Thcentric.zip.

@bowring this is the end of the Perm1-based comparisons; I don't intend to add anything further. I will make a separate "non-issue" for the results of the combined Perm2/Perm4 comparisons (I will handle them together because they are so similar), so feel free to close this one whenever you're ready.

bowring commented 5 years ago

Great news!