API Reference

SCEBasis

Material + basis: crystal structure, symmetry, and SALC basis. The heavy part is the SALC construction; an SCEBasis can be persisted and reused. Cluster is a construction step, not a stored field — it is computed inside the constructor to build salcbasis, then discarded.

Fields

• structure::Structure: Crystal structure information.
• symmetry::Symmetry: Symmetry operations.
• salcbasis::SALCBasis: SALC basis functions.
• isotropy::Bool: Whether the basis was built with the isotropic restriction (only Lf = 0 terms). Provenance metadata: it is not derivable from salcbasis alone without inspecting every basis function.
• salc_fingerprint::UInt64: Structural fingerprint of salcbasis computed via SALCBases.salc_fingerprint(salcbasis). Used by the basis-identity check between an SCEModel / SCEFit and an SCEDataset so that a basis reloaded from disk still matches the in-memory original. Stable across Magesty.save / Magesty.load because the recipe hashes only integer structural identifiers; see SALCBases.salc_fingerprint for the included / excluded fields.

Examples

# Build from a TOML input file (most common path):
basis = SCEBasis("input.toml")

# The SALC construction is the expensive step; persist for reuse:
Magesty.save(basis, "basis.xml")
basis2 = Magesty.load(SCEBasis, "basis.xml")

SCEDataset

A SCEBasis paired with training data: the spin configurations and the unweighted design matrices and observation vectors derived from them.

X_E and X_T are stored unweighted. The torque weight is applied later at fit time, so a weight sweep reuses one SCEDataset without rebuilding design matrices.

Fields

• basis::SCEBasis: The SCE basis the design matrices were built from.
• spinconfigs::Vector{SpinConfig}: Training spin configurations.
• X_E::Matrix{Float64}: Energy design matrix, one row per configuration, one column per SALC. No bias column — the reference energy j0 is recovered analytically at fit time and is not a column of X_E. Unweighted.
• X_T::Matrix{Float64}: Torque design matrix, no bias column, a 3 * num_atoms block of rows per configuration. Unweighted.
• y_E::Vector{Float64}: Observed energies, length n_configs.
• y_T::Vector{Float64}: Observed torques, flattened, length 3 * num_atoms * n_configs.

Examples

# Pair a basis with training data read from an EMBSET file:
dataset = SCEDataset(basis, "EMBSET")

# Indexing yields a new SCEDataset — handy for train/test splits:
train = dataset[1:80]
test  = dataset[81:end]

# Concatenate datasets that share the same basis (e.g., for cross-validation):
combined = vcat(train, test)

SCEFit

A fitted SCE regression. Holds the dataset it was fit on, the fitted coefficients, the estimator and torque weight used, and the residuals of the augmented (weighted) least-squares system.

SCEFit <: StatsAPI.RegressionModel. The response-block-independent verbs coef, intercept, nobs, dof are defined for it.

Fields

• dataset::SCEDataset: The training dataset.
• j0::Float64: Fitted reference energy (bias term).
• jphi::Vector{Float64}: Fitted SCE coefficients.
• estimator::AbstractEstimator: The estimator used (OLS / Ridge).
• torque_weight::Float64: Torque weight in [0, 1] used at fit time.
• residuals::Vector{Float64}: Residuals of the augmented weighted least-squares system (energy rows stacked above flattened torque rows). These carry the torque_weight scaling and are not in physical units; use rmse_energy(f) / rmse_torque(f) / r2_energy(f) / r2_torque(f) for interpretable in-sample errors.

Examples

# Fit (see `fit(SCEFit, ...)` for full options):
f = fit(SCEFit, dataset, Ridge(lambda = 1e-4))

# Inspect:
coef(f)         # SCE coefficient vector
intercept(f)    # fitted reference energy j0 (eV)
rmse_torque(f)  # in-sample torque RMSE
r2_torque(f)    # in-sample torque R²

# Convert to the persistable predictor:
model = SCEModel(f)

SCEModel

A fitted SCE model: a SCEBasis together with the fitted reference energy j0 and SCE coefficients jphi. This is the lightweight, persistable predictor — predict_energy / predict_torque and the evaluation verbs accept it. Build one from a SCEFit via SCEModel(f).

Fields

• basis::SCEBasis: Structure, symmetry, and SALC basis the fit used.
• j0::Float64: Reference energy (bias term) in eV.
• jphi::Vector{Float64}: SCE coefficients.

The constructor checks that length(jphi) matches the number of SALCs in basis — there is one coefficient per SALC.

Examples

# From a fit:
model = SCEModel(f)

# Predict on any dataset (e.g., the held-out split):
ŷ_E = predict_energy(model, test)
ŷ_T = predict_torque(model, test)

# Persist and reload (round-trip is byte-stable):
Magesty.save(model, "model.xml")
model2 = Magesty.load(SCEModel, "model.xml")

fit(::Type{SCEFit}, dataset::SCEDataset, estimator::AbstractEstimator;
    torque_weight::Real = 1.0, verbosity::Bool = true) -> SCEFit

Fit SCE coefficients on dataset with estimator, returning a SCEFit.

torque_weight in [0, 1] sets the convex combination of the per-sample energy and torque mean squared errors that the augmented least-squares problem minimizes:

loss(jphi, j0) = (1 - torque_weight) * MSE_energy + torque_weight * MSE_torque

with

MSE_energy = (1 / n_E) * Σ_{i=1..n_E}            (y_E[i] - ŷ_E[i])^2
MSE_torque = (1 / n_T) * Σ_{k=1..n_T}            (y_T[k] - ŷ_T[k])^2

where n_E is the number of configurations, n_T = 3 * num_atoms * n_E is the total number of torque components, y_E / y_T are the observed energies and flattened torques, and ŷ_E / ŷ_T are the predictions of the model parameterised by (j0, jphi). Dividing each block by its sample count makes the two terms commensurate, so the convex combination is meaningful regardless of how many torque components a system has.

Limiting cases:

The default 1.0 is chosen on physical grounds: the SCE coefficients jphi are best determined by torque residuals, which carry the per-atom directional information that drives the response of the spin model. Energies enter the fit through the closed-form reference-energy recovery j0 = mean(y_E - X_E * jphi), so a torque-only fit still yields a usable j0. Set torque_weight < 1 only when the energies carry information that the torques do not — typically because the dataset is energy-rich and torque-poor.

Both MSEs are computed on the raw (unscaled) energy and torque residuals — the energy unit of the DFT input (typically eV) for the energy term, and eV per unit spin direction for the torque term. The two terms therefore live on the same physical scale (energy squared), so the default and the limiting cases are meaningful out of the box.

The design matrices stored in dataset are unweighted, so a torque_weight sweep reuses one SCEDataset without rebuilding them.

This is the StatsAPI fit verb; using Magesty re-exports it.

Arguments

• ::Type{SCEFit}: Target type — dispatch tag, written literally.
• dataset::SCEDataset: Training dataset (basis + design matrices + observations).
• estimator::AbstractEstimator: Regression estimator (OLS() or Ridge(lambda=...)).
• torque_weight::Real = 1.0: Convex weight described above.
• verbosity::Bool = true: Whether to print a summary of the fit (estimator, sizes, j0, in-sample RSS / RMSE / R² on energy and torque blocks, elapsed time) to stdout after solving.

Returns

• SCEFit: The fitted model. Inspect with coef, intercept, r2_energy, rmse_torque, …; persist with Magesty.save(SCEModel(f), path).

Examples

# Default: torque-only fit with Ridge regularization.
f = fit(SCEFit, dataset, Ridge(lambda = 1e-4))

# Energy-only fit, ignoring torques entirely.
f0 = fit(SCEFit, dataset, OLS(); torque_weight = 0.0)

