diff --git a/pwiz_tools/Osprey/Osprey.ML/GradientBoostedTrees.cs b/pwiz_tools/Osprey/Osprey.ML/GradientBoostedTrees.cs
index 2c8e1d8c12..67a2f992d2 100644
--- a/pwiz_tools/Osprey/Osprey.ML/GradientBoostedTrees.cs
+++ b/pwiz_tools/Osprey/Osprey.ML/GradientBoostedTrees.cs
@@ -21,31 +21,59 @@
* limitations under the License.
*/
-// Pure-managed gradient-boosted decision trees for binary classification, used
-// as a non-linear alternative to the linear Percolator SVM for FDR scoring.
+// Pure-managed gradient-boosted decision trees, used as a non-linear alternative to the
+// linear Percolator SVM for FDR scoring (binary logistic) and as the m/z calibration
+// model for MARS (squared error).
//
// Second-order (Newton) boosting with the XGBoost regularized objective
-// (Chen & Guestrin 2016): logistic loss, per-leaf L2 (lambda) + L1 (alpha)
-// penalties, minimum split gain (gamma), minimum child hessian, row/column
-// subsampling, and shrinkage. Histogram split finding over quantile-binned
-// features. No native dependencies (builds on net472 + net8.0).
+// (Chen & Guestrin 2016): per-leaf L2 (lambda) + L1 (alpha) penalties, minimum split
+// gain (gamma), minimum child hessian, row/column subsampling, and shrinkage.
+// Histogram split finding over quantile-binned features. No native dependencies
+// (builds on net472 + net8.0).
//
-// The model output is a raw log-odds margin; the caller ranks by it exactly as
-// with the SVM discriminant (target-decoy competition, q-values, PEP).
+// Everything except the base score and the per-round gradient is loss-agnostic: quantile
+// binning, histogram split finding, the L1 soft-threshold and L2 leaf weight, subsampling
+// and the flat node arrays all apply unchanged to either objective.
+//
+// The model output is the raw additive margin with no link function. For logistic that is
+// a log-odds the caller ranks by exactly as with the SVM discriminant (target-decoy
+// competition, q-values, PEP); for squared error it is the prediction itself.
using System;
using System.Collections.Generic;
+using System.Threading.Tasks;
namespace pwiz.Osprey.ML
{
+ /// Loss function optimized by .
+ public enum GbtObjective
+ {
+ /// Binary logistic loss. Base score is the log-odds of the weighted
+ /// positive fraction; h = p(1-p)w, which never exceeds 0.25w.
+ LogisticBinary,
+
+ /// Squared error. Base score is the weighted mean of y; g = (f-y)w and
+ /// h = w, so an unweighted hessian sum is exactly a sample count.
+ SquaredError
+ }
+
/// Hyper-parameters for . Defaults are a
/// conservative, regularized setting matching the validated Python XGBoost run.
public sealed class GbtParams
{
+ /// Loss function. Defaults to the binary logistic loss used by the FDR path.
+ public GbtObjective Objective = GbtObjective.LogisticBinary;
public int NTrees = 200;
public int MaxDepth = 6;
public double LearningRate = 0.1;
- /// Minimum summed hessian (Σ p(1-p)) per leaf; blocks leaves that fit a handful of points.
+ /// Minimum summed hessian per leaf; blocks leaves that fit a handful of
+ /// points. NOTE that the hessian means different things under the two objectives:
+ /// under it is p(1-p), at most 0.25 and
+ /// shrinking as the model sharpens, so a summed hessian is far below the sample
+ /// count; under it is the sample weight, so
+ /// with unit weights a threshold of 1.0 means exactly one sample. Carry the
+ /// hyper-parameters over from the reference run rather than assuming these defaults
+ /// transfer between objectives.
public double MinChildWeight = 1.0;
/// Row subsample fraction per tree (stochastic boosting).
public double Subsample = 0.8;
@@ -57,35 +85,90 @@ public sealed class GbtParams
public double RegLambda = 1.0;
/// L1 penalty on leaf weights (alpha).
public double RegAlpha = 0.0;
- /// Histogram bins per feature (<= 255 so bin indices fit a byte).
+ /// Histogram bins per feature. CLAMPED to [2, 255] at the start of training,
+ /// because a bin index has to fit a byte. Note that XGBoost's own max_bin
+ /// default is 256, so a value transcribed from a Python configuration trains with 255
+ /// here rather than the 256 it asks for.
public int MaxBins = 64;
/// Seed for the row/column subsampling PRNG. Drives
/// -- see the determinism note on
/// .
public ulong Seed = 42;
+ /// Threads used to accumulate histograms. Parallelism is applied ACROSS
+ /// FEATURES only, so every histogram is still summed in ascending row order by a
+ /// single thread and the trained model is bit-identical at any value. Defaults to
+ /// 1, which keeps the FDR path exactly sequential; raise it only for training sets
+ /// large enough that the accumulation dominates (millions of rows).
+ public int MaxDegreeOfParallelism = 1;
+ }
+
+ ///
+ /// A trained ensemble reduced to plain arrays, so callers that need to persist a model
+ /// can serialize it without reflecting over private state. Round-trips exactly: the
+ /// arrays ARE the model, and rebuilds
+ /// a scorer that returns bit-identical margins.
+ ///
+ /// Internal nodes have Feature in [0, FeatureCount) and branch on Threshold
+ /// (value <= Threshold goes to Left), with Left and Right both greater than the
+ /// node's own index and different from each other. Leaves have Feature == -1, carry
+ /// -1 in BOTH Left and Right, and contribute Leaf, already scaled by the learning
+ /// rate. TreeRoot holds the node index each tree starts at.
+ ///
+ /// A writer that omits Left/Right for leaves rather than emitting -1 will be rejected
+ /// on load; those two fields are part of the contract, not an implementation detail.
+ ///
+ public sealed class GbtModelData
+ {
+ public int[] Feature;
+ public double[] Threshold;
+ public int[] Left;
+ public int[] Right;
+ public double[] Leaf;
+ public int[] TreeRoot;
+ public double BaseScore;
+
+ /// Feature-vector width the model was trained on. Lets the load bounds-check
+ /// every split feature, so a corrupted index fails there rather than as an
+ /// index-out-of-range inside .
+ public int FeatureCount;
+
+ /// Objective the model was trained under. Recorded because the node arrays
+ /// alone cannot distinguish a squared-error margin from a log-odds one, and feeding a
+ /// reloaded regression margin to q-value or PEP estimation would be silently wrong.
+ public GbtObjective Objective;
}
///
- /// Gradient-boosted decision trees (Newton boosting, logistic loss) with L1/L2
- /// leaf regularization. Trained via ; scored via
- /// , which returns a raw log-odds margin.
+ /// Gradient-boosted decision trees (Newton boosting) with L1/L2 leaf regularization.
+ /// Trained via for binary
+ /// classification or for
+ /// regression; scored via , which returns the raw additive
+ /// margin.
///
/// DETERMINISTIC by construction, to the same standard as the linear SVM it stands
/// in for: identical input produces a bit-identical model and bit-identical scores,
- /// on every target framework. The pieces that guarantee it:
+ /// on every target framework and at any .
+ /// The pieces that guarantee it:
///
/// - subsampling draws from -- the same seeded,
/// bit-exact-by-definition PRNG the rest of Osprey.ML uses -- NOT
/// System.Random, whose seeded sequence is a framework implementation detail
/// (this builds net472 AND net8.0, so a divergence there would silently train two
/// different models from one source);
- /// - every float accumulation (histograms, leaf gradients/hessians) runs
- /// single-threaded in a fixed row order, so no summation order can drift;
+ /// - every float accumulation (histograms, leaf gradients/hessians) runs in a fixed
+ /// row order. Histogram work may be spread across threads, but only ACROSS FEATURES:
+ /// one thread owns a feature's histogram and walks the node's rows in ascending order,
+ /// so no summation order can drift with the thread count;
+ /// - row partitioning is stable, so each child sees its rows in the same relative
+ /// order they had in the parent;
+ /// - split selection scans features and bins in ascending order and takes a new
+ /// best only on a strict improvement, so ties resolve to the lowest (feature, bin);
/// - the one Array.Sort is over a primitive array read only by quantile
/// index, where equal values are interchangeable.
///
- /// Callers may train folds in parallel: each call owns its PRNG
- /// and touches no shared state. is pure and thread-safe.
+ /// Callers may train folds in parallel: each
+ /// call owns its PRNG and touches no shared state. is pure
+ /// and thread-safe.
///
public sealed class GradientBoostedTrees
{
@@ -99,33 +182,24 @@ public sealed class GradientBoostedTrees
private readonly double[] _leaf;
private readonly int[] _treeRoot;
private readonly double _baseScore;
+ private readonly int _featureCount;
+ private readonly GbtObjective _objective;
private GradientBoostedTrees(int[] feature, double[] threshold, int[] left, int[] right,
- double[] leaf, int[] treeRoot, double baseScore)
+ double[] leaf, int[] treeRoot, double baseScore, int featureCount, GbtObjective objective)
{
_feature = feature; _threshold = threshold; _left = left; _right = right;
_leaf = leaf; _treeRoot = treeRoot; _baseScore = baseScore;
+ _featureCount = featureCount; _objective = objective;
}
- /// Raw log-odds margin for one feature vector.
- public double ScoreSingle(double[] x)
- {
- double f = _baseScore;
- for (int t = 0; t < _treeRoot.Length; t++)
- {
- int node = _treeRoot[t];
- while (_feature[node] >= 0)
- node = x[_feature[node]] <= _threshold[node] ? _left[node] : _right[node];
- f += _leaf[node];
- }
- return f;
- }
+ /// Feature-vector width this model expects.
+ public int FeatureCount { get { return _featureCount; } }
- private static double Sigmoid(double z)
- {
- if (z >= 0) { double e = Math.Exp(-z); return 1.0 / (1.0 + e); }
- double ez = Math.Exp(z); return ez / (1.0 + ez);
- }
+ /// Objective this model was trained under. A margin from
+ /// is a prediction, not a log-odds, and must
+ /// not be handed to q-value or PEP estimation.
+ public GbtObjective Objective { get { return _objective; } }
///
/// Train on (rows = samples, cols = features) with binary
@@ -134,15 +208,118 @@ private static double Sigmoid(double z)
///
public static GradientBoostedTrees Train(double[][] x, bool[] isDecoy, GbtParams p, double[] sampleWeight = null)
{
+ if (isDecoy == null)
+ throw new ArgumentNullException(nameof(isDecoy));
+ if (p == null)
+ throw new ArgumentNullException(nameof(p));
+
+ // This overload promises a log-odds margin from binary labels, and callers rank
+ // by it. A GbtParams instance carried over from a regression call would otherwise
+ // quietly fit squared error to 0/1 targets and return something that is not a
+ // log-odds at all, which q-value and PEP estimation downstream would not survive.
+ if (p.Objective != GbtObjective.LogisticBinary)
+ {
+ throw new ArgumentException(string.Format(
+ @"GradientBoostedTrees.Train: the binary-label overload requires GbtObjective.LogisticBinary, not {0}. Use the continuous-target overload for regression.",
+ p.Objective));
+ }
+
+ var y = new double[isDecoy.Length];
+ for (int i = 0; i < isDecoy.Length; i++)
+ y[i] = isDecoy[i] ? 0.0 : 1.0;
+ return Train(x, y, p, sampleWeight);
+ }
+
+ ///
+ /// Train on (rows = samples, cols = features) against a
+ /// continuous target , using the loss named by
+ /// . Optional per-sample weights.
+ ///
+ public static GradientBoostedTrees Train(double[][] x, double[] y, GbtParams p, double[] sampleWeight = null)
+ {
+ if (x == null)
+ throw new ArgumentNullException(nameof(x));
+ if (p == null)
+ throw new ArgumentNullException(nameof(p));
int n = x.Length;
if (n == 0) throw new ArgumentException(@"GradientBoostedTrees.Train: empty training set");
+ if (y == null || y.Length != n)
+ throw new ArgumentException(@"GradientBoostedTrees.Train: target length must match the row count");
+
+ // Caught here rather than as an index-out-of-range partway through boosting,
+ // which would leave the caller guessing which array was the wrong length.
+ if (sampleWeight != null && sampleWeight.Length != n)
+ {
+ throw new ArgumentException(
+ @"GradientBoostedTrees.Train: sample weight length must match the row count");
+ }
+
+ // Under squared error the hessian IS the weight, with no positive floor of the
+ // kind the logistic branch applies. A negative weight can then drive a node's
+ // summed hessian to exactly -RegLambda and divide by zero in LeafValue, and one
+ // NaN leaf poisons every later round through the margin update. RegLambda is
+ // settable to 0 from the environment, so this is reachable in production.
+ if (sampleWeight != null)
+ {
+ for (int i = 0; i < n; i++)
+ {
+ if (double.IsNaN(sampleWeight[i]) || double.IsInfinity(sampleWeight[i]) || sampleWeight[i] < 0)
+ {
+ throw new ArgumentException(string.Format(
+ @"GradientBoostedTrees.Train: sample weight {0} is {1}; weights must be finite and non-negative",
+ i, sampleWeight[i]));
+ }
+ }
+ }
+
+ switch (p.Objective)
+ {
+ case GbtObjective.LogisticBinary:
+ // The logistic gradient is sigmoid(f) - y, which only means anything for a
+ // y in [0, 1]. Leaving Objective at its default and passing a continuous
+ // target through this overload would otherwise train a finite-looking but
+ // meaningless model with no diagnostic at all.
+ for (int i = 0; i < n; i++)
+ {
+ if (double.IsNaN(y[i]) || y[i] < 0.0 || y[i] > 1.0)
+ {
+ throw new ArgumentException(string.Format(
+ @"GradientBoostedTrees.Train: target {0} is {1}; GbtObjective.LogisticBinary requires targets in [0, 1]. Set Objective to SquaredError for regression.",
+ i, y[i]));
+ }
+ }
+
+ break;
+
+ case GbtObjective.SquaredError:
+ for (int i = 0; i < n; i++)
+ {
+ if (double.IsNaN(y[i]) || double.IsInfinity(y[i]))
+ {
+ throw new ArgumentException(string.Format(
+ @"GradientBoostedTrees.Train: target {0} is {1}; targets must be finite", i, y[i]));
+ }
+ }
+
+ break;
+
+ default:
+ // An out-of-range cast would otherwise fall through to the logistic branch
+ // and train the wrong loss silently.
+ throw new ArgumentException(string.Format(
+ @"GradientBoostedTrees.Train: unknown objective {0}", p.Objective));
+ }
+
int nFeat = x[0].Length;
int maxBins = Math.Max(2, Math.Min(255, p.MaxBins));
+ if ((long)n * nFeat > int.MaxValue)
+ throw new ArgumentException(@"GradientBoostedTrees.Train: training matrix exceeds the addressable bin array size");
// --- 1. Quantile bin edges per feature; precompute byte bin indices ---
+ // Bins are stored column-major so histogram accumulation for one feature walks
+ // a contiguous run instead of striding across one small array per row.
var cuts = new double[nFeat][];
- var bin = new byte[n][];
- for (int i = 0; i < n; i++) bin[i] = new byte[nFeat];
+ var bin = new byte[n * nFeat];
var col = new double[n];
for (int j = 0; j < nFeat; j++)
{
@@ -153,22 +330,29 @@ public static GradientBoostedTrees Train(double[][] x, bool[] isDecoy, GbtParams
}
cuts[j] = QuantileCuts(col, maxBins);
var cj = cuts[j];
+ int colStart = j * n;
for (int i = 0; i < n; i++)
- bin[i][j] = (byte)BinOf(cj, x[i][j]);
+ bin[colStart + i] = (byte)BinOf(cj, x[i][j]);
}
- // --- 2. Labels, weights, base score ---
- var y = new double[n];
+ // --- 2. Weights and base score ---
var w = sampleWeight;
double pos = 0, tot = 0;
for (int i = 0; i < n; i++)
{
- y[i] = isDecoy[i] ? 0.0 : 1.0;
double wi = w != null ? w[i] : 1.0;
pos += y[i] * wi; tot += wi;
}
- double frac = tot > 0 ? Math.Min(Math.Max(pos / tot, 1e-6), 1 - 1e-6) : 0.5;
- double baseScore = Math.Log(frac / (1 - frac));
+ double baseScore;
+ if (p.Objective == GbtObjective.SquaredError)
+ {
+ baseScore = tot > 0 ? pos / tot : 0.0;
+ }
+ else
+ {
+ double frac = tot > 0 ? Math.Min(Math.Max(pos / tot, 1e-6), 1 - 1e-6) : 0.5;
+ baseScore = Math.Log(frac / (1 - frac));
+ }
var f = new double[n];
for (int i = 0; i < n; i++) f[i] = baseScore;
@@ -185,15 +369,29 @@ public static GradientBoostedTrees Train(double[][] x, bool[] isDecoy, GbtParams
var allFeat = new int[nFeat];
for (int j = 0; j < nFeat; j++) allFeat[j] = j;
+ var workspace = new TreeWorkspace(n, nColUse, maxBins, p, bin, g, h);
+
// --- 3. Boosting rounds ---
for (int t = 0; t < p.NTrees; t++)
{
- for (int i = 0; i < n; i++)
+ if (p.Objective == GbtObjective.SquaredError)
{
- double pi = Sigmoid(f[i]);
- double wi = w != null ? w[i] : 1.0;
- g[i] = (pi - y[i]) * wi;
- h[i] = Math.Max(pi * (1 - pi) * wi, 1e-6);
+ for (int i = 0; i < n; i++)
+ {
+ double wi = w != null ? w[i] : 1.0;
+ g[i] = (f[i] - y[i]) * wi;
+ h[i] = wi;
+ }
+ }
+ else
+ {
+ for (int i = 0; i < n; i++)
+ {
+ double pi = Sigmoid(f[i]);
+ double wi = w != null ? w[i] : 1.0;
+ g[i] = (pi - y[i]) * wi;
+ h[i] = Math.Max(pi * (1 - pi) * wi, 1e-6);
+ }
}
// Row subsample (paired grouping is enforced upstream in fold assignment).
@@ -201,11 +399,14 @@ public static GradientBoostedTrees Train(double[][] x, bool[] isDecoy, GbtParams
// Column subsample for this tree.
var feats = SampleColumns(allFeat, nColUse, rng);
- int root = BuildTree(rows, 0, bin, cuts, feats, g, h, p, maxBins,
+ workspace.Reset(rows, feats);
+ int root = BuildTree(workspace, 0, rows.Length, 0, cuts, p,
nodesFeature, nodesThresh, nodesLeft, nodesRight, nodesLeaf);
treeRoots.Add(root);
- // Update margins for ALL samples with the new tree.
+ // Update margins for ALL samples with the new tree. The walk compares raw
+ // feature values, not bins, so a NaN feature takes the right branch here
+ // even though binning maps it to bin 0.
for (int i = 0; i < n; i++)
{
int node = root;
@@ -217,40 +418,259 @@ public static GradientBoostedTrees Train(double[][] x, bool[] isDecoy, GbtParams
return new GradientBoostedTrees(nodesFeature.ToArray(), nodesThresh.ToArray(),
nodesLeft.ToArray(), nodesRight.ToArray(), nodesLeaf.ToArray(),
- treeRoots.ToArray(), baseScore);
+ treeRoots.ToArray(), baseScore, nFeat, p.Objective);
+ }
+
+ /// Raw additive margin for one feature vector: a log-odds under
+ /// , the prediction itself under
+ /// .
+ public double ScoreSingle(double[] x)
+ {
+ if (x == null)
+ throw new ArgumentNullException(nameof(x));
+
+ // A short vector would otherwise read past the caller's array only for whichever
+ // features the traversal happens to touch, so the failure would depend on the
+ // data rather than on the mistake.
+ if (x.Length < _featureCount)
+ {
+ throw new ArgumentException(string.Format(
+ @"GradientBoostedTrees.ScoreSingle: model expects {0} features, got {1}",
+ _featureCount, x.Length));
+ }
+
+ double f = _baseScore;
+ for (int t = 0; t < _treeRoot.Length; t++)
+ {
+ int node = _treeRoot[t];
+ while (_feature[node] >= 0)
+ node = x[_feature[node]] <= _threshold[node] ? _left[node] : _right[node];
+ f += _leaf[node];
+ }
+ return f;
}
- // Recursively build one tree; appends nodes to the shared flat lists, returns the
- // node index of this subtree's root.
- private static int BuildTree(int[] rows, int depth, byte[][] bin, double[][] cuts,
- int[] feats, double[] g, double[] h, GbtParams p, int maxBins,
+ /// Flatten this model to plain arrays for persistence. The arrays are
+ /// copies; mutating them does not affect this instance.
+ public GbtModelData ToModelData()
+ {
+ return new GbtModelData
+ {
+ Feature = (int[])_feature.Clone(),
+ Threshold = (double[])_threshold.Clone(),
+ Left = (int[])_left.Clone(),
+ Right = (int[])_right.Clone(),
+ Leaf = (double[])_leaf.Clone(),
+ TreeRoot = (int[])_treeRoot.Clone(),
+ BaseScore = _baseScore,
+ FeatureCount = _featureCount,
+ Objective = _objective
+ };
+ }
+
+ /// Rebuild a scorer from persisted arrays. Validates the node graph rather
+ /// than trusting it: a truncated or hand-edited model file would otherwise surface
+ /// as an index-out-of-range deep inside scoring, or worse, as silently wrong
+ /// scores.
+ public static GradientBoostedTrees FromModelData(GbtModelData data)
+ {
+ if (data == null)
+ throw new ArgumentNullException(nameof(data));
+ if (data.Feature == null || data.Threshold == null || data.Left == null ||
+ data.Right == null || data.Leaf == null || data.TreeRoot == null)
+ {
+ throw new ArgumentException(@"GradientBoostedTrees.FromModelData: incomplete model data");
+ }
+
+ int nodes = data.Feature.Length;
+ if (data.Threshold.Length != nodes || data.Left.Length != nodes ||
+ data.Right.Length != nodes || data.Leaf.Length != nodes)
+ {
+ throw new ArgumentException(@"GradientBoostedTrees.FromModelData: node arrays must be the same length");
+ }
+ if (data.TreeRoot.Length == 0)
+ throw new ArgumentException(@"GradientBoostedTrees.FromModelData: model has no trees");
+ if (data.FeatureCount <= 0)
+ throw new ArgumentException(@"GradientBoostedTrees.FromModelData: feature count must be positive");
+
+ for (int i = 0; i < nodes; i++)
+ {
+ // Bounds-checking the split feature here is what keeps a corrupted index from
+ // surfacing as an index-out-of-range inside ScoreSingle instead.
+ if (data.Feature[i] >= data.FeatureCount)
+ {
+ throw new ArgumentException(string.Format(
+ @"GradientBoostedTrees.FromModelData: node {0} splits on feature {1}, outside the {2} the model was trained on",
+ i, data.Feature[i], data.FeatureCount));
+ }
+
+ if (data.Feature[i] < 0)
+ {
+ // A leaf owns no children. Rejecting stale indices here stops a partial
+ // edit from leaving a node that scores as a leaf but still points somewhere.
+ if (data.Left[i] != -1 || data.Right[i] != -1)
+ {
+ throw new ArgumentException(string.Format(
+ @"GradientBoostedTrees.FromModelData: leaf {0} carries a child index", i));
+ }
+
+ continue;
+ }
+
+ if (data.Left[i] < 0 || data.Left[i] >= nodes || data.Right[i] < 0 || data.Right[i] >= nodes)
+ {
+ throw new ArgumentException(string.Format(
+ @"GradientBoostedTrees.FromModelData: node {0} has a child index outside the node array", i));
+ }
+
+ // BuildTree appends a node before recursing into its children, so a child
+ // index always exceeds its parent's and the two differ. Requiring that on load
+ // is what rules out a cycle: without it, corrupted data makes ScoreSingle spin
+ // forever instead of failing.
+ if (data.Left[i] <= i || data.Right[i] <= i || data.Left[i] == data.Right[i])
+ {
+ throw new ArgumentException(string.Format(
+ @"GradientBoostedTrees.FromModelData: node {0} has child indices {1} and {2} that do not both increase, which would make scoring cycle",
+ i, data.Left[i], data.Right[i]));
+ }
+ }
+ for (int t = 0; t < data.TreeRoot.Length; t++)
+ {
+ if (data.TreeRoot[t] < 0 || data.TreeRoot[t] >= nodes)
+ throw new ArgumentException(string.Format(
+ @"GradientBoostedTrees.FromModelData: tree {0} root is outside the node array", t));
+ }
+
+ return new GradientBoostedTrees((int[])data.Feature.Clone(), (double[])data.Threshold.Clone(),
+ (int[])data.Left.Clone(), (int[])data.Right.Clone(), (double[])data.Leaf.Clone(),
+ (int[])data.TreeRoot.Clone(), data.BaseScore, data.FeatureCount, data.Objective);
+ }
+
+ private static double Sigmoid(double z)
+ {
+ if (z >= 0) { double e = Math.Exp(-z); return 1.0 / (1.0 + e); }
+ double ez = Math.Exp(z); return ez / (1.0 + ez);
+ }
+
+ // Per-tree scratch that outlives the recursion: the row index permutation being
+ // partitioned in place, the pooled per-depth histograms, and the partition buffer.
+ // Pooling matters at scale. A fresh double[maxBins] pair per feature per node is
+ // millions of short-lived arrays on a multi-million-row training set.
+ private sealed class TreeWorkspace
+ {
+ // Only these two change per boosting round; everything else is fixed for the
+ // whole Train call and is therefore set once, in the constructor.
+ public int[] Rows;
+ public int[] Feats;
+
+ public readonly byte[] Bin;
+ public readonly int RowStride;
+ public readonly double[] G;
+ public readonly double[] H;
+ public readonly int MaxBins;
+ public readonly int MaxDegreeOfParallelism;
+
+ private readonly double[][] _gradHist;
+ private readonly double[][] _hessHist;
+ private readonly int[] _partition;
+ private readonly ParallelOptions _parallelOptions;
+
+ public TreeWorkspace(int n, int nColUse, int maxBins, GbtParams p,
+ byte[] bin, double[] g, double[] h)
+ {
+ MaxBins = maxBins;
+ MaxDegreeOfParallelism = Math.Max(1, p.MaxDegreeOfParallelism);
+ Bin = bin;
+ RowStride = n;
+ G = g;
+ H = h;
+ _partition = new int[n];
+
+ // Allocated once rather than at every node. A depth-6 tree has up to 63
+ // internal nodes, so per-node allocation would be ~12,600 throwaway objects
+ // per fold per boosting run.
+ _parallelOptions = MaxDegreeOfParallelism > 1
+ ? new ParallelOptions { MaxDegreeOfParallelism = MaxDegreeOfParallelism }
+ : null;
+
+ // One histogram buffer per depth that can still split: a node's histogram is
+ // dead once its split is chosen, so the two children share the next level's
+ // buffer in turn, and a node at MaxDepth is always a leaf and never builds one.
+ // Width is nColUse, not nFeat, because only the sampled columns are indexed.
+ int levels = Math.Max(1, p.MaxDepth);
+ _gradHist = new double[levels][];
+ _hessHist = new double[levels][];
+ for (int d = 0; d < levels; d++)
+ {
+ _gradHist[d] = new double[nColUse * maxBins];
+ _hessHist[d] = new double[nColUse * maxBins];
+ }
+ }
+
+ public ParallelOptions ParallelOptions
+ {
+ get { return _parallelOptions; }
+ }
+
+ public void Reset(int[] rows, int[] feats)
+ {
+ Rows = rows; Feats = feats;
+ }
+
+ // Indexed directly: AccumulateHistograms is only reached for a node that can
+ // still split, so depth is always below MaxDepth. Clamping instead would let an
+ // out-of-contract depth quietly alias a parent's buffer rather than throw.
+ public double[] GradHist(int depth)
+ {
+ return _gradHist[depth];
+ }
+
+ public double[] HessHist(int depth)
+ {
+ return _hessHist[depth];
+ }
+
+ public int[] PartitionBuffer
+ {
+ get { return _partition; }
+ }
+ }
+
+ // Recursively build one tree over rows [start, start + count) of the workspace's
+ // row permutation; appends nodes to the shared flat lists and returns the node
+ // index of this subtree's root.
+ private static int BuildTree(TreeWorkspace ws, int start, int count, int depth,
+ double[][] cuts, GbtParams p,
List nFeat, List nThr, List nLeft, List nRight, List nLeaf)
{
- double G = 0, H = 0;
- for (int r = 0; r < rows.Length; r++) { int i = rows[r]; G += g[i]; H += h[i]; }
+ int maxBins = ws.MaxBins;
+ var rows = ws.Rows;
+ var g = ws.G;
+ var h = ws.H;
+ double gSum = 0, hSum = 0;
+ for (int r = start; r < start + count; r++) { int i = rows[r]; gSum += g[i]; hSum += h[i]; }
- bool leaf = depth >= p.MaxDepth || rows.Length < 2 || H < 2 * p.MinChildWeight;
+ bool leaf = depth >= p.MaxDepth || count < 2 || hSum < 2 * p.MinChildWeight;
int bestFeat = -1, bestBin = -1;
double bestGain = p.Gamma; // require gain strictly above gamma
if (!leaf)
{
- double parentTerm = G * G / (H + p.RegLambda);
+ var hg = ws.GradHist(depth);
+ var hh = ws.HessHist(depth);
+ var feats = ws.Feats;
+ AccumulateHistograms(ws, start, count, depth, hg, hh, feats);
+
+ double parentTerm = gSum * gSum / (hSum + p.RegLambda);
for (int fi = 0; fi < feats.Length; fi++)
{
int j = feats[fi];
- var hg = new double[maxBins];
- var hh = new double[maxBins];
- for (int r = 0; r < rows.Length; r++)
- {
- int i = rows[r]; int b = bin[i][j];
- hg[b] += g[i]; hh[b] += h[i];
- }
+ int histStart = fi * maxBins;
double gl = 0, hl = 0;
for (int b = 0; b < maxBins - 1; b++)
{
- gl += hg[b]; hl += hh[b];
+ gl += hg[histStart + b]; hl += hh[histStart + b];
if (hl < 1e-12 && gl == 0) continue;
- double gr = G - gl, hr = H - hl;
+ double gr = gSum - gl, hr = hSum - hl;
if (hl < p.MinChildWeight || hr < p.MinChildWeight) continue;
double gain = 0.5 * (gl * gl / (hl + p.RegLambda) + gr * gr / (hr + p.RegLambda) - parentTerm) - p.Gamma;
if (gain > bestGain) { bestGain = gain; bestFeat = j; bestBin = b; }
@@ -263,26 +683,97 @@ private static int BuildTree(int[] rows, int depth, byte[][] bin, double[][] cut
{
int idx = nFeat.Count;
nFeat.Add(-1); nThr.Add(0); nLeft.Add(-1); nRight.Add(-1);
- nLeaf.Add(LeafValue(G, H, p));
+ nLeaf.Add(LeafValue(gSum, hSum, p));
return idx;
}
- // Partition rows by the chosen bin threshold.
- var left = new List(); var right = new List();
- for (int r = 0; r < rows.Length; r++)
- {
- int i = rows[r];
- if (bin[i][bestFeat] <= bestBin) left.Add(i); else right.Add(i);
- }
+ int leftCount = Partition(ws, start, count, bestFeat, bestBin);
int self = nFeat.Count;
nFeat.Add(bestFeat); nThr.Add(cuts[bestFeat][bestBin]); nLeft.Add(-1); nRight.Add(-1); nLeaf.Add(0);
- int lc = BuildTree(left.ToArray(), depth + 1, bin, cuts, feats, g, h, p, maxBins, nFeat, nThr, nLeft, nRight, nLeaf);
- int rc = BuildTree(right.ToArray(), depth + 1, bin, cuts, feats, g, h, p, maxBins, nFeat, nThr, nLeft, nRight, nLeaf);
+ int lc = BuildTree(ws, start, leftCount, depth + 1, cuts, p, nFeat, nThr, nLeft, nRight, nLeaf);
+ int rc = BuildTree(ws, start + leftCount, count - leftCount, depth + 1, cuts, p, nFeat, nThr, nLeft, nRight, nLeaf);
nLeft[self] = lc; nRight[self] = rc;
return self;
}
+ // Fill this depth's pooled histogram with the node's gradient and hessian sums per
+ // (sampled feature, bin). One thread owns a feature and walks the node's rows in
+ // ascending order, so the sums do not depend on the thread count.
+ // Smallest node worth handing to the scheduler, in row-by-feature accumulation steps.
+ // Node population halves at every level, so most nodes in a depth-6 tree are far too
+ // small to repay a parallel dispatch; the million-row cost this exists for lives in
+ // the handful of wide nodes near the root.
+ private const long PARALLEL_WORK_THRESHOLD = 1L << 16;
+
+ private static void AccumulateHistograms(TreeWorkspace ws, int start, int count, int depth,
+ double[] hg, double[] hh, int[] feats)
+ {
+ int used = feats.Length * ws.MaxBins;
+ Array.Clear(hg, 0, used);
+ Array.Clear(hh, 0, used);
+
+ // NOTE on the choice of Parallel.For over this assembly's own OspreyParallel.For:
+ // OspreyParallel allocates dedicated Threads per call, which is right at fold
+ // granularity (a handful of calls per run) but not here, where a call happens at
+ // every internal node. Thread creation would swamp the work it parallelizes.
+ // The trade-off is that this path uses the shared ThreadPool, so it must not be
+ // turned on inside fold-parallel training; MaxDegreeOfParallelism defaults to 1
+ // and the FDR path never enters this branch.
+ if (ws.MaxDegreeOfParallelism <= 1 ||
+ (long)count * feats.Length < PARALLEL_WORK_THRESHOLD)
+ {
+ for (int fi = 0; fi < feats.Length; fi++)
+ AccumulateFeature(ws, start, count, fi, hg, hh);
+ return;
+ }
+
+ Parallel.For(0, feats.Length, ws.ParallelOptions,
+ fi => AccumulateFeature(ws, start, count, fi, hg, hh));
+ }
+
+ private static void AccumulateFeature(TreeWorkspace ws, int start, int count, int fi,
+ double[] hg, double[] hh)
+ {
+ var rows = ws.Rows;
+ var bin = ws.Bin;
+ var g = ws.G;
+ var h = ws.H;
+ int colStart = ws.Feats[fi] * ws.RowStride;
+ int histStart = fi * ws.MaxBins;
+ for (int r = start; r < start + count; r++)
+ {
+ int i = rows[r];
+ int b = bin[colStart + i];
+ hg[histStart + b] += g[i];
+ hh[histStart + b] += h[i];
+ }
+ }
+
+ // Stable in-place partition of rows [start, start + count) around the chosen bin.
+ // Rows keep their relative order on both sides, so each child sees exactly the row
+ // sequence the previous list-building implementation produced.
+ private static int Partition(TreeWorkspace ws, int start, int count, int bestFeat, int bestBin)
+ {
+ var rows = ws.Rows;
+ var bin = ws.Bin;
+ var right = ws.PartitionBuffer;
+ int colStart = bestFeat * ws.RowStride;
+
+ int leftCount = 0, rightCount = 0;
+ for (int r = start; r < start + count; r++)
+ {
+ int i = rows[r];
+ if (bin[colStart + i] <= bestBin)
+ rows[start + leftCount++] = i;
+ else
+ right[rightCount++] = i;
+ }
+
+ Array.Copy(right, 0, rows, start + leftCount, rightCount);
+ return leftCount;
+ }
+
// Optimal leaf weight with L1 soft-threshold + L2 shrinkage, times learning rate.
private static double LeafValue(double g, double h, GbtParams p)
{
diff --git a/pwiz_tools/Osprey/Osprey.Test/MLTest.cs b/pwiz_tools/Osprey/Osprey.Test/MLTest.cs
index 014b20b2ee..fdc1917fcf 100644
--- a/pwiz_tools/Osprey/Osprey.Test/MLTest.cs
+++ b/pwiz_tools/Osprey/Osprey.Test/MLTest.cs
@@ -499,6 +499,295 @@ public void TestGbtSingleClassNoCrash()
Assert.IsFalse(double.IsNaN(s) || double.IsInfinity(s));
}
+ ///
+ /// The squared-error objective: it must fit a continuous target, be invariant to
+ /// the histogram thread count, and, the part that protects the FDR path, leave
+ /// the logistic objective producing exactly the scores it produced before
+ /// regression existed. Everything but the base score and the per-round gradient is
+ /// shared between the two objectives, so a change to split finding, binning, or row
+ /// partitioning made for the benefit of regression would silently move every
+ /// Percolator-replacement q-value.
+ ///
+ [TestMethod]
+ public void TestGbtSquaredErrorObjective()
+ {
+ AssertLogisticGoldenUnchanged(1);
+
+ // The class documents bit-identity "at any MaxDegreeOfParallelism", and the FDR
+ // path is what that promise protects. Asserting it only under squared error
+ // would leave the logistic path unverified against the parallel branch.
+ AssertLogisticGoldenUnchanged(4);
+
+ AssertRegressionFitsContinuousTarget();
+ AssertRegressionThreadInvariant();
+ }
+
+ // Owned by the golden assertion rather than shared with the other GBT tests: these
+ // constants were captured under exactly these parameters, so a future tweak to
+ // FastGbt to speed the suite up must not fail here blaming the boosting code.
+ private static GbtParams GoldenGbt()
+ {
+ return new GbtParams { NTrees = 40, MaxDepth = 4, Seed = 42 };
+ }
+
+ ///
+ /// Golden logistic scores, captured from the implementation that predates the
+ /// objective switch (pwiz commit dd9e84581). Exact equality, not a tolerance: the
+ /// guarantee being protected is bit-identity, and a tolerance would let real
+ /// numerical drift through.
+ ///
+ private static void AssertLogisticGoldenUnchanged(int threads)
+ {
+ double[][] x = GoldenFeatures();
+ var isDecoy = new bool[x.Length];
+ for (int i = 0; i < x.Length; i++)
+ isDecoy[i] = x[i][0] + x[i][3] < 2.5;
+
+ var parameters = GoldenGbt();
+ parameters.MaxDegreeOfParallelism = threads;
+ var model = GradientBoostedTrees.Train(x, isDecoy, parameters);
+
+ var expected = new[]
+ {
+ -2.4677282578605486,
+ 4.097252275267094,
+ -0.8225901814943258,
+ 3.5439334614396674,
+ 2.610757110470939,
+ 2.677520123143121,
+ 2.2090870478970692,
+ 1.7584981803881718,
+ 1.7584981803881718,
+ 1.4031884801190766
+ };
+ double[][] probes = GoldenProbes();
+ for (int k = 0; k < probes.Length; k++)
+ {
+ Assert.AreEqual(expected[k], model.ScoreSingle(probes[k]), string.Format(
+ @"logistic score for probe {0} at {1} thread(s) moved: the shared boosting code is no longer bit-identical",
+ k, threads));
+ }
+ }
+
+ ///
+ /// Regression must actually reduce error against the weighted-mean baseline the
+ /// model starts from. A model that returned its base score everywhere would land
+ /// exactly on that baseline, so the margin here is the capacity check.
+ ///
+ private static void AssertRegressionFitsContinuousTarget()
+ {
+ double[][] x = GoldenFeatures();
+ double[] y = RegressionTarget(x);
+
+ var model = GradientBoostedTrees.Train(x, y, RegressionGbt());
+
+ double mean = 0;
+ for (int i = 0; i < y.Length; i++)
+ mean += y[i];
+ mean /= y.Length;
+
+ double sse = 0, sseBaseline = 0;
+ for (int i = 0; i < x.Length; i++)
+ {
+ double residual = y[i] - model.ScoreSingle(x[i]);
+ sse += residual * residual;
+ sseBaseline += (y[i] - mean) * (y[i] - mean);
+ }
+
+ Assert.IsFalse(double.IsNaN(sse) || double.IsInfinity(sse), @"regression residuals must be finite");
+ Assert.IsTrue(sse < 0.05 * sseBaseline, string.Format(
+ @"squared error {0} must be well below the weighted-mean baseline {1}", sse, sseBaseline));
+ }
+
+ ///
+ /// Histogram accumulation may be spread across threads, but only across features,
+ /// so the model must not depend on how many are used. This is the property that
+ /// makes the parallel path safe to turn on for large training sets.
+ ///
+ private static void AssertRegressionThreadInvariant()
+ {
+ double[][] x = GoldenFeatures();
+ double[] y = RegressionTarget(x);
+
+ var sequential = GradientBoostedTrees.Train(x, y, RegressionGbt());
+ var parallelParams = RegressionGbt();
+ parallelParams.MaxDegreeOfParallelism = 4;
+ var concurrent = GradientBoostedTrees.Train(x, y, parallelParams);
+
+ double[][] probes = GoldenProbes();
+ for (int k = 0; k < probes.Length; k++)
+ {
+ Assert.AreEqual(sequential.ScoreSingle(probes[k]), concurrent.ScoreSingle(probes[k]), string.Format(
+ @"regression score for probe {0} depends on the histogram thread count", k));
+ }
+ }
+
+ // Carries hyper-parameters from a reference XGBoost regression run rather than the
+ // logistic defaults: under squared error the hessian is the sample weight, so a
+ // MinChildWeight of 1.0 means one sample, not the several it means under logistic.
+ private static GbtParams RegressionGbt()
+ {
+ return new GbtParams
+ {
+ Objective = GbtObjective.SquaredError,
+ NTrees = 60,
+ MaxDepth = 6,
+ LearningRate = 0.1,
+ MinChildWeight = 1.0,
+ Subsample = 1.0,
+ ColSample = 1.0,
+ MaxBins = 256,
+ RegLambda = 1.0,
+ Seed = 42
+ };
+ }
+
+ private static double[] RegressionTarget(double[][] x)
+ {
+ var y = new double[x.Length];
+ for (int i = 0; i < x.Length; i++)
+ y[i] = (0.31 * x[i][0]) - (0.12 * x[i][3]) + (0.02 * x[i][2] * x[i][2]);
+ return y;
+ }
+
+ private static double[][] GoldenFeatures()
+ {
+ var x = new double[96][];
+ for (int i = 0; i < x.Length; i++)
+ x[i] = new[] { (i % 7) * 0.5, (i % 3) * 1.3, i * 0.1, ((i * 13) % 11) * 0.4 };
+ return x;
+ }
+
+ ///
+ /// A model flattened to and rebuilt must score
+ /// identically. Persistence that changed a margin would move q-values on reload.
+ /// The validation on the way back in matters just as much: a truncated or
+ /// hand-edited model file must fail at load rather than score silently wrong.
+ ///
+ [TestMethod]
+ public void TestGbtModelDataRoundTrip()
+ {
+ double[][] x = GoldenFeatures();
+ var isDecoy = new bool[x.Length];
+ for (int i = 0; i < x.Length; i++)
+ isDecoy[i] = x[i][0] + x[i][3] < 2.5;
+
+ var model = GradientBoostedTrees.Train(x, isDecoy, FastGbt());
+ var reloaded = GradientBoostedTrees.FromModelData(model.ToModelData());
+
+ double[][] probes = GoldenProbes();
+ for (int k = 0; k < probes.Length; k++)
+ {
+ Assert.AreEqual(model.ScoreSingle(probes[k]), reloaded.ScoreSingle(probes[k]), string.Format(
+ @"round-tripped model scores differ at probe {0}", k));
+ }
+
+ // Mutating the snapshot must not reach back into the model it came from. The
+ // reference score has to be taken BEFORE the mutation: comparing two values that
+ // are both recomputed afterwards passes even when ToModelData hands out the live
+ // internal arrays, which is the bug this is here to catch.
+ double beforeMutation = model.ScoreSingle(probes[0]);
+ var data = model.ToModelData();
+ data.Leaf[0] = 12345.0;
+ data.Feature[0] = 0;
+ Assert.AreEqual(beforeMutation, model.ScoreSingle(probes[0]));
+
+ var truncated = model.ToModelData();
+ truncated.Leaf = new double[truncated.Leaf.Length - 1];
+ Assert.ThrowsException(() => GradientBoostedTrees.FromModelData(truncated));
+
+ var dangling = model.ToModelData();
+ for (int i = 0; i < dangling.Feature.Length; i++)
+ {
+ if (dangling.Feature[i] < 0)
+ continue;
+ dangling.Left[i] = dangling.Feature.Length + 10;
+ break;
+ }
+ Assert.ThrowsException(() => GradientBoostedTrees.FromModelData(dangling));
+
+ var noTrees = model.ToModelData();
+ noTrees.TreeRoot = new int[0];
+ Assert.ThrowsException(() => GradientBoostedTrees.FromModelData(noTrees));
+
+ // A child index that does not increase is a cycle. Range checks alone let it
+ // through, and ScoreSingle would then spin forever rather than fail.
+ var cycle = model.ToModelData();
+ for (int i = 0; i < cycle.Feature.Length; i++)
+ {
+ if (cycle.Feature[i] < 0)
+ continue;
+ cycle.Left[i] = i;
+ break;
+ }
+ Assert.ThrowsException(() => GradientBoostedTrees.FromModelData(cycle));
+
+ var leafWithChild = model.ToModelData();
+ for (int i = 0; i < leafWithChild.Feature.Length; i++)
+ {
+ if (leafWithChild.Feature[i] >= 0)
+ continue;
+ leafWithChild.Left[i] = 0;
+ break;
+ }
+ Assert.ThrowsException(() => GradientBoostedTrees.FromModelData(leafWithChild));
+ }
+
+ ///
+ /// Argument validation on both Train overloads. Every case here would otherwise
+ /// surface as a NullReferenceException or an index-out-of-range partway through
+ /// boosting, pointing at boosting internals rather than at the bad argument.
+ ///
+ [TestMethod]
+ public void TestGbtTrainArgumentValidation()
+ {
+ double[][] x = GoldenFeatures();
+ var isDecoy = new bool[x.Length];
+ double[] y = RegressionTarget(x);
+
+ Assert.ThrowsException(
+ () => GradientBoostedTrees.Train(null, y, RegressionGbt()));
+ Assert.ThrowsException(
+ () => GradientBoostedTrees.Train(x, y, null));
+ Assert.ThrowsException(
+ () => GradientBoostedTrees.Train(x, (bool[])null, FastGbt()));
+ Assert.ThrowsException(
+ () => GradientBoostedTrees.Train(x, isDecoy, null));
+
+ Assert.ThrowsException(
+ () => GradientBoostedTrees.Train(x, new double[x.Length - 1], RegressionGbt()));
+ Assert.ThrowsException(
+ () => GradientBoostedTrees.Train(x, y, RegressionGbt(), new double[x.Length - 1]));
+
+ // The binary-label overload returns a log-odds margin that callers rank by, so it
+ // must refuse a params instance carried over from a regression call rather than
+ // silently fitting squared error to 0/1 labels.
+ var regressionParams = RegressionGbt();
+ Assert.ThrowsException(
+ () => GradientBoostedTrees.Train(x, isDecoy, regressionParams));
+
+ // The same instance is fine once the objective matches the entry point.
+ regressionParams.Objective = GbtObjective.LogisticBinary;
+ Assert.IsNotNull(GradientBoostedTrees.Train(x, isDecoy, regressionParams));
+ }
+
+ private static double[][] GoldenProbes()
+ {
+ return new[]
+ {
+ new[] { 0.0, 0.0, 0.0, 0.0 },
+ new[] { 3.0, 2.6, 9.5, 4.0 },
+ new[] { 1.5, 1.3, 6.3, 0.8 },
+ new[] { 2.5, 0.0, 1.1, 2.4 },
+ new[] { 0.5, 2.6, 5.0, 3.2 },
+ new[] { 2.0, 1.3, 2.2, 1.6 },
+ new[] { 1.0, 0.0, 7.7, 2.0 },
+ new[] { 3.0, 1.3, 0.4, 0.4 },
+ new[] { 100.0, -100.0, 0.0, 0.0 },
+ new[] { -5.0, 5.0, 50.0, 12.0 }
+ };
+ }
+
#endregion
#region FeatureStandardizer Tests
diff --git a/pwiz_tools/Osprey/docs/16-determinism.md b/pwiz_tools/Osprey/docs/16-determinism.md
index dbe813a0a1..5752290022 100644
--- a/pwiz_tools/Osprey/docs/16-determinism.md
+++ b/pwiz_tools/Osprey/docs/16-determinism.md
@@ -166,6 +166,32 @@ port also avoids nondeterministic parallel float reduction: the SVM training and
scoring reductions are per-row and per-lane deterministic, not a thread-order
`Parallel` sum.
+### The one parallel float accumulation, and why it is still deterministic
+
+`GradientBoostedTrees.AccumulateHistograms`
+(`Osprey.ML/GradientBoostedTrees.cs`) is the single place in the port that sums
+floats across threads, gated behind `GbtParams.MaxDegreeOfParallelism`. It is
+deterministic only because of one non-obvious invariant, which has to survive
+any future edit:
+
+> **Work is partitioned ACROSS FEATURES, never across rows.** One thread owns
+> one feature's histogram for the whole node and walks that node's rows in
+> ascending order, so every `hg[histStart + b] += g[i]` for a given bin happens
+> in the same sequence no matter how many threads run. No two threads ever touch
+> the same histogram slot, and there is no cross-thread reduction step.
+
+Partitioning over ROWS instead - the obvious "optimization" if the invariant is
+not understood - would make each bin's summation order depend on where the range
+boundaries fell, and would silently break the `1e-9` cross-impl parity gate that
+`regression.ps1` enforces. Sibling-subtraction histogram construction (deriving
+a child's histogram as parent minus its sibling) is off limits for the same
+reason: it changes the arithmetic, not merely its order.
+
+`MaxDegreeOfParallelism` defaults to **1**, so the FDR path stays sequential
+unless a caller opts in. `MLTest.TestGbtSquaredErrorObjective` asserts the ten
+golden logistic scores at both 1 and 4 threads, which is what keeps the claim
+above enforced rather than aspirational.
+
## Step 8 — Canonical entry order across Parquet and process boundaries
The deterministic entry order established during scoring must be reproduced after