# Reuse one dataset across a `torque_weight` sweep without rebuilding
# the design matrices.
for w in (0.0, 0.25, 0.5, 0.75, 1.0)
    println(w, "  ", rmse_torque(fit(SCEFit, dataset, OLS();
                                     torque_weight = w)))
end

refit(fit::SCEFit, estimator::AbstractEstimator = OLS();
      threshold::Real = 0.0, verbosity::Bool = true) -> SCEFit

Post-selection refit on the basis support of fit.

A basis index j enters the support iff abs(coef(fit)[j]) * norm(X[:, j]) > threshold, where X is the weighted, energy-centered design matrix used by fit (reassembled from fit.dataset with the input fit's torque_weight). The criterion is scale-equivariant: multiplying column j by any positive constant leaves the product invariant, because the OLS coefficient scales inversely while the column norm scales linearly. This cancels both the per-cluster (4π)^(N/2) factor and any other column-scale convention, so a single criterion serves both selector families: L1 estimators (Lasso, AdaptiveLasso) have exact-zero coefficients that fall out at the default threshold = 0.0, while AdaptiveRidge users pass a positive threshold to drop near-negligible bases.

The support columns are resolved on the same weighted problem with estimator; dropped bases stay at 0.0, so jphi keeps its full length and SALC ordering — a refit SCEFit round-trips through SCEModel and Magesty.save like any other fit. j0 is recovered from the uncentered energy data via extract_j0_jphi, identical to fit.

A PrecomputedPilot-backed estimator (including an AdaptiveLasso whose pilot is a PrecomputedPilot, such as AdaptiveLasso(::SCEFit)) is rejected upfront: the fixed pilot vector has the original column count, not the refit support length. A precomputed pilot is also meaningless once a support has been chosen.

Arguments

• fit::SCEFit: The fitted model whose coefficient support is reused. fit.dataset and fit.torque_weight are reused verbatim, so the refit minimizes the same weighted objective as the original fit.
• estimator::AbstractEstimator = OLS(): Estimator for the support resolve. Defaults to OLS() — classic post-selection debiasing. Pass Ridge(lambda = small) when the selected support is still rank-deficient or near-collinear.
• threshold::Real = 0.0: Scaled-magnitude cutoff. Bases with abs(coef(fit)[j]) * norm(X[:, j]) <= threshold are dropped. Must be >= 0. The cutoff is in the scale of the row-whitened augmented residual block (energy rows scaled by sqrt((1-torque_weight)/n_E), torque rows by sqrt(torque_weight/n_T)) rather than raw eV, so pick it relative to the largest scaled magnitude in the input fit, e.g. maximum(abs.(coef(fit)) .* [norm(@view X[:, j]) for j in axes(X, 2)]).
• verbosity::Bool = true: Whether to print the standard fit summary, exactly as fit does.

Returns

• SCEFit: A fresh fit; dropped bases carry coefficient 0.0. The returned estimator field is the refit estimator; the selection estimator is not recorded.

Examples

Both snippets assume dataset::SCEDataset is already in scope (built from an SCEBasis plus spin configurations).

# L1 selection + OLS debiasing. The default threshold = 0.0 keeps the
# Lasso's exact-zero support verbatim.
fit_lasso = fit(SCEFit, dataset, Lasso(lambda = 1e-3))
fit_db    = refit(fit_lasso, OLS())

# AdaptiveRidge fit: pass a positive threshold to drop near-zero bases,
# then refit with Ridge against residual collinearity.
fit_ar    = fit(SCEFit, dataset, AdaptiveRidge(lambda = 1e-2))
fit_db    = refit(fit_ar, Ridge(lambda = 1e-6); threshold = 1e-4)

coef(f::SCEFit) -> Vector{Float64}
coef(m::SCEModel) -> Vector{Float64}

The fitted SCE coefficients jphi — a single set shared by the energy and torque models.

intercept(f::SCEFit) -> Float64
intercept(m::SCEModel) -> Float64

The fitted reference energy j0 (bias term), in the energy unit of the DFT input (typically eV). The returned value is exactly the j0 field of the underlying SCEFit / SCEModel, kept separate from the SCE coefficients jphi; XML persistence stores it as the j0 attribute of the <JPhi> block. intercept is a Magesty-native verb adopted because StatsAPI has no intercept concept; the two names refer to the same quantity.

Returns

• Float64: the reference energy j0 (the bias term), in the energy unit of the DFT input.

nobs(f::SCEFit) -> Int

Number of observations (spin configurations) the fit was trained on, i.e. the energy-block observation count n_configs.

Returns

• Int: the number of training spin configurations.

dof(f::SCEFit) -> Int

Degrees of freedom consumed by the fit: length(coef(f)) + 1 — the SCE coefficients plus the intercept.

Returns

• Int: the number of SCE coefficients plus one for the intercept.

predict_energy(model::SCEModel, spin_directions::AbstractMatrix{<:Real}) -> Float64
predict_energy(model::SCEModel, sc::SpinConfig) -> Float64
predict_energy(model::SCEModel, sd_list::AbstractVector{<:AbstractMatrix{<:Real}}) -> Vector{Float64}
predict_energy(model::SCEModel, configs::AbstractVector{SpinConfig}) -> Vector{Float64}
predict_energy(model::SCEModel, dataset::SCEDataset) -> Vector{Float64}

f::SCEFit may be passed in place of model; the SCEFit overloads delegate through SCEModel(f).

Predict SCE energies for one or more spin configurations.

Arguments

• model::SCEModel or f::SCEFit: Trained predictor. SCEFit inputs delegate through SCEModel(f).
• spin_directions::AbstractMatrix{<:Real}: Spin direction matrix of size 3 × num_atoms (rows = x, y, z). Columns must be unit vectors (‖·‖ = 1 within 1e-6); a non-unit or non-finite column throws an ArgumentError.
• sc::SpinConfig: Single spin configuration; equivalent to passing sc.spin_directions.
• sd_list::AbstractVector{<:AbstractMatrix{<:Real}}: Sequence of spin direction matrices. The number of atoms must match across entries and the predictor's SCEBasis; the matrices themselves can be views or freshly allocated.
• configs::AbstractVector{SpinConfig}: Sequence of spin configurations (e.g. the output of read_embset); only spin_directions is read, the other fields are ignored.
• dataset::SCEDataset: Batch evaluation reusing the dataset's stored energy design matrix. Must share the predictor's SCEBasis (same (l, m, site) column ordering as the fitted coefficients).

Returns

Return type depends on the second positional argument:

predict_torque(model::SCEModel, spin_directions::AbstractMatrix{<:Real}) -> Matrix{Float64}
predict_torque(model::SCEModel, sc::SpinConfig) -> Matrix{Float64}
predict_torque(model::SCEModel, sd_list::AbstractVector{<:AbstractMatrix{<:Real}}) -> Vector{Matrix{Float64}}
predict_torque(model::SCEModel, configs::AbstractVector{SpinConfig}) -> Vector{Matrix{Float64}}
predict_torque(model::SCEModel, dataset::SCEDataset) -> Vector{Matrix{Float64}}

f::SCEFit may be passed in place of model; the SCEFit overloads delegate through SCEModel(f).

Predict per-atom SCE torques for one or more spin configurations.

Arguments

• model::SCEModel or f::SCEFit: Trained predictor. SCEFit inputs delegate through SCEModel(f).
• spin_directions::AbstractMatrix{<:Real}: Spin direction matrix of size 3 × num_atoms (rows = x, y, z). Columns must be unit vectors (‖·‖ = 1 within 1e-6); a non-unit or non-finite column throws an ArgumentError.
• sc::SpinConfig: Single spin configuration; equivalent to passing sc.spin_directions.
• sd_list::AbstractVector{<:AbstractMatrix{<:Real}}: Sequence of spin direction matrices. The number of atoms must match across entries and the predictor's SCEBasis.
• configs::AbstractVector{SpinConfig}: Sequence of spin configurations (e.g. the output of read_embset); only spin_directions is read.
• dataset::SCEDataset: Batch evaluation reusing the dataset's stored torque design matrix. Must share the predictor's SCEBasis (same (l, m, site) column ordering as the fitted coefficients).

Returns

Return type depends on the second positional argument:

r2_energy(predictor, data) -> Float64
r2_energy(f::SCEFit) -> Float64

Coefficient of determination (R²) of the SCE energy predictions.

Call as r2_energy(f) to evaluate in-sample on the training dataset embedded in f; call as r2_energy(predictor, data) to evaluate on a different dataset — typically a held-out test set.

Arguments

• predictor::Union{SCEModel, SCEFit}: Trained predictor.
• data: Evaluation data — SCEDataset, AbstractVector{SpinConfig}, or an EMBSET file path (AbstractString).
• f::SCEFit (single-argument form): Evaluates f in-sample on its own training dataset.

Returns

• Float64: R² over the observed and predicted energies (eV).

r2_torque(predictor, data) -> Float64
r2_torque(f::SCEFit) -> Float64

Coefficient of determination (R²) of the SCE torque predictions, computed over the flattened torque components.

Call as r2_torque(f) to evaluate in-sample on the training dataset embedded in f; call as r2_torque(predictor, data) to evaluate on a different dataset — typically a held-out test set.

Arguments

• predictor::Union{SCEModel, SCEFit}: Trained predictor.
• data: Evaluation data — SCEDataset, AbstractVector{SpinConfig}, or an EMBSET file path (AbstractString).
• f::SCEFit (single-argument form): Evaluates f in-sample on its own training dataset.

Returns

• Float64: R² over the flattened observed and predicted torques (eV), length 3 * num_atoms * n_configs.

rmse_energy(predictor, data) -> Float64
rmse_energy(f::SCEFit) -> Float64

Root mean squared error of the SCE energy predictions (eV).

Call as rmse_energy(f) to evaluate in-sample on the training dataset embedded in f; call as rmse_energy(predictor, data) to evaluate on a different dataset — typically a held-out test set.

Arguments

• predictor::Union{SCEModel, SCEFit}: Trained predictor.
• data: Evaluation data — SCEDataset, AbstractVector{SpinConfig}, or an EMBSET file path (AbstractString).
• f::SCEFit (single-argument form): Evaluates f in-sample on its own training dataset.

Returns

• Float64: sqrt(mean((observed - predicted).^2)) over the energies (eV).

rmse_torque(predictor, data) -> Float64
rmse_torque(f::SCEFit) -> Float64

Root mean squared error of the SCE torque predictions (eV), over the flattened torque components.

Call as rmse_torque(f) to evaluate in-sample on the training dataset embedded in f; call as rmse_torque(predictor, data) to evaluate on a different dataset — typically a held-out test set.

Arguments

• predictor::Union{SCEModel, SCEFit}: Trained predictor.
• data: Evaluation data — SCEDataset, AbstractVector{SpinConfig}, or an EMBSET file path (AbstractString).
• f::SCEFit (single-argument form): Evaluates f in-sample on its own training dataset.

Returns

• Float64: sqrt(mean((observed - predicted).^2)) over the flattened torques (eV), length 3 * num_atoms * n_configs.

rss_energy(predictor, data) -> Float64
rss_energy(f::SCEFit) -> Float64

Residual sum of squares of the SCE energy predictions.

Call as rss_energy(f) to evaluate in-sample on the training dataset embedded in f; call as rss_energy(predictor, data) to evaluate on a different dataset — typically a held-out test set.

Arguments

• predictor::Union{SCEModel, SCEFit}: Trained predictor.
• data: Evaluation data — SCEDataset, AbstractVector{SpinConfig}, or an EMBSET file path (AbstractString).
• f::SCEFit (single-argument form): Evaluates f in-sample on its own training dataset.

Returns

• Float64: sum((observed - predicted).^2) over the energies (eV²).

rss_torque(predictor, data) -> Float64
rss_torque(f::SCEFit) -> Float64

Residual sum of squares of the SCE torque predictions, over the flattened torque components.

Call as rss_torque(f) to evaluate in-sample on the training dataset embedded in f; call as rss_torque(predictor, data) to evaluate on a different dataset — typically a held-out test set.

Arguments

• predictor::Union{SCEModel, SCEFit}: Trained predictor.
• data: Evaluation data — SCEDataset, AbstractVector{SpinConfig}, or an EMBSET file path (AbstractString).
• f::SCEFit (single-argument form): Evaluates f in-sample on its own training dataset.

Returns

• Float64: sum((observed - predicted).^2) over the flattened torques (eV²), length 3 * num_atoms * n_configs.

residuals_energy(predictor, data) -> Vector{Float64}
residuals_energy(f::SCEFit) -> Vector{Float64}

Per-configuration energy residuals observed - predicted (eV).

Call as residuals_energy(f) to evaluate in-sample on the training dataset embedded in f; call as residuals_energy(predictor, data) to evaluate on a different dataset — typically a held-out test set.

Arguments

• predictor::Union{SCEModel, SCEFit}: Trained predictor.
• data: Evaluation data — SCEDataset, AbstractVector{SpinConfig}, or an EMBSET file path (AbstractString).
• f::SCEFit (single-argument form): Evaluates f in-sample on its own training dataset.

Returns

• Vector{Float64} of length n_configs, in input order.

residuals_torque(predictor, data) -> Vector{Float64}
residuals_torque(f::SCEFit) -> Vector{Float64}

Flattened torque residuals observed - predicted (eV).

Call as residuals_torque(f) to evaluate in-sample on the training dataset embedded in f; call as residuals_torque(predictor, data) to evaluate on a different dataset — typically a held-out test set.

Arguments

• predictor::Union{SCEModel, SCEFit}: Trained predictor.
• data: Evaluation data — SCEDataset, AbstractVector{SpinConfig}, or an EMBSET file path (AbstractString).
• f::SCEFit (single-argument form): Evaluates f in-sample on its own training dataset.

Returns

• Vector{Float64} of length 3 * num_atoms * n_configs, flattened over (component, atom, configuration).

write_energies(f::SCEFit, filename::AbstractString = "energy_list.txt")
write_energies(predictor::Union{SCEModel, SCEFit},
               data::Union{SCEDataset, AbstractVector{SpinConfig}, AbstractString},
               filename::AbstractString)

Write observed (DFT) versus predicted (SCE) energies to a text file for fit-quality inspection. The output is consumed by the FitCheck_energy.py visualization script under tools/.

The two-argument form evaluates f against its own training dataset; the lone string argument is always the output path. The three-argument form evaluates any predictor against an explicit dataset, e.g. a held-out validation or test set, and requires an explicit output path so the data and output arguments cannot be confused.

Arguments

• f::SCEFit: A fitted model; evaluated against f.dataset.
• predictor::Union{SCEModel, SCEFit}: The predictor producing SCE energies.
• data::Union{SCEDataset, AbstractVector{SpinConfig}, AbstractString}: The configurations supplying observed energies. An AbstractString is read as an EMBSET file path.
• filename::AbstractString: Output path. Defaults to "energy_list.txt" only in the two-argument SCEFit form.

Returns

• nothing. The file at filename is created (overwritten if it exists).

The file has a two-line comment header followed by one row per configuration: data_index DFT_Energy SCE_Energy. Energies are in the unit of the training data (typically eV).

Examples

f = fit(SCEFit, dataset, Ridge(lambda = 1e-4))
write_energies(f)                                       # -> energy_list.txt
write_energies(f, "train_E.txt")                        # training set
write_energies(f, test_dataset, "test_E.txt")           # held-out set
write_energies(SCEModel(f), "EMBSET", "E.txt")          # model + EMBSET path

write_torques(f::SCEFit, filename::AbstractString = "torque_list.txt")
write_torques(predictor::Union{SCEModel, SCEFit},
              data::Union{SCEDataset, AbstractVector{SpinConfig}, AbstractString},
              filename::AbstractString)

Write observed (DFT) versus predicted (SCE) per-atom torques to a text file for fit-quality inspection. The output is consumed by the FitCheck_torque.py visualization script under tools/.

The two-argument form evaluates f against its own training dataset; the lone string argument is always the output path. The three-argument form evaluates any predictor against an explicit dataset, e.g. a held-out validation or test set, and requires an explicit output path so the data and output arguments cannot be confused.

Arguments

• f::SCEFit: A fitted model; evaluated against f.dataset.
• predictor::Union{SCEModel, SCEFit}: The predictor producing SCE torques.
• data::Union{SCEDataset, AbstractVector{SpinConfig}, AbstractString}: The configurations supplying observed torques. An AbstractString is read as an EMBSET file path.
• filename::AbstractString: Output path. Defaults to "torque_list.txt" only in the two-argument SCEFit form.

Returns

• nothing. The file at filename is created (overwritten if it exists).

The file has a two-line comment header followed by a # data index: N block per configuration. Each block holds one row per atom: atom_index element DFT_torque_{x,y,z} SCE_torque_{x,y,z}. Torque components are in the unit of the training data (typically eV).

Examples

f = fit(SCEFit, dataset, Ridge(lambda = 1e-4))
write_torques(f)                                # -> torque_list.txt
write_torques(f, "train_T.txt")                 # training set
write_torques(f, test_dataset, "test_T.txt")    # held-out set
write_torques(SCEModel(f), "EMBSET", "T.txt")   # model + EMBSET path

gcv(f::SCEFit) -> Float64

Combined energy+torque generalized cross-validation score for the fitted model f, evaluated on its training dataset and the weighted objective f was fit with (same torque_weight and estimator).

GCV estimates the out-of-sample prediction error from a single fit via the hat matrix H (ŷ = H y):

GCV = (‖r‖² / N) / (1 − tr(H)/N)²,

where r is the augmented weighted residual, tr(H) the effective degrees of freedom, and N the number of live rows — energy plus torque, minus any block zeroed by the weighting (torque_weight = 1 drops the energy block, torque_weight = 0 drops the torque block). The eliminated reference energy j0 counts one degree of freedom only when the energy block is live (torque_weight < 1). The score is in the weighted-objective unit, not eV²; compare scores (e.g. across penalties or data sizes), not the absolute magnitude.

Defined only for linear estimators (OLS, Ridge, AdaptiveRidge).

Arguments

• f::SCEFit: A fitted model whose estimator is linear.

Returns

• Float64: The GCV score, or NaN if the model is numerically saturated (tr(H) ≥ N).

Throws

• ArgumentError if f.estimator is non-linear (ElasticNet / Lasso / AdaptiveLasso).

Examples

f = fit(SCEFit, dataset, Ridge(lambda = 1e-4))
gcv(f)

gcv_r2(f::SCEFit) -> Float64

GCV-based predictive R² for the fitted model f: the gcv score normalized against the null-model mean square msy = ‖y‖² / N, namely

R²_gcv = 1 − GCV / msy.

The null model is β = 0 on the weighted, energy-centered augmented system (energy predicted at its mean, torque predicted as zero), so R²_gcv measures the cross-validated variance explained on a fixed scale: 1 is a perfect fit, 0 matches the null model, and a negative value means the fit predicts worse than the null (over-parameterized / too little data). Unlike the raw gcv score — which is in the weighted-objective unit and only meaningful in relative comparison — this value can be read in isolation.

Defined only for linear estimators (OLS, Ridge, AdaptiveRidge).

Arguments

• f::SCEFit: A fitted model whose estimator is linear.

Returns

• Float64: The predictive R², or NaN if the model is numerically saturated (tr(H) ≥ N).

Throws

• ArgumentError if f.estimator is non-linear (ElasticNet / Lasso / AdaptiveLasso).

Examples

f = fit(SCEFit, dataset, Ridge(lambda = 1e-4))
gcv_r2(f)    # ~1 good, ~0 no better than the mean / zero-torque model

gcv_lambda(dataset::SCEDataset, lambdas::AbstractVector{<:Real};
           torque_weight::Real = 1.0) -> GCVLambdaPath

Ridge GCV penalty sweep: compute the combined energy+torque GCV score for every lambda and report the minimizer. A single SVD of the weighted, energy-centered design matrix serves the whole path, so passing a fine lambdas grid is cheap.

Arguments

• dataset::SCEDataset: The training data (design matrices built once).
• lambdas::AbstractVector{<:Real}: Non-negative ridge penalties to evaluate.
• torque_weight::Real = 1.0: Convex energy/torque weight in [0, 1], as in fit.

Returns

• GCVLambdaPath: Per-lambda GCV, predictive R² (gcv_r2), and effective dof, plus lambda_best.

Throws

• ArgumentError if lambdas is empty, contains a negative value, torque_weight is outside [0, 1], or every penalty gives a non-finite GCV.

Examples

path = gcv_lambda(dataset, 10.0 .^ (-6:0.5:0))
f    = fit(SCEFit, dataset, Ridge(lambda = path.lambda_best))

gcv_learning_curve(dataset::SCEDataset, estimator::AbstractEstimator = OLS();
                   sizes::AbstractVector{<:Integer} = <auto grid>,
                   repeats::Integer = 5, seed::Integer = 0,
                   torque_weight::Real = 1.0) -> GCVSizeCurve

Data-sufficiency GCV learning curve. At each training-set size, draw repeats random config subsets, fit estimator to each, and average their combined GCV scores; a curve that flattens with size indicates enough training data.

Subsets are drawn from the full dataset with a seeded RNG (reproducible). Each draw reuses the prebuilt design matrices via row slicing — the heavy design-matrix construction is not repeated. A draw that yields a rank-deficient OLS solve or a saturated model (tr(H) ≥ N) contributes NaN and is dropped from that size's statistics (with a warning); if every draw at a size fails, the size reports NaN.

Arguments

• dataset::SCEDataset: The full training data.
• estimator::AbstractEstimator = OLS(): Linear estimator fit on each subset.
• sizes::AbstractVector{<:Integer}: Training-set sizes. Defaults to six points spanning max(p + 2, 10) to length(dataset) (p = number of SALCs).
• repeats::Integer = 5: Random draws averaged per size.
• seed::Integer = 0: RNG seed.
• torque_weight::Real = 1.0: Convex energy/torque weight in [0, 1], as in fit.

Returns

• GCVSizeCurve: sizes, gcv_mean, gcv_std, the predictive-R² summary (gcv_r2_mean, gcv_r2_std), and the sweep settings.

Throws

• ArgumentError if estimator is non-linear, repeats < 1, torque_weight is outside [0, 1], or any size is outside 1:length(dataset).

Examples

curve = gcv_learning_curve(dataset, Ridge(lambda = 1e-4); repeats = 8)
curve.sizes        # training-set sizes
curve.gcv_mean     # mean GCV at each size — look for a plateau
curve.gcv_r2_mean  # mean predictive R² — ~1 good, plateau ⇒ enough data
curve.gcv_std      # spread across draws — large spread suggests more data needed

GCVLambdaPath

Result of a ridge GCV penalty sweep (gcv_lambda). The combined energy+torque GCV score is evaluated at each penalty lambda from a single SVD of the weighted design matrix.

Fields

• lambdas::Vector{Float64}: The penalty values swept, in input order.
• gcv_scores::Vector{Float64}: GCV score at each lambda (NaN where the model is numerically saturated).
• gcv_r2::Vector{Float64}: GCV-based predictive R² at each lambda, 1 - gcv / msy (1 perfect, 0 matches the null model, < 0 worse than null). Interpretable on a fixed scale, unlike the raw gcv_scores.
• dof::Vector{Float64}: Effective degrees of freedom tr(H) at each lambda.
• lambda_best::Float64: The lambda minimizing the GCV score.
• torque_weight::Float64: The torque weight the sweep used.

GCVSizeCurve

Result of a data-sufficiency GCV sweep (gcv_learning_curve). At each training-set size, repeats random config subsets are fit and scored; the mean and standard deviation across draws are reported, so a flattening curve signals that enough data is present.

Fields

• sizes::Vector{Int}: Training-set sizes, ascending.
• gcv_mean::Vector{Float64}: Mean GCV over the random draws at each size.
• gcv_std::Vector{Float64}: Standard deviation over the draws at each size.
• gcv_r2_mean::Vector{Float64}: Mean GCV-based predictive R² over the draws at each size (1 - gcv / msy; 1 perfect, 0 matches the null model). Read on a fixed scale, unlike the raw gcv_mean.
• gcv_r2_std::Vector{Float64}: Standard deviation of the predictive R² over the draws at each size.
• repeats::Int: Random draws per size.
• seed::Int: RNG seed used for reproducibility.
• estimator::AbstractEstimator: The estimator fit on each subset.
• torque_weight::Float64: The torque weight the sweep used.

write_gcv_lambda(path::GCVLambdaPath,
                 filename::AbstractString = "gcv_lambda.txt")

Write a ridge GCV penalty sweep (gcv_lambda) to a text file for inspection or plotting. The output is consumed by the FitCheck_gcv_lambda.py script under tools/.

The file has a comment header followed by one row per penalty: lambda gcv gcv_r2 dof. The GCV score is in the weighted-objective unit (not eV²), gcv_r2 is the GCV-based predictive R² (1 - gcv / msy; 1 perfect, 0 matches the null model) and reads on a fixed scale, and dof is the effective degrees of freedom tr(H). The selected lambda_best and the torque weight are recorded in the header.

Arguments

• path::GCVLambdaPath: The sweep result.
• filename::AbstractString: Output path. Defaults to "gcv_lambda.txt".

Returns

• nothing. The file at filename is created (overwritten if it exists). Any filesystem error is logged and re-thrown.

Examples

path = gcv_lambda(dataset, 10.0 .^ (-6:0.5:0))
write_gcv_lambda(path, "gcv_lambda.txt")

write_gcv_learning_curve(curve::GCVSizeCurve,
                         filename::AbstractString = "gcv_learning_curve.txt")

Write a data-sufficiency GCV learning curve (gcv_learning_curve) to a text file for inspection or plotting. The output is consumed by the FitCheck_gcv_learning_curve.py script under tools/.

The file has a comment header followed by one row per training-set size: size gcv_mean gcv_std gcv_r2_mean gcv_r2_std, where the means and standard deviations are taken over the random subset draws at that size. The GCV score is in the weighted-objective unit (not eV²); gcv_r2 is the GCV-based predictive R² (1 - gcv / msy), which reads on a fixed scale (1 perfect, 0 null). The estimator, torque weight, repeats, and seed are recorded in the header.

Arguments

• curve::GCVSizeCurve: The sweep result.
• filename::AbstractString: Output path. Defaults to "gcv_learning_curve.txt".

Returns

• nothing. The file at filename is created (overwritten if it exists). Any filesystem error is logged and re-thrown.

Examples

curve = gcv_learning_curve(dataset, Ridge(lambda = 1e-4); repeats = 8)
write_gcv_learning_curve(curve, "gcv_learning_curve.txt")

vasp_to_extxyz(vasprun; oszicar=nothing, output=nothing) -> String

Convert a VASP run to extended XYZ (extxyz) format and return the extxyz text.

When only vasprun is given, the extxyz carries structure, forces, stress, and energies. Passing oszicar additionally writes per-atom magnetic moments (magmom_smoothed, magmom_raw) and the constraint field (constr_field).

The header also records soc (spin-orbit coupling flag) and, when the electronic-convergence status can be determined, converged (T/F for the electronic SCF loop). A non-converged run is still written, with converged=F and a warning. The converged key is omitted when the status is indeterminate (no NELM, single-shot NELM = 1, EDIFF <= 0, or no SCF steps).

Arguments

• vasprun::AbstractString: path to vasprun.xml.

Keyword arguments

• oszicar::Union{AbstractString, Nothing} = nothing: path to OSZICAR; when given, the magnetic-moment and constraint-field columns are added.
• output::Union{AbstractString, Nothing} = nothing: when given, the extxyz text is also written to this file (.extxyz is appended if the name does not already end with it).

Returns

• String: the full extxyz text.

Throws

• VaspParseError: if vasprun (or oszicar) cannot be parsed. The error message names the offending file path.

Examples

text = vasp_to_extxyz("vasprun.xml")
vasp_to_extxyz("vasprun.xml"; oszicar = "OSZICAR", output = "frame.extxyz")

poscar_to_toml(poscar; output=nothing) -> String

Convert a VASP POSCAR structure file to a Magesty input TOML configuration and return the TOML text.

The generated configuration is a starting point for an SCE input file: it fills [general], [symmetry], [interaction], and [structure] from the POSCAR, with placeholder interaction settings (lmax = 0, cutoff = -1) meant to be edited before use.

Arguments

• poscar::AbstractString: path to a POSCAR structure file.

Keyword arguments

• output::Union{AbstractString, Nothing} = nothing: when given, the TOML text is also written to this file (.toml is appended if the name does not already end with it).

Returns

• String: the full TOML text.

Examples

text = poscar_to_toml("POSCAR")
poscar_to_toml("POSCAR"; output = "input.toml")

oszicar_to_embset(oszicars; saxis=[0.0, 0.0, 1.0], energy_kind="f", mint=false, output=nothing) -> String

Convert one or more VASP OSZICAR files to the EMBSET training-data format and return the EMBSET text.

Each OSZICAR contributes one configuration block: the final-step energy, the per-atom magnetic moments, and the per-atom constraining field. The magnetic moments and fields are rotated by the saxis quantization-axis rotation Rz(alpha) * Ry(beta).

Arguments

• oszicars::AbstractVector{<:AbstractString}: paths to the OSZICAR files; each becomes one configuration, numbered in the given order.

Keyword arguments

• saxis::AbstractVector{<:Real} = [0.0, 0.0, 1.0]: quantization axis.
• energy_kind::AbstractString = "f": "f" for the free energy, "e0" for energy(sigma->0).
• mint::Bool = false: when true, read the magnetic moment from the M_int columns; otherwise from MW_int.
• output::Union{AbstractString, Nothing} = nothing: when given, the EMBSET text is also written to this file.

Returns

• String: the full EMBSET text.

Examples

text = oszicar_to_embset(["run1/OSZICAR", "run2/OSZICAR"])
oszicar_to_embset(["OSZICAR"]; saxis = [1.0, 0.0, 0.0], output = "EMBSET")

sample_mfa_incar(incar_path; variable, start, stop, num_points,
                 num_samples=1, randomize=false, fix="", uniform_atoms="",
                 outdir=".", prefix="sample") -> Vector{String}

Sample thermally conditioned spin configurations from a VASP INCAR and write each one to its own INCAR file.

The initial spin matrix is read from MAGMOM (or M_CONSTR if MAGMOM is absent). For every value in an evenly spaced sweep of the control variable, num_samples configurations are drawn with the Mean-Field Approximation sampler (per-atom directions from a von Mises-Fisher distribution, magnitudes preserved). Each output file copies all keys from the input INCAR and sets both MAGMOM and M_CONSTR to the sampled configuration.

Arguments

• incar_path::AbstractString: path to the input INCAR.

Keyword arguments

• variable::AbstractString (required): control variable, "tau" (scaled temperature T/Tc, expected in (0, 1]) or "m" (magnetization, expected in [0, 1)). Values outside the range are clamped to the corresponding ordered/disordered limit.
• start::Real, stop::Real, num_points::Integer (required): the sweep values are range(start, stop; length = num_points).
• num_samples::Integer = 1: configurations drawn per sweep value.
• randomize::Bool = false: apply a Haar-uniform random global rotation (quantization-axis randomization) to each drawn configuration. Uniform over all of SO(3), so the sampled orientations are isotropic regardless of the direction the reference MAGMOM happens to be written along.
• fix::AbstractString = "": 1-based atom indices kept at their input directions (rotated by the same global rotation when randomize), e.g. "1-10,12,20-22".
• uniform_atoms::AbstractString = "": 1-based atom indices whose direction is redrawn uniformly on the sphere instead of from the vMF distribution (same index syntax). These carry no mean-field alignment — fully isotropic for every sweep value, independent of variable (the disordered κ → 0 limit), unlike a default atom (partially aligned) or a fix atom (frozen). Magnitudes are preserved; if an index is also in fix, fix takes precedence.
• outdir::AbstractString = ".": directory for the output files (created if needed).
• prefix::AbstractString = "sample": output file-name prefix. The magesty vasp mfa command does not expose this; it always uses "sample".

Returns

• Vector{String}: the written file paths, in sweep order. Files are named joinpath(outdir, "<prefix>-NN.INCAR") where NN runs from 1 to num_points * num_samples in (point, sample) order, zero-padded to the width of that count.

Examples

# 3 temperatures from 0.1 to 0.3, two samples each -> 6 INCAR files.
sample_mfa_incar("INCAR"; variable="tau", start=0.1, stop=0.3,
                 num_points=3, num_samples=2, outdir="samples")

sce_to_sunny(model::SCEModel; spin, g=2, mode=:auto, scaling=:auto, output=nothing, placement=:auto) -> String

Export a fitted spin-cluster-expansion model to a runnable Sunny.jl script that computes a linear spin-wave-theory (LSWT) magnon dispersion, and return the script text.

The lowest-order SALCs are converted to a conventional spin Hamiltonian: two-site l₁ = l₂ = 1 terms become 3×3 bilinear exchange matrices (Heisenberg, Dzyaloshinskii–Moriya, and anisotropic symmetric parts together), and single-site l = 2 terms become single-ion anisotropy. Higher-order SALCs (higher-l pairs, three-body and beyond) cannot be represented in Sunny and are skipped, with a warning listing what was dropped. The reference energy j0 and spin-independent terms are dropped (Sunny carries no constant energy term); the dispersion is unaffected. Energies are in the unit of the fit (typically eV).

Physical spin

The SCE couplings are fit with unit spin directions, so they absorb the spin magnitude (J_SCE = J_phys·S²). The classical energy is therefore independent of the spin length, but the magnon dispersion scales as ħω ∝ 1/S for a fixed energy landscape. You must pass the physical effective spin S_eff = m/(g μ_B) (the local-moment magnitude); using s = 1 would inflate the dispersion by a factor ~S (for MnTe, S = 5/2 ⇒ ~2.5× too high).

scaling selects how S_eff is encoded, since Sunny's Moment only accepts spin lengths that are exact multiples of 1/2:

• :moment — put S_eff directly into Moment. Each bilinear bond is rescaled by 1/(s_i s_j) and each single-ion term by a mode-dependent factor, so energy(sys) still reproduces predict_energy(model, …) - j0 and the dispersion is physical. Requires S_eff to be a half-integer.
• :coupling — keep Moment at a fixed half-integer placeholder (s₀ = 1) and let the couplings carry S_eff (J = M/(s₀·√(S_i S_j)), single-ion 1/(s₀ S_i)). Accepts any positive real S_eff (itinerant / non-half-integer moments). The dispersion is invariant under an overall spin scale, so this still yields the physical magnon dispersion, but the represented energy landscape is rescaled — energy(sys) is then not the SCE energy. Exact for a uniform S_eff; for a non-uniform S_eff the off-diagonal exchange stays exact while the on-site (Larmor) term is approximate (a warning is emitted).

:auto (default) uses :moment when every magnetic S_eff is a half-integer and :coupling otherwise, so half-integer moments keep the exact (energy-preserving) behavior while itinerant moments still produce a physical dispersion.

Arguments

• model::SCEModel: a fitted model (e.g. from SCEModel(fit) or Magesty.load(SCEModel, path)).

Keyword arguments

• spin::Union{Real, AbstractDict} (required): the effective spin length S_eff = m/(g μ_B). A scalar applies to every magnetic species; a Dict(species => S_eff) sets it per species (every magnetic species must be present). Must be positive. With scaling = :moment it must also be a half-integer; scaling = :coupling (or :auto) accepts any positive real. Omitting it is an error.
• g::Union{Real, AbstractDict} = 2: the g-factor passed to Moment; scalar or per-species Dict. It does not affect the bare dispersion (only an external field or neutron intensities would use it).
• mode::Symbol = :auto: Sunny system mode. :auto selects :dipole when every magnetic spin is a half-integer, otherwise :dipole_uncorrected (the classical limit, appropriate for non-half-integer / itinerant moments). :dipole applies the quantum single-ion renormalization s(2s-1)/2 (undefined for s ≤ 1/2); :dipole_uncorrected uses the classical s². Note: mode = :dipole combined with scaling = :coupling and single-ion anisotropy is rejected — the quantum renormalization cannot ride the placeholder Moment (use :dipole_uncorrected).
• scaling::Symbol = :auto: how the physical S_eff is encoded (see Physical spin). :moment puts S_eff in Moment (energy-preserving, half-integer only); :coupling keeps Moment at a placeholder s₀ = 1 and rescales the couplings (dispersion-only, any positive real S_eff); :auto picks :moment for all-half-integer magnetic spins and :coupling otherwise.
• output::Union{AbstractString, Nothing} = nothing: when given, the script is also written to this file (.jl is appended if absent).
• placement::Symbol = :auto: :primitive maps interactions onto the chemical primitive cell for an unfolded dispersion; :explicit keeps the training supercell (the dispersion is folded into the supercell Brillouin zone) but is exact for any model. :auto chooses :primitive when the model is cleanly unfoldable (interaction range below half the supercell) and :explicit otherwise.

Returns

• String: the full Sunny.jl script.

Throws

• ArgumentError: if spin is omitted; if a spin / g Dict is missing a magnetic species; if any magnetic spin is not positive; if scaling = :moment (or :auto resolving to it) is used with a non-half-integer spin; if mode, scaling, or placement is invalid; if a model with single-ion anisotropy uses mode = :dipole with s ≤ 1/2; or if scaling = :coupling is combined with mode = :dipole and single-ion anisotropy.

Examples

model = Magesty.load(SCEModel, "model.xml")
# MnTe: Mn²⁺ has S = 5/2 (half-integer ⇒ :moment route, energy-preserving).
script = sce_to_sunny(model; spin = 5//2, output = "lswt.jl")
# Itinerant moment, e.g. Fe 2.2 μB ⇒ S_eff = 1.1 (:auto picks the :coupling route).
script = sce_to_sunny(model; spin = 1.1)
# Per-species:
script = sce_to_sunny(model; spin = Dict("Mn" => 5//2, "Fe" => 1.1))

save(obj::SCEBasis, path::AbstractString)
save(obj::SCEModel, path::AbstractString)
save(f::SCEFit, path::AbstractString)

Write obj to an XML file at path. An SCEBasis is written as structure, symmetry parameters, and SALC basis; an SCEModel adds the fitted reference energy and SCE coefficients in a <JPhi> block. An SCEFit is serialized as the corresponding SCEModel(f) — the fit-time dataset and estimator are not persisted.

The path must end in .xml; any other extension is an error. Use load to read the file back.

Arguments

• obj::Union{SCEBasis, SCEModel}: The object to serialize.
• f::SCEFit: Trained fit; saved as SCEModel(f).
• path::AbstractString: Output file path; must end in .xml.

Returns

• nothing. The XML file is written as a side effect.

Throws

• ArgumentError if path does not end in .xml.

Examples

Magesty.save(basis, "basis.xml")
Magesty.save(model, "model.xml")

load(::Type{SCEBasis}, path::AbstractString) -> SCEBasis
load(::Type{SCEModel}, path::AbstractString) -> SCEModel

Read an SCEBasis or SCEModel from the XML file at path. structure is parsed from the file, symmetry is recomputed from structure and the stored tolerance_sym, and the SALC basis is reconstructed from the stored SALC data (the expensive SALC computation is skipped).

load(SCEModel, path) additionally reads the <JPhi> block; load(SCEBasis, path) accepts an SCEModel XML as well and simply ignores the <JPhi> block.

The path must end in .xml; any other extension is an error.

Arguments

• T::Type: SCEBasis or SCEModel.
• path::AbstractString: Input file path; must end in .xml.

Returns

• SCEBasis or SCEModel: the deserialized object, matching the requested type T.

Throws

• ArgumentError if path does not end in .xml, if a required schema attribute is missing, or (for SCEModel) if the <JPhi> block is absent.

Examples

basis = Magesty.load(SCEBasis, "basis.xml")
model = Magesty.load(SCEModel, "model.xml")

AbstractEstimator

Abstract type for SCE coefficient estimation methods. Concrete subtypes are passed to fit(SCEFit, dataset, estimator; torque_weight) and carry the estimator's hyperparameters in their fields; estimators without hyperparameters are zero-field singleton types (OLS), and those with hyperparameters expose them through keyword constructors (Ridge).

OLS()

Ordinary least-squares estimator (no regularization). OLS is a zero-field singleton because the OLS solver has no hyperparameter to tune — there is no lambda-like knob, by design. Construct as OLS(); contrast with Ridge(lambda = ...).

The solve goes through Cholesky on the normal equations: cholesky(Symmetric(X'X)) \ (X'y). Cholesky is the fastest stable factorization for symmetric positive-definite systems and keeps memory to the num_salcs × num_salcs Gram matrix.

If the design matrix is rank-deficient or numerically near-collinear, OLS throws ArgumentError whose message explains the cause and recommends Ridge(lambda = ε). An unregularized fit on such data is physically meaningless (the SCE coefficients are not identifiable), so the right answer is an explicit error rather than a silent fallback.

Examples

est = OLS()
f   = fit(SCEFit, dataset, OLS(); torque_weight = 0.3)

Ridge(; lambda::Real = 0.0)

L2-regularized least-squares (ridge) estimator. The penalty applies uniformly to every SCE coefficient; the bias term j0 does not need to be excluded explicitly because it is eliminated analytically before the solve (see assemble_weighted_problem / extract_j0_jphi).

The solve is cholesky(Symmetric(X'X + lambda * I)) \ (X'y) — the same Cholesky-on-normal-equations route as OLS, with the lambda * I shift guaranteeing strict positive-definiteness for lambda > 0 (so Cholesky cannot fail here). When lambda ≈ 0, the call delegates to the OLS solver and inherits its PosDefException → ArgumentError behavior.

Unlike OLS, Ridge carries one hyperparameter, lambda, and is therefore a regular struct rather than a singleton. Constructing Ridge(lambda = 0.0) is exactly equivalent to OLS numerically — the extra type exists so that estimator sweeps can iterate over OLS() and Ridge(lambda = λ) without special-casing.

Fields

• lambda::Float64: Regularization strength λ ∈ [0, ∞). λ = 0 reduces to OLS; larger values shrink jphi toward zero.

Examples

# Default keyword form (lambda = 0.0 -> equivalent to OLS).
est = Ridge()

# Typical regularized fit.
est = Ridge(lambda = 1e-4)

ElasticNet(; alpha::Real, lambda::Real, standardize::Bool = true)

Elastic-Net estimator backed by GLMNet.jl. Covers Lasso (alpha = 1), GLMNet-style L2 (alpha = 0), and honest Elastic Net (mixed norm). standardize = true divides each column of X by its empirical standard deviation before the solve. The per-cluster (4π)^(N/2) basis normalization already puts the columns on an equal footing in the sphere-averaged (population) sense — that is what makes the fitted coefficients map onto conventional spin-model parameters — so standardization corrects only the residual per-column scale that finite, non-uniform sampling leaves behind. L1 and mixed-norm selection are sensitive to that residual scale (it shifts which clusters enter the active set), so the default is true.

The bias term j0 is eliminated analytically inside assemble_weighted_problem (the energy block of X is mean-centered before the solve) and re-fit afterward by extract_j0_jphi from the un-scaled energy residual. This estimator therefore calls GLMNet with intercept = false; adding GLMNet's own intercept on top of the already-centered system would re-introduce a uniform offset across both energy and torque rows and bias jphi.

Fields

• alpha::Float64: Mixing parameter, 0 ≤ alpha ≤ 1. alpha = 1 is pure Lasso, alpha = 0 is pure L2 (GLMNet's coordinate-descent variant), in between is Elastic Net.
• lambda::Float64: Penalty strength, λ ≥ 0. λ = 0 reduces to OLS up to GLMNet's coordinate-descent precision.
• standardize::Bool: Forwarded to GLMNet, which divides each column by its empirical standard deviation. Default true: it equalizes the residual per-column scale left by finite sampling, to which L1 / mixed-norm selection is sensitive. (The (4π)^(N/2) basis factor already handles the population-level normalization.)

Examples

est = ElasticNet(alpha = 0.5, lambda = 1e-3)

# Lasso(λ = ...) is a convenience function returning ElasticNet(alpha = 1, ...).
est = Lasso(lambda = 1e-3)

Lasso(; lambda::Real, standardize::Bool = true) -> ElasticNet

Convenience function (not a type) returning ElasticNet(alpha = 1.0, lambda = lambda, standardize = standardize).

There is no separate Lasso struct: Lasso(lambda = ...) isa ElasticNet is true, and code that needs to detect the Lasso case should test e isa ElasticNet && e.alpha == 1.0. A function (rather than a const alias) prevents Lasso(alpha = 0.5, ...) from parsing, which would contradict the intended meaning "α = 1 only".

Examples

est = Lasso(lambda = 1e-3)
est isa ElasticNet  # true
est.alpha           # 1.0

AdaptiveLasso(; pilot::AbstractEstimator = OLS(),
                lambda::Real,
                gamma::Real = 1.0,
                epsilon::Real = eps(Float64),
                standardize::Bool = true)

One-shot Adaptive Lasso (Zou 2006). Runs pilot on (X, y) to obtain beta_pilot, then solves the weighted-L1 Lasso

min_b ||y - X * b||^2 / (2 n) + lambda * sum_j w_j * |b_j|

with w_j = 1 / max(|beta_pilot[j]|, epsilon)^gamma. Backed by GLMNet.jl with intercept = false; the energy block of X is already mean-centered upstream by assemble_weighted_problem and j0 is recovered downstream by extract_j0_jphi, the same post-processing OLS, Ridge, and ElasticNet use.

Defaults match the ALAMODE adaptive-LASSO recipe (OLS pilot, gamma = 1). gamma = 0 reduces to plain Lasso, which the test suite exploits as a correctness anchor.

Note: GLMNet internally rescales penalty_factor so the supplied weights sum to nvars. The user-supplied lambda therefore interacts with the rescaled weights, and matched-lambda comparison against plain Lasso is not apples-to-apples once gamma > 0. Per-side lambda tuning is the recommended pattern.

Fields

• pilot::AbstractEstimator: First-stage estimator producing beta_pilot. Default OLS() (Zou 2006 verbatim; matches ALAMODE). For SCE designs that are rank-deficient or near-collinear (e.g. num_salcs >= num_spinconfigs at torque_weight = 0), use pilot = Ridge(lambda = small) instead – the OLS minimum-norm solution populates null-space directions with ~1e-10 noise that is not clipped by eps(Float64) and miscalibrates the adaptive weights.
• lambda::Float64: Final L1 penalty strength, lambda >= 0.
• gamma::Float64: Weight exponent, gamma >= 0. Default 1.0. gamma = 0 reduces to plain Lasso.
• epsilon::Float64: Floor on |beta_pilot[j]| before reciprocation, epsilon > 0. Default eps(Float64). Prevents penalty_factor = Inf when a pilot coefficient is numerically zero (e.g. when pilot is itself a Lasso).
• standardize::Bool: Forwarded to GLMNet. Default true for consistency with ElasticNet / Lasso; on SCE designs this is partially redundant with the adaptive reweighting but does not hurt.

Examples

# Default: OLS pilot, gamma = 1.0 (ALAMODE-style).
est = AdaptiveLasso(lambda = 1e-3)

# Recommended for rank-deficient designs.
est = AdaptiveLasso(pilot = Ridge(lambda = 1e-4), lambda = 1e-3)

# Sanity check: gamma = 0 reduces to plain Lasso.
est = AdaptiveLasso(lambda = 1e-3, gamma = 0.0)

# Reuse a previously fitted SCEFit / SCEModel as the pilot, skipping
# the pilot regression. See `PrecomputedPilot` for the underlying
# adapter and `AdaptiveLasso(::SCEFit; ...)` /
# `AdaptiveLasso(::SCEModel; ...)` for the convenience constructors.
est = AdaptiveLasso(pilot = PrecomputedPilot(coef(prior_fit)), lambda = 1e-3)

PrecomputedPilot(beta::AbstractVector{<:Real})

Estimator adapter that returns a fixed coefficient vector from solve_coefficients, ignoring the supplied (X, y) except for a length check against size(X, 2). Designed as an AdaptiveLasso.pilot choice: lets the adaptive call reuse coefficients from a previous fit instead of running a fresh pilot regression.

The input vector is copied at construction (enforced by the inner constructor), so later mutation of the caller's storage does not leak into PrecomputedPilot.beta.

Fields

• beta::Vector{Float64}: Pilot coefficient vector. Named beta rather than coef to avoid visual collision with the StatsAPI.coef function that Magesty extends.

Examples

# Reuse an existing fit's coefficients as the AdaptiveLasso pilot.
est = AdaptiveLasso(pilot = PrecomputedPilot(coef(fit)), lambda = 1e-3)

AdaptiveRidge(; lambda::Real,
                epsilon::Real = 1e-8,
                max_iter::Integer = 50,
                tol::Real = 1e-6)

Iterative Adaptive Ridge estimator (Frommlet & Nuel 2016). Approximates an L0-penalized fit by repeatedly refitting a per-coefficient weighted ridge problem

min_b  ||y - X b||^2 + lambda * sum_j w_j * b_j^2

and updating the weights w_j = 1 / (b_j^2 + epsilon) between iterations. Iteration zero is a plain ridge solve (uniform weights); each subsequent step rebuilds the weights from the current coefficients. Large coefficients receive a light penalty and small ones a heavy penalty, so iterating drives the small coefficients toward zero – an L0 approximation.

Each weighted ridge subproblem is solved analytically via the closed form b = (X'X + lambda * Diagonal(w)) \ (X'y), the same analytic family as Ridge. Unlike ElasticNet / AdaptiveLasso no GLMNet call is involved, so there is no standardize keyword: the penalty acts directly on the coefficients, where the per-cluster (4π)^(N/2) basis normalization already places them on the conventional spin-model scale, so a single epsilon is a roughly uniform magnitude floor across clusters. Because the reweighting approximates L0 (it selects), residual per-column scale from finite, non-uniform sampling can still influence which coefficients survive; rescale the design upstream if that matters under strongly non-uniform sampling.

The iteration stops when the relative infinity-norm change in the coefficient vector drops below tol, or after max_iter reweighting steps, whichever comes first.

Fields

• lambda::Float64: Ridge penalty strength, lambda >= 0. lambda = 0 reduces to OLS – the penalty term vanishes and the iteration is a no-op.
• epsilon::Float64: Floor added to b_j^2 before reciprocation, epsilon > 0. Default 1e-8. Keeps the weights finite when a coefficient is numerically zero, and sets the scale below which a coefficient is treated as negligible.
• max_iter::Int: Maximum number of reweighting iterations, max_iter >= 1. Default 50.
• tol::Float64: Convergence threshold on the relative infinity-norm coefficient change, tol > 0. Default 1e-6.

Examples

# Default iterative Adaptive Ridge.
est = AdaptiveRidge(lambda = 1e-3)

# Tighter convergence, more iterations.
est = AdaptiveRidge(lambda = 1e-3, tol = 1e-8, max_iter = 200)

read_embset(filepath::AbstractString) -> Vector{SpinConfig}

Read spin configurations from an EMBSET file.

Arguments

• filepath::AbstractString: Path to the EMBSET file

Returns

• Vector{SpinConfig}: Array of spin configurations

Throws

• ErrorException if the file format is invalid
• ArgumentError if the file does not exist

SpinConfig

A single noncollinear spin configuration plus the DFT observables that go with it. Typically produced by read_embset, then handed to SCEDataset, predict_energy, predict_torque, or the evaluation verbs (r2_energy, rmse_torque, …).

Fields

• energy::Float64: DFT total energy of the configuration [eV].
• magmom_size::Vector{Float64}: Magnitude of the magnetic moment on each atom [μB]. Length is `numatoms`.
• spin_directions::Matrix{Float64}: Unit-vector spin directions, laid out as 3 × num_atoms (rows = x, y, z; each column has unit norm).
• local_magfield::Matrix{Float64}: Local constraining magnetic field at each atom [eV / μB], same `3 × numatomslayout asspindirections. The values come from VASP'slambda*MWperpblock, which has dimensions of energy per magnetic moment so thatE = −m · B` is in eV — not the Tesla a conventional magnetic field would carry.
• local_magfield_vertical::Matrix{Float64}: Component of local_magfield perpendicular to the spin direction at each atom [eV / μ_B]. Computed at construction; the parallel component is dropped because classical-spin torque depends only on the perpendicular field.
• torques::Matrix{Float64}: Per-atom torque vectors τᵢ = −mᵢ (eᵢ × Bᵢ) (where mᵢ is the moment magnitude, eᵢ the unit-vector spin direction, Bᵢ the local field), laid out as 3 × num_atoms. These are the observables compared against predict_torque during fitting and evaluation.

Constructors

• SpinConfig(energy, magmom_size, spin_directions, local_magfield; atol_unit_norm = 1e-6) — local_magfield_vertical and torques are computed from the other four fields. atol_unit_norm is the absolute tolerance applied to ‖spin_directions[:, i]‖ - 1 for every atom column.

Throws

• ArgumentError if spin_directions and local_magfield do not have matching 3 × num_atoms shapes.
• ArgumentError if any entry of magmom_size is negative.
• ArgumentError if any column of spin_directions deviates from unit norm by more than atol_unit_norm (NaN columns are also rejected).

Examples

# Two atoms: one along +z, one in-plane along +x.
energy          = -12.34
magmom_size     = [1.5, 1.5]
spin_directions = [0.0 1.0; 0.0 0.0; 1.0 0.0]   # 3 × 2 (columns = atoms)
local_magfield  = [0.0 0.1; 0.2 0.0; 0.0 0.0]   # 3 × 2, eV / μ_B
sc = SpinConfig(energy, magmom_size, spin_directions, local_magfield)
sc.torques    # per-atom τᵢ = −mᵢ (eᵢ × Bᵢ), 3 × 2