pub fn build_spherical_spline_basis(
data: ArrayView2<'_, f64>,
spec: &SphericalSplineBasisSpec,
) -> Result<BasisBuildResult, BasisError> {
if matches!(spec.method, SphereMethod::Harmonic) {
return build_spherical_harmonic_basis(data, spec);
}
validate_lat_lon_matrix(data, "spherical spline", spec.radians)?;
if !(1..=4).contains(&spec.penalty_order) {
crate::bail_invalid_basis!(
"spherical spline penalty_order must be one of 1, 2, 3, 4; got {}",
spec.penalty_order
);
}
let centers = match realized_center_strategy(&spec.center_strategy) {
CenterStrategy::FarthestPoint { num_centers } => {
select_spherical_farthest_point_centers(data, *num_centers, spec.radians)?
}
_ => select_centers_by_strategy(data, &spec.center_strategy)?,
};
validate_lat_lon_matrix(centers.view(), "spherical spline centers", spec.radians)?;
if centers.nrows() < 2 {
return Err(BasisError::InsufficientColumnsForConstraint {
found: centers.nrows(),
});
}
let raw_penalty = spherical_wahba_kernel_matrix_with_kind(
centers.view(),
centers.view(),
spec.penalty_order,
spec.radians,
spec.wahba_kernel,
)?;
let z = match &spec.identifiability {
SphericalSplineIdentifiability::FrozenTransform { transform } => {
if transform.nrows() != centers.nrows() {
crate::bail_dim_basis!(
"frozen spherical identifiability transform mismatch: {} centers but transform has {} rows",
centers.nrows(),
transform.nrows()
);
}
transform.clone()
}
SphericalSplineIdentifiability::CenterSumToZero => {
let weights = sphere_area_weights(centers.view(), spec.radians);
weighted_coefficient_sum_to_zero_transform(weights.view())?
}
};
let penalty = z.t().dot(&raw_penalty).dot(&z);
let gpu_raw_design = try_build_truncated_sphere_design_gpu(
data,
centers.view(),
spec.wahba_kernel,
spec.penalty_order,
spec.radians,
);
let sphere_auto_chunk = if gpu_raw_design.is_some() {
None
} else {
auto_streaming_chunk_size_for_dense(data.nrows(), z.ncols())
};
let design = if let Some(raw_design) = gpu_raw_design {
DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(raw_design.dot(&z)))
} else if let Some(chunk) = sphere_auto_chunk {
log::info!(
"Sphere basis auto-streaming evaluator: n={} p={} chunk_size={}",
data.nrows(),
z.ncols(),
chunk,
);
let op = StreamingSphereEvaluator::new(
Arc::new(data.as_standard_layout().to_owned()),
Arc::new(centers.clone()),
spec.penalty_order,
spec.radians,
spec.wahba_kernel,
Some(Arc::new(z.clone())),
Some(chunk),
)
.map_err(BasisError::InvalidInput)?;
DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(Arc::new(op)))
} else {
let raw_design = spherical_wahba_kernel_matrix_with_kind(
data,
centers.view(),
spec.penalty_order,
spec.radians,
spec.wahba_kernel,
)?;
DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(raw_design.dot(&z)))
};
let (penalty_norm, c_primary) = normalize_penalty(&((&penalty + &penalty.t()) * 0.5));
let mut candidates = vec![PenaltyCandidate {
matrix: penalty_norm,
nullspace_dim_hint: 0,
source: PenaltySource::Primary,
normalization_scale: c_primary,
kronecker_factors: None,
op: None,
}];
if spec.double_penalty {
let ridge = Array2::<f64>::eye(design.ncols());
let (ridge_norm, c_ridge) = normalize_penalty(&ridge);
candidates.push(PenaltyCandidate {
matrix: ridge_norm,
nullspace_dim_hint: 0,
source: PenaltySource::DoublePenaltyNullspace,
normalization_scale: c_ridge,
kronecker_factors: None,
op: None,
});
}
let (penalties, nullspace_dims, penaltyinfo, null_eigenvectors, ops) =
filter_active_penalty_candidates_with_ops(candidates)?;
Ok(BasisBuildResult {
design,
penalties,
nullspace_dims,
penaltyinfo,
metadata: BasisMetadata::Sphere {
centers,
penalty_order: spec.penalty_order,
method: SphereMethod::Wahba,
max_degree: None,
wahba_kernel: spec.wahba_kernel,
constraint_transform: Some(z),
},
kronecker_factored: None,
ops,
null_eigenvectors,
joint_null_rotation: None,
})
}
fn precompute_harmonic_norms(max_degree: usize) -> Vec<f64> {
let l_cap = max_degree + 1;
let sqrt2 = std::f64::consts::SQRT_2;
let mut out = vec![0.0_f64; l_cap * l_cap];
for l in 0..=max_degree {
let mut ratio = 1.0_f64;
let base = ((2 * l + 1) as f64) / (4.0 * std::f64::consts::PI);
out[l * l_cap] = base.sqrt(); for m in 1..=l {
ratio /= ((l - m + 1) * (l + m)) as f64;
out[l * l_cap + m] = sqrt2 * (base * ratio).sqrt();
}
}
out
}
fn fill_real_spherical_harmonics_row(
lat: f64,
lon: f64,
max_degree: usize,
p_buf: &mut [f64],
norms: &[f64],
mut row: ndarray::ArrayViewMut1<'_, f64>,
) {
let l_cap = max_degree + 1;
assert_eq!(p_buf.len(), l_cap * l_cap);
assert_eq!(norms.len(), l_cap * l_cap);
let x = lat.sin();
let somx2 = (1.0 - x * x).max(0.0).sqrt();
for slot in p_buf.iter_mut() {
*slot = 0.0;
}
let idx = |l: usize, m: usize| l * l_cap + m;
p_buf[idx(0, 0)] = 1.0;
for m in 1..=max_degree {
p_buf[idx(m, m)] = -((2 * m - 1) as f64) * somx2 * p_buf[idx(m - 1, m - 1)];
}
for m in 0..max_degree {
p_buf[idx(m + 1, m)] = ((2 * m + 1) as f64) * x * p_buf[idx(m, m)];
}
for m in 0..=max_degree {
for l in (m + 2)..=max_degree {
p_buf[idx(l, m)] = (((2 * l - 1) as f64) * x * p_buf[idx(l - 1, m)]
- ((l + m - 1) as f64) * p_buf[idx(l - 2, m)])
/ ((l - m) as f64);
}
}
let (sin1, cos1) = lon.sin_cos();
let mut sin_buf = [0.0_f64; 33];
let mut cos_buf = [0.0_f64; 33];
sin_buf[0] = 0.0;
cos_buf[0] = 1.0;
if max_degree >= 1 {
sin_buf[1] = sin1;
cos_buf[1] = cos1;
}
let two_cos1 = 2.0 * cos1;
for m in 2..=max_degree {
sin_buf[m] = two_cos1 * sin_buf[m - 1] - sin_buf[m - 2];
cos_buf[m] = two_cos1 * cos_buf[m - 1] - cos_buf[m - 2];
}
let mut col = 0usize;
for l in 1..=max_degree {
for m_pos in (1..=l).rev() {
row[col] = norms[idx(l, m_pos)] * p_buf[idx(l, m_pos)] * sin_buf[m_pos];
col += 1;
}
row[col] = norms[idx(l, 0)] * p_buf[idx(l, 0)];
col += 1;
for m in 1..=l {
row[col] = norms[idx(l, m)] * p_buf[idx(l, m)] * cos_buf[m];
col += 1;
}
}
}
pub fn default_spherical_harmonic_degree(n_rows: usize) -> usize {
let target_cols = ((n_rows as f64) * 0.25).min(50.0).max(3.0);
let mut l = 1usize;
while (l as f64) * (l as f64 + 2.0) < target_cols && l < 12 {
l += 1;
}
l.max(2)
}
fn build_spherical_harmonic_basis(
data: ArrayView2<'_, f64>,
spec: &SphericalSplineBasisSpec,
) -> Result<BasisBuildResult, BasisError> {
validate_lat_lon_matrix(data, "spherical-harmonic", spec.radians)?;
let n = data.nrows();
let l_max = spec
.max_degree
.unwrap_or_else(|| default_spherical_harmonic_degree(n));
if l_max < 1 {
crate::bail_invalid_basis!("spherical-harmonic max_degree must be >= 1");
}
if l_max > 32 {
crate::bail_invalid_basis!("spherical-harmonic max_degree {l_max} too large; cap is 32");
}
if !(1..=4).contains(&spec.penalty_order) {
crate::bail_invalid_basis!(
"spherical-harmonic penalty_order must be one of 1, 2, 3, 4; got {}",
spec.penalty_order
);
}
let p = l_max * (l_max + 2);
let to_rad = if spec.radians {
1.0
} else {
std::f64::consts::PI / 180.0
};
let norms = precompute_harmonic_norms(l_max);
let l_cap = l_max + 1;
let mut design = Array2::<f64>::zeros((n, p));
{
let mut row_blocks = design
.axis_chunks_iter_mut(ndarray::Axis(0), 1024)
.collect::<Vec<_>>();
let chunks_iter = row_blocks.par_iter_mut().enumerate();
let chunk_size = 1024usize;
chunks_iter.try_for_each(|(chunk_idx, block)| -> Result<(), BasisError> {
let mut p_buf = vec![0.0_f64; l_cap * l_cap];
let row_offset = chunk_idx * chunk_size;
for (local_i, mut out_row) in block.outer_iter_mut().enumerate() {
let i = row_offset + local_i;
let lat_raw = data[(i, 0)] * to_rad;
let lat = lat_raw.clamp(-std::f64::consts::FRAC_PI_2, std::f64::consts::FRAC_PI_2);
let lon = data[(i, 1)] * to_rad;
fill_real_spherical_harmonics_row(
lat,
lon,
l_max,
p_buf.as_mut_slice(),
norms.as_slice(),
out_row.view_mut(),
);
}
Ok(())
})?;
}
let mut penalty = Array2::<f64>::zeros((p, p));
let mut col = 0usize;
for l in 1..=l_max {
let eig = (l as f64 * (l as f64 + 1.0)).powi(spec.penalty_order as i32);
for _ in 0..(2 * l + 1) {
penalty[[col, col]] = eig;
col += 1;
}
}
let mut candidates = vec![PenaltyCandidate {
matrix: penalty,
nullspace_dim_hint: 0,
source: PenaltySource::Primary,
normalization_scale: 1.0,
kronecker_factors: None,
op: None,
}];
if spec.double_penalty {
let ridge = Array2::<f64>::eye(p);
let (ridge_norm, c_ridge) = normalize_penalty(&ridge);
candidates.push(PenaltyCandidate {
matrix: ridge_norm,
nullspace_dim_hint: 0,
source: PenaltySource::DoublePenaltyNullspace,
normalization_scale: c_ridge,
kronecker_factors: None,
op: None,
});
}
let (penalties, nullspace_dims, penaltyinfo, null_eigenvectors, ops) =
filter_active_penalty_candidates_with_ops(candidates)?;
Ok(BasisBuildResult {
design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(design)),
penalties,
nullspace_dims,
penaltyinfo,
metadata: BasisMetadata::Sphere {
centers: Array2::<f64>::zeros((0, 2)),
penalty_order: spec.penalty_order,
method: SphereMethod::Harmonic,
max_degree: Some(l_max),
wahba_kernel: spec.wahba_kernel,
constraint_transform: None,
},
kronecker_factored: None,
ops,
null_eigenvectors,
joint_null_rotation: None,
})
}
pub fn build_matern_basis(
data: ArrayView2<'_, f64>,
spec: &MaternBasisSpec,
) -> Result<BasisBuildResult, BasisError> {
let mut workspace = BasisWorkspace::default();
build_matern_basiswithworkspace(data, spec, &mut workspace)
}
pub fn build_matern_basis_literal_aniso(
data: ArrayView2<'_, f64>,
spec: &MaternBasisSpec,
) -> Result<BasisBuildResult, BasisError> {
let mut workspace = BasisWorkspace::default();
build_matern_basis_seeded(data, spec, &mut workspace, AnisoSeedMode::Literal)
}
pub fn build_matern_basiswithworkspace(
data: ArrayView2<'_, f64>,
spec: &MaternBasisSpec,
workspace: &mut BasisWorkspace,
) -> Result<BasisBuildResult, BasisError> {
build_matern_basis_seeded(data, spec, workspace, AnisoSeedMode::AutoSeedFromGeometry)
}
fn build_matern_basis_seeded(
data: ArrayView2<'_, f64>,
spec: &MaternBasisSpec,
workspace: &mut BasisWorkspace,
aniso_seed_mode: AnisoSeedMode,
) -> Result<BasisBuildResult, BasisError> {
let selected_centers = select_centers_by_strategy(data, &spec.center_strategy)?;
let original_centers = if matches!(
spec.identifiability,
MaternIdentifiability::FrozenTransform { .. }
) {
selected_centers
} else {
let reduce_aniso = resolve_matern_forward_aniso(
aniso_seed_mode,
selected_centers.view(),
spec.aniso_log_scales.as_deref(),
);
matern_rank_reduce_centers(
data,
&selected_centers,
spec.length_scale,
spec.nu,
reduce_aniso.as_deref(),
)?
};
let centers = expand_periodic_centers(&original_centers, spec.periodic.as_deref())?;
let aniso = resolve_matern_forward_aniso(
aniso_seed_mode,
centers.view(),
spec.aniso_log_scales.as_deref(),
);
let z_opt = matern_identifiability_transform(centers.view(), &spec.identifiability)?;
let identifiability_transform = z_opt.clone();
let full_transform = z_opt.as_ref().map(|z| {
if spec.include_intercept {
append_intercept_to_transform(z)
} else {
z.clone()
}
});
let frozen_nullspace_shrinkage_survived = match &spec.identifiability {
MaternIdentifiability::FrozenTransform {
nullspace_shrinkage_survived,
..
} => *nullspace_shrinkage_survived,
_ => None,
};
let mut realized_nullspace_shrinkage_survived = false;
let design_cols =
z_opt.as_ref().map_or(centers.nrows(), Array2::ncols) + usize::from(spec.include_intercept);
let dense_bytes = dense_design_bytes(data.nrows(), design_cols);
let matern_auto_chunk = auto_streaming_chunk_size_for_dense(data.nrows(), design_cols);
let use_streaming = matern_auto_chunk.is_some();
let use_lazy = !use_streaming
&& should_use_lazy_spatial_design(data.nrows(), design_cols, workspace.policy());
let (design, candidates) = if let Some(chunk) = matern_auto_chunk {
log::info!(
"Matérn basis auto-streaming evaluator: n={} p={} chunk_size={}",
data.nrows(),
design_cols,
chunk,
);
let shared_data = shared_owned_data_matrix(data, &workspace.cache);
let op = StreamingMaternEvaluator::new(
shared_data,
Arc::new(centers.clone()),
spec.length_scale,
spec.nu,
aniso.clone(),
z_opt.as_ref().map(|z| Arc::new(z.clone())),
spec.include_intercept,
Some(chunk),
)
.map_err(BasisError::InvalidInput)?;
let design = DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(Arc::new(op)));
let candidates = if spec.double_penalty {
let penalty_kernel = build_matern_kernel_penalty(
centers.view(),
spec.length_scale,
spec.nu,
spec.include_intercept,
aniso.as_deref(),
)?;
let primary = project_penalty_matrix(&penalty_kernel, full_transform.as_ref());
let (candidates, survived) = matern_double_penalty_candidates_with_decision(
&primary,
frozen_nullspace_shrinkage_survived,
)?;
realized_nullspace_shrinkage_survived = survived;
candidates
} else {
build_matern_operator_penalty_candidates(
centers.view(),
spec.length_scale,
spec.nu,
spec.include_intercept,
z_opt.as_ref(),
aniso.as_deref(),
)?
};
(design, candidates)
} else if use_lazy {
log::info!(
"Matérn basis switching to lazy chunked design: n={} p={} ({:.1} MiB dense)",
data.nrows(),
design_cols,
dense_bytes as f64 / (1024.0 * 1024.0),
);
let shared_data = shared_owned_data_matrix(data, &workspace.cache);
let d = data.ncols();
let length_scale = spec.length_scale;
let nu = spec.nu;
let poly_basis = if spec.include_intercept {
Some(Arc::new(Array2::<f64>::ones((data.nrows(), 1))))
} else {
None
};
let design = if let Some(eta) = aniso.as_ref() {
let metric_weights = eta.iter().map(|&v| (2.0 * v).exp()).collect::<Vec<_>>();
let kernel = move |data_row: &[f64], center_row: &[f64]| -> f64 {
let mut q = 0.0f64;
for axis in 0..data_row.len() {
let delta = data_row[axis] - center_row[axis];
q += metric_weights[axis] * delta * delta;
}
matern_kernel_from_distance(q.sqrt(), length_scale, nu)
.expect("validated Matérn inputs should not fail")
};
let op = ChunkedKernelDesignOperator::new(
shared_data.clone(),
Arc::new(centers.clone()),
kernel,
z_opt.as_ref().map(|z| Arc::new(z.clone())),
poly_basis.clone(),
)
.map_err(BasisError::InvalidInput)?;
DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(Arc::new(op)))
} else {
let kernel = move |data_row: &[f64], center_row: &[f64]| -> f64 {
let r = stable_euclidean_norm((0..d).map(|axis| data_row[axis] - center_row[axis]));
matern_kernel_from_distance(r, length_scale, nu)
.expect("validated Matérn inputs should not fail")
};
let op = ChunkedKernelDesignOperator::new(
shared_data,
Arc::new(centers.clone()),
kernel,
z_opt.as_ref().map(|z| Arc::new(z.clone())),
poly_basis,
)
.map_err(BasisError::InvalidInput)?;
DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(Arc::new(op)))
};
let candidates = if spec.double_penalty {
let penalty_kernel = build_matern_kernel_penalty(
centers.view(),
spec.length_scale,
spec.nu,
spec.include_intercept,
aniso.as_deref(),
)?;
let primary = project_penalty_matrix(&penalty_kernel, full_transform.as_ref());
let (candidates, survived) = matern_double_penalty_candidates_with_decision(
&primary,
frozen_nullspace_shrinkage_survived,
)?;
realized_nullspace_shrinkage_survived = survived;
candidates
} else {
build_matern_operator_penalty_candidates(
centers.view(),
spec.length_scale,
spec.nu,
spec.include_intercept,
z_opt.as_ref(),
aniso.as_deref(),
)?
};
(design, candidates)
} else {
let m = create_matern_spline_basiswithworkspace(
data,
centers.view(),
spec.length_scale,
spec.nu,
spec.include_intercept,
aniso.as_deref(),
workspace,
)?;
let design = if let Some(transform) = full_transform.as_ref() {
DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(fast_ab(
&m.basis, transform,
)))
} else {
DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(m.basis.clone()))
};
let candidates = if spec.double_penalty {
let (candidates, survived) = build_matern_double_penalty_candidates(
&m,
full_transform.as_ref(),
frozen_nullspace_shrinkage_survived,
)?;
realized_nullspace_shrinkage_survived = survived;
candidates
} else {
build_matern_operator_penalty_candidates(
centers.view(),
spec.length_scale,
spec.nu,
spec.include_intercept,
z_opt.as_ref(),
aniso.as_deref(),
)?
};
(design, candidates)
};
let (penalties, nullspace_dims, penaltyinfo, null_eigenvectors, ops) =
filter_active_penalty_candidates_with_ops(candidates)?;
Ok(BasisBuildResult {
design,
penalties,
nullspace_dims,
penaltyinfo,
metadata: BasisMetadata::Matern {
centers: original_centers,
length_scale: spec.length_scale,
periodic: spec.periodic.clone(),
nu: spec.nu,
include_intercept: spec.include_intercept,
identifiability_transform,
input_scales: None,
aniso_log_scales: aniso,
nullspace_shrinkage_survived: realized_nullspace_shrinkage_survived,
},
kronecker_factored: None,
ops,
null_eigenvectors,
joint_null_rotation: None,
})
}
#[inline(always)]
fn eval_polywith_derivatives(coeffs: &[f64], a: f64) -> (f64, f64, f64) {
let mut p = 0.0;
let mut p1 = 0.0;
let mut p2 = 0.0;
for (i, &c) in coeffs.iter().enumerate() {
p += c * a.powi(i as i32);
if i >= 1 {
p1 += (i as f64) * c * a.powi((i - 1) as i32);
}
if i >= 2 {
p2 += (i as f64) * ((i - 1) as f64) * c * a.powi((i - 2) as i32);
}
}
(p, p1, p2)
}
#[inline(always)]
fn maternvalue_psi_triplet(
r: f64,
length_scale: f64,
nu: MaternNu,
) -> Result<(f64, f64, f64), BasisError> {
validate_matern_inputs(r, length_scale)?;
let kappa = 1.0 / length_scale;
let (s, p): (f64, &[f64]) = match nu {
MaternNu::Half => (kappa, &[1.0]),
MaternNu::ThreeHalves => (3.0_f64.sqrt() * kappa, &[1.0, 1.0]),
MaternNu::FiveHalves => (5.0_f64.sqrt() * kappa, &[1.0, 1.0, 1.0 / 3.0]),
MaternNu::SevenHalves => (7.0_f64.sqrt() * kappa, &[1.0, 1.0, 2.0 / 5.0, 1.0 / 15.0]),
MaternNu::NineHalves => (
9.0_f64.sqrt() * kappa,
&[1.0, 1.0, 3.0 / 7.0, 2.0 / 21.0, 1.0 / 105.0],
),
};
let a = s * r;
if a > 700.0 {
return Ok((0.0, 0.0, 0.0));
}
let e = (-a).exp();
let (p0, p1, p2) = eval_polywith_derivatives(p, a);
let value = e * p0;
let value_psi = e * a * (p1 - p0);
let value_psi_psi = e * (a * (p1 - p0) + a * a * (p2 - 2.0 * p1 + p0));
Ok((value, value_psi, value_psi_psi))
}
#[inline(always)]
fn exp_poly_scaled_s2_psi_triplet(s: f64, a: f64, coeffs: &[f64], scalar: f64) -> (f64, f64, f64) {
if a > 700.0 {
return (0.0, 0.0, 0.0);
}
let e = (-a).exp();
let (p0, p1, p2) = eval_polywith_derivatives(coeffs, a);
let d = p1 - p0;
let y = scalar * s * s * e * p0;
let y_psi = scalar * s * s * e * (2.0 * p0 + a * d);
let y_psi_psi = scalar * s * s * e * (4.0 * p0 + 5.0 * a * d + a * a * (p2 - 2.0 * p1 + p0));
(y, y_psi, y_psi_psi)
}
#[inline(always)]
fn matern_operator_psi_triplet(
r: f64,
length_scale: f64,
nu: MaternNu,
dimension: usize,
) -> Result<
(
f64, // phi
f64, // phi_psi
f64, // phi_psi_psi
f64, // phi_r_over_r
f64, // derivative of phi_r_over_r with respect to psi
f64, // second derivative of phi_r_over_r with respect to psi
f64, // lap
f64, // lap_psi
f64, // lap_psi_psi
),
BasisError,
> {
let (phi, phi_psi, phi_psi_psi) = maternvalue_psi_triplet(r, length_scale, nu)?;
let kappa = 1.0 / length_scale;
let d = dimension as f64;
let (s, q, rr): (f64, &[f64], &[f64]) = match nu {
MaternNu::Half => (kappa, &[1.0], &[1.0]),
MaternNu::ThreeHalves => (3.0_f64.sqrt() * kappa, &[1.0], &[-1.0, 1.0]),
MaternNu::FiveHalves => (
5.0_f64.sqrt() * kappa,
&[1.0 / 3.0, 1.0 / 3.0],
&[-1.0 / 3.0, -1.0 / 3.0, 1.0 / 3.0],
),
MaternNu::SevenHalves => (
7.0_f64.sqrt() * kappa,
&[1.0 / 5.0, 1.0 / 5.0, 1.0 / 15.0],
&[-1.0 / 5.0, -1.0 / 5.0, 0.0, 1.0 / 15.0],
),
MaternNu::NineHalves => (
9.0_f64.sqrt() * kappa,
&[1.0 / 7.0, 1.0 / 7.0, 2.0 / 35.0, 1.0 / 105.0],
&[
-1.0 / 7.0,
-1.0 / 7.0,
-1.0 / 35.0,
2.0 / 105.0,
1.0 / 105.0,
],
),
};
let a = s * r;
let (phi_rr, phi_rr_psi, phi_rr_psi_psi) = exp_poly_scaled_s2_psi_triplet(s, a, rr, 1.0);
let (ratio, ratio_psi, ratio_psi_psi) = if matches!(nu, MaternNu::Half) {
let r_eff = r.max(1e-12);
let e_eff = (-a).exp();
let g = -(s / r_eff) * e_eff;
let g_psi = -(s / r_eff) * e_eff * (1.0 - a);
let g_psi_psi = -(s / r_eff) * e_eff * (1.0 - 3.0 * a + a * a);
(g, g_psi, g_psi_psi)
} else {
exp_poly_scaled_s2_psi_triplet(s, a, q, -1.0)
};
let lap = phi_rr + (d - 1.0) * ratio;
let lap_psi = phi_rr_psi + (d - 1.0) * ratio_psi;
let lap_psi_psi = phi_rr_psi_psi + (d - 1.0) * ratio_psi_psi;
if !phi.is_finite()
|| !phi_psi.is_finite()
|| !phi_psi_psi.is_finite()
|| !ratio.is_finite()
|| !ratio_psi.is_finite()
|| !ratio_psi_psi.is_finite()
|| !lap.is_finite()
|| !lap_psi.is_finite()
|| !lap_psi_psi.is_finite()
{
crate::bail_invalid_basis!(
"non-finite Matérn psi-derivative operator values at r={r}, length_scale={length_scale}, nu={nu:?}"
);
}
Ok((
phi,
phi_psi,
phi_psi_psi,
ratio,
ratio_psi,
ratio_psi_psi,
lap,
lap_psi,
lap_psi_psi,
))
}
fn gram_and_psi_derivatives_from_operator(
d: &Array2<f64>,
d_psi: &Array2<f64>,
d_psi_psi: &Array2<f64>,
) -> (Array2<f64>, Array2<f64>, Array2<f64>) {
let s_raw = symmetrize(&fast_ata(d));
let s_raw_psi = symmetrize(&(d_psi.t().dot(d) + d.t().dot(d_psi)));
let s_raw_psi_psi =
symmetrize(&(d_psi_psi.t().dot(d) + d.t().dot(d_psi_psi) + 2.0 * d_psi.t().dot(d_psi)));
(s_raw, s_raw_psi, s_raw_psi_psi)
}
fn gram_cross_psi_derivative_from_operator(
d: &Array2<f64>,
d_a: &Array2<f64>,
d_b: &Array2<f64>,
d_ab: &Array2<f64>,
) -> Array2<f64> {
symmetrize(&(d_ab.t().dot(d) + d.t().dot(d_ab) + d_a.t().dot(d_b) + d_b.t().dot(d_a)))
}
fn normalize_penalty_cross_psi_derivative(
s: &Array2<f64>,
s_a: &Array2<f64>,
s_b: &Array2<f64>,
s_ab: &Array2<f64>,
c: f64,
) -> Array2<f64> {
if !c.is_finite() || c <= 1e-12 {
return Array2::<f64>::zeros(s.raw_dim());
}
let c2 = c * c;
let c3 = c2 * c;
let a_val = trace_of_product(s, s_a);
let c_a = a_val / c;
let b_val = trace_of_product(s, s_b);
let c_b = b_val / c;
let cross_val = trace_of_product(s_a, s_b) + trace_of_product(s, s_ab);
let c_ab = cross_val / c - c_a * c_b / c;
let coeff_s = 2.0 * c_a * c_b / c3 - c_ab / c2;
s_ab.mapv(|v| v / c) - s_b.mapv(|v| c_a / c2 * v) - s_a.mapv(|v| c_b / c2 * v)
+ s.mapv(|v| coeff_s * v)
}
#[inline(always)]
fn trace_of_product(a: &Array2<f64>, b: &Array2<f64>) -> f64 {
a.t().dot(b).diag().sum()
}
fn normalize_penaltywith_psi_derivatives(
s: &Array2<f64>,
s_psi: &Array2<f64>,
s_psi_psi: &Array2<f64>,
) -> (Array2<f64>, Array2<f64>, Array2<f64>, f64) {
let fro2 = trace_of_product(s, s);
let c = fro2.sqrt();
if !c.is_finite() || c <= 1e-12 {
return (
s.clone(),
Array2::<f64>::zeros(s.raw_dim()),
Array2::<f64>::zeros(s.raw_dim()),
1.0,
);
}
let a = trace_of_product(s, s_psi);
let b = trace_of_product(s_psi, s_psi) + trace_of_product(s, s_psi_psi);
let c_psi = a / c;
let c_psi_psi = b / c - (a * a) / (c * c * c);
let s_tilde = s.mapv(|v| v / c);
let s_tilde_psi = s_psi.mapv(|v| v / c) - s.mapv(|v| (c_psi / (c * c)) * v);
let s_tilde_psi_psi = s_psi_psi.mapv(|v| v / c) - s_psi.mapv(|v| 2.0 * c_psi / (c * c) * v)
+ s.mapv(|v| ((2.0 * c_psi * c_psi) / (c * c * c) - c_psi_psi / (c * c)) * v);
(s_tilde, s_tilde_psi, s_tilde_psi_psi, c)
}
fn build_matern_operator_penalty_psi_derivatives(
centers: ArrayView2<'_, f64>,
length_scale: f64,
nu: MaternNu,
include_intercept: bool,
z_opt: Option<&Array2<f64>>,
aniso_log_scales: Option<&[f64]>,
) -> Result<(Vec<Array2<f64>>, Vec<Array2<f64>>), BasisError> {
let p = centers.nrows();
let d = centers.ncols();
let mut d0_raw = Array2::<f64>::zeros((p, p));
let mut d1_raw = Array2::<f64>::zeros((p * d, p));
let mut d2_raw = Array2::<f64>::zeros((p * d * d, p));
let mut d0_raw_psi = Array2::<f64>::zeros((p, p));
let mut d1_raw_psi = Array2::<f64>::zeros((p * d, p));
let mut d2_raw_psi = Array2::<f64>::zeros((p * d * d, p));
let mut d0_raw_psi_psi = Array2::<f64>::zeros((p, p));
let mut d1_raw_psi_psi = Array2::<f64>::zeros((p * d, p));
let mut d2_raw_psi_psi = Array2::<f64>::zeros((p * d * d, p));
let metric_weights = aniso_log_scales
.map(centered_aniso_metric_weights)
.unwrap_or_else(|| vec![1.0; d]);
for k in 0..p {
for j in 0..p {
let (r, _s_vec) = if let Some(eta) = aniso_log_scales {
aniso_distance_and_components(
centers.row(k).as_slice().unwrap(),
centers.row(j).as_slice().unwrap(),
eta,
)
} else {
(
stable_euclidean_norm((0..d).map(|c| centers[[k, c]] - centers[[j, c]])),
(0..d)
.map(|c| {
let h = centers[[k, c]] - centers[[j, c]];
h * h
})
.collect(),
)
};
let (
phi,
phi_psi,
phi_psi_psi,
ratio,
ratio_psi,
ratio_psi_psi,
_lap,
_lap_psi,
_lap_psi_psi,
) = matern_operator_psi_triplet(r, length_scale, nu, d)?;
let (_, _q_shape, t, t_r, t_rr) =
matern_aniso_extended_radial_scalars(r, length_scale, nu)?;
let q = ratio;
let q_psi = ratio_psi;
let q_psi_psi = ratio_psi_psi;
let t_psi = 4.0 * t + r * t_r;
let t_psi_psi = 16.0 * t + 9.0 * r * t_r + r * r * t_rr;
d0_raw[[k, j]] = phi;
d0_raw_psi[[k, j]] = phi_psi;
d0_raw_psi_psi[[k, j]] = phi_psi_psi;
for axis in 0..d {
let delta = centers[[k, axis]] - centers[[j, axis]];
let axis_scale = metric_weights[axis];
let row = k * d + axis;
d1_raw[[row, j]] = ratio * axis_scale * delta;
d1_raw_psi[[row, j]] = ratio_psi * axis_scale * delta;
d1_raw_psi_psi[[row, j]] = ratio_psi_psi * axis_scale * delta;
}
for b in 0..d {
let h_b = centers[[k, b]] - centers[[j, b]];
let w_b = metric_weights[b];
for c in 0..d {
let h_c = centers[[k, c]] - centers[[j, c]];
let w_c = metric_weights[c];
let row = (k * d + b) * d + c;
d2_raw[[row, j]] = hessian_operator_entry(q, t, h_b, h_c, w_b, w_c, b, c);
d2_raw_psi[[row, j]] =
hessian_operator_entry(q_psi, t_psi, h_b, h_c, w_b, w_c, b, c);
d2_raw_psi_psi[[row, j]] =
hessian_operator_entry(q_psi_psi, t_psi_psi, h_b, h_c, w_b, w_c, b, c);
}
}
}
}
let project = |mat: Array2<f64>| {
if let Some(z) = z_opt {
fast_ab(&mat, z)
} else {
mat
}
};
let d0_kernel = project(d0_raw);
let d0_kernel_psi = project(d0_raw_psi);
let d0_kernel_psi_psi = project(d0_raw_psi_psi);
let d1_kernel = project(d1_raw);
let d1_kernel_psi = project(d1_raw_psi);
let d1_kernel_psi_psi = project(d1_raw_psi_psi);
let d2_kernel = project(d2_raw);
let d2_kernel_psi = project(d2_raw_psi);
let d2_kernel_psi_psi = project(d2_raw_psi_psi);
let kernel_cols = d0_kernel.ncols();
let total_cols = kernel_cols + usize::from(include_intercept);
let mut d0 = Array2::<f64>::zeros((p, total_cols));
let mut d1 = Array2::<f64>::zeros((p * d, total_cols));
let mut d2 = Array2::<f64>::zeros((p * d * d, total_cols));
let mut d0_psi = Array2::<f64>::zeros((p, total_cols));
let mut d1_psi = Array2::<f64>::zeros((p * d, total_cols));
let mut d2_psi = Array2::<f64>::zeros((p * d * d, total_cols));
let mut d0_psi_psi = Array2::<f64>::zeros((p, total_cols));
let mut d1_psi_psi = Array2::<f64>::zeros((p * d, total_cols));
let mut d2_psi_psi = Array2::<f64>::zeros((p * d * d, total_cols));
d0.slice_mut(s![.., 0..kernel_cols]).assign(&d0_kernel);
d1.slice_mut(s![.., 0..kernel_cols]).assign(&d1_kernel);
d2.slice_mut(s![.., 0..kernel_cols]).assign(&d2_kernel);
d0_psi
.slice_mut(s![.., 0..kernel_cols])
.assign(&d0_kernel_psi);
d1_psi
.slice_mut(s![.., 0..kernel_cols])
.assign(&d1_kernel_psi);
d2_psi
.slice_mut(s![.., 0..kernel_cols])
.assign(&d2_kernel_psi);
d0_psi_psi
.slice_mut(s![.., 0..kernel_cols])
.assign(&d0_kernel_psi_psi);
d1_psi_psi
.slice_mut(s![.., 0..kernel_cols])
.assign(&d1_kernel_psi_psi);
d2_psi_psi
.slice_mut(s![.., 0..kernel_cols])
.assign(&d2_kernel_psi_psi);
if include_intercept {
d0.column_mut(kernel_cols).fill(1.0);
}
let (s0, s0_psi, s0_psi_psi) =
gram_and_psi_derivatives_from_operator(&d0, &d0_psi, &d0_psi_psi);
let (s1, s1_psi, s1_psi_psi) =
gram_and_psi_derivatives_from_operator(&d1, &d1_psi, &d1_psi_psi);
let (s2, s2_psi, s2_psi_psi) =
gram_and_psi_derivatives_from_operator(&d2, &d2_psi, &d2_psi_psi);
let (s0_norm, s0_norm_psi, s0_norm_psi_psi, c0) =
normalize_penaltywith_psi_derivatives(&s0, &s0_psi, &s0_psi_psi);
let (s1_norm, s1_norm_psi, s1_norm_psi_psi, c1) =
normalize_penaltywith_psi_derivatives(&s1, &s1_psi, &s1_psi_psi);
let (s2_norm, s2_norm_psi, s2_norm_psi_psi, c2) =
normalize_penaltywith_psi_derivatives(&s2, &s2_psi, &s2_psi_psi);
let matern_spec = DuchonOperatorPenaltySpec::matern_for_smoothness(nu, d);
let mut candidates = Vec::with_capacity(3);
for (spec_gate, source, matrix, normalization_scale) in [
(&matern_spec.mass, PenaltySource::OperatorMass, s0_norm, c0),
(
&matern_spec.tension,
PenaltySource::OperatorTension,
s1_norm,
c1,
),
(
&matern_spec.stiffness,
PenaltySource::OperatorStiffness,
s2_norm,
c2,
),
] {
if !matches!(spec_gate, OperatorPenaltySpec::Active { .. }) {
continue;
}
candidates.push(PenaltyCandidate {
matrix,
nullspace_dim_hint: 0,
source,
normalization_scale,
kronecker_factors: None,
op: None,
});
}
let (_, _, penaltyinfo) = filter_active_penalty_candidates(candidates)?;
let penalties_derivative = active_operator_penalty_derivatives(
&penaltyinfo,
&[s0_norm_psi, s1_norm_psi, s2_norm_psi],
"Matérn",
)?;
let penaltiessecond_derivative = active_operator_penalty_derivatives(
&penaltyinfo,
&[s0_norm_psi_psi, s1_norm_psi_psi, s2_norm_psi_psi],
"Matérn",
)?;
Ok((penalties_derivative, penaltiessecond_derivative))
}
struct DuchonRawPenaltyPsiDerivativeBlocks {
d0: Array2<f64>,
d1: Array2<f64>,
d2: Array2<f64>,
d0_psi: Array2<f64>,
d1_psi: Array2<f64>,
d2_psi: Array2<f64>,
d0_psi_psi: Array2<f64>,
d1_psi_psi: Array2<f64>,
d2_psi_psi: Array2<f64>,
}
impl DuchonRawPenaltyPsiDerivativeBlocks {
fn zeros(p: usize, d: usize, cols: usize) -> Self {
Self {
d0: Array2::<f64>::zeros((p, cols)),
d1: Array2::<f64>::zeros((p * d, cols)),
d2: Array2::<f64>::zeros((p * d * d, cols)),
d0_psi: Array2::<f64>::zeros((p, cols)),
d1_psi: Array2::<f64>::zeros((p * d, cols)),
d2_psi: Array2::<f64>::zeros((p * d * d, cols)),
d0_psi_psi: Array2::<f64>::zeros((p, cols)),
d1_psi_psi: Array2::<f64>::zeros((p * d, cols)),
d2_psi_psi: Array2::<f64>::zeros((p * d * d, cols)),
}
}
fn add_assign(&mut self, rhs: &Self) {
self.d0 += &rhs.d0;
self.d1 += &rhs.d1;
self.d2 += &rhs.d2;
self.d0_psi += &rhs.d0_psi;
self.d1_psi += &rhs.d1_psi;
self.d2_psi += &rhs.d2_psi;
self.d0_psi_psi += &rhs.d0_psi_psi;
self.d1_psi_psi += &rhs.d1_psi_psi;
self.d2_psi_psi += &rhs.d2_psi_psi;
}
}
fn build_duchon_operator_penalty_psi_derivatives(
collocation_points: ArrayView2<'_, f64>,
centers: ArrayView2<'_, f64>,
spec: &DuchonBasisSpec,
identifiability_transform: Option<&Array2<f64>>,
workspace: &mut BasisWorkspace,
) -> Result<(Vec<PenaltySource>, Vec<Array2<f64>>, Vec<Array2<f64>>), BasisError> {
let length_scale = spec.length_scale.ok_or_else(|| {
BasisError::InvalidInput(
"exact Duchon log-kappa derivatives require hybrid Duchon with length_scale"
.to_string(),
)
})?;
let effective_nullspace_order = duchon_effective_nullspace_order(centers, spec.nullspace_order);
let p_order = duchon_p_from_nullspace_order(effective_nullspace_order);
let s_order = spec.power_as_usize();
let dim = centers.ncols();
let two_pps = 2.0 * (p_order as f64 + spec.power);
let mut effective_operator_penalties = spec.operator_penalties.clone();
if two_pps <= dim as f64 + 1.0 {
effective_operator_penalties.tension = OperatorPenaltySpec::Disabled;
}
if two_pps <= dim as f64 + 2.0 {
effective_operator_penalties.stiffness = OperatorPenaltySpec::Disabled;
}
let max_derivative_order =
duchon_max_active_operator_derivative_order(&effective_operator_penalties);
if max_derivative_order == 0
&& !matches!(
effective_operator_penalties.mass,
OperatorPenaltySpec::Active { .. }
)
{
return Ok((Vec::new(), Vec::new(), Vec::new()));
}
validate_duchon_collocation_orders(
Some(length_scale),
p_order,
s_order as f64,
dim,
max_derivative_order,
)?;
let coeffs = duchon_partial_fraction_coeffs(p_order, s_order, 1.0 / length_scale);
let z_kernel =
kernel_constraint_nullspace(centers, effective_nullspace_order, &mut workspace.cache)?;
let n_basis = centers.nrows();
if collocation_points.ncols() != dim {
crate::bail_dim_basis!(
"Duchon psi-derivative collocation dim {} != centers dim {dim}",
collocation_points.ncols()
);
}
let p_colloc = collocation_points.nrows();
let d = dim;
let kernel_cols = z_kernel.ncols();
let aniso = spec.aniso_log_scales.as_deref();
if let Some(eta) = aniso
&& eta.len() != d
{
crate::bail_dim_basis!(
"Duchon anisotropy dimension mismatch: got {}, expected {d}",
eta.len()
);
}
let metric_weights: Option<Vec<f64>> = aniso.map(centered_aniso_metric_weights);
let need_d1 = max_derivative_order >= 1;
let need_d2 = max_derivative_order >= 2;
let chunk_count = rayon::current_num_threads().max(1);
let chunk_size = p_colloc.div_ceil(chunk_count).max(1);
let chunks: Vec<(usize, usize)> = (0..p_colloc)
.step_by(chunk_size)
.map(|start| (start, (start + chunk_size).min(p_colloc)))
.collect();
let partial_blocks = chunks
.into_par_iter()
.map(
|(start, end)| -> Result<DuchonRawPenaltyPsiDerivativeBlocks, BasisError> {
let mut local =
DuchonRawPenaltyPsiDerivativeBlocks::zeros(p_colloc, d, kernel_cols);
for i in start..end {
for j in 0..n_basis {
let r = if let Some(eta) = aniso {
let row_i: Vec<f64> =
(0..d).map(|a| collocation_points[[i, a]]).collect();
let row_j: Vec<f64> = (0..d).map(|a| centers[[j, a]]).collect();
let (r, _) = aniso_distance_and_components(&row_i, &row_j, eta);
r
} else {
stable_euclidean_norm(
(0..d)
.map(|axis| collocation_points[[i, axis]] - centers[[j, axis]]),
)
};
let core = duchon_radial_core_psi_triplet(
r,
length_scale,
p_order,
s_order,
d,
&coeffs,
)?;
for col in 0..kernel_cols {
let z_jc = z_kernel[[j, col]];
local.d0[[i, col]] += core.phi.value * z_jc;
local.d0_psi[[i, col]] += core.phi.psi * z_jc;
local.d0_psi_psi[[i, col]] += core.phi.psi_psi * z_jc;
}
if !need_d1 && !need_d2 {
continue;
}
if r > 1e-10 {
let jets =
duchon_radial_jets(r, length_scale, p_order, s_order, d, &coeffs)?;
let q = jets.q;
let (q_psi, q_psi_psi) =
duchon_q_psi_triplet_from_jets(&jets, p_order, s_order, d, r);
let t_exponent = duchon_scaling_exponent(p_order, s_order, d) + 4.0;
let (t_psi, t_psi_psi) = scaled_log_kappa_derivatives(
jets.t, jets.t_r, jets.t_rr, t_exponent, r,
);
if need_d1 {
for axis in 0..d {
let delta = collocation_points[[i, axis]] - centers[[j, axis]];
let axis_scale = metric_weights
.as_ref()
.map(|weights| weights[axis])
.unwrap_or(1.0);
let row = i * d + axis;
for col in 0..kernel_cols {
let z_jc = z_kernel[[j, col]];
local.d1[[row, col]] += q * axis_scale * delta * z_jc;
local.d1_psi[[row, col]] +=
q_psi * axis_scale * delta * z_jc;
local.d1_psi_psi[[row, col]] +=
q_psi_psi * axis_scale * delta * z_jc;
}
}
}
if need_d2 {
for col in 0..kernel_cols {
let z_jc = z_kernel[[j, col]];
for axis_b in 0..d {
let h_b =
collocation_points[[i, axis_b]] - centers[[j, axis_b]];
let w_b = metric_weights
.as_ref()
.map(|weights| weights[axis_b])
.unwrap_or(1.0);
for axis_c in 0..d {
let h_c = collocation_points[[i, axis_c]]
- centers[[j, axis_c]];
let w_c = metric_weights
.as_ref()
.map(|weights| weights[axis_c])
.unwrap_or(1.0);
let row = (i * d + axis_b) * d + axis_c;
local.d2[[row, col]] += hessian_operator_entry(
q, jets.t, h_b, h_c, w_b, w_c, axis_b, axis_c,
) * z_jc;
local.d2_psi[[row, col]] += hessian_operator_entry(
q_psi, t_psi, h_b, h_c, w_b, w_c, axis_b, axis_c,
) * z_jc;
local.d2_psi_psi[[row, col]] += hessian_operator_entry(
q_psi_psi, t_psi_psi, h_b, h_c, w_b, w_c, axis_b,
axis_c,
) * z_jc;
}
}
}
}
} else if need_d2 {
let (phi_rr, phi_rr_psi, phi_rr_psi_psi) =
duchonphi_rr_collision_psi_triplet(
length_scale,
p_order,
s_order,
d,
&coeffs,
)?;
for col in 0..kernel_cols {
let z_jc = z_kernel[[j, col]];
for axis in 0..d {
let w_axis = metric_weights
.as_ref()
.map(|weights| weights[axis])
.unwrap_or(1.0);
let row = (i * d + axis) * d + axis;
local.d2[[row, col]] += w_axis * phi_rr * z_jc;
local.d2_psi[[row, col]] += w_axis * phi_rr_psi * z_jc;
local.d2_psi_psi[[row, col]] += w_axis * phi_rr_psi_psi * z_jc;
}
}
}
}
}
Ok(local)
},
)
.collect::<Result<Vec<_>, BasisError>>()?;
let mut raw = DuchonRawPenaltyPsiDerivativeBlocks::zeros(p_colloc, d, kernel_cols);
for partial in &partial_blocks {
raw.add_assign(partial);
}
let d0_raw = raw.d0;
let d1_raw = raw.d1;
let d2_raw = raw.d2;
let d0_raw_psi = raw.d0_psi;
let d1_raw_psi = raw.d1_psi;
let d2_raw_psi = raw.d2_psi;
let d0_raw_psi_psi = raw.d0_psi_psi;
let d1_raw_psi_psi = raw.d1_psi_psi;
let d2_raw_psi_psi = raw.d2_psi_psi;
let poly = polynomial_block_from_order(centers, effective_nullspace_order);
let kernel_cols = d0_raw.ncols();
let poly_cols = poly.ncols();
let total_cols = kernel_cols + poly_cols;
let mut d0 = Array2::<f64>::zeros((p_colloc, total_cols));
let mut d1 = Array2::<f64>::zeros((p_colloc * d, total_cols));
let mut d2 = Array2::<f64>::zeros((p_colloc * d * d, total_cols));
let mut d0_psi = Array2::<f64>::zeros((p_colloc, total_cols));
let mut d1_psi = Array2::<f64>::zeros((p_colloc * d, total_cols));
let mut d2_psi = Array2::<f64>::zeros((p_colloc * d * d, total_cols));
let mut d0_psi_psi = Array2::<f64>::zeros((p_colloc, total_cols));
let mut d1_psi_psi = Array2::<f64>::zeros((p_colloc * d, total_cols));
let mut d2_psi_psi = Array2::<f64>::zeros((p_colloc * d * d, total_cols));
d0.slice_mut(s![.., 0..kernel_cols]).assign(&d0_raw);
d1.slice_mut(s![.., 0..kernel_cols]).assign(&d1_raw);
d2.slice_mut(s![.., 0..kernel_cols]).assign(&d2_raw);
d0_psi.slice_mut(s![.., 0..kernel_cols]).assign(&d0_raw_psi);
d1_psi.slice_mut(s![.., 0..kernel_cols]).assign(&d1_raw_psi);
d2_psi.slice_mut(s![.., 0..kernel_cols]).assign(&d2_raw_psi);
d0_psi_psi
.slice_mut(s![.., 0..kernel_cols])
.assign(&d0_raw_psi_psi);
d1_psi_psi
.slice_mut(s![.., 0..kernel_cols])
.assign(&d1_raw_psi_psi);
d2_psi_psi
.slice_mut(s![.., 0..kernel_cols])
.assign(&d2_raw_psi_psi);
let project = |mat: Array2<f64>| {
if let Some(z) = identifiability_transform {
fast_ab(&mat, z)
} else {
mat
}
};
let d0 = project(d0);
let d1 = project(d1);
let d2 = project(d2);
let d0_psi = project(d0_psi);
let d1_psi = project(d1_psi);
let d2_psi = project(d2_psi);
let d0_psi_psi = project(d0_psi_psi);
let d1_psi_psi = project(d1_psi_psi);
let d2_psi_psi = project(d2_psi_psi);
let (s0, s0_psi, s0_psi_psi) =
centered_operator_gram_and_psi_derivatives(&d0, &d0_psi, &d0_psi_psi);
let (mut s1, mut s1_psi, mut s1_psi_psi) =
gram_and_psi_derivatives_from_operator(&d1, &d1_psi, &d1_psi_psi);
let (mut s2, mut s2_psi, mut s2_psi_psi) =
gram_and_psi_derivatives_from_operator(&d2, &d2_psi, &d2_psi_psi);
let kappa = 1.0 / length_scale.max(1e-300);
let aniso = spec.aniso_log_scales.as_deref();
if duchon_closed_form_operator_penalty_converges(1, p_order, s_order as f64, d) {
let (cf_s, cf_s_psi, cf_s_psi_psi) = closed_form_psi_derivatives_in_total_basis(
centers,
1,
p_order,
s_order,
kappa,
aniso,
Some(&z_kernel),
poly_cols,
identifiability_transform,
);
s1 = cf_s;
s1_psi = cf_s_psi;
s1_psi_psi = cf_s_psi_psi;
}
if duchon_closed_form_operator_penalty_converges(2, p_order, s_order as f64, d) {
let (cf_s, cf_s_psi, cf_s_psi_psi) = closed_form_psi_derivatives_in_total_basis(
centers,
2,
p_order,
s_order,
kappa,
aniso,
Some(&z_kernel),
poly_cols,
identifiability_transform,
);
s2 = cf_s;
s2_psi = cf_s_psi;
s2_psi_psi = cf_s_psi_psi;
}
let (s0_norm, s0_norm_psi, s0_norm_psi_psi, c0) =
normalize_penaltywith_psi_derivatives(&s0, &s0_psi, &s0_psi_psi);
let (s1_norm, s1_norm_psi, s1_norm_psi_psi, c1) =
normalize_penaltywith_psi_derivatives(&s1, &s1_psi, &s1_psi_psi);
let (s2_norm, s2_norm_psi, s2_norm_psi_psi, c2) =
normalize_penaltywith_psi_derivatives(&s2, &s2_psi, &s2_psi_psi);
let candidates = vec![
PenaltyCandidate {
matrix: s0_norm,
nullspace_dim_hint: 0,
source: PenaltySource::OperatorMass,
normalization_scale: c0,
kronecker_factors: None,
op: None,
},
PenaltyCandidate {
matrix: s1_norm,
nullspace_dim_hint: 0,
source: PenaltySource::OperatorTension,
normalization_scale: c1,
kronecker_factors: None,
op: None,
},
PenaltyCandidate {
matrix: s2_norm,
nullspace_dim_hint: 0,
source: PenaltySource::OperatorStiffness,
normalization_scale: c2,
kronecker_factors: None,
op: None,
},
];
let candidates = operator_penalty_candidates_from_derivative_candidates(
candidates,
&effective_operator_penalties,
);
let first_derivs = vec![s0_norm_psi, s1_norm_psi, s2_norm_psi];
let second_derivs = vec![s0_norm_psi_psi, s1_norm_psi_psi, s2_norm_psi_psi];
let (_, _, penaltyinfo) = filter_active_penalty_candidates(candidates)?;
let active_sources = penaltyinfo
.iter()
.filter(|info| info.active)
.map(|info| info.source.clone())
.collect::<Vec<_>>();
let penalties_derivative =
active_operator_penalty_derivatives(&penaltyinfo, &first_derivs, "Duchon")?;
let penaltiessecond_derivative =
active_operator_penalty_derivatives(&penaltyinfo, &second_derivs, "Duchon")?;
Ok((
active_sources,
penalties_derivative,
penaltiessecond_derivative,
))
}
fn operator_penalty_candidates_from_derivative_candidates(
candidates: Vec<PenaltyCandidate>,
spec: &DuchonOperatorPenaltySpec,
) -> Vec<PenaltyCandidate> {
candidates
.into_iter()
.filter(|candidate| match candidate.source {
PenaltySource::OperatorMass => matches!(spec.mass, OperatorPenaltySpec::Active { .. }),
PenaltySource::OperatorTension => {
matches!(spec.tension, OperatorPenaltySpec::Active { .. })
}
PenaltySource::OperatorStiffness => {
matches!(spec.stiffness, OperatorPenaltySpec::Active { .. })
}
_ => true,
})
.collect()
}
fn build_duchon_native_penalty_psi_derivatives(
centers: ArrayView2<'_, f64>,
spec: &DuchonBasisSpec,
identifiability_transform: Option<&Array2<f64>>,
workspace: &mut BasisWorkspace,
) -> Result<(Vec<PenaltySource>, Vec<Array2<f64>>, Vec<Array2<f64>>), BasisError> {
let length_scale = spec.length_scale.ok_or_else(|| {
BasisError::InvalidInput(
"exact Duchon native penalty log-kappa derivatives require hybrid Duchon with length_scale"
.to_string(),
)
})?;
let effective_nullspace_order = duchon_effective_nullspace_order(centers, spec.nullspace_order);
let p_order = duchon_p_from_nullspace_order(effective_nullspace_order);
let s_order = spec.power_as_usize();
let dim = centers.ncols();
let z = kernel_constraint_nullspace(centers, effective_nullspace_order, &mut workspace.cache)?;
let kernel_cols = z.ncols();
let poly_cols = polynomial_block_from_order(centers, effective_nullspace_order).ncols();
let total_cols = kernel_cols + poly_cols;
let coeffs = duchon_partial_fraction_coeffs(p_order, s_order, 1.0 / length_scale.max(1e-300));
let kernel_amp = duchon_kernel_amplification(
centers,
Some(length_scale),
p_order,
s_order,
dim,
spec.aniso_log_scales.as_deref(),
Some(&coeffs),
None,
);
let axis_scales = spec.aniso_log_scales.as_deref().map(aniso_axis_scales);
let n_centers = centers.nrows();
let mut kernel = Array2::<f64>::zeros((n_centers, n_centers));
let mut kernel_psi = Array2::<f64>::zeros((n_centers, n_centers));
let mut kernel_psi_psi = Array2::<f64>::zeros((n_centers, n_centers));
for i in 0..n_centers {
for j in i..n_centers {
let r = if let Some(scales) = axis_scales.as_deref() {
aniso_distance_rows_with_scales(centers, i, centers, j, scales)
} else {
euclidean_distance_rows(centers, i, centers, j)
};
let core =
duchon_radial_core_psi_triplet(r, length_scale, p_order, s_order, dim, &coeffs)?;
kernel[[i, j]] = core.phi.value;
kernel[[j, i]] = core.phi.value;
kernel_psi[[i, j]] = core.phi.psi;
kernel_psi[[j, i]] = core.phi.psi;
kernel_psi_psi[[i, j]] = core.phi.psi_psi;
kernel_psi_psi[[j, i]] = core.phi.psi_psi;
}
}
let amp2 = kernel_amp * kernel_amp;
let project_kernel = |k: &Array2<f64>| fast_ab(&fast_atb(&z, k), &z).mapv(|v| v * amp2);
let omega = project_kernel(&kernel);
let omega_psi = project_kernel(&kernel_psi);
let omega_psi_psi = project_kernel(&kernel_psi_psi);
let embed = |block: Array2<f64>| {
let mut out = Array2::<f64>::zeros((total_cols, total_cols));
out.slice_mut(s![..kernel_cols, ..kernel_cols])
.assign(&block);
symmetrize(&project_penalty_matrix(&out, identifiability_transform))
};
let primary = embed(omega);
let primary_psi = embed(omega_psi);
let primary_psi_psi = embed(omega_psi_psi);
let (_, primary_psi_norm, primary_psi_psi_norm, _) =
normalize_penaltywith_psi_derivatives(&primary, &primary_psi, &primary_psi_psi);
let candidates = duchon_native_penalty_candidates(
centers,
spec.length_scale,
spec.power,
effective_nullspace_order,
spec.aniso_log_scales.as_deref(),
&z,
identifiability_transform,
poly_cols,
)?;
let (_, _, penaltyinfo) = filter_active_penalty_candidates(candidates)?;
let mut sources = Vec::new();
let mut first = Vec::new();
let mut second = Vec::new();
for info in penaltyinfo.iter().filter(|info| info.active) {
sources.push(info.source.clone());
match info.source {
PenaltySource::Primary => {
first.push(primary_psi_norm.clone());
second.push(primary_psi_psi_norm.clone());
}
PenaltySource::DoublePenaltyNullspace => {
first.push(Array2::<f64>::zeros(primary_psi_norm.raw_dim()));
second.push(Array2::<f64>::zeros(primary_psi_psi_norm.raw_dim()));
}
ref other => {
crate::bail_invalid_basis!(
"unexpected Duchon native penalty source in derivative path: {other:?}"
);
}
}
}
Ok((sources, first, second))
}
fn prepare_duchon_derivative_contextwithworkspace(
data: ArrayView2<'_, f64>,
spec: &DuchonBasisSpec,
workspace: &mut BasisWorkspace,
) -> Result<(Array2<f64>, Option<Array2<f64>>), BasisError> {
let original_centers = select_centers_by_strategy(data, &spec.center_strategy)?;
let centers = expand_periodic_centers(&original_centers, spec.periodic.as_deref())?;
assert_spatial_centers_below_large_scale_cap(data.ncols(), centers.view())?;
let raw_design = build_duchon_basis_designwithworkspace(
data,
centers.view(),
spec.length_scale,
spec.power,
spec.nullspace_order,
spec.aniso_log_scales.as_deref(),
workspace,
)?;
let identifiability_transform = spatial_identifiability_transform_from_design(
data,
raw_design.basis.view(),
&spec.identifiability,
"Duchon",
)?;
Ok((centers, identifiability_transform))
}
fn prepare_periodic_duchon_centers_1d(
centers: Array2<f64>,
) -> Result<(Array2<f64>, f64, f64), BasisError> {
prepare_periodic_duchon_centers_1d_with_period(centers, None)
}
fn prepare_periodic_duchon_centers_1d_with_period(
centers: Array2<f64>,
explicit_period: Option<f64>,
) -> Result<(Array2<f64>, f64, f64), BasisError> {
if centers.ncols() != 1 {
crate::bail_invalid_basis!(
"periodic Duchon smooths currently require exactly one covariate"
);
}
let left = centers
.column(0)
.iter()
.fold(f64::INFINITY, |a, &b| a.min(b));
let right = centers
.column(0)
.iter()
.fold(f64::NEG_INFINITY, |a, &b| a.max(b));
if !left.is_finite() || !right.is_finite() || left >= right {
return Err(BasisError::InvalidRange(left, right));
}
let span = right - left;
let period = match explicit_period {
Some(p) => {
if !p.is_finite() || p <= 0.0 {
crate::bail_invalid_basis!(
"periodic Duchon period must be finite and positive; got {p}"
);
}
if p < span - 1.0e-10 * span.max(1.0) {
crate::bail_invalid_basis!(
"periodic Duchon period ({p}) is smaller than the center span ({span}); \
every center must lie within a single period"
);
}
p
}
None => span,
};
let centers = collapse_periodic_endpoint(centers, left, period);
Ok((centers, left, period))
}
fn fill_periodic_duchon_kernel_psi_matrices(
rows: ArrayView2<'_, f64>,
centers: ArrayView2<'_, f64>,
left: f64,
period: f64,
length_scale: f64,
p_order: usize,
s_order: usize,
coeffs: &DuchonPartialFractionCoeffs,
) -> Result<(Array2<f64>, Array2<f64>, Array2<f64>), BasisError> {
let n = rows.nrows();
let k = centers.nrows();
let mut kernel = Array2::<f64>::zeros((n, k));
let mut kernel_psi = Array2::<f64>::zeros((n, k));
let mut kernel_psi_psi = Array2::<f64>::zeros((n, k));
for i in 0..n {
let x = wrap_to_period(rows[[i, 0]], left, period);
for j in 0..k {
let r = periodic_distance_1d(x, centers[[j, 0]], period);
let core =
duchon_radial_core_psi_triplet(r, length_scale, p_order, s_order, 1, coeffs)?;
kernel[[i, j]] = core.phi.value;
kernel_psi[[i, j]] = core.phi.psi;
kernel_psi_psi[[i, j]] = core.phi.psi_psi;
}
}
Ok((kernel, kernel_psi, kernel_psi_psi))
}
fn periodic_duchon_identifiability_transformwithworkspace(
data: ArrayView2<'_, f64>,
spec: &DuchonBasisSpec,
centers: Array2<f64>,
workspace: &mut BasisWorkspace,
) -> Result<Option<Array2<f64>>, BasisError> {
let built = build_periodic_duchon_basis_1d(data, spec, centers, workspace)?;
match built.metadata {
BasisMetadata::Duchon {
identifiability_transform,
..
} => Ok(identifiability_transform),
other => Err(BasisError::InvalidInput(format!(
"periodic Duchon builder must return Duchon metadata, got {:?}",
std::mem::discriminant(&other)
))),
}
}
fn build_periodic_duchon_basis_log_kappa_derivativeswithworkspace(
data: ArrayView2<'_, f64>,
spec: &DuchonBasisSpec,
workspace: &mut BasisWorkspace,
) -> Result<BasisPsiDerivativeBundle, BasisError> {
if data.ncols() != 1 {
crate::bail_invalid_basis!(
"periodic Duchon log-kappa derivatives require exactly one covariate"
);
}
let length_scale = spec.length_scale.ok_or_else(|| {
BasisError::InvalidInput(
"periodic Duchon log-kappa derivatives require hybrid Duchon with length_scale"
.to_string(),
)
})?;
let centers = select_centers_by_strategy(data, &spec.center_strategy)?;
assert_spatial_centers_below_large_scale_cap(data.ncols(), centers.view())?;
let (centers, left, period) = prepare_periodic_duchon_centers_1d(centers)?;
let effective_nullspace_order = DuchonNullspaceOrder::Zero;
let p_order = duchon_p_from_nullspace_order(effective_nullspace_order);
let s_order = spec.power_as_usize();
validate_duchon_kernel_orders(Some(length_scale), p_order, s_order as f64, 1)?;
let coeffs = duchon_partial_fraction_coeffs(p_order, s_order, 1.0 / length_scale.max(1e-300));
let z_kernel = kernel_constraint_nullspace(
centers.view(),
effective_nullspace_order,
&mut workspace.cache,
)?;
let identifiability_transform = periodic_duchon_identifiability_transformwithworkspace(
data,
spec,
centers.clone(),
workspace,
)?;
let (_data_kernel, data_kernel_psi, data_kernel_psi_psi) =
fill_periodic_duchon_kernel_psi_matrices(
data,
centers.view(),
left,
period,
length_scale,
p_order,
s_order,
&coeffs,
)?;
let kernel_amp = duchon_kernel_amplification(
centers.view(),
Some(length_scale),
p_order,
s_order,
1,
None,
Some(&coeffs),
None,
);
let kernel_cols = z_kernel.ncols();
let total_cols = kernel_cols + 1;
let mut design_first = Array2::<f64>::zeros((data.nrows(), total_cols));
let mut design_second = Array2::<f64>::zeros((data.nrows(), total_cols));
design_first
.slice_mut(s![.., 0..kernel_cols])
.assign(&(fast_ab(&data_kernel_psi, &z_kernel) * kernel_amp));
design_second
.slice_mut(s![.., 0..kernel_cols])
.assign(&(fast_ab(&data_kernel_psi_psi, &z_kernel) * kernel_amp));
if let Some(transform) = identifiability_transform.as_ref() {
design_first = fast_ab(&design_first, transform);
design_second = fast_ab(&design_second, transform);
}
let (center_kernel, center_kernel_psi, center_kernel_psi_psi) =
fill_periodic_duchon_kernel_psi_matrices(
centers.view(),
centers.view(),
left,
period,
length_scale,
p_order,
s_order,
&coeffs,
)?;
let omega = fast_ab(&fast_atb(&z_kernel, ¢er_kernel), &z_kernel);
let omega_psi = fast_ab(&fast_atb(&z_kernel, ¢er_kernel_psi), &z_kernel);
let omega_psi_psi = fast_ab(&fast_atb(&z_kernel, ¢er_kernel_psi_psi), &z_kernel);
let mut penalty = Array2::<f64>::zeros((total_cols, total_cols));
let mut penalty_psi = Array2::<f64>::zeros((total_cols, total_cols));
let mut penalty_psi_psi = Array2::<f64>::zeros((total_cols, total_cols));
penalty
.slice_mut(s![0..kernel_cols, 0..kernel_cols])
.assign(&omega);
penalty_psi
.slice_mut(s![0..kernel_cols, 0..kernel_cols])
.assign(&omega_psi);
penalty_psi_psi
.slice_mut(s![0..kernel_cols, 0..kernel_cols])
.assign(&omega_psi_psi);
if let Some(transform) = identifiability_transform.as_ref() {
penalty = fast_ab(&fast_atb(transform, &penalty), transform);
penalty_psi = fast_ab(&fast_atb(transform, &penalty_psi), transform);
penalty_psi_psi = fast_ab(&fast_atb(transform, &penalty_psi_psi), transform);
}
let (penalty_norm, penalty_norm_psi, penalty_norm_psi_psi, normalization_scale) =
normalize_penaltywith_psi_derivatives(
&symmetrize(&penalty),
&symmetrize(&penalty_psi),
&symmetrize(&penalty_psi_psi),
);
let (_, _, penaltyinfo) = filter_active_penalty_candidates(vec![PenaltyCandidate {
matrix: penalty_norm,
nullspace_dim_hint: 1,
source: PenaltySource::Primary,
normalization_scale,
kronecker_factors: None,
op: None,
}])?;
let mut penalties_derivative = Vec::new();
let mut penaltiessecond_derivative = Vec::new();
for info in penaltyinfo.iter().filter(|info| info.active) {
match info.source {
PenaltySource::Primary => {
penalties_derivative.push(penalty_norm_psi.clone());
penaltiessecond_derivative.push(penalty_norm_psi_psi.clone());
}
ref other => {
crate::bail_invalid_basis!(
"unexpected periodic Duchon penalty source in derivative path: {other:?}"
);
}
}
}
Ok(BasisPsiDerivativeBundle {
first: BasisPsiDerivativeResult {
design_derivative: design_first,
penalties_derivative,
implicit_operator: None,
},
second: BasisPsiSecondDerivativeResult {
designsecond_derivative: design_second,
penaltiessecond_derivative,
implicit_operator: None,
},
implicit_operator: None,
})
}
fn build_matern_design_psi_derivatives(
data: ArrayView2<'_, f64>,
centers: ArrayView2<'_, f64>,
length_scale: f64,
nu: MaternNu,
include_intercept: bool,
z_opt: Option<&Array2<f64>>,
aniso_log_scales: Option<&[f64]>,
) -> Result<ScalarDesignPsiDerivatives, BasisError> {
let k = centers.nrows();
let kernel_cols = z_opt.map(|z| z.ncols()).unwrap_or(k);
let total_cols = kernel_cols + usize::from(include_intercept);
build_scalar_design_psi_derivatives_shared(
data,
centers,
aniso_log_scales,
total_cols,
z_opt.cloned(),
None,
usize::from(include_intercept),
RadialScalarKind::Matern { length_scale, nu },
0.0,
)
}
fn build_matern_double_penalty_primarywith_psi_derivatives(
centers: ArrayView2<'_, f64>,
length_scale: f64,
nu: MaternNu,
include_intercept: bool,
z_opt: Option<&Array2<f64>>,
aniso_log_scales: Option<&[f64]>,
) -> Result<(Array2<f64>, Array2<f64>, Array2<f64>, f64), BasisError> {
let k = centers.nrows();
let kernel_cols = z_opt.map(|z| z.ncols()).unwrap_or(k);
let total_cols = kernel_cols + usize::from(include_intercept);
let mut kernel = Array2::<f64>::zeros((k, k));
let mut kernel_psi = Array2::<f64>::zeros((k, k));
let mut kernel_psi_psi = Array2::<f64>::zeros((k, k));
for i in 0..k {
for j in i..k {
let r = if let Some(eta) = aniso_log_scales {
aniso_distance(
centers.row(i).as_slice().unwrap(),
centers.row(j).as_slice().unwrap(),
eta,
)
} else {
stable_euclidean_norm(
(0..centers.ncols()).map(|axis| centers[[i, axis]] - centers[[j, axis]]),
)
};
let value = matern_kernel_from_distance(r, length_scale, nu)?;
let d1 = matern_kernel_log_kappa_derivative_from_distance(r, length_scale, nu)?;
let d2 = matern_kernel_log_kappasecond_derivative_from_distance(r, length_scale, nu)?;
kernel[[i, j]] = value;
kernel[[j, i]] = value;
kernel_psi[[i, j]] = d1;
kernel_psi[[j, i]] = d1;
kernel_psi_psi[[i, j]] = d2;
kernel_psi_psi[[j, i]] = d2;
}
}
let (kernel, kernel_psi, kernel_psi_psi) = if let Some(z) = z_opt {
let zt_s = z.t().dot(&kernel);
let zt_d1 = z.t().dot(&kernel_psi);
let zt_d2 = z.t().dot(&kernel_psi_psi);
(zt_s.dot(z), zt_d1.dot(z), zt_d2.dot(z))
} else {
(kernel, kernel_psi, kernel_psi_psi)
};
let mut s = Array2::<f64>::zeros((total_cols, total_cols));
let mut s_psi = Array2::<f64>::zeros((total_cols, total_cols));
let mut s_psi_psi = Array2::<f64>::zeros((total_cols, total_cols));
s.slice_mut(s![0..kernel_cols, 0..kernel_cols])
.assign(&kernel);
s_psi
.slice_mut(s![0..kernel_cols, 0..kernel_cols])
.assign(&kernel_psi);
s_psi_psi
.slice_mut(s![0..kernel_cols, 0..kernel_cols])
.assign(&kernel_psi_psi);
let (s_norm, s_norm_psi, s_norm_psi_psi, c) =
normalize_penaltywith_psi_derivatives(&s, &s_psi, &s_psi_psi);
Ok((s_norm, s_norm_psi, s_norm_psi_psi, c))
}
fn active_matern_double_penalty_derivatives(
penaltyinfo: &[PenaltyInfo],
primary_derivative: &Array2<f64>,
) -> Result<Vec<Array2<f64>>, BasisError> {
penaltyinfo
.iter()
.filter(|info| info.active)
.map(|info| match &info.source {
PenaltySource::Primary => Ok(primary_derivative.clone()),
PenaltySource::DoublePenaltyNullspace => {
Ok(Array2::<f64>::zeros(primary_derivative.raw_dim()))
}
other => Err(BasisError::InvalidInput(format!(
"unexpected Matérn penalty source in double-penalty path: {other:?}"
))),
})
.collect()
}
pub fn build_matern_basis_log_kappa_derivative(
data: ArrayView2<'_, f64>,
spec: &MaternBasisSpec,
) -> Result<BasisPsiDerivativeResult, BasisError> {
let mut workspace = BasisWorkspace::default();
build_matern_basis_log_kappa_derivativewithworkspace(data, spec, &mut workspace)
}
pub fn build_matern_basis_log_kappa_derivativewithworkspace(
data: ArrayView2<'_, f64>,
spec: &MaternBasisSpec,
workspace: &mut BasisWorkspace,
) -> Result<BasisPsiDerivativeResult, BasisError> {
let mut bundle = build_matern_basis_log_kappa_derivativeswithworkspace(data, spec, workspace)?;
bundle.first.implicit_operator = bundle.implicit_operator;
Ok(bundle.first)
}
pub fn build_matern_basis_log_kappa_derivatives(
data: ArrayView2<'_, f64>,
spec: &MaternBasisSpec,
) -> Result<BasisPsiDerivativeBundle, BasisError> {
let mut workspace = BasisWorkspace::default();
build_matern_basis_log_kappa_derivativeswithworkspace(data, spec, &mut workspace)
}
pub fn build_matern_basis_log_kappa_derivativeswithworkspace(
data: ArrayView2<'_, f64>,
spec: &MaternBasisSpec,
workspace: &mut BasisWorkspace,
) -> Result<BasisPsiDerivativeBundle, BasisError> {
let centers = select_centers_by_strategy(data, &spec.center_strategy)?;
let z_opt = matern_identifiability_transform(centers.view(), &spec.identifiability)?;
let aniso = spec.aniso_log_scales.as_deref();
let design_derivatives = build_matern_design_psi_derivatives(
data,
centers.view(),
spec.length_scale,
spec.nu,
spec.include_intercept,
z_opt.as_ref(),
aniso,
)?;
let (penalties_derivative, penaltiessecond_derivative) = if spec.double_penalty {
let base = build_matern_basiswithworkspace(data, spec, workspace)?;
let (_, primary_derivative, primarysecond_derivative, _) =
build_matern_double_penalty_primarywith_psi_derivatives(
centers.view(),
spec.length_scale,
spec.nu,
spec.include_intercept,
z_opt.as_ref(),
aniso,
)?;
(
active_matern_double_penalty_derivatives(&base.penaltyinfo, &primary_derivative)?,
active_matern_double_penalty_derivatives(&base.penaltyinfo, &primarysecond_derivative)?,
)
} else {
build_matern_operator_penalty_psi_derivatives(
centers.view(),
spec.length_scale,
spec.nu,
spec.include_intercept,
z_opt.as_ref(),
aniso,
)?
};
Ok(BasisPsiDerivativeBundle {
first: BasisPsiDerivativeResult {
design_derivative: design_derivatives.design_first,
penalties_derivative,
implicit_operator: None,
},
second: BasisPsiSecondDerivativeResult {
designsecond_derivative: design_derivatives.design_second_diag,
penaltiessecond_derivative,
implicit_operator: None,
},
implicit_operator: design_derivatives.implicit_operator,
})
}
pub fn build_matern_basis_log_kappasecond_derivative(
data: ArrayView2<'_, f64>,
spec: &MaternBasisSpec,
) -> Result<BasisPsiSecondDerivativeResult, BasisError> {
let mut workspace = BasisWorkspace::default();
build_matern_basis_log_kappasecond_derivativewithworkspace(data, spec, &mut workspace)
}
pub fn build_matern_basis_log_kappasecond_derivativewithworkspace(
data: ArrayView2<'_, f64>,
spec: &MaternBasisSpec,
workspace: &mut BasisWorkspace,
) -> Result<BasisPsiSecondDerivativeResult, BasisError> {
let mut bundle = build_matern_basis_log_kappa_derivativeswithworkspace(data, spec, workspace)?;
bundle.second.implicit_operator = bundle.implicit_operator;
Ok(bundle.second)
}
fn build_matern_design_psi_aniso_derivatives(
data: ArrayView2<'_, f64>,
centers: ArrayView2<'_, f64>,
length_scale: f64,
nu: MaternNu,
eta: &[f64],
include_intercept: bool,
z_opt: Option<&Array2<f64>>,
) -> Result<AnisoBasisPsiDerivatives, BasisError> {
let k = centers.nrows();
let p_constrained = z_opt.map(|z| z.ncols()).unwrap_or(k);
let n_poly = usize::from(include_intercept);
let p_smooth = p_constrained + n_poly;
build_aniso_design_psi_derivatives_shared(
data,
centers,
eta,
p_smooth,
z_opt.cloned(),
None,
n_poly,
RadialScalarKind::Matern { length_scale, nu },
)
}
fn build_matern_aniso_primary_raw_derivative_matrices(
centers: ArrayView2<'_, f64>,
eta: &[f64],
length_scale: f64,
nu: MaternNu,
) -> Result<(Vec<Array2<f64>>, Vec<Array2<f64>>), BasisError> {
let k = centers.nrows();
let dim = centers.ncols();
let row_blocks: Result<Vec<_>, BasisError> = (0..k)
.into_par_iter()
.map(|i| {
let ci: Vec<f64> = (0..dim).map(|a| centers[[i, a]]).collect();
let mut first_by_axis: Vec<Vec<f64>> =
(0..dim).map(|_| Vec::with_capacity(k - i)).collect();
let mut second_diag_by_axis: Vec<Vec<f64>> =
(0..dim).map(|_| Vec::with_capacity(k - i)).collect();
for j in i..k {
let cj: Vec<f64> = (0..dim).map(|a| centers[[j, a]]).collect();
let (r, s_vec) = aniso_distance_and_components(&ci, &cj, eta);
let (_, q, t, _, _) = matern_aniso_extended_radial_scalars(r, length_scale, nu)?;
for a in 0..dim {
let s_a = s_vec[a];
first_by_axis[a].push(q * s_a);
second_diag_by_axis[a].push(2.0 * q * s_a + t * s_a * s_a);
}
}
Ok((first_by_axis, second_diag_by_axis))
})
.collect();
let row_blocks = row_blocks?;
let mut raw_first = vec![Array2::<f64>::zeros((k, k)); dim];
let mut raw_second_diag = vec![Array2::<f64>::zeros((k, k)); dim];
for (i, (first_by_axis, second_diag_by_axis)) in row_blocks.into_iter().enumerate() {
for (offset, j) in (i..k).enumerate() {
for a in 0..dim {
let d1 = first_by_axis[a][offset];
let d2 = second_diag_by_axis[a][offset];
raw_first[a][[i, j]] = d1;
raw_first[a][[j, i]] = d1;
raw_second_diag[a][[i, j]] = d2;
raw_second_diag[a][[j, i]] = d2;
}
}
}
Ok((raw_first, raw_second_diag))
}
fn build_matern_aniso_raw_cross_derivative_matrix(
centers: ArrayView2<'_, f64>,
eta: &[f64],
length_scale: f64,
nu: MaternNu,
axis_a: usize,
axis_b: usize,
) -> Result<Array2<f64>, BasisError> {
let k = centers.nrows();
let dim = centers.ncols();
let row_blocks: Result<Vec<_>, BasisError> = (0..k)
.into_par_iter()
.map(|i| {
let ci: Vec<f64> = (0..dim).map(|ax| centers[[i, ax]]).collect();
let mut values = Vec::with_capacity(k - i);
for j in i..k {
let cj: Vec<f64> = (0..dim).map(|ax| centers[[j, ax]]).collect();
let (r, s_vec) = aniso_distance_and_components(&ci, &cj, eta);
let (_, _, t_val, _, _) =
matern_aniso_extended_radial_scalars(r, length_scale, nu)?;
values.push(t_val * s_vec[axis_a] * s_vec[axis_b]);
}
Ok(values)
})
.collect();
let row_blocks = row_blocks?;
let mut raw_cross = Array2::<f64>::zeros((k, k));
for (i, values) in row_blocks.into_iter().enumerate() {
for (offset, j) in (i..k).enumerate() {
let value = values[offset];
raw_cross[[i, j]] = value;
raw_cross[[j, i]] = value;
}
}
Ok(raw_cross)
}
pub fn build_matern_basis_log_kappa_aniso_derivatives(
data: ArrayView2<'_, f64>,
spec: &MaternBasisSpec,
) -> Result<AnisoBasisPsiDerivatives, BasisError> {
let eta = spec.aniso_log_scales.as_deref().ok_or_else(|| {
BasisError::InvalidInput("aniso derivatives require aniso_log_scales to be set".to_string())
})?;
let dim = data.ncols();
if eta.len() != dim {
crate::bail_dim_basis!(
"aniso_log_scales length {} != data dimension {dim}",
eta.len()
);
}
let centers = select_centers_by_strategy(data, &spec.center_strategy)?;
let z_opt = matern_identifiability_transform(centers.view(), &spec.identifiability)?;
let mut result = build_matern_design_psi_aniso_derivatives(
data,
centers.view(),
spec.length_scale,
spec.nu,
eta,
spec.include_intercept,
z_opt.as_ref(),
)?;
if spec.double_penalty {
let k = centers.nrows();
let kernel_cols = z_opt.as_ref().map(|z| z.ncols()).unwrap_or(k);
let total_cols = kernel_cols + usize::from(spec.include_intercept);
let mut primary_first = vec![Array2::<f64>::zeros((total_cols, total_cols)); dim];
let mut primary_second_diag = vec![Array2::<f64>::zeros((total_cols, total_cols)); dim];
let (mut raw_first, mut raw_second_diag) =
build_matern_aniso_primary_raw_derivative_matrices(
centers.view(),
eta,
spec.length_scale,
spec.nu,
)?;
for a in 0..dim {
let projected_first = if let Some(z) = z_opt.as_ref() {
z.t().dot(&raw_first[a]).dot(z)
} else {
std::mem::take(&mut raw_first[a])
};
let projected_second = if let Some(z) = z_opt.as_ref() {
z.t().dot(&raw_second_diag[a]).dot(z)
} else {
std::mem::take(&mut raw_second_diag[a])
};
primary_first[a]
.slice_mut(s![0..kernel_cols, 0..kernel_cols])
.assign(&projected_first);
primary_second_diag[a]
.slice_mut(s![0..kernel_cols, 0..kernel_cols])
.assign(&projected_second);
}
let mut dp_cross_pairs: Vec<(usize, usize)> = Vec::new();
for a in 0..dim {
for b in (a + 1)..dim {
dp_cross_pairs.push((a, b));
}
}
let base = build_matern_basiswithworkspace(data, spec, &mut BasisWorkspace::default())?;
result.penalties_first = Vec::with_capacity(dim);
result.penalties_second_diag = Vec::with_capacity(dim);
for a in 0..dim {
let pf =
active_matern_double_penalty_derivatives(&base.penaltyinfo, &primary_first[a])?;
let ps = active_matern_double_penalty_derivatives(
&base.penaltyinfo,
&primary_second_diag[a],
)?;
result.penalties_first.push(pf);
result.penalties_second_diag.push(ps);
}
result.penalties_cross_pairs = dp_cross_pairs;
let centers_owned = centers.to_owned();
let eta_owned = eta.to_vec();
let z_owned = z_opt.clone();
let penaltyinfo = base.penaltyinfo.clone();
let length_scale = spec.length_scale;
let nu = spec.nu;
result.penalties_cross_provider = Some(AnisoPenaltyCrossProvider::new(
move |axis_a: usize, axis_b: usize| {
let (a, b) = if axis_a < axis_b {
(axis_a, axis_b)
} else {
(axis_b, axis_a)
};
if a == b || b >= eta_owned.len() {
return Ok(Vec::new());
}
let raw_cross = build_matern_aniso_raw_cross_derivative_matrix(
centers_owned.view(),
&eta_owned,
length_scale,
nu,
a,
b,
)?;
let projected: Array2<f64> = if let Some(z) = z_owned.as_ref() {
z.t().dot(&raw_cross).dot(z)
} else {
raw_cross
};
let mut padded = Array2::<f64>::zeros((total_cols, total_cols));
padded
.slice_mut(s![0..kernel_cols, 0..kernel_cols])
.assign(&projected);
active_matern_double_penalty_derivatives(&penaltyinfo, &padded)
},
));
} else {
let (per_axis, cross_pairs, cross_provider) =
build_matern_operator_penalty_aniso_derivatives(
centers.view(),
spec.length_scale,
spec.nu,
spec.include_intercept,
z_opt.as_ref(),
eta,
)?;
result.penalties_first = Vec::with_capacity(dim);
result.penalties_second_diag = Vec::with_capacity(dim);
for (pen_first, pen_second) in per_axis {
result.penalties_first.push(pen_first);
result.penalties_second_diag.push(pen_second);
}
result.penalties_cross_pairs = cross_pairs;
result.penalties_cross_provider = Some(cross_provider);
}
Ok(result)
}
fn duchon_coeff_exponents(p_order: usize, s_order: usize, m_or_n: usize) -> f64 {
-2.0 * (p_order + s_order - m_or_n) as f64
}
#[inline(always)]
fn duchon_scaling_exponent(p_order: usize, s_order: usize, k_dim: usize) -> f64 {
k_dim as f64 - 2.0 * (p_order + s_order) as f64
}
#[derive(Clone, Copy)]
struct DuchonMaternDerivativeTerm {
coeff: f64,
kappa_power: usize,
r_power: f64,
bessel_order: f64,
}
#[derive(Clone, Copy, Debug, Default)]
struct PsiTriplet {
value: f64,
psi: f64,
psi_psi: f64,
}
#[derive(Clone, Copy, Debug, Default)]
struct DuchonRadialCore {
phi: PsiTriplet,
}
#[derive(Clone, Copy, Debug, Default)]
struct DuchonRadialJets {
phi: f64,
phi_r: f64,
phi_rr: f64,
phi_rrr: f64,
q: f64,
q_r: f64,
q_rr: f64,
lap: f64,
lap_r: f64,
lap_rr: f64,
t: f64,
t_r: f64,
t_rr: f64,
}
#[derive(Clone, Copy, Debug, Default)]
struct DuchonRegularizedOperatorCore {
q: f64,
t: f64,
t_r: f64,
t_rr: f64,
}
#[inline(always)]
fn duchon_operator_jets_from_primary_core(
core: DuchonRegularizedOperatorCore,
r: f64,
d: f64,
) -> DuchonRadialJets {
let r2 = r * r;
let mut out = DuchonRadialJets {
q: core.q,
t: core.t,
t_r: core.t_r,
t_rr: core.t_rr,
..DuchonRadialJets::default()
};
out.q_r = r * out.t;
out.q_rr = out.t + r * out.t_r;
out.lap = d * out.q + r2 * out.t;
out.lap_r = (d + 2.0) * r * out.t + r2 * out.t_r;
out.lap_rr = (d + 2.0) * out.t + (d + 4.0) * r * out.t_r + r2 * out.t_rr;
out.phi_r = r * out.q;
out.phi_rr = out.q + r2 * out.t;
out.phi_rrr = 3.0 * r * out.t + r2 * out.t_r;
assert!(
((out.phi_rr - (out.q + r * out.q_r)).abs()) <= 1e-10 * out.phi_rr.abs().max(1.0),
"radial scalar identity failed: phi_rr != q + r*q_r, phi_rr={}, q={}, r={}, q_r={}",
out.phi_rr,
out.q,
r,
out.q_r
);
assert!(
((out.phi_rr - (out.q + r2 * out.t)).abs()) <= 1e-10 * out.phi_rr.abs().max(1.0),
"radial scalar identity failed: phi_rr != q + r2*t, phi_rr={}, q={}, r2={}, t={}",
out.phi_rr,
out.q,
r2,
out.t
);
assert!(
((out.lap - (d * out.q + r2 * out.t)).abs()) <= 1e-10 * out.lap.abs().max(1.0),
"radial scalar identity failed: lap != d*q + r2*t, lap={}, d={}, q={}, r2={}, t={}",
out.lap,
d,
out.q,
r2,
out.t
);
out
}
#[inline(always)]
fn scaled_log_kappa_derivatives(
value: f64,
radial_first: f64,
radialsecond: f64,
exponent: f64,
r: f64,
) -> (f64, f64) {
let first = exponent * value + r * radial_first;
let second = exponent * exponent * value
+ (2.0 * exponent + 1.0) * r * radial_first
+ r * r * radialsecond;
(first, second)
}
#[inline(always)]
fn duchon_q_psi_triplet_from_jets(
jets: &DuchonRadialJets,
p_order: usize,
s_order: usize,
k_dim: usize,
r: f64,
) -> (f64, f64) {
scaled_log_kappa_derivatives(
jets.q,
jets.q_r,
jets.q_rr,
duchon_operator_scaling_exponent(p_order, s_order, k_dim),
r,
)
}
#[inline(always)]
fn duchon_operator_scaling_exponent(p_order: usize, s_order: usize, k_dim: usize) -> f64 {
duchon_scaling_exponent(p_order, s_order, k_dim) + 2.0
}
fn duchon_regularized_operator_core(
r_eval: f64,
kappa: f64,
k_dim: usize,
coeffs: &DuchonPartialFractionCoeffs,
) -> Result<DuchonRegularizedOperatorCore, BasisError> {
let mut q_sum = KahanSum::default();
let mut t_sum = KahanSum::default();
let mut t_r_sum = KahanSum::default();
let mut t_rr_sum = KahanSum::default();
for (m, coeff) in coeffs.a.iter().enumerate().skip(1) {
if *coeff == 0.0 {
continue;
}
let (q, t, t_r, t_rr) = duchon_polyharmonic_operator_block_jets(r_eval, m as f64, k_dim)?;
q_sum.add(coeff * q);
t_sum.add(coeff * t);
t_r_sum.add(coeff * t_r);
t_rr_sum.add(coeff * t_rr);
}
let max_ladder_steps = coeffs
.b
.iter()
.enumerate()
.skip(1)
.filter(|(_, coeff)| **coeff != 0.0)
.map(|(n, _)| duchon_matern_block_max_ladder_steps(n, k_dim))
.max();
if let Some(max_ladder_steps) = max_ladder_steps {
let ladder =
BesselKLadder::build(kappa * r_eval, !k_dim.is_multiple_of(2), max_ladder_steps);
for (n, coeff) in coeffs.b.iter().enumerate().skip(1) {
if *coeff == 0.0 {
continue;
}
let (q, t, t_r, t_rr) =
duchon_matern_operator_block_jets_with_ladder(r_eval, kappa, n, k_dim, &ladder)?;
q_sum.add(coeff * q);
t_sum.add(coeff * t);
t_r_sum.add(coeff * t_r);
t_rr_sum.add(coeff * t_rr);
}
}
Ok(DuchonRegularizedOperatorCore {
q: q_sum.sum(),
t: t_sum.sum(),
t_r: t_r_sum.sum(),
t_rr: t_rr_sum.sum(),
})
}
#[inline(always)]
fn duchon_collision_taylor_operator_core(
r: f64,
phi_rr_collision: f64,
t_collision: f64,
t_rr_collision: f64,
) -> DuchonRegularizedOperatorCore {
let r2 = r * r;
let r4 = r2 * r2;
DuchonRegularizedOperatorCore {
q: phi_rr_collision + 0.5 * t_collision * r2 + 0.125 * t_rr_collision * r4,
t: t_collision + 0.5 * t_rr_collision * r2,
t_r: t_rr_collision * r,
t_rr: t_rr_collision,
}
}
fn duchon_radial_jets(
r: f64,
length_scale: f64,
p_order: usize,
s_order: usize,
k_dim: usize,
coeffs: &DuchonPartialFractionCoeffs,
) -> Result<DuchonRadialJets, BasisError> {
let kappa = 1.0 / length_scale.max(1e-300);
let r_floor = DUCHON_DERIVATIVE_R_FLOOR_REL * length_scale.max(1e-8);
let collision_taylor_radius = DUCHON_COLLISION_TAYLOR_REL * length_scale.max(1e-8);
let r_eval = r.max(r_floor);
let d = k_dim as f64;
let phi = duchon_matern_kernel_general_from_distance(
r,
Some(length_scale),
p_order,
s_order,
k_dim,
Some(coeffs),
)?;
if !phi.is_finite() {
crate::bail_invalid_basis!(
"non-finite Duchon radial kernel value at r={r}, length_scale={length_scale}, p={p_order}, s={s_order}, dim={k_dim}"
);
}
let generic_jets = duchon_operator_jets_from_primary_core(
duchon_regularized_operator_core(r_eval, kappa, k_dim, coeffs)?,
r_eval,
d,
);
let mut out = DuchonRadialJets {
phi,
..generic_jets
};
let smoothness_order = 2 * (p_order + s_order);
let collision_q_exists = smoothness_order > k_dim + 2;
let collision_t_exists = smoothness_order > k_dim + 4;
let collision_t_rr_exists = smoothness_order > k_dim + 6;
if r <= collision_taylor_radius.max(r_floor) && collision_t_exists {
let (analytic_phi_rr, _, _) =
duchonphi_rr_collision_psi_triplet(length_scale, p_order, s_order, k_dim, coeffs)?;
let analytic_t_collision =
duchon_phi_rrrr_collision(length_scale, p_order, s_order, k_dim, coeffs)? / 3.0;
let analytic_t_rr_collision = if collision_t_rr_exists {
duchon_phi_rrrrrr_collision(length_scale, p_order, s_order, k_dim, coeffs)? / 15.0
} else {
0.0
};
let collision_jets = duchon_operator_jets_from_primary_core(
duchon_collision_taylor_operator_core(
r,
analytic_phi_rr,
analytic_t_collision,
analytic_t_rr_collision,
),
r,
d,
);
out = DuchonRadialJets {
phi: out.phi,
..collision_jets
};
} else if r < r_floor && collision_q_exists {
let (analytic_phi_rr, _, _) =
duchonphi_rr_collision_psi_triplet(length_scale, p_order, s_order, k_dim, coeffs)?;
out.phi_r = analytic_phi_rr * r;
out.phi_rr = analytic_phi_rr;
out.q = analytic_phi_rr;
out.q_r = 0.0;
out.lap = d * analytic_phi_rr;
out.lap_r = 0.0;
}
if !out.phi_r.is_finite()
|| !out.phi_rr.is_finite()
|| !out.phi_rrr.is_finite()
|| !out.q.is_finite()
|| !out.q_r.is_finite()
|| !out.q_rr.is_finite()
|| !out.lap.is_finite()
|| !out.lap_r.is_finite()
|| !out.lap_rr.is_finite()
|| !out.t.is_finite()
|| !out.t_r.is_finite()
|| !out.t_rr.is_finite()
{
crate::bail_invalid_basis!(
"non-finite Duchon radial jets at r={r}, length_scale={length_scale}, p={p_order}, s={s_order}, dim={k_dim}"
);
}
Ok(out)
}
fn duchon_radial_core_psi_triplet(
r: f64,
length_scale: f64,
p_order: usize,
s_order: usize,
k_dim: usize,
coeffs: &DuchonPartialFractionCoeffs,
) -> Result<DuchonRadialCore, BasisError> {
let delta = duchon_scaling_exponent(p_order, s_order, k_dim);
let jets = duchon_radial_jets(r, length_scale, p_order, s_order, k_dim, coeffs)?;
let phi = jets.phi;
let (phi_psi, phi_psi_psi) =
scaled_log_kappa_derivatives(phi, jets.phi_r, jets.phi_rr, delta, r);
if r > 1e-10 {
assert!(
((delta * phi + r * jets.phi_r) - phi_psi).abs() < 1e-7_f64.max(1e-7_f64 * phi.abs())
);
return Ok(DuchonRadialCore {
phi: PsiTriplet {
value: phi,
psi: phi_psi,
psi_psi: phi_psi_psi,
},
});
}
Ok(DuchonRadialCore {
phi: PsiTriplet {
value: phi,
psi: phi_psi,
psi_psi: phi_psi_psi,
},
})
}
fn duchonphi_rr_collision_psi_triplet(
length_scale: f64,
p_order: usize,
s_order: usize,
k_dim: usize,
coeffs: &DuchonPartialFractionCoeffs,
) -> Result<(f64, f64, f64), BasisError> {
duchon_phi_even_derivative_collision_psi_triplet(
length_scale,
p_order,
s_order,
k_dim,
coeffs,
1,
)
}
const EULER_MASCHERONI: f64 = 0.577_215_664_901_532_9;
#[inline(always)]
fn digamma_pos_int(n: usize) -> f64 {
assert!(n >= 1, "digamma_pos_int requires n >= 1: n={n}");
let mut h = 0.0_f64;
for j in 1..n {
h += 1.0 / j as f64;
}
-EULER_MASCHERONI + h
}
fn duchon_matern_block_taylor_r2j(
kappa: f64,
n_order: usize,
k_dim: usize,
j: usize,
) -> (f64, f64) {
let n = n_order as f64;
let k_half = 0.5 * k_dim as f64;
let nu = n - k_half;
let c = kappa.powf(k_half - n)
/ ((2.0 * std::f64::consts::PI).powf(k_half) * 2.0_f64.powf(n - 1.0) * gamma_lanczos(n));
if k_dim.is_multiple_of(2) {
let nu_int = n_order as i64 - (k_dim as i64) / 2;
duchon_matern_block_taylor_r2j_integer_nu(kappa, c, nu_int, j)
} else {
duchon_matern_block_taylor_r2j_half_integer_nu(kappa, c, nu, j)
}
}
#[inline(always)]
fn psi_power_triplet(value: f64, exponent: f64) -> (f64, f64, f64) {
(value, exponent * value, exponent * exponent * value)
}
#[inline(always)]
fn psi_power_log_triplet(base: f64, exponent: f64, log_kappa_half: f64) -> (f64, f64, f64) {
(
base * log_kappa_half,
base * (exponent * log_kappa_half + 1.0),
base * (exponent * exponent * log_kappa_half + 2.0 * exponent),
)
}
#[inline(always)]
fn add_triplet(dst: &mut (f64, f64, f64), inc: (f64, f64, f64)) {
dst.0 += inc.0;
dst.1 += inc.1;
dst.2 += inc.2;
}
fn duchon_matern_block_taylor_r2j_triplet(
kappa: f64,
n_order: usize,
k_dim: usize,
j: usize,
) -> ((f64, f64, f64), (f64, f64, f64)) {
let n = n_order as f64;
let k_half = 0.5 * k_dim as f64;
let nu = n - k_half;
let c_const = 1.0
/ ((2.0 * std::f64::consts::PI).powf(k_half) * 2.0_f64.powf(n - 1.0) * gamma_lanczos(n));
let c_exp = k_half - n;
let mut pure = (0.0, 0.0, 0.0);
let mut log_part = (0.0, 0.0, 0.0);
let log_kappa_half = (0.5 * kappa).ln();
if k_dim.is_multiple_of(2) {
let nu_int = n_order as i64 - (k_dim as i64) / 2;
let mu = nu_int.unsigned_abs() as usize;
let sign_mu = if mu.is_multiple_of(2) { 1.0 } else { -1.0 };
if nu_int >= 0 {
let nu_usize = nu_int as usize;
if j < nu_usize {
let sign = if j.is_multiple_of(2) { 1.0 } else { -1.0 };
let power = 2 * j as i32 - nu_usize as i32;
let coeff = 0.5 * sign * gamma_lanczos((nu_usize - j) as f64)
/ gamma_lanczos((j + 1) as f64)
* 2.0_f64.powi(-power);
let exponent = c_exp + power as f64;
let value = c_const * coeff * kappa.powf(exponent);
add_triplet(&mut pure, psi_power_triplet(value, exponent));
}
if j >= nu_usize {
let k = j - nu_usize;
let inv_fac = 1.0
/ (gamma_lanczos((k + 1) as f64) * gamma_lanczos((nu_usize + k + 1) as f64));
let power = (2 * k + nu_usize) as i32;
let exponent = c_exp + power as f64;
let kp_base = c_const * kappa.powf(exponent) * 2.0_f64.powi(-power);
let log_base = -sign_mu * kp_base * inv_fac;
add_triplet(&mut log_part, psi_power_triplet(log_base, exponent));
add_triplet(
&mut pure,
psi_power_log_triplet(log_base, exponent, log_kappa_half),
);
let psi_sum = digamma_pos_int(k + 1) + digamma_pos_int(nu_usize + k + 1);
let digamma_base = sign_mu * 0.5 * kp_base * inv_fac * psi_sum;
add_triplet(&mut pure, psi_power_triplet(digamma_base, exponent));
}
} else {
let k = j;
let inv_fac =
1.0 / (gamma_lanczos((k + 1) as f64) * gamma_lanczos((mu + k + 1) as f64));
let power = (mu + 2 * k) as i32;
let exponent = c_exp + power as f64;
let kp_base = c_const * kappa.powf(exponent) * 2.0_f64.powi(-power);
let log_base = -sign_mu * kp_base * inv_fac;
add_triplet(&mut log_part, psi_power_triplet(log_base, exponent));
add_triplet(
&mut pure,
psi_power_log_triplet(log_base, exponent, log_kappa_half),
);
let psi_sum = digamma_pos_int(k + 1) + digamma_pos_int(mu + k + 1);
let digamma_base = sign_mu * 0.5 * kp_base * inv_fac * psi_sum;
add_triplet(&mut pure, psi_power_triplet(digamma_base, exponent));
}
} else {
let nu_abs = nu.abs();
let l = (2.0 * nu_abs - 1.0).round().max(0.0) as usize;
let prefactor_const = (std::f64::consts::PI / 2.0).sqrt();
let prefactor_exp = -0.5;
let target = 2 * j;
for i in 0..=l {
let c_i = gamma_lanczos((l + i + 1) as f64)
/ (gamma_lanczos((i + 1) as f64) * gamma_lanczos((l - i + 1) as f64));
let p_f64 = nu - 0.5 - i as f64;
let p_round = p_f64.round() as i64;
if (p_f64 - p_round as f64).abs() > 1e-12 {
continue;
}
let q_needed = target as i64 - p_round;
if q_needed < 0 {
continue;
}
let q = q_needed as usize;
let sign = if q.is_multiple_of(2) { 1.0 } else { -1.0 };
let exponent = c_exp + prefactor_exp - i as f64 + q as f64;
let value = c_const * prefactor_const * c_i * 2.0_f64.powi(-(i as i32)) * sign
/ gamma_lanczos((q + 1) as f64)
* kappa.powf(exponent);
add_triplet(&mut pure, psi_power_triplet(value, exponent));
}
}
(pure, log_part)
}
fn duchon_matern_block_taylor_r2j_integer_nu(
kappa: f64,
c: f64,
nu_int: i64,
j: usize,
) -> (f64, f64) {
let mu = nu_int.unsigned_abs() as usize;
let kappa_half = 0.5 * kappa;
if nu_int >= 0 {
let nu = nu_int as usize;
let mut pure = 0.0;
let mut log_part = 0.0;
if j < nu {
let sign = if j.is_multiple_of(2) { 1.0 } else { -1.0 };
let coeff = sign * gamma_lanczos((nu - j) as f64) / gamma_lanczos((j + 1) as f64)
* kappa_half.powi(2 * j as i32 - nu as i32)
* 0.5;
pure += coeff;
}
if j >= nu {
let k = j - nu;
let inv_fac =
1.0 / (gamma_lanczos((k + 1) as f64) * gamma_lanczos((nu + k + 1) as f64));
let kp = kappa_half.powi(2 * k as i32 + nu as i32);
let sign_mu = if mu.is_multiple_of(2) { 1.0 } else { -1.0 };
log_part += -sign_mu * kp * inv_fac;
let psi_sum = digamma_pos_int(k + 1) + digamma_pos_int(nu + k + 1);
pure += -sign_mu * kp * inv_fac * kappa_half.ln();
pure += sign_mu * 0.5 * kp * inv_fac * psi_sum;
}
(c * pure, c * log_part)
} else {
let k = j;
let inv_fac = 1.0 / (gamma_lanczos((k + 1) as f64) * gamma_lanczos((mu + k + 1) as f64));
let kp = kappa_half.powi(mu as i32 + 2 * k as i32);
let sign_mu = if mu.is_multiple_of(2) { 1.0 } else { -1.0 };
let log_part = -sign_mu * kp * inv_fac;
let psi_sum = digamma_pos_int(k + 1) + digamma_pos_int(mu + k + 1);
let pure =
-sign_mu * kp * inv_fac * kappa_half.ln() + sign_mu * 0.5 * kp * inv_fac * psi_sum;
(c * pure, c * log_part)
}
}
fn duchon_matern_block_taylor_r2j_half_integer_nu(
kappa: f64,
c: f64,
nu: f64,
j: usize,
) -> (f64, f64) {
let nu_abs = nu.abs();
let l = (2.0 * nu_abs - 1.0).round().max(0.0) as usize;
let prefactor = (std::f64::consts::PI / (2.0 * kappa)).sqrt();
let target = 2 * j;
let mut pure = 0.0;
for i in 0..=l {
let c_i = gamma_lanczos((l + i + 1) as f64)
/ (gamma_lanczos((i + 1) as f64) * gamma_lanczos((l - i + 1) as f64));
let inv_2kappa_i = (2.0 * kappa).powi(-(i as i32));
let p_f64 = nu - 0.5 - i as f64;
let p_round = p_f64.round() as i64;
if (p_f64 - p_round as f64).abs() > 1e-12 {
continue;
}
let q_needed = target as i64 - p_round;
if q_needed < 0 {
continue;
}
let q = q_needed as usize;
let exp_coeff = (-kappa).powi(q as i32) / gamma_lanczos((q + 1) as f64);
pure += c_i * inv_2kappa_i * exp_coeff;
}
(c * prefactor * pure, 0.0) }
fn duchon_polyharmonic_block_taylor_r2j(m: usize, k_dim: usize, j: usize) -> (f64, f64) {
let k_half = 0.5 * k_dim as f64;
let alpha = 2 * m as i64 - k_dim as i64;
if alpha != 2 * j as i64 {
return (0.0, 0.0);
}
if k_dim.is_multiple_of(2) && m >= k_dim / 2 {
let c = polyharmonic_log_sign(m, k_dim)
/ (2.0_f64.powi((2 * m - 1) as i32)
* std::f64::consts::PI.powf(k_half)
* gamma_lanczos(m as f64)
* gamma_lanczos((m - k_dim / 2 + 1) as f64));
(0.0, c)
} else {
let c = gamma_lanczos(k_half - m as f64)
/ (4.0_f64.powi(m as i32)
* std::f64::consts::PI.powf(k_half)
* gamma_lanczos(m as f64));
(c, 0.0)
}
}
fn duchon_phi_even_derivative_collision(
length_scale: f64,
p_order: usize,
s_order: usize,
k_dim: usize,
coeffs: &DuchonPartialFractionCoeffs,
j: usize,
) -> Result<f64, BasisError> {
let smoothness_order = 2 * (p_order + s_order);
let required = k_dim + 2 * j;
if smoothness_order <= required {
let min_power = (required / 2 + 1).saturating_sub(p_order);
crate::bail_invalid_basis!(
"Duchon collision derivative phi^({}) requires 2*(p+s) > dimension+{}; got 2*(p+s)={}, dimension={}, p={}, s={}. \
This path needs the {}-order radial-kernel derivative at the origin, which is finite only for a smoother spline: raise power to >= {} (or reduce the joint smooth's dimension).",
2 * j,
2 * j,
smoothness_order,
k_dim,
p_order,
s_order,
2 * j,
min_power
);
}
let kappa = 1.0 / length_scale.max(1e-300);
let mut total_pure = KahanSum::default();
let mut total_log = KahanSum::default();
let mut total_log_abs_scale = KahanSum::default();
for (m, &a_m) in coeffs.a.iter().enumerate().skip(1) {
if a_m == 0.0 {
continue;
}
let (pure, log) = duchon_polyharmonic_block_taylor_r2j(m, k_dim, j);
total_pure.add(a_m * pure);
total_log.add(a_m * log);
total_log_abs_scale.add((a_m * log).abs());
}
for (n, &b_n) in coeffs.b.iter().enumerate().skip(1) {
if b_n == 0.0 {
continue;
}
let (pure, log) = duchon_matern_block_taylor_r2j(kappa, n, k_dim, j);
total_pure.add(b_n * pure);
total_log.add(b_n * log);
total_log_abs_scale.add((b_n * log).abs());
}
let total_pure = total_pure.sum();
let total_log = total_log.sum();
let total_log_abs_scale = total_log_abs_scale.sum();
let log_cancel_tol = 1e-10 * total_log_abs_scale.max(total_pure.abs()).max(1e-30);
if total_log.abs() > log_cancel_tol {
crate::bail_invalid_basis!(
"Duchon Taylor a_{} log-coefficient did not cancel: log={total_log:.6e}, pure={total_pure:.6e}; \
log_abs_scale={total_log_abs_scale:.6e}, tol={log_cancel_tol:.6e}; p={p_order}, s={s_order}, d={k_dim}",
2 * j
);
}
let factorial_2j = gamma_lanczos((2 * j + 1) as f64);
Ok(factorial_2j * total_pure)
}
fn duchon_phi_even_derivative_collision_psi_triplet(
length_scale: f64,
p_order: usize,
s_order: usize,
k_dim: usize,
coeffs: &DuchonPartialFractionCoeffs,
j: usize,
) -> Result<(f64, f64, f64), BasisError> {
let smoothness_order = 2 * (p_order + s_order);
let required = k_dim + 2 * j;
if smoothness_order <= required {
let min_power = (required / 2 + 1).saturating_sub(p_order);
crate::bail_invalid_basis!(
"Duchon collision derivative phi^({}) psi triplet requires 2*(p+s) > dimension+{}; got 2*(p+s)={}, dimension={}, p={}, s={}. \
The exact two-block / transformation-normal path needs analytic length-scale derivatives of the kernel, which are finite only for a smoother spline: raise power to >= {} (or reduce the joint smooth's dimension).",
2 * j,
2 * j,
smoothness_order,
k_dim,
p_order,
s_order,
min_power
);
}
let kappa = 1.0 / length_scale.max(1e-300);
let mut value = KahanSum::default();
let mut psi = KahanSum::default();
let mut psi_psi = KahanSum::default();
let mut log_value = KahanSum::default();
let mut log_psi = KahanSum::default();
let mut log_psi_psi = KahanSum::default();
let mut log_abs_scale = KahanSum::default();
for (m, &a_m) in coeffs.a.iter().enumerate().skip(1) {
if a_m == 0.0 {
continue;
}
let alpha_m = duchon_coeff_exponents(p_order, s_order, m);
let (pure, log) = duchon_polyharmonic_block_taylor_r2j(m, k_dim, j);
value.add(a_m * pure);
psi.add(alpha_m * a_m * pure);
psi_psi.add(alpha_m * alpha_m * a_m * pure);
log_value.add(a_m * log);
log_psi.add(alpha_m * a_m * log);
log_psi_psi.add(alpha_m * alpha_m * a_m * log);
log_abs_scale.add((a_m * log).abs());
log_abs_scale.add((alpha_m * a_m * log).abs());
log_abs_scale.add((alpha_m * alpha_m * a_m * log).abs());
}
for (n, &b_n) in coeffs.b.iter().enumerate().skip(1) {
if b_n == 0.0 {
continue;
}
let beta_n = duchon_coeff_exponents(p_order, s_order, n);
let (pure, log) = duchon_matern_block_taylor_r2j_triplet(kappa, n, k_dim, j);
value.add(b_n * pure.0);
psi.add(beta_n * b_n * pure.0 + b_n * pure.1);
psi_psi.add(beta_n * beta_n * b_n * pure.0 + 2.0 * beta_n * b_n * pure.1 + b_n * pure.2);
log_value.add(b_n * log.0);
log_psi.add(beta_n * b_n * log.0 + b_n * log.1);
log_psi_psi.add(beta_n * beta_n * b_n * log.0 + 2.0 * beta_n * b_n * log.1 + b_n * log.2);
let log_v = b_n * log.0;
let log_p = beta_n * b_n * log.0 + b_n * log.1;
let log_pp = beta_n * beta_n * b_n * log.0 + 2.0 * beta_n * b_n * log.1 + b_n * log.2;
log_abs_scale.add(log_v.abs());
log_abs_scale.add(log_p.abs());
log_abs_scale.add(log_pp.abs());
}
let value = value.sum();
let psi = psi.sum();
let psi_psi = psi_psi.sum();
let log_value = log_value.sum();
let log_psi = log_psi.sum();
let log_psi_psi = log_psi_psi.sum();
let log_abs_scale = log_abs_scale.sum();
let scale = value.abs().max(psi.abs()).max(psi_psi.abs()).max(1e-30);
let log_cancel_tol = 1e-10 * log_abs_scale.max(scale);
if log_value.abs().max(log_psi.abs()).max(log_psi_psi.abs()) > log_cancel_tol {
crate::bail_invalid_basis!(
"Duchon Taylor a_{} log-coefficient derivative did not cancel: \
log=({log_value:.6e}, {log_psi:.6e}, {log_psi_psi:.6e}), \
value=({value:.6e}, {psi:.6e}, {psi_psi:.6e}), log_abs_scale={log_abs_scale:.6e}, tol={log_cancel_tol:.6e}; \
p={p_order}, s={s_order}, d={k_dim}",
2 * j
);
}
let factorial_2j = gamma_lanczos((2 * j + 1) as f64);
Ok((
factorial_2j * value,
factorial_2j * psi,
factorial_2j * psi_psi,
))
}
fn duchon_phi_rrrr_collision(
length_scale: f64,
p_order: usize,
s_order: usize,
k_dim: usize,
coeffs: &DuchonPartialFractionCoeffs,
) -> Result<f64, BasisError> {
duchon_phi_even_derivative_collision(length_scale, p_order, s_order, k_dim, coeffs, 2)
}
fn duchon_phi_rrrrrr_collision(
length_scale: f64,
p_order: usize,
s_order: usize,
k_dim: usize,
coeffs: &DuchonPartialFractionCoeffs,
) -> Result<f64, BasisError> {
duchon_phi_even_derivative_collision(length_scale, p_order, s_order, k_dim, coeffs, 3)
}
fn build_duchon_design_psi_derivativeswithworkspace(
data: ArrayView2<'_, f64>,
centers: ArrayView2<'_, f64>,
spec: &DuchonBasisSpec,
identifiability_transform: Option<&Array2<f64>>,
workspace: &mut BasisWorkspace,
) -> Result<ScalarDesignPsiDerivatives, BasisError> {
let length_scale = spec.length_scale.ok_or_else(|| {
BasisError::InvalidInput(
"exact Duchon log-kappa derivatives require hybrid Duchon with length_scale"
.to_string(),
)
})?;
let effective_nullspace_order = duchon_effective_nullspace_order(centers, spec.nullspace_order);
let p_order = duchon_p_from_nullspace_order(effective_nullspace_order);
let s_order = spec.power_as_usize();
let kappa = 1.0 / length_scale;
let coeffs = duchon_partial_fraction_coeffs(p_order, s_order, kappa);
let z_kernel =
kernel_constraint_nullspace(centers, effective_nullspace_order, &mut workspace.cache)?;
let poly_cols = polynomial_block_from_order(data, effective_nullspace_order).ncols();
let p_padded = z_kernel.ncols() + poly_cols;
if let Some(zf) = identifiability_transform
&& p_padded != zf.nrows()
{
crate::bail_dim_basis!(
"Duchon identifiability transform mismatch in design derivatives: local cols={}, transform rows={}",
p_padded,
zf.nrows()
);
}
let p_final = identifiability_transform
.map(|zf| zf.ncols())
.unwrap_or(p_padded);
build_scalar_design_psi_derivatives_shared(
data,
centers,
spec.aniso_log_scales.as_deref(),
p_final,
Some(z_kernel),
identifiability_transform.cloned(),
poly_cols,
RadialScalarKind::Duchon {
length_scale,
p_order,
s_order,
dim: data.ncols(),
coeffs,
},
duchon_scaling_exponent(p_order, s_order, data.ncols()),
)
}
pub fn build_duchon_basis_log_kappa_derivative(
data: ArrayView2<'_, f64>,
spec: &DuchonBasisSpec,
) -> Result<BasisPsiDerivativeResult, BasisError> {
let mut workspace = BasisWorkspace::default();
build_duchon_basis_log_kappa_derivativewithworkspace(data, spec, &mut workspace)
}
pub fn build_duchon_basis_log_kappa_derivativewithworkspace(
data: ArrayView2<'_, f64>,
spec: &DuchonBasisSpec,
workspace: &mut BasisWorkspace,
) -> Result<BasisPsiDerivativeResult, BasisError> {
let mut bundle = build_duchon_basis_log_kappa_derivativeswithworkspace(data, spec, workspace)?;
bundle.first.implicit_operator = bundle.implicit_operator;
Ok(bundle.first)
}
pub fn build_duchon_basis_log_kappa_derivatives(
data: ArrayView2<'_, f64>,
spec: &DuchonBasisSpec,
) -> Result<BasisPsiDerivativeBundle, BasisError> {
let mut workspace = BasisWorkspace::default();
build_duchon_basis_log_kappa_derivativeswithworkspace(data, spec, &mut workspace)
}
fn duchon_operator_penalties_requested(spec: &DuchonOperatorPenaltySpec) -> bool {
matches!(spec.mass, OperatorPenaltySpec::Active { .. })
|| matches!(spec.tension, OperatorPenaltySpec::Active { .. })
|| matches!(spec.stiffness, OperatorPenaltySpec::Active { .. })
}
pub fn build_duchon_basis_log_kappa_derivativeswithworkspace(
data: ArrayView2<'_, f64>,
spec: &DuchonBasisSpec,
workspace: &mut BasisWorkspace,
) -> Result<BasisPsiDerivativeBundle, BasisError> {
if spec.periodic.is_some() {
return build_periodic_duchon_basis_log_kappa_derivativeswithworkspace(
data, spec, workspace,
);
}
let (centers, identifiability_transform) =
prepare_duchon_derivative_contextwithworkspace(data, spec, workspace)?;
let operator_collocation_points =
if duchon_operator_penalties_requested(&spec.operator_penalties) {
let m = (DUCHON_COLLOCATION_OVERSAMPLE * centers.nrows()).min(data.nrows());
Some(select_thin_plate_knots(data, m)?)
} else {
None
};
build_duchon_basis_log_kappa_derivativeswith_collocationwithworkspace(
data,
spec,
centers.view(),
identifiability_transform.as_ref(),
operator_collocation_points
.as_ref()
.map(|points| points.view()),
workspace,
)
}
pub(crate) fn build_duchon_basis_log_kappa_derivativeswith_collocationwithworkspace(
data: ArrayView2<'_, f64>,
spec: &DuchonBasisSpec,
centers: ArrayView2<'_, f64>,
identifiability_transform: Option<&Array2<f64>>,
operator_collocation_points: Option<ArrayView2<'_, f64>>,
workspace: &mut BasisWorkspace,
) -> Result<BasisPsiDerivativeBundle, BasisError> {
let design_derivatives = build_duchon_design_psi_derivativeswithworkspace(
data,
centers,
spec,
identifiability_transform,
workspace,
)?;
let (native_sources, native_first, native_second) =
build_duchon_native_penalty_psi_derivatives(
centers,
spec,
identifiability_transform,
workspace,
)?;
let (operator_sources, operator_first, operator_second) = if duchon_operator_penalties_requested(
&spec.operator_penalties,
) {
let Some(collocation_points) = operator_collocation_points else {
crate::bail_invalid_basis!(
"Duchon log-kappa operator penalty derivatives require realized collocation points"
);
};
build_duchon_operator_penalty_psi_derivatives(
collocation_points,
centers,
spec,
identifiability_transform,
workspace,
)?
} else {
(Vec::new(), Vec::new(), Vec::new())
};
let mut penalties_derivative = Vec::with_capacity(native_first.len() + operator_first.len());
penalties_derivative.extend(native_first);
penalties_derivative.extend(operator_first);
let mut penaltiessecond_derivative =
Vec::with_capacity(native_second.len() + operator_second.len());
penaltiessecond_derivative.extend(native_second);
penaltiessecond_derivative.extend(operator_second);
let expected_derivative_count = native_sources.len() + operator_sources.len();
if penalties_derivative.len() != expected_derivative_count {
crate::bail_invalid_basis!(
"Duchon penalty derivative count mismatch: assembled {}, expected {} from active penalty sources",
penalties_derivative.len(),
expected_derivative_count
);
}
Ok(BasisPsiDerivativeBundle {
first: BasisPsiDerivativeResult {
design_derivative: design_derivatives.design_first,
penalties_derivative,
implicit_operator: None,
},
second: BasisPsiSecondDerivativeResult {
designsecond_derivative: design_derivatives.design_second_diag,
penaltiessecond_derivative,
implicit_operator: None,
},
implicit_operator: design_derivatives.implicit_operator,
})
}
pub fn build_duchon_basis_log_kappasecond_derivative(
data: ArrayView2<'_, f64>,
spec: &DuchonBasisSpec,
) -> Result<BasisPsiSecondDerivativeResult, BasisError> {
let mut workspace = BasisWorkspace::default();
build_duchon_basis_log_kappasecond_derivativewithworkspace(data, spec, &mut workspace)
}
pub fn build_duchon_basis_log_kappasecond_derivativewithworkspace(
data: ArrayView2<'_, f64>,
spec: &DuchonBasisSpec,
workspace: &mut BasisWorkspace,
) -> Result<BasisPsiSecondDerivativeResult, BasisError> {
let mut bundle = build_duchon_basis_log_kappa_derivativeswithworkspace(data, spec, workspace)?;
bundle.second.implicit_operator = bundle.implicit_operator;
Ok(bundle.second)
}
fn duchon_kernel_amplification(
centers: ArrayView2<'_, f64>,
length_scale: Option<f64>,
p_order: usize,
s_order: usize,
d: usize,
aniso_log_scales: Option<&[f64]>,
coeffs: Option<&DuchonPartialFractionCoeffs>,
pure_poly_coeff: Option<&PolyharmonicBlockCoeff>,
) -> f64 {
let k = centers.nrows();
if k == 0 {
return 1.0;
}
let axis_scales = aniso_log_scales.map(aniso_axis_scales);
let mut max_abs = 0.0_f64;
for i in 0..k {
for j in i..k {
let r = if let Some(scales) = axis_scales.as_deref() {
aniso_distance_rows_with_scales(centers, i, centers, j, scales)
} else {
euclidean_distance_rows(centers, i, centers, j)
};
let val = if let Some(ppc) = pure_poly_coeff {
ppc.eval(r)
} else {
match duchon_matern_kernel_general_from_distance(
r,
length_scale,
p_order,
s_order,
d,
coeffs,
) {
Ok(v) => v,
Err(_) => continue,
}
};
if val.abs() > max_abs {
max_abs = val.abs();
}
}
}
if max_abs > 0.0 && max_abs < 1e-10 {
1.0 / max_abs
} else {
1.0
}
}
pub fn duchon_pure_kernel_amplification(
centers: ArrayView2<'_, f64>,
order: DuchonNullspaceOrder,
power: f64,
) -> f64 {
let dim = centers.ncols();
if dim == 0 || centers.nrows() == 0 {
return 1.0;
}
let effective_order = duchon_effective_nullspace_order(centers, order);
let p_order = duchon_p_from_nullspace_order(effective_order);
let s_order: f64 = power;
let pure_poly_coeff =
PolyharmonicBlockCoeff::new(pure_duchon_block_order(p_order, s_order), dim);
duchon_kernel_amplification(
centers,
None,
p_order,
duchon_power_to_usize(s_order),
dim,
None,
None,
Some(&pure_poly_coeff),
)
}
fn build_duchon_basis_designwithworkspace(
data: ArrayView2<'_, f64>,
centers: ArrayView2<'_, f64>,
length_scale: Option<f64>,
power: f64,
nullspace_order: DuchonNullspaceOrder,
aniso_log_scales: Option<&[f64]>,
workspace: &mut BasisWorkspace,
) -> Result<DuchonBasisDesign, BasisError> {
let n = data.nrows();
let d = data.ncols();
let k = centers.nrows();
if d == 0 {
crate::bail_invalid_basis!("Duchon basis requires at least one covariate dimension");
}
if k == 0 {
crate::bail_invalid_basis!("Duchon basis requires at least one center");
}
if centers.ncols() != d {
crate::bail_dim_basis!(
"Duchon basis dimension mismatch: data has {d} columns, centers have {}",
centers.ncols()
);
}
if data.iter().any(|v| !v.is_finite()) || centers.iter().any(|v| !v.is_finite()) {
crate::bail_invalid_basis!("Duchon basis requires finite data and center values");
}
let nullspace_order = duchon_effective_nullspace_order(centers, nullspace_order);
let p_order = duchon_p_from_nullspace_order(nullspace_order);
let s_order: f64 = power;
let validation_power = if length_scale.is_some() {
duchon_power_to_usize(s_order) as f64
} else {
s_order
};
validate_duchon_kernel_orders(length_scale, p_order, validation_power, d)?;
let poly_block = polynomial_block_from_order(data, nullspace_order);
let z = kernel_constraint_nullspace(centers, nullspace_order, &mut workspace.cache)?;
let coeffs = length_scale.map(|ls| {
duchon_partial_fraction_coeffs(
p_order,
duchon_power_to_usize(s_order),
1.0 / ls.max(1e-300),
)
});
let warn_bounds = match (length_scale, aniso_log_scales) {
(Some(_), Some(eta)) => {
let y_centers = points_in_aniso_y_space(centers, eta);
pairwise_distance_bounds(y_centers.view())
}
(Some(_), None) => pairwise_distance_bounds(centers),
(None, _) => None,
};
if let (Some(length_scale), Some((r_min, r_max))) = (length_scale, warn_bounds) {
let kappa = 1.0 / length_scale.max(1e-300);
let kappa_lo = 1e-2 / r_max;
let kappa_hi = 1e2 / r_min;
if kappa < kappa_lo || kappa > kappa_hi {
log::debug!(
"Duchon κ={} is outside recommended range [{}, {}] derived from centers (r_min={}, r_max={}); numerical conditioning may degrade",
kappa,
kappa_lo,
kappa_hi,
r_min,
r_max
);
}
}
let kernel_cols = z.ncols();
let poly_cols = poly_block.ncols();
let total_cols = kernel_cols + poly_cols;
let pure_poly_coeff = if length_scale.is_none() {
Some(PolyharmonicBlockCoeff::new(
(pure_duchon_block_order(p_order, s_order)) as f64,
d,
))
} else {
None
};
let axis_scales = aniso_log_scales.map(aniso_axis_scales);
let kernel_amp = duchon_kernel_amplification(
centers,
length_scale,
p_order,
duchon_power_to_usize(s_order),
d,
aniso_log_scales,
coeffs.as_ref(),
pure_poly_coeff.as_ref(),
);
let hybrid_kind = match (length_scale, coeffs.as_ref()) {
(Some(ls), Some(c)) if pure_poly_coeff.is_none() => Some(RadialScalarKind::Duchon {
length_scale: ls,
p_order,
s_order: duchon_power_to_usize(s_order),
dim: d,
coeffs: c.clone(),
}),
_ => None,
};
let value_profile = hybrid_kind.as_ref().and_then(|kind| {
if n.saturating_mul(k) < RADIAL_PROFILE_MIN_PAIRS {
return None;
}
let (r_lo, r_hi) = (0..n)
.into_par_iter()
.map(|i| {
let mut lo = f64::INFINITY;
let mut hi = 0.0_f64;
for j in 0..k {
let r = if let Some(scales) = axis_scales.as_deref() {
aniso_distance_rows_with_scales(data, i, centers, j, scales)
} else {
euclidean_distance_rows(data, i, centers, j)
};
if r > 0.0 {
lo = lo.min(r);
hi = hi.max(r);
}
}
(lo, hi)
})
.reduce(
|| (f64::INFINITY, 0.0_f64),
|a, b| (a.0.min(b.0), a.1.max(b.1)),
);
if r_lo.is_finite() && r_hi > r_lo {
radial_profile::RadialProfile::build(kind, r_lo, r_hi)
} else {
None
}
});
let mut basis = Array2::<f64>::zeros((n, total_cols));
let chunk_size = 1024.min(n);
let basis_result: Result<(), BasisError> = basis
.axis_chunks_iter_mut(Axis(0), chunk_size)
.into_par_iter()
.enumerate()
.try_for_each(|(ci, mut chunk)| {
let mut kernel_row = vec![0.0; k];
let chunk_start = ci * chunk_size;
for local_i in 0..chunk.nrows() {
let i = chunk_start + local_i;
for j in 0..k {
let r = if let Some(scales) = axis_scales.as_deref() {
aniso_distance_rows_with_scales(data, i, centers, j, scales)
} else {
euclidean_distance_rows(data, i, centers, j)
};
let raw = if let Some(ref ppc) = pure_poly_coeff {
ppc.eval(r)
} else if let (Some(profile), Some(kind)) =
(value_profile.as_ref(), hybrid_kind.as_ref())
{
profile.eval_or_exact(kind, r)?.0
} else {
duchon_matern_kernel_general_from_distance(
r,
length_scale,
p_order,
duchon_power_to_usize(s_order),
d,
coeffs.as_ref(),
)?
};
kernel_row[j] = raw * kernel_amp;
}
let mut row = chunk.row_mut(local_i);
row.slice_mut(s![..kernel_cols]).fill(0.0);
for j in 0..k {
let kv = kernel_row[j];
if kv != 0.0 {
let z_row = z.row(j);
for col in 0..kernel_cols {
row[col] += kv * z_row[col];
}
}
}
}
Ok(())
});
basis_result?;
if poly_cols > 0 {
basis.slice_mut(s![.., kernel_cols..]).assign(&poly_block);
}
Ok(DuchonBasisDesign { basis })
}
fn build_cyclic_duchon_basis_1dwithworkspace(
data: ArrayView2<'_, f64>,
spec: &DuchonBasisSpec,
start: f64,
end: f64,
) -> Result<BasisBuildResult, BasisError> {
if data.ncols() != 1 {
crate::bail_invalid_basis!("cyclic Duchon smooths currently require exactly one covariate");
}
if end <= start {
return Err(BasisError::InvalidRange(start, end));
}
let centers = select_centers_by_strategy(data, &spec.center_strategy)?;
if centers.ncols() != 1 {
crate::bail_dim_basis!(
"cyclic Duchon centers must have one column, got {}",
centers.ncols()
);
}
let k = centers.nrows();
let s_order_usize = spec.power_as_usize();
if k <= s_order_usize.max(1) {
crate::bail_invalid_basis!(
"cyclic Duchon basis requires more centers ({k}) than power ({})",
spec.power
);
}
let period = end - start;
let p_order = duchon_p_from_nullspace_order(DuchonNullspaceOrder::Zero);
let validation_power = if spec.length_scale.is_some() {
s_order_usize as f64
} else {
spec.power
};
validate_duchon_kernel_orders(spec.length_scale, p_order, validation_power, 1)?;
let coeffs = spec
.length_scale
.map(|ls| duchon_partial_fraction_coeffs(p_order, s_order_usize, 1.0 / ls.max(1e-300)));
let pure_poly_coeff = if spec.length_scale.is_none() {
Some(PolyharmonicBlockCoeff::new(
pure_duchon_block_order(p_order, spec.power),
1,
))
} else {
None
};
let kernel_amp = duchon_kernel_amplification(
centers.view(),
spec.length_scale,
p_order,
duchon_power_to_usize(spec.power),
1,
None,
coeffs.as_ref(),
pure_poly_coeff.as_ref(),
);
let mut basis = Array2::<f64>::zeros((data.nrows(), k + 1));
for i in 0..data.nrows() {
let x = wrap_to_period(data[[i, 0]], start, period);
for j in 0..k {
let c = wrap_to_period(centers[[j, 0]], start, period);
let r = cyclic_distance_1d(x, c, period);
let raw = if let Some(ref ppc) = pure_poly_coeff {
ppc.eval(r)
} else {
duchon_matern_kernel_general_from_distance(
r,
spec.length_scale,
p_order,
s_order_usize,
1,
coeffs.as_ref(),
)?
};
basis[[i, j]] = raw * kernel_amp;
}
basis[[i, k]] = 1.0;
}
let identifiability_transform = spatial_identifiability_transform_from_design(
data,
basis.view(),
&spec.identifiability,
"cyclic Duchon",
)?;
let (design_matrix, transform_for_penalty) = if let Some(z) = identifiability_transform.as_ref()
{
(fast_ab(&basis, z), Some(z))
} else {
(basis, None)
};
let mut s_kernel = Array2::<f64>::zeros((k + 1, k + 1));
let s_cyclic = create_cyclic_difference_penalty_matrix(k, s_order_usize.max(1).min(k - 1))?;
s_kernel.slice_mut(s![..k, ..k]).assign(&s_cyclic);
let s_final = if let Some(z) = transform_for_penalty {
fast_ab(&fast_atb(z, &s_kernel), z)
} else {
s_kernel
};
let candidates = vec![PenaltyCandidate {
matrix: s_final,
nullspace_dim_hint: 1,
source: PenaltySource::Primary,
normalization_scale: 1.0,
kronecker_factors: None,
op: None,
}];
let (penalties, nullspace_dims, penaltyinfo, null_eigenvectors, ops) =
filter_active_penalty_candidates_with_ops(candidates)?;
Ok(BasisBuildResult {
design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(design_matrix)),
penalties,
nullspace_dims,
penaltyinfo,
metadata: BasisMetadata::Duchon {
centers,
length_scale: spec.length_scale,
periodic: Some(vec![Some(period)]),
power: spec.power,
nullspace_order: DuchonNullspaceOrder::Zero,
identifiability_transform,
input_scales: None,
aniso_log_scales: None,
operator_collocation_points: None,
},
kronecker_factored: None,
ops,
null_eigenvectors,
joint_null_rotation: None,
})
}
pub fn build_duchon_basis(
data: ArrayView2<'_, f64>,
spec: &DuchonBasisSpec,
) -> Result<BasisBuildResult, BasisError> {
let mut workspace = BasisWorkspace::default();
build_duchon_basiswithworkspace(data, spec, &mut workspace)
}
pub fn create_duchon_basis_1d_derivative_dense(
t: ArrayView1<'_, f64>,
centers: ArrayView1<'_, f64>,
power: f64,
nullspace_order: DuchonNullspaceOrder,
periodic: bool,
period: Option<f64>,
order: usize,
) -> Result<Array2<f64>, BasisError> {
if order > 2 {
crate::bail_invalid_basis!(
"Duchon basis derivative supports orders 0, 1, and 2; got order={order}"
);
}
if t.is_empty() || centers.is_empty() {
crate::bail_invalid_basis!("Duchon basis derivative requires non-empty t and centers");
}
if t.iter().any(|v| !v.is_finite()) || centers.iter().any(|v| !v.is_finite()) {
crate::bail_invalid_basis!("Duchon basis derivative requires finite t and center values");
}
if !periodic && period.is_some() {
crate::bail_invalid_basis!(
"Duchon basis derivative period is only valid when periodic=true"
);
}
let data = t.to_owned().insert_axis(Axis(1));
let center_matrix = centers.to_owned().insert_axis(Axis(1));
let mut workspace = BasisWorkspace::default();
let user_m = duchon_p_from_nullspace_order(nullspace_order);
let effective_order = if periodic {
DuchonNullspaceOrder::Zero
} else {
duchon_effective_nullspace_order(center_matrix.view(), nullspace_order)
};
let p_order = duchon_p_from_nullspace_order(effective_order);
let s_order = duchon_power_to_usize(power);
validate_duchon_kernel_orders(None, p_order, s_order as f64, 1)?;
if periodic {
let (collapsed_centers, left, resolved_period) =
prepare_periodic_duchon_centers_1d_with_period(center_matrix, period)?;
let z = kernel_constraint_nullspace(
collapsed_centers.view(),
effective_order,
&mut workspace.cache,
)?;
let kernel_cols = z.ncols();
let k_centers = collapsed_centers.nrows();
let centers_col0: Vec<f64> = collapsed_centers.column(0).to_vec();
let mut raw_kernel = Array2::<f64>::zeros((t.len(), k_centers));
for i in 0..t.len() {
let x = wrap_to_period(t[i], left, resolved_period);
for j in 0..k_centers {
let mut delta = (x - centers_col0[j]).rem_euclid(resolved_period);
if delta > 0.5 * resolved_period {
delta -= resolved_period;
}
let r = delta.abs();
let sign = if delta > 0.0 {
1.0
} else if delta < 0.0 {
-1.0
} else {
0.0
};
let (phi, phi_r, phi_rr) =
periodic_duchon_kernel_bernoulli_triplet(r, user_m, resolved_period)?;
raw_kernel[[i, j]] = match order {
0 => phi,
1 => phi_r * sign,
2 => phi_rr,
other => {
crate::bail_invalid_basis!(
"Duchon basis derivative supports orders 0, 1, and 2; got order={other}"
);
}
};
}
}
let mut basis = Array2::<f64>::zeros((t.len(), kernel_cols + 1));
let design_kernel = fast_ab(&raw_kernel, &z);
basis
.slice_mut(s![.., 0..kernel_cols])
.assign(&design_kernel);
if order == 0 {
basis.column_mut(kernel_cols).fill(1.0);
}
return Ok(basis);
}
let z =
kernel_constraint_nullspace(center_matrix.view(), effective_order, &mut workspace.cache)?;
let kernel_cols = z.ncols();
let poly_cols = polynomial_block_from_order(data.view(), effective_order).ncols();
let pure_coeff =
PolyharmonicBlockCoeff::new((pure_duchon_block_order(p_order, s_order as f64)) as f64, 1);
let kernel_amp = duchon_kernel_amplification(
center_matrix.view(),
None,
p_order,
s_order,
1,
None,
None,
Some(&pure_coeff),
);
let mut raw_kernel = Array2::<f64>::zeros((t.len(), centers.len()));
for i in 0..t.len() {
let x = t[i];
for j in 0..centers.len() {
let delta = x - centers[j];
let r = delta.abs();
let sign = if delta > 0.0 {
1.0
} else if delta < 0.0 {
-1.0
} else {
0.0
};
let (phi, phi_r, phi_rr) =
duchon_kernel_radial_triplet(r, None, p_order, s_order as f64, 1, None)?;
raw_kernel[[i, j]] = match order {
0 => phi,
1 => phi_r * sign,
2 => phi_rr,
other => {
crate::bail_invalid_basis!(
"Duchon basis derivative supports orders 0, 1, and 2; got order={other}"
);
}
} * kernel_amp;
}
}
let mut basis = Array2::<f64>::zeros((t.len(), kernel_cols + poly_cols));
let design_kernel = fast_ab(&raw_kernel, &z);
basis
.slice_mut(s![.., 0..kernel_cols])
.assign(&design_kernel);
fill_duchon_1d_polynomial_derivative(&mut basis, kernel_cols, t, effective_order, order);
Ok(basis)
}
struct DuchonRadialJetsNd {
phi_r: Array2<f64>,
phi_rr: Array2<f64>,
phi_rrr: Array2<f64>,
}
fn duchon_effective_length_scale(length_scale: Option<f64>, centers: ArrayView2<'_, f64>) -> f64 {
length_scale.unwrap_or_else(|| {
let n_centers = centers.nrows();
let dim = centers.ncols();
let mut acc = 0.0_f64;
let mut cnt = 0usize;
for i in 0..n_centers.min(8) {
for j in (i + 1)..n_centers.min(8) {
let mut r2 = 0.0_f64;
for a in 0..dim {
let dv = centers[[i, a]] - centers[[j, a]];
r2 += dv * dv;
}
acc += r2.sqrt();
cnt += 1;
}
}
if cnt == 0 || acc <= 0.0 {
1.0
} else {
acc / cnt as f64
}
})
}
fn duchon_radial_jets_nd(
max_order: usize,
t: ArrayView2<'_, f64>,
centers: ArrayView2<'_, f64>,
length_scale: Option<f64>,
nullspace_order: DuchonNullspaceOrder,
power: usize,
caller: &str,
) -> Result<DuchonRadialJetsNd, BasisError> {
let n_rows = t.nrows();
let n_centers = centers.nrows();
let dim = centers.ncols();
if dim == 0 {
crate::bail_invalid_basis!("{caller}: centers must have at least one column");
}
if t.ncols() != dim {
crate::bail_invalid_basis!(
"{caller}: t has {} cols but centers have {}",
t.ncols(),
dim
);
}
assert!(
(1..=3).contains(&max_order),
"duchon_radial_jets_nd supports radial orders 1..=3; got {max_order}"
);
let effective_order = duchon_effective_nullspace_order(centers, nullspace_order);
let p_order = duchon_p_from_nullspace_order(effective_order);
let (jet_p_order, s_order) = match length_scale {
Some(_) => (p_order, power),
None => (p_order + power, 0usize),
};
let kappa = length_scale.map(|l| 1.0 / l.max(1e-300)).unwrap_or(0.0);
let coeffs = duchon_partial_fraction_coeffs(jet_p_order, s_order, kappa);
let effective_length_scale = duchon_effective_length_scale(length_scale, centers);
let mut phi_r = Array2::<f64>::zeros((n_rows, n_centers));
let mut phi_rr = if max_order >= 2 {
Array2::<f64>::zeros((n_rows, n_centers))
} else {
Array2::<f64>::zeros((0, 0))
};
let mut phi_rrr = if max_order >= 3 {
Array2::<f64>::zeros((n_rows, n_centers))
} else {
Array2::<f64>::zeros((0, 0))
};
for n in 0..n_rows {
for k in 0..n_centers {
let mut r2 = 0.0_f64;
for a in 0..dim {
let dv = t[[n, a]] - centers[[k, a]];
r2 += dv * dv;
}
let r = r2.sqrt();
let jets = duchon_radial_jets(
r,
effective_length_scale,
jet_p_order,
s_order,
dim,
&coeffs,
)?;
phi_r[[n, k]] = jets.phi_r;
if max_order >= 2 {
phi_rr[[n, k]] = jets.phi_rr;
}
if max_order >= 3 {
phi_rrr[[n, k]] = jets.phi_rrr;
}
}
}
Ok(DuchonRadialJetsNd {
phi_r,
phi_rr,
phi_rrr,
})
}
pub(crate) fn radial_basis_cartesian_derivative(
order: usize,
t: ArrayView2<'_, f64>,
centers: ArrayView2<'_, f64>,
coeffs: ArrayView2<'_, f64>,
length_scale: Option<f64>,
nullspace_order: DuchonNullspaceOrder,
power: usize,
) -> Result<Array2<f64>, BasisError> {
assert!(
order == 2 || order == 3,
"radial_basis_cartesian_derivative supports Cartesian orders 2 and 3; got {order}"
);
let n_rows = t.nrows();
let n_centers = centers.nrows();
let d = centers.ncols();
let p_out = coeffs.ncols();
assert_eq!(
coeffs.nrows(),
n_centers,
"radial_basis_cartesian_derivative: coeffs has {} rows but centers have {n_centers}",
coeffs.nrows()
);
let jets = duchon_radial_jets_nd(
order,
t,
centers,
length_scale,
nullspace_order,
power,
"radial_basis_cartesian_derivative",
)?;
let d_pow = d.pow(order as u32);
let mut out = Array2::<f64>::zeros((n_rows, p_out * d_pow));
for n in 0..n_rows {
for k in 0..n_centers {
let mut r2 = 0.0_f64;
for a in 0..d {
let delta = t[[n, a]] - centers[[k, a]];
r2 += delta * delta;
}
let r = r2.sqrt();
match order {
2 => {
for a in 0..d {
for c in 0..d {
let basis_hess = if r == 0.0 {
if a == c { jets.phi_rr[[n, k]] } else { 0.0 }
} else {
let inv_r = 1.0 / r;
let u_a = (t[[n, a]] - centers[[k, a]]) * inv_r;
let u_c = (t[[n, c]] - centers[[k, c]]) * inv_r;
let q = jets.phi_r[[n, k]] * inv_r;
let eye = if a == c { 1.0 } else { 0.0 };
q * eye + (jets.phi_rr[[n, k]] - q) * u_a * u_c
};
if basis_hess == 0.0 {
continue;
}
let m = a * d + c;
for i in 0..p_out {
out[[n, i * d_pow + m]] += coeffs[[k, i]] * basis_hess;
}
}
}
}
_ => {
if r == 0.0 {
continue;
}
let inv_r = 1.0 / r;
let q = jets.phi_r[[n, k]] * inv_r;
let b_coef = (jets.phi_rr[[n, k]] - q) * inv_r;
let a_coef = jets.phi_rrr[[n, k]] - 3.0 * b_coef;
for a in 0..d {
let u_a = (t[[n, a]] - centers[[k, a]]) * inv_r;
for c in 0..d {
let u_c = (t[[n, c]] - centers[[k, c]]) * inv_r;
for e in 0..d {
let u_e = (t[[n, e]] - centers[[k, e]]) * inv_r;
let eye_ac = if a == c { 1.0 } else { 0.0 };
let eye_ae = if a == e { 1.0 } else { 0.0 };
let eye_ce = if c == e { 1.0 } else { 0.0 };
let basis_third = a_coef * u_a * u_c * u_e
+ b_coef * (eye_ac * u_e + eye_ae * u_c + eye_ce * u_a);
if basis_third == 0.0 {
continue;
}
let m = (a * d + c) * d + e;
for i in 0..p_out {
out[[n, i * d_pow + m]] += coeffs[[k, i]] * basis_third;
}
}
}
}
}
}
}
}
Ok(out)
}
pub fn duchon_radial_first_derivative_nd(
t: ArrayView2<'_, f64>,
centers: ArrayView2<'_, f64>,
length_scale: Option<f64>,
nullspace_order: DuchonNullspaceOrder,
power: usize,
) -> Result<Array2<f64>, BasisError> {
Ok(duchon_radial_jets_nd(
1,
t,
centers,
length_scale,
nullspace_order,
power,
"duchon_radial_first_derivative_nd",
)?
.phi_r)
}
pub fn duchon_radial_second_derivative_nd(
t: ArrayView2<'_, f64>,
centers: ArrayView2<'_, f64>,
length_scale: Option<f64>,
nullspace_order: DuchonNullspaceOrder,
power: usize,
) -> Result<Array2<f64>, BasisError> {
Ok(duchon_radial_jets_nd(
2,
t,
centers,
length_scale,
nullspace_order,
power,
"duchon_radial_second_derivative_nd",
)?
.phi_rr)
}
pub fn duchon_radial_third_derivative_nd(
t: ArrayView2<'_, f64>,
centers: ArrayView2<'_, f64>,
length_scale: Option<f64>,
nullspace_order: DuchonNullspaceOrder,
power: usize,
) -> Result<Array2<f64>, BasisError> {
Ok(duchon_radial_jets_nd(
3,
t,
centers,
length_scale,
nullspace_order,
power,
"duchon_radial_third_derivative_nd",
)?
.phi_rrr)
}
fn fill_duchon_1d_polynomial_derivative(
basis: &mut Array2<f64>,
start_col: usize,
t: ArrayView1<'_, f64>,
nullspace_order: DuchonNullspaceOrder,
derivative_order: usize,
) {
let exponents: Vec<usize> = match nullspace_order {
DuchonNullspaceOrder::Zero => vec![0],
DuchonNullspaceOrder::Linear => vec![0, 1],
DuchonNullspaceOrder::Degree(degree) => monomial_exponents(1, degree)
.into_iter()
.map(|exponent| exponent[0])
.collect(),
};
for (offset, exponent) in exponents.into_iter().enumerate() {
if exponent < derivative_order {
continue;
}
let coefficient =
(0..derivative_order).fold(1.0, |acc, step| acc * (exponent - step) as f64);
let remaining = exponent - derivative_order;
for row in 0..t.len() {
basis[[row, start_col + offset]] = if remaining == 0 {
coefficient
} else {
coefficient * t[row].powi(remaining as i32)
};
}
}
}
pub fn duchon_polynomial_first_derivative_nd(
t: ArrayView2<'_, f64>,
nullspace_order: DuchonNullspaceOrder,
) -> Array3<f64> {
let n_rows = t.nrows();
let dim = t.ncols();
let max_degree = match nullspace_order {
DuchonNullspaceOrder::Zero => 0usize,
DuchonNullspaceOrder::Linear => 1usize,
DuchonNullspaceOrder::Degree(k) => k,
};
let exponents = monomial_exponents(dim, max_degree);
let n_poly = exponents.len();
let mut out = Array3::<f64>::zeros((n_rows, n_poly, dim));
if dim == 0 || n_poly == 0 {
return out;
}
for (col, alpha) in exponents.iter().enumerate() {
for axis in 0..dim {
let a_axis = alpha[axis];
if a_axis == 0 {
continue;
}
for row in 0..n_rows {
let mut value = a_axis as f64;
for a in 0..dim {
let exp_a = if a == axis { a_axis - 1 } else { alpha[a] };
if exp_a != 0 {
value *= t[[row, a]].powi(exp_a as i32);
}
}
out[[row, col, axis]] = value;
}
}
}
out
}
pub fn radial_input_location_jet_nd(
t: ArrayView2<'_, f64>,
centers: ArrayView2<'_, f64>,
phi_r: ArrayView2<'_, f64>,
) -> Result<Array3<f64>, BasisError> {
let n_rows = t.nrows();
let n_centers = centers.nrows();
let dim = t.ncols();
if phi_r.dim() != (n_rows, n_centers) {
crate::bail_dim_basis!(
"radial_input_location_jet_nd: phi_r shape {:?} != ({n_rows}, {n_centers})",
phi_r.dim()
);
}
if centers.ncols() != dim {
crate::bail_dim_basis!(
"radial_input_location_jet_nd: t has {dim} cols but centers have {}",
centers.ncols()
);
}
let mut out = Array3::<f64>::zeros((n_rows, n_centers, dim));
for n in 0..n_rows {
for k in 0..n_centers {
let mut r2 = 0.0_f64;
for a in 0..dim {
let delta = t[[n, a]] - centers[[k, a]];
r2 += delta * delta;
}
let r = r2.sqrt();
if r <= 1.0e-12 {
continue;
}
let scale = phi_r[[n, k]] / r;
for a in 0..dim {
out[[n, k, a]] = scale * (t[[n, a]] - centers[[k, a]]);
}
}
}
Ok(out)
}
pub fn duchon_sae_atom_basis_with_jet(
t: ArrayView2<'_, f64>,
centers: ArrayView2<'_, f64>,
nullspace_order: DuchonNullspaceOrder,
) -> Result<(Array2<f64>, Array3<f64>), BasisError> {
let dim = centers.ncols();
if dim == 0 {
crate::bail_invalid_basis!(
"duchon_sae_atom_basis_with_jet: centers must have at least one column"
);
}
if t.ncols() != dim {
crate::bail_dim_basis!(
"duchon_sae_atom_basis_with_jet: t has {} cols but centers have {dim}",
t.ncols()
);
}
let effective_order = duchon_effective_nullspace_order(centers, nullspace_order);
let p_order = duchon_p_from_nullspace_order(effective_order);
let s_order: f64 = 0.0;
let poly_block_centers = polynomial_block_from_order(centers, effective_order);
let z = kernel_constraint_nullspace_from_matrix(poly_block_centers.view())?;
let n_kernel = z.ncols();
let pure_poly_coeff =
PolyharmonicBlockCoeff::new(pure_duchon_block_order(p_order, s_order), dim);
let kernel_amp = duchon_kernel_amplification(
centers,
None,
p_order,
duchon_power_to_usize(s_order),
dim,
None,
None,
Some(&pure_poly_coeff),
);
let n_rows = t.nrows();
let n_centers = centers.nrows();
let mut radial_value = Array2::<f64>::zeros((n_rows, n_centers));
for n in 0..n_rows {
for k in 0..n_centers {
let r = euclidean_distance_rows(t, n, centers, k);
radial_value[[n, k]] = pure_poly_coeff.eval(r) * kernel_amp;
}
}
let kernel_design = fast_ab(&radial_value, &z);
let poly_design = polynomial_block_from_order(t, effective_order);
let n_poly = poly_design.ncols();
let mut phi = Array2::<f64>::zeros((n_rows, n_kernel + n_poly));
phi.slice_mut(s![.., ..n_kernel]).assign(&kernel_design);
if n_poly > 0 {
phi.slice_mut(s![.., n_kernel..]).assign(&poly_design);
}
let radial_first = duchon_radial_first_derivative_nd(t, centers, None, effective_order, 0)?;
let radial_jet = radial_input_location_jet_nd(t, centers, radial_first.view())?;
let poly_jet = duchon_polynomial_first_derivative_nd(t, effective_order);
if poly_jet.shape()[1] != n_poly {
crate::bail_dim_basis!(
"duchon_sae_atom_basis_with_jet: polynomial jet has {} columns but design has {n_poly}",
poly_jet.shape()[1]
);
}
let mut jet = Array3::<f64>::zeros((n_rows, n_kernel + n_poly, dim));
for axis in 0..dim {
let projected = radial_jet.index_axis(Axis(2), axis).dot(&z);
let mut block = jet.slice_mut(s![.., ..n_kernel, axis]);
block.assign(&projected);
block *= kernel_amp;
}
jet.slice_mut(s![.., n_kernel.., ..]).assign(&poly_jet);
Ok((phi, jet))
}
pub fn duchon_sae_atom_second_jet(
t: ArrayView2<'_, f64>,
centers: ArrayView2<'_, f64>,
nullspace_order: DuchonNullspaceOrder,
) -> Result<Array4<f64>, BasisError> {
let dim = centers.ncols();
if dim == 0 {
crate::bail_invalid_basis!(
"duchon_sae_atom_second_jet: centers must have at least one column"
);
}
if t.ncols() != dim {
crate::bail_dim_basis!(
"duchon_sae_atom_second_jet: t has {} cols but centers have {dim}",
t.ncols()
);
}
let effective_order = duchon_effective_nullspace_order(centers, nullspace_order);
let p_order = duchon_p_from_nullspace_order(effective_order);
let s_order: f64 = 0.0;
let poly_block_centers = polynomial_block_from_order(centers, effective_order);
let z = kernel_constraint_nullspace_from_matrix(poly_block_centers.view())?;
let n_kernel = z.ncols();
let pure_poly_coeff =
PolyharmonicBlockCoeff::new(pure_duchon_block_order(p_order, s_order), dim);
let kernel_amp = duchon_kernel_amplification(
centers,
None,
p_order,
duchon_power_to_usize(s_order),
dim,
None,
None,
Some(&pure_poly_coeff),
);
let n_rows = t.nrows();
let poly_block_t_cols = polynomial_block_from_order(t, effective_order).ncols();
let coeffs = &z * kernel_amp;
let flat =
radial_basis_cartesian_derivative(2, t, centers, coeffs.view(), None, effective_order, 0)?;
let mut out = Array4::<f64>::zeros((n_rows, n_kernel + poly_block_t_cols, dim, dim));
for n in 0..n_rows {
for i in 0..n_kernel {
for a in 0..dim {
for c in 0..dim {
out[[n, i, a, c]] = flat[[n, i * dim * dim + a * dim + c]];
}
}
}
}
let exponents = monomial_exponents(
dim,
match effective_order {
DuchonNullspaceOrder::Zero => 0,
DuchonNullspaceOrder::Linear => 1,
DuchonNullspaceOrder::Degree(k) => k,
},
);
if exponents.len() != poly_block_t_cols {
crate::bail_dim_basis!(
"duchon_sae_atom_second_jet: monomial count {} != polynomial block columns {poly_block_t_cols}",
exponents.len()
);
}
for (col, alpha) in exponents.iter().enumerate() {
for a in 0..dim {
for c in 0..dim {
for row in 0..n_rows {
let coeff_a = alpha[a];
if coeff_a == 0 {
continue;
}
let coeff_c = if a == c {
alpha[c].saturating_sub(1)
} else {
alpha[c]
};
if a != c && coeff_c == 0 {
continue;
}
let lead = (coeff_a as f64) * (coeff_c as f64);
if lead == 0.0 {
continue;
}
let mut value = lead;
for axis in 0..dim {
let mut exp = alpha[axis];
if axis == a {
exp = exp.saturating_sub(1);
}
if axis == c {
exp = exp.saturating_sub(1);
}
if exp != 0 {
value *= t[[row, axis]].powi(exp as i32);
}
}
out[[row, n_kernel + col, a, c]] = value;
}
}
}
}
Ok(out)
}
pub fn duchon_sae_atom_third_jet(
t: ArrayView2<'_, f64>,
centers: ArrayView2<'_, f64>,
nullspace_order: DuchonNullspaceOrder,
) -> Result<Array5<f64>, BasisError> {
let dim = centers.ncols();
if dim == 0 {
crate::bail_invalid_basis!(
"duchon_sae_atom_third_jet: centers must have at least one column"
);
}
if t.ncols() != dim {
crate::bail_dim_basis!(
"duchon_sae_atom_third_jet: t has {} cols but centers have {dim}",
t.ncols()
);
}
let effective_order = duchon_effective_nullspace_order(centers, nullspace_order);
let p_order = duchon_p_from_nullspace_order(effective_order);
let s_order: f64 = 0.0;
let poly_block_centers = polynomial_block_from_order(centers, effective_order);
let z = kernel_constraint_nullspace_from_matrix(poly_block_centers.view())?;
let n_kernel = z.ncols();
let pure_poly_coeff =
PolyharmonicBlockCoeff::new(pure_duchon_block_order(p_order, s_order), dim);
let kernel_amp = duchon_kernel_amplification(
centers,
None,
p_order,
duchon_power_to_usize(s_order),
dim,
None,
None,
Some(&pure_poly_coeff),
);
let n_rows = t.nrows();
let poly_block_t_cols = polynomial_block_from_order(t, effective_order).ncols();
let coeffs = &z * kernel_amp;
let flat =
radial_basis_cartesian_derivative(3, t, centers, coeffs.view(), None, effective_order, 0)?;
let d_pow = dim * dim * dim;
let mut out = Array5::<f64>::zeros((n_rows, n_kernel + poly_block_t_cols, dim, dim, dim));
for n in 0..n_rows {
for i in 0..n_kernel {
for a in 0..dim {
for c in 0..dim {
for e in 0..dim {
out[[n, i, a, c, e]] = flat[[n, i * d_pow + ((a * dim) + c) * dim + e]];
}
}
}
}
}
let exponents = monomial_exponents(
dim,
match effective_order {
DuchonNullspaceOrder::Zero => 0,
DuchonNullspaceOrder::Linear => 1,
DuchonNullspaceOrder::Degree(k) => k,
},
);
if exponents.len() != poly_block_t_cols {
crate::bail_dim_basis!(
"duchon_sae_atom_third_jet: monomial count {} != polynomial block columns {poly_block_t_cols}",
exponents.len()
);
}
let falling = |alpha: usize, k: usize| -> f64 {
let mut acc = 1.0_f64;
for j in 0..k {
acc *= (alpha as f64) - (j as f64);
}
acc
};
for (col, alpha) in exponents.iter().enumerate() {
for a in 0..dim {
if alpha[a] == 0 {
continue;
}
for c in 0..dim {
for e in 0..dim {
let mut order = vec![0usize; dim];
order[a] += 1;
order[c] += 1;
order[e] += 1;
if (0..dim).any(|axis| order[axis] > alpha[axis]) {
continue;
}
let mut lead = 1.0_f64;
for axis in 0..dim {
lead *= falling(alpha[axis], order[axis]);
}
if lead == 0.0 {
continue;
}
for row in 0..n_rows {
let mut value = lead;
for axis in 0..dim {
let exp = alpha[axis] - order[axis];
if exp != 0 {
value *= t[[row, axis]].powi(exp as i32);
}
}
out[[row, n_kernel + col, a, c, e]] = value;
}
}
}
}
}
Ok(out)
}
pub fn build_duchon_basis_design_and_jets(
t: ArrayView2<'_, f64>,
centers: ArrayView2<'_, f64>,
length_scale: Option<f64>,
power: f64,
nullspace_order: DuchonNullspaceOrder,
periodic_per_axis: &[bool],
periods: &[f64],
) -> Result<(Array2<f64>, Array3<f64>, Array4<f64>), BasisError> {
let dim = centers.ncols();
if dim == 0 {
crate::bail_invalid_basis!(
"build_duchon_basis_design_and_jets: centers must have at least one column"
);
}
if t.ncols() != dim {
crate::bail_dim_basis!(
"build_duchon_basis_design_and_jets: t has {} cols but centers have {dim}",
t.ncols()
);
}
if periodic_per_axis.len() != dim {
crate::bail_invalid_basis!(
"build_duchon_basis_design_and_jets: periodic_per_axis must have length {dim}, got {}",
periodic_per_axis.len()
);
}
if periods.len() != dim {
crate::bail_invalid_basis!(
"build_duchon_basis_design_and_jets: periods must have length {dim}, got {}",
periods.len()
);
}
let any_periodic = periodic_per_axis.iter().any(|&p| p);
let n_rows = t.nrows();
let n_centers = centers.nrows();
if any_periodic {
if length_scale.is_some() {
crate::bail_invalid_basis!(
"mixed-periodicity Duchon basis currently only supports the pure polyharmonic spectrum (length_scale=None)"
);
}
if power != 0.0 {
crate::bail_invalid_basis!(
"mixed-periodicity Duchon basis currently requires power = 0 (pure polyharmonic); got power={power}"
);
}
for (j, (&per, &period)) in periodic_per_axis.iter().zip(periods.iter()).enumerate() {
if per && !(period.is_finite() && period > 0.0) {
crate::bail_invalid_basis!(
"axis {j} is periodic but period={period} is not finite & positive"
);
}
}
}
let effective_order = duchon_effective_nullspace_order(centers, nullspace_order);
let kernel_nullspace_order = if any_periodic {
DuchonNullspaceOrder::Zero
} else {
effective_order
};
let p_order = duchon_p_from_nullspace_order(kernel_nullspace_order);
let kernel_m = if any_periodic {
duchon_p_from_nullspace_order(nullspace_order)
} else {
p_order
};
let s_order_int = duchon_power_to_usize(power);
let s_order_f = power;
let poly_block_centers = polynomial_block_from_order(centers, kernel_nullspace_order);
let z = kernel_constraint_nullspace_from_matrix(poly_block_centers.view())?;
let n_kernel = z.ncols();
let coeffs = length_scale
.map(|ls| duchon_partial_fraction_coeffs(p_order, s_order_int, 1.0 / ls.max(1e-300)));
let pure_poly_coeff = if length_scale.is_none() {
Some(PolyharmonicBlockCoeff::new(
pure_duchon_block_order(kernel_m, s_order_f),
dim,
))
} else {
None
};
let kernel_amp = if any_periodic {
1.0
} else {
duchon_kernel_amplification(
centers,
length_scale,
p_order,
s_order_int,
dim,
None,
coeffs.as_ref(),
pure_poly_coeff.as_ref(),
)
};
let mut radial_value = Array2::<f64>::zeros((n_rows, n_centers));
let mut radial_first = Array2::<f64>::zeros((n_rows, n_centers));
let mut radial_second = Array2::<f64>::zeros((n_rows, n_centers));
let mut delta = Array3::<f64>::zeros((n_rows, n_centers, dim));
let mut metric_d1 = Array3::<f64>::zeros((n_rows, n_centers, dim));
let mut metric_d2 = Array3::<f64>::zeros((n_rows, n_centers, dim));
let m_pure = pure_duchon_block_order(kernel_m, s_order_f);
let pi = std::f64::consts::PI;
for n in 0..n_rows {
for k in 0..n_centers {
let mut r2 = 0.0_f64;
for a in 0..dim {
let raw = t[[n, a]] - centers[[k, a]];
let (d_a, d1_a, d2_a) = if periodic_per_axis[a] {
let p = periods[a];
let theta = pi * raw / p;
((p / pi) * theta.sin(), theta.cos(), -(pi / p) * theta.sin())
} else {
(raw, 1.0, 0.0)
};
delta[[n, k, a]] = d_a;
metric_d1[[n, k, a]] = d1_a;
metric_d2[[n, k, a]] = d2_a;
r2 += d_a * d_a;
}
let r = r2.sqrt();
let (phi, phi_r, phi_rr) = if pure_poly_coeff.is_some() {
polyharmonic_kernel_triplet(r, m_pure, dim)?
} else {
let jets = duchon_radial_jets(
r,
length_scale.expect("hybrid Duchon requires length_scale"),
p_order,
s_order_int,
dim,
coeffs.as_ref().expect("hybrid Duchon requires coeffs"),
)?;
(jets.phi, jets.phi_r, jets.phi_rr)
};
radial_value[[n, k]] = phi * kernel_amp;
radial_first[[n, k]] = phi_r;
radial_second[[n, k]] = phi_rr;
}
}
let kernel_design = fast_ab(&radial_value, &z);
let poly_design = polynomial_block_from_order(t, kernel_nullspace_order);
let n_poly = poly_design.ncols();
let mut phi_design = Array2::<f64>::zeros((n_rows, n_kernel + n_poly));
phi_design
.slice_mut(s![.., ..n_kernel])
.assign(&kernel_design);
if n_poly > 0 {
phi_design
.slice_mut(s![.., n_kernel..])
.assign(&poly_design);
}
let mut radial_jet = Array3::<f64>::zeros((n_rows, n_centers, dim));
let mut radial_hess = Array4::<f64>::zeros((n_rows, n_centers, dim, dim));
for n in 0..n_rows {
for k in 0..n_centers {
let mut r2 = 0.0_f64;
for a in 0..dim {
let da = delta[[n, k, a]];
r2 += da * da;
}
let r = r2.sqrt();
let phi_r = radial_first[[n, k]];
let phi_rr = radial_second[[n, k]];
if r <= 1.0e-12 {
for a in 0..dim {
let d1a = metric_d1[[n, k, a]];
radial_hess[[n, k, a, a]] = phi_rr * kernel_amp * d1a * d1a;
}
continue;
}
let q = phi_r / r;
let s_scalar = (phi_rr - q) / r2;
for a in 0..dim {
let ga = delta[[n, k, a]] * metric_d1[[n, k, a]];
radial_jet[[n, k, a]] = q * ga * kernel_amp;
}
for a in 0..dim {
let ga = delta[[n, k, a]] * metric_d1[[n, k, a]];
for c in 0..dim {
let gc = delta[[n, k, c]] * metric_d1[[n, k, c]];
let mut value = s_scalar * ga * gc;
if a == c {
let d1a = metric_d1[[n, k, a]];
let curvature = d1a * d1a + delta[[n, k, a]] * metric_d2[[n, k, a]];
value += q * curvature;
}
radial_hess[[n, k, a, c]] = value * kernel_amp;
}
}
}
}
let poly_jet = duchon_polynomial_first_derivative_nd(t, kernel_nullspace_order);
if poly_jet.shape()[1] != n_poly {
crate::bail_dim_basis!(
"build_duchon_basis_design_and_jets: polynomial jet has {} columns but design has {n_poly}",
poly_jet.shape()[1]
);
}
let mut jet = Array3::<f64>::zeros((n_rows, n_kernel + n_poly, dim));
for axis in 0..dim {
let projected = radial_jet.index_axis(Axis(2), axis).dot(&z);
jet.slice_mut(s![.., ..n_kernel, axis]).assign(&projected);
}
if n_poly > 0 {
jet.slice_mut(s![.., n_kernel.., ..]).assign(&poly_jet);
}
let mut hess = Array4::<f64>::zeros((n_rows, n_kernel + n_poly, dim, dim));
for a in 0..dim {
for c in 0..dim {
let slab = radial_hess.slice(s![.., .., a, c]);
let projected = slab.dot(&z);
hess.slice_mut(s![.., ..n_kernel, a, c]).assign(&projected);
}
}
if n_poly > 0 {
let max_degree = match kernel_nullspace_order {
DuchonNullspaceOrder::Zero => 0,
DuchonNullspaceOrder::Linear => 1,
DuchonNullspaceOrder::Degree(deg) => deg,
};
let exponents = monomial_exponents(dim, max_degree);
if exponents.len() != n_poly {
crate::bail_dim_basis!(
"build_duchon_basis_design_and_jets: monomial count {} != polynomial block columns {n_poly}",
exponents.len()
);
}
for (col, alpha) in exponents.iter().enumerate() {
for a in 0..dim {
let coeff_a = alpha[a];
if coeff_a == 0 {
continue;
}
for c in 0..dim {
let coeff_c = if a == c {
alpha[c].saturating_sub(1)
} else {
alpha[c]
};
if a != c && coeff_c == 0 {
continue;
}
let lead = (coeff_a as f64) * (coeff_c as f64);
if lead == 0.0 {
continue;
}
for row in 0..n_rows {
let mut value = lead;
for axis in 0..dim {
let mut exp = alpha[axis];
if axis == a {
exp = exp.saturating_sub(1);
}
if axis == c {
exp = exp.saturating_sub(1);
}
if exp != 0 {
value *= t[[row, axis]].powi(exp as i32);
}
}
hess[[row, n_kernel + col, a, c]] = value;
}
}
}
}
}
Ok((phi_design, jet, hess))
}
pub fn matern_radial_first_derivative_nd(
t: ArrayView2<'_, f64>,
centers: ArrayView2<'_, f64>,
length_scale: f64,
nu: MaternNu,
) -> Result<Array2<f64>, BasisError> {
let n_rows = t.nrows();
let n_centers = centers.nrows();
let dim = centers.ncols();
if dim == 0 {
crate::bail_invalid_basis!(
"matern_radial_first_derivative_nd: centers must have at least one column".into(),
);
}
if t.ncols() != dim {
crate::bail_invalid_basis!(
"matern_radial_first_derivative_nd: t has {} cols but centers have {}",
t.ncols(),
dim
);
}
let mut out = Array2::<f64>::zeros((n_rows, n_centers));
for n in 0..n_rows {
for k in 0..n_centers {
let mut r2 = 0.0_f64;
for a in 0..dim {
let dv = t[[n, a]] - centers[[k, a]];
r2 += dv * dv;
}
let r = r2.sqrt();
let (_phi, phi_r, _phi_rr, _ratio) =
matern_kernel_radial_tripletwith_safe_ratio(r, length_scale, nu)?;
out[[n, k]] = phi_r;
}
}
Ok(out)
}
pub fn matern_radial_second_derivative_nd(
t: ArrayView2<'_, f64>,
centers: ArrayView2<'_, f64>,
length_scale: f64,
nu: MaternNu,
) -> Result<Array2<f64>, BasisError> {
let n_rows = t.nrows();
let n_centers = centers.nrows();
let dim = centers.ncols();
if dim == 0 {
crate::bail_invalid_basis!(
"matern_radial_second_derivative_nd: centers must have at least one column".into(),
);
}
if t.ncols() != dim {
crate::bail_invalid_basis!(
"matern_radial_second_derivative_nd: t has {} cols but centers have {}",
t.ncols(),
dim
);
}
let mut out = Array2::<f64>::zeros((n_rows, n_centers));
for n in 0..n_rows {
for k in 0..n_centers {
let mut r2 = 0.0_f64;
for a in 0..dim {
let dv = t[[n, a]] - centers[[k, a]];
r2 += dv * dv;
}
let r = r2.sqrt();
let (_phi, _phi_r, phi_rr, _ratio) =
matern_kernel_radial_tripletwith_safe_ratio(r, length_scale, nu)?;
out[[n, k]] = phi_rr;
}
}
Ok(out)
}
fn matern_metric_weights(dim: usize, aniso: Option<&[f64]>) -> Vec<f64> {
match centered_aniso_contrasts(aniso) {
Some(psi) => psi.iter().map(|&v| (2.0 * v).exp()).collect(),
None => vec![1.0; dim],
}
}
pub fn matern_input_location_jet_nd(
t: ArrayView2<'_, f64>,
centers: ArrayView2<'_, f64>,
length_scale: f64,
nu: MaternNu,
aniso_log_scales: Option<&[f64]>,
) -> Result<Array3<f64>, BasisError> {
let n_rows = t.nrows();
let n_centers = centers.nrows();
let dim = centers.ncols();
if dim == 0 {
crate::bail_invalid_basis!(
"matern_input_location_jet_nd: centers must have at least one column".into(),
);
}
if t.ncols() != dim {
crate::bail_invalid_basis!(
"matern_input_location_jet_nd: t has {} cols but centers have {}",
t.ncols(),
dim
);
}
if let Some(eta) = aniso_log_scales
&& eta.len() != dim
{
crate::bail_invalid_basis!(
"matern_input_location_jet_nd: aniso_log_scales len {} != dim {}",
eta.len(),
dim
);
}
let weights = matern_metric_weights(dim, aniso_log_scales);
let mut out = Array3::<f64>::zeros((n_rows, n_centers, dim));
for n in 0..n_rows {
for k in 0..n_centers {
let mut r2 = 0.0_f64;
for a in 0..dim {
let h = t[[n, a]] - centers[[k, a]];
r2 += weights[a] * h * h;
}
let r = r2.sqrt();
if r <= 0.0 {
continue;
}
let (_phi, phi_r, _phi_rr, _ratio) =
matern_kernel_radial_tripletwith_safe_ratio(r, length_scale, nu)?;
let scale = phi_r / r;
for a in 0..dim {
let h = t[[n, a]] - centers[[k, a]];
out[[n, k, a]] = scale * weights[a] * h;
}
}
}
Ok(out)
}
pub fn matern_input_location_hessian_nd(
t: ArrayView2<'_, f64>,
centers: ArrayView2<'_, f64>,
length_scale: f64,
nu: MaternNu,
aniso_log_scales: Option<&[f64]>,
) -> Result<Array4<f64>, BasisError> {
let n_rows = t.nrows();
let n_centers = centers.nrows();
let dim = centers.ncols();
if dim == 0 {
crate::bail_invalid_basis!(
"matern_input_location_hessian_nd: centers must have at least one column".into(),
);
}
if t.ncols() != dim {
crate::bail_invalid_basis!(
"matern_input_location_hessian_nd: t has {} cols but centers have {}",
t.ncols(),
dim
);
}
if let Some(eta) = aniso_log_scales
&& eta.len() != dim
{
crate::bail_invalid_basis!(
"matern_input_location_hessian_nd: aniso_log_scales len {} != dim {}",
eta.len(),
dim
);
}
let weights = matern_metric_weights(dim, aniso_log_scales);
let mut out = Array4::<f64>::zeros((n_rows, n_centers, dim, dim));
for n in 0..n_rows {
for k in 0..n_centers {
let mut r2 = 0.0_f64;
for a in 0..dim {
let h = t[[n, a]] - centers[[k, a]];
r2 += weights[a] * h * h;
}
let r = r2.sqrt();
let (_phi, phi_r, phi_rr, phi_r_over_r) =
matern_kernel_radial_tripletwith_safe_ratio(r, length_scale, nu)?;
if r <= 0.0 {
for a in 0..dim {
out[[n, k, a, a]] = phi_r_over_r * weights[a];
}
continue;
}
let q = phi_r / r; let inv_r2 = 1.0 / r2;
for a in 0..dim {
let ga = weights[a] * (t[[n, a]] - centers[[k, a]]); for c in 0..dim {
let gc = weights[c] * (t[[n, c]] - centers[[k, c]]);
let mut value = phi_rr * (ga / r) * (gc / r);
value -= q * ga * gc * inv_r2;
if a == c {
value += q * weights[a];
}
out[[n, k, a, c]] = value;
}
}
}
}
Ok(out)
}
pub fn sphere_first_derivative_nd(
points: ArrayView2<'_, f64>,
centers: ArrayView2<'_, f64>,
penalty_order: usize,
project_to_tangent: bool,
) -> Result<Array3<f64>, BasisError> {
let n_rows = points.nrows();
let n_centers = centers.nrows();
let dim = points.ncols();
if !(1..=4).contains(&penalty_order) {
crate::bail_invalid_basis!(
"sphere_first_derivative_nd: penalty_order must be in 1..=4; got {penalty_order}"
);
}
if dim == 0 {
crate::bail_invalid_basis!(
"sphere_first_derivative_nd: points must have at least one column".into(),
);
}
if centers.ncols() != dim {
crate::bail_invalid_basis!(
"sphere_first_derivative_nd: points have dim {} but centers have dim {}",
dim,
centers.ncols()
);
}
let tangent_projector =
project_to_tangent.then_some(crate::terms::latent_coord::LatentManifold::Sphere { dim });
let mut out = Array3::<f64>::zeros((n_rows, n_centers, dim));
let mut ambient = Array1::<f64>::zeros(dim);
for n in 0..n_rows {
for k in 0..n_centers {
let mut cos_g = 0.0_f64;
for a in 0..dim {
cos_g += points[[n, a]] * centers[[k, a]];
}
let dk = wahba_sphere_kernel_sobolev_derivative_dcos(cos_g, penalty_order);
for a in 0..dim {
ambient[a] = dk * centers[[k, a]];
}
if let Some(manifold) = &tangent_projector {
let tangent = manifold.project_to_tangent(points.row(n), ambient.view());
for a in 0..dim {
out[[n, k, a]] = tangent[a];
}
} else {
for a in 0..dim {
out[[n, k, a]] = ambient[a];
}
}
}
}
Ok(out)
}
fn spherical_wahba_kernel_jet_with_kind(
data: ArrayView2<'_, f64>,
centers: ArrayView2<'_, f64>,
penalty_order: usize,
radians: bool,
kernel: SphereWahbaKernel,
) -> Result<Array3<f64>, BasisError> {
validate_lat_lon_matrix(data, "spherical spline jet data", radians)?;
validate_lat_lon_matrix(centers, "spherical spline jet centers", radians)?;
if !(1..=4).contains(&penalty_order) {
crate::bail_invalid_basis!(
"spherical spline jet penalty_order must be one of 1, 2, 3, 4; got {penalty_order}"
);
}
let n = data.nrows();
let k = centers.nrows();
let deg = if radians {
1.0
} else {
std::f64::consts::PI / 180.0
};
let mut sin_lat_c = Vec::<f64>::with_capacity(k);
let mut cos_lat_c = Vec::<f64>::with_capacity(k);
let mut sin_lon_c = Vec::<f64>::with_capacity(k);
let mut cos_lon_c = Vec::<f64>::with_capacity(k);
for c in centers.outer_iter() {
let (s_lat, c_lat) = (c[0] * deg).sin_cos();
let (s_lon, c_lon) = (c[1] * deg).sin_cos();
sin_lat_c.push(s_lat);
cos_lat_c.push(c_lat);
sin_lon_c.push(s_lon);
cos_lon_c.push(c_lon);
}
let mut out = Array3::<f64>::zeros((n, k, 2));
let err_flag = std::sync::atomic::AtomicBool::new(false);
out.axis_chunks_iter_mut(ndarray::Axis(0), 256)
.into_par_iter()
.enumerate()
.for_each(|(chunk_idx, mut block)| {
let row_offset = chunk_idx * 256;
for (local_i, mut out_row) in block.outer_iter_mut().enumerate() {
let i = row_offset + local_i;
let (sin_lat, cos_lat) = (data[(i, 0)] * deg).sin_cos();
let (sin_lon, cos_lon) = (data[(i, 1)] * deg).sin_cos();
for j in 0..k {
let dlon_cos = cos_lon * cos_lon_c[j] + sin_lon * sin_lon_c[j];
let dlon_sin = sin_lon * cos_lon_c[j] - cos_lon * sin_lon_c[j];
let cos_gamma = sin_lat * sin_lat_c[j] + cos_lat * cos_lat_c[j] * dlon_cos;
let dk =
wahba_sphere_kernel_derivative_dcos_kind(cos_gamma, penalty_order, kernel);
let dcos_dphi = cos_lat * sin_lat_c[j] - sin_lat * cos_lat_c[j] * dlon_cos;
let dcos_dpsi = -cos_lat * cos_lat_c[j] * dlon_sin;
let dphi = dk * dcos_dphi * deg;
let dpsi = dk * dcos_dpsi * deg;
if !dphi.is_finite() || !dpsi.is_finite() {
err_flag.store(true, std::sync::atomic::Ordering::Relaxed);
return;
}
out_row[[j, 0]] = dphi;
out_row[[j, 1]] = dpsi;
}
}
});
if err_flag.load(std::sync::atomic::Ordering::Relaxed) {
crate::bail_invalid_basis!("spherical spline kernel jet produced a non-finite value");
}
Ok(out)
}
fn apply_identifiability_to_jet(raw_jet: &Array3<f64>, z: &Array2<f64>) -> Array3<f64> {
let n = raw_jet.shape()[0];
let k = raw_jet.shape()[1];
let kp = z.ncols();
assert_eq!(
z.nrows(),
k,
"apply_identifiability_to_jet: identifiability transform rows ({}) must match raw jet basis dim ({})",
z.nrows(),
k
);
let mut out = Array3::<f64>::zeros((n, kp, 2));
for axis in 0..2 {
let raw_axis: ndarray::ArrayView2<'_, f64> = raw_jet.index_axis(ndarray::Axis(2), axis);
let projected = raw_axis.dot(z);
out.slice_mut(ndarray::s![.., .., axis]).assign(&projected);
}
out
}
fn spherical_harmonic_jet(
data: ArrayView2<'_, f64>,
max_degree: usize,
radians: bool,
) -> Result<Array3<f64>, BasisError> {
validate_lat_lon_matrix(data, "spherical-harmonic jet", radians)?;
if max_degree < 1 {
crate::bail_invalid_basis!("spherical-harmonic jet max_degree must be >= 1");
}
if max_degree > 32 {
crate::bail_invalid_basis!(
"spherical-harmonic jet max_degree {max_degree} too large; cap is 32"
);
}
let n = data.nrows();
let p = max_degree * (max_degree + 2);
let deg = if radians {
1.0
} else {
std::f64::consts::PI / 180.0
};
let norms = precompute_harmonic_norms(max_degree);
let l_cap = max_degree + 1;
let mut out = Array3::<f64>::zeros((n, p, 2));
let idx = |l: usize, m: usize| l * l_cap + m;
{
let mut row_blocks = out
.axis_chunks_iter_mut(ndarray::Axis(0), 1024)
.collect::<Vec<_>>();
let chunk_size = 1024usize;
row_blocks
.par_iter_mut()
.enumerate()
.for_each(|(chunk_idx, block)| {
let mut p_buf = vec![0.0_f64; l_cap * l_cap];
let row_offset = chunk_idx * chunk_size;
for (local_i, mut out_row) in block.outer_iter_mut().enumerate() {
let i = row_offset + local_i;
let lat_raw = data[(i, 0)] * deg;
let lat =
lat_raw.clamp(-std::f64::consts::FRAC_PI_2, std::f64::consts::FRAC_PI_2);
let lon = data[(i, 1)] * deg;
let cos_lat = lat.cos();
let x = lat.sin();
let somx2 = (1.0 - x * x).max(0.0).sqrt();
let one_minus_x2 = (1.0 - x * x).max(f64::EPSILON);
for slot in p_buf.iter_mut() {
*slot = 0.0;
}
p_buf[idx(0, 0)] = 1.0;
for m in 1..=max_degree {
p_buf[idx(m, m)] = -((2 * m - 1) as f64) * somx2 * p_buf[idx(m - 1, m - 1)];
}
for m in 0..max_degree {
p_buf[idx(m + 1, m)] = ((2 * m + 1) as f64) * x * p_buf[idx(m, m)];
}
for m in 0..=max_degree {
for l in (m + 2)..=max_degree {
p_buf[idx(l, m)] = (((2 * l - 1) as f64) * x * p_buf[idx(l - 1, m)]
- ((l + m - 1) as f64) * p_buf[idx(l - 2, m)])
/ ((l - m) as f64);
}
}
let dp = |l: usize, m: usize| -> f64 {
let p_lm1 = if l >= 1 { p_buf[idx(l - 1, m)] } else { 0.0 };
(-(l as f64) * x * p_buf[idx(l, m)] + ((l + m) as f64) * p_lm1)
/ one_minus_x2
};
let (sin1, cos1) = lon.sin_cos();
let mut sin_buf = [0.0_f64; 33];
let mut cos_buf = [0.0_f64; 33];
sin_buf[0] = 0.0;
cos_buf[0] = 1.0;
if max_degree >= 1 {
sin_buf[1] = sin1;
cos_buf[1] = cos1;
}
let two_cos1 = 2.0 * cos1;
for m in 2..=max_degree {
sin_buf[m] = two_cos1 * sin_buf[m - 1] - sin_buf[m - 2];
cos_buf[m] = two_cos1 * cos_buf[m - 1] - cos_buf[m - 2];
}
let mut col = 0usize;
for l in 1..=max_degree {
for m_pos in (1..=l).rev() {
let nlm = norms[idx(l, m_pos)];
let mf = m_pos as f64;
out_row[[col, 0]] = nlm * sin_buf[m_pos] * dp(l, m_pos) * cos_lat * deg;
out_row[[col, 1]] =
nlm * mf * cos_buf[m_pos] * p_buf[idx(l, m_pos)] * deg;
col += 1;
}
let nl0 = norms[idx(l, 0)];
out_row[[col, 0]] = nl0 * dp(l, 0) * cos_lat * deg;
out_row[[col, 1]] = 0.0;
col += 1;
for m in 1..=l {
let nlm = norms[idx(l, m)];
let mf = m as f64;
out_row[[col, 0]] = nlm * cos_buf[m] * dp(l, m) * cos_lat * deg;
out_row[[col, 1]] = -nlm * mf * sin_buf[m] * p_buf[idx(l, m)] * deg;
col += 1;
}
}
}
});
}
if out.iter().any(|v| !v.is_finite()) {
crate::bail_invalid_basis!("spherical-harmonic jet produced a non-finite value");
}
Ok(out)
}
pub fn spherical_spline_design_jet(
data: ArrayView2<'_, f64>,
spec: &SphericalSplineBasisSpec,
) -> Result<Array3<f64>, BasisError> {
if matches!(spec.method, SphereMethod::Harmonic) {
let l_max = spec
.max_degree
.unwrap_or_else(|| default_spherical_harmonic_degree(data.nrows()));
if !(1..=4).contains(&spec.penalty_order) {
crate::bail_invalid_basis!(
"spherical-harmonic jet penalty_order must be one of 1, 2, 3, 4; got {}",
spec.penalty_order
);
}
return spherical_harmonic_jet(data, l_max, spec.radians);
}
validate_lat_lon_matrix(data, "spherical spline jet", spec.radians)?;
let centers = match realized_center_strategy(&spec.center_strategy) {
CenterStrategy::FarthestPoint { num_centers } => {
select_spherical_farthest_point_centers(data, *num_centers, spec.radians)?
}
_ => select_centers_by_strategy(data, &spec.center_strategy)?,
};
validate_lat_lon_matrix(centers.view(), "spherical spline jet centers", spec.radians)?;
if centers.nrows() < 2 {
return Err(BasisError::InsufficientColumnsForConstraint {
found: centers.nrows(),
});
}
let z = match &spec.identifiability {
SphericalSplineIdentifiability::FrozenTransform { transform } => {
if transform.nrows() != centers.nrows() {
crate::bail_dim_basis!(
"frozen spherical identifiability transform mismatch: {} centers but transform has {} rows",
centers.nrows(),
transform.nrows()
);
}
transform.clone()
}
SphericalSplineIdentifiability::CenterSumToZero => {
let weights = sphere_area_weights(centers.view(), spec.radians);
weighted_coefficient_sum_to_zero_transform(weights.view())?
}
};
let raw_jet = spherical_wahba_kernel_jet_with_kind(
data,
centers.view(),
spec.penalty_order,
spec.radians,
spec.wahba_kernel,
)?;
Ok(apply_identifiability_to_jet(&raw_jet, &z))
}
pub fn periodic_bspline_first_derivative_nd(
t: ArrayView2<'_, f64>,
data_range: (f64, f64),
degree: usize,
num_basis: usize,
) -> Result<Array3<f64>, BasisError> {
if t.ncols() != 1 {
crate::bail_invalid_basis!(
"periodic_bspline_first_derivative_nd: t must have exactly 1 column; got {}",
t.ncols()
);
}
if degree == 0 {
crate::bail_invalid_basis!("periodic_bspline_first_derivative_nd requires degree >= 1");
}
if num_basis < degree + 1 {
crate::bail_invalid_basis!(
"periodic_bspline_first_derivative_nd requires num_basis >= degree + 1 (got num_basis={num_basis}, degree={degree})"
);
}
let (start, end) = data_range;
if !(start.is_finite() && end.is_finite()) || end <= start {
crate::bail_invalid_basis!(
"periodic_bspline_first_derivative_nd: data_range must be finite and ordered, got {data_range:?}"
);
}
let period = end - start;
let n_rows = t.nrows();
let t_col = t.column(0);
let mut phi = vec![0.0_f64; num_basis];
let mut dphi = vec![0.0_f64; num_basis];
let mut out = Array3::<f64>::zeros((n_rows, num_basis, 1));
for row in 0..n_rows {
let xi = t_col[row];
if !xi.is_finite() {
crate::bail_invalid_basis!(
"periodic_bspline_first_derivative_nd: non-finite latent at row {row}"
);
}
let rowsum =
fill_periodic_bspline_unnormalized_value_row(xi, start, period, degree, &mut phi);
if !rowsum.is_finite() || rowsum <= 0.0 {
crate::bail_invalid_basis!(
"periodic_bspline_first_derivative_nd: non-positive rowsum at row {row}: {rowsum}"
);
}
let rowsum_derivative =
fill_periodic_bspline_unnormalized_derivative_row(xi, start, period, degree, &mut dphi);
if !rowsum_derivative.is_finite() {
crate::bail_invalid_basis!(
"periodic_bspline_first_derivative_nd: non-finite rowsum derivative at row {row}: {rowsum_derivative}"
);
}
let rowsum_squared = rowsum * rowsum;
for i in 0..num_basis {
out[[row, i, 0]] = dphi[i] / rowsum - phi[i] * rowsum_derivative / rowsum_squared;
}
}
Ok(out)
}
pub fn bspline_tensor_first_derivative(
t: ArrayView2<'_, f64>,
knots_per_axis: &[ArrayView1<'_, f64>],
degrees: &[usize],
) -> Result<Array3<f64>, BasisError> {
let n_axes = t.ncols();
if knots_per_axis.len() != n_axes || degrees.len() != n_axes {
crate::bail_invalid_basis!(
"bspline_tensor_first_derivative: t has {n_axes} axes but received \
{} knot vectors and {} degrees",
knots_per_axis.len(),
degrees.len(),
);
}
if n_axes == 0 {
crate::bail_invalid_basis!(
"bspline_tensor_first_derivative: t must have at least one axis".into(),
);
}
let n_rows = t.nrows();
let mut k_per_axis = Vec::<usize>::with_capacity(n_axes);
let mut total = 1usize;
for a in 0..n_axes {
let k = knots_per_axis[a]
.len()
.checked_sub(degrees[a] + 1)
.ok_or_else(|| {
BasisError::InvalidInput(format!(
"bspline_tensor_first_derivative: axis {a} knot vector too short \
for degree {}",
degrees[a]
))
})?;
k_per_axis.push(k);
total = total.checked_mul(k).ok_or_else(|| {
BasisError::InvalidInput(
"bspline_tensor_first_derivative: tensor-product basis size overflow".into(),
)
})?;
}
let mut out = Array3::<f64>::zeros((n_rows, total, n_axes));
let mut values_per_axis: Vec<Vec<f64>> = k_per_axis.iter().map(|&k| vec![0.0; k]).collect();
let mut derivs_per_axis: Vec<Vec<f64>> = k_per_axis.iter().map(|&k| vec![0.0; k]).collect();
let mut value_scratch_per_axis: Vec<internal::BsplineScratch> = degrees
.iter()
.map(|&d| internal::BsplineScratch::new(d))
.collect();
let mut lower_basis_per_axis: Vec<Vec<f64>> = knots_per_axis
.iter()
.zip(degrees.iter())
.map(|(knots, &d)| vec![0.0; knots.len().saturating_sub(d)])
.collect();
let mut lower_scratch_per_axis: Vec<internal::BsplineScratch> = degrees
.iter()
.map(|&d| internal::BsplineScratch::new(d.saturating_sub(1)))
.collect();
let mut idx = vec![0usize; n_axes];
let mut prefix = vec![1.0; n_axes + 1];
let mut suffix = vec![1.0; n_axes + 1];
for n in 0..n_rows {
for a in 0..n_axes {
internal::evaluate_splines_at_point_into(
t[[n, a]],
degrees[a],
knots_per_axis[a],
&mut values_per_axis[a],
&mut value_scratch_per_axis[a],
);
evaluate_bspline_derivative_scalar_into(
t[[n, a]],
knots_per_axis[a],
degrees[a],
&mut derivs_per_axis[a],
&mut lower_basis_per_axis[a],
&mut lower_scratch_per_axis[a],
)?;
}
for k in 0..total {
let mut rem = k;
for a in (0..n_axes).rev() {
idx[a] = rem % k_per_axis[a];
rem /= k_per_axis[a];
}
prefix[0] = 1.0;
for a in 0..n_axes {
prefix[a + 1] = prefix[a] * values_per_axis[a][idx[a]];
}
suffix[n_axes] = 1.0;
for a in (0..n_axes).rev() {
suffix[a] = suffix[a + 1] * values_per_axis[a][idx[a]];
}
for axis in 0..n_axes {
let leave_one_out = prefix[axis] * suffix[axis + 1];
out[[n, k, axis]] = derivs_per_axis[axis][idx[axis]] * leave_one_out;
}
}
}
Ok(out)
}
#[inline]
fn periodic_distance_1d(x: f64, c: f64, period: f64) -> f64 {
let dx = (x - c).rem_euclid(period).abs();
dx.min(period - dx).abs()
}
fn even_bernoulli_polynomial(degree: usize, t: f64) -> Result<f64, BasisError> {
let t2 = t * t;
match degree {
2 => Ok(t2 - t + 1.0 / 6.0),
4 => Ok(t2 * t2 - 2.0 * t2 * t + t2 - 1.0 / 30.0),
6 => {
let t4 = t2 * t2;
let t6 = t4 * t2;
Ok(t6 - 3.0 * t4 * t + 2.5 * t4 - 0.5 * t2 + 1.0 / 42.0)
}
8 => {
let t4 = t2 * t2;
let t6 = t4 * t2;
let t8 = t4 * t4;
Ok(
t8 - 4.0 * t6 * t + (14.0 / 3.0) * t6 - (7.0 / 3.0) * t4 + (2.0 / 3.0) * t2
- 1.0 / 30.0,
)
}
other => Err(BasisError::InvalidInput(format!(
"periodic Duchon Bernoulli kernel only implemented for B_{{2m}} with m ∈ {{1, 2, 3, 4}}; got degree {other}"
))),
}
}
fn periodic_duchon_kernel_bernoulli(r: f64, m: usize, period: f64) -> Result<f64, BasisError> {
if !period.is_finite() || period <= 0.0 {
crate::bail_invalid_basis!(
"periodic Duchon kernel requires positive finite period; got {period}"
);
}
if m == 0 {
crate::bail_invalid_basis!("periodic Duchon order m must be at least 1");
}
let t = (r / period).rem_euclid(1.0);
let sign = if m % 2 == 1 { 1.0 } else { -1.0 };
Ok(sign * even_bernoulli_polynomial(2 * m, t)?)
}
fn even_bernoulli_polynomial_derivatives(degree: usize, s: f64) -> Result<(f64, f64), BasisError> {
let s2 = s * s;
match degree {
2 => Ok((2.0 * s - 1.0, 2.0)),
4 => {
let d1 = 4.0 * s2 * s - 6.0 * s2 + 2.0 * s;
let d2 = 12.0 * s2 - 12.0 * s + 2.0;
Ok((d1, d2))
}
6 => {
let s3 = s2 * s;
let s4 = s2 * s2;
let s5 = s4 * s;
let d1 = 6.0 * s5 - 15.0 * s4 + 10.0 * s3 - s;
let d2 = 30.0 * s4 - 60.0 * s3 + 30.0 * s2 - 1.0;
Ok((d1, d2))
}
8 => {
let s3 = s2 * s;
let s4 = s2 * s2;
let s5 = s4 * s;
let s6 = s4 * s2;
let s7 = s6 * s;
let d1 = 8.0 * s7 - 28.0 * s6 + 28.0 * s5 - (28.0 / 3.0) * s3 + (4.0 / 3.0) * s;
let d2 = 56.0 * s6 - 168.0 * s5 + 140.0 * s4 - 28.0 * s2 + 4.0 / 3.0;
Ok((d1, d2))
}
other => Err(BasisError::InvalidInput(format!(
"periodic Duchon Bernoulli kernel derivative only implemented for B_{{2m}} with m ∈ {{1, 2, 3, 4}}; got degree {other}"
))),
}
}
fn periodic_duchon_kernel_bernoulli_triplet(
r: f64,
m: usize,
period: f64,
) -> Result<(f64, f64, f64), BasisError> {
if !period.is_finite() || period <= 0.0 {
crate::bail_invalid_basis!(
"periodic Duchon kernel requires positive finite period; got {period}"
);
}
if m == 0 {
crate::bail_invalid_basis!("periodic Duchon order m must be at least 1");
}
let s = (r / period).rem_euclid(1.0);
let sign = if m % 2 == 1 { 1.0 } else { -1.0 };
let phi = sign * even_bernoulli_polynomial(2 * m, s)?;
let (b1, b2) = even_bernoulli_polynomial_derivatives(2 * m, s)?;
let dphi_dr = sign * b1 / period;
let d2phi_dr2 = sign * b2 / (period * period);
Ok((phi, dphi_dr, d2phi_dr2))
}
fn collapse_periodic_endpoint(centers: Array2<f64>, left: f64, period: f64) -> Array2<f64> {
if period <= 0.0 || !period.is_finite() {
return centers;
}
let tol = period.max(1.0) * 1.0e-10;
let col = centers.column(0);
let n_rows = col.len();
let mut seen_left = false;
let keep: Vec<usize> = (0..n_rows)
.filter(|&i| {
if periodic_distance_1d(col[i], left, period) <= tol {
if seen_left {
return false;
}
seen_left = true;
}
true
})
.collect();
if keep.len() == n_rows {
return centers;
}
let mut trimmed = Array2::<f64>::zeros((keep.len(), centers.ncols()));
for (out_row, &src_row) in keep.iter().enumerate() {
for c in 0..centers.ncols() {
trimmed[[out_row, c]] = centers[[src_row, c]];
}
}
trimmed
}
fn build_periodic_duchon_basis_1d(
data: ArrayView2<'_, f64>,
spec: &DuchonBasisSpec,
centers: Array2<f64>,
workspace: &mut BasisWorkspace,
) -> Result<BasisBuildResult, BasisError> {
if data.ncols() != 1 {
crate::bail_invalid_basis!(
"periodic Duchon smooths currently require exactly one covariate"
);
}
let explicit_period = spec
.periodic
.as_ref()
.and_then(|axes| axes.first().copied().flatten());
let (centers, left, period) =
prepare_periodic_duchon_centers_1d_with_period(centers, explicit_period)?;
let user_m = duchon_p_from_nullspace_order(spec.nullspace_order);
let effective_nullspace_order = DuchonNullspaceOrder::Zero;
let p_order = duchon_p_from_nullspace_order(effective_nullspace_order);
let s_order = spec.power_as_usize();
validate_duchon_kernel_orders(spec.length_scale, p_order, s_order as f64, 1)?;
let z = kernel_constraint_nullspace(
centers.view(),
effective_nullspace_order,
&mut workspace.cache,
)?;
let kernel_cols = z.ncols();
let mut basis = Array2::<f64>::zeros((data.nrows(), kernel_cols + 1));
let coeffs = spec
.length_scale
.map(|ls| duchon_partial_fraction_coeffs(p_order, s_order, 1.0 / ls.max(1e-300)));
let pure_poly_coeff = if spec.length_scale.is_none() {
Some(PolyharmonicBlockCoeff::new(
(pure_duchon_block_order(p_order, s_order as f64)) as f64,
1,
))
} else {
None
};
let kernel_amp = duchon_kernel_amplification(
centers.view(),
spec.length_scale,
p_order,
s_order,
1,
None,
coeffs.as_ref(),
pure_poly_coeff.as_ref(),
);
let centers_col0: Vec<f64> = centers.column(0).to_vec();
let n_data = data.nrows();
let k_centers = centers_col0.len();
let len_scale = spec.length_scale;
let mut raw_kernel = Array2::<f64>::zeros((n_data, k_centers));
let err_flag = std::sync::atomic::AtomicBool::new(false);
let amp = kernel_amp;
if pure_poly_coeff.is_some() {
raw_kernel
.axis_chunks_iter_mut(ndarray::Axis(0), 1024)
.into_par_iter()
.enumerate()
.for_each(|(chunk_idx, mut block)| {
let row_offset = chunk_idx * 1024;
for (local_i, mut out_row) in block.outer_iter_mut().enumerate() {
let i = row_offset + local_i;
let x = wrap_to_period(data[[i, 0]], left, period);
for j in 0..k_centers {
let r = periodic_distance_1d(x, centers_col0[j], period);
match periodic_duchon_kernel_bernoulli(r, user_m, period) {
Ok(v) => out_row[j] = v,
Err(_) => {
err_flag.store(true, std::sync::atomic::Ordering::Relaxed);
return;
}
}
}
}
});
} else {
raw_kernel
.axis_chunks_iter_mut(ndarray::Axis(0), 1024)
.into_par_iter()
.enumerate()
.for_each(|(chunk_idx, mut block)| {
let row_offset = chunk_idx * 1024;
for (local_i, mut out_row) in block.outer_iter_mut().enumerate() {
let i = row_offset + local_i;
let x = wrap_to_period(data[[i, 0]], left, period);
for j in 0..k_centers {
let r = periodic_distance_1d(x, centers_col0[j], period);
match duchon_matern_kernel_general_from_distance(
r,
len_scale,
p_order,
s_order,
1,
coeffs.as_ref(),
) {
Ok(v) => out_row[j] = v * amp,
Err(_) => {
err_flag.store(true, std::sync::atomic::Ordering::Relaxed);
return;
}
}
}
}
});
}
if err_flag.load(std::sync::atomic::Ordering::Relaxed) {
crate::bail_invalid_basis!("periodic Duchon kernel evaluation produced a non-finite value");
}
let design_kernel = fast_ab(&raw_kernel, &z);
basis
.slice_mut(s![.., 0..kernel_cols])
.assign(&design_kernel);
basis.column_mut(kernel_cols).fill(1.0);
let mut center_kernel = Array2::<f64>::zeros((centers.nrows(), centers.nrows()));
fill_symmetric_from_row_kernel(&mut center_kernel, |i, j| {
let r = periodic_distance_1d(centers[[i, 0]], centers[[j, 0]], period);
if pure_poly_coeff.is_some() {
periodic_duchon_kernel_bernoulli(r, user_m, period)
} else {
Ok(duchon_matern_kernel_general_from_distance(
r,
spec.length_scale,
p_order,
s_order,
1,
coeffs.as_ref(),
)? * kernel_amp)
}
})?;
let omega = fast_ab(&fast_atb(&z, ¢er_kernel), &z);
let mut penalty = Array2::<f64>::zeros((basis.ncols(), basis.ncols()));
penalty
.slice_mut(s![0..kernel_cols, 0..kernel_cols])
.assign(&omega);
let base_design = DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(basis));
let identifiability_transform = spatial_identifiability_transform_from_design_matrix(
data,
&base_design,
&spec.identifiability,
"periodic Duchon",
)?;
let (design, primary) = if let Some(transform) = identifiability_transform.as_ref() {
let design = wrap_dense_design_with_transform(base_design, transform, "periodic Duchon")?;
let transformed = fast_ab(&fast_atb(transform, &penalty), transform);
(design, transformed)
} else {
(base_design, penalty)
};
let candidates = vec![normalize_penalty_candidate(
primary,
1,
PenaltySource::Primary,
)];
let (penalties, nullspace_dims, penaltyinfo, null_eigenvectors, ops) =
filter_active_penalty_candidates_with_ops(candidates)?;
Ok(BasisBuildResult {
design,
penalties,
nullspace_dims,
penaltyinfo,
ops,
null_eigenvectors,
joint_null_rotation: None,
metadata: BasisMetadata::Duchon {
centers,
length_scale: spec.length_scale,
periodic: Some(vec![Some(period)]),
power: spec.power,
nullspace_order: effective_nullspace_order,
identifiability_transform,
input_scales: None,
aniso_log_scales: None,
operator_collocation_points: None,
},
kronecker_factored: None,
})
}
#[inline]
fn duchon_mixed_periodicity_distance(
x: ArrayView1<'_, f64>,
y: ArrayView1<'_, f64>,
periodic_per_axis: &[bool],
periods: &[f64],
) -> f64 {
let d = x.len();
assert_eq!(d, y.len());
assert_eq!(d, periodic_per_axis.len());
assert_eq!(d, periods.len());
let mut acc = 0.0_f64;
for j in 0..d {
let delta = if periodic_per_axis[j] {
let p = periods[j];
(p / std::f64::consts::PI) * (std::f64::consts::PI * (x[j] - y[j]) / p).sin()
} else {
x[j] - y[j]
};
acc += delta * delta;
}
acc.sqrt()
}
fn build_duchon_basis_mixed_periodicity(
data: ArrayView2<'_, f64>,
spec: &DuchonBasisSpec,
centers: Array2<f64>,
periodic_per_axis: &[bool],
periods: &[f64],
workspace: &mut BasisWorkspace,
) -> Result<BasisBuildResult, BasisError> {
let d = data.ncols();
if d == 0 {
crate::bail_invalid_basis!("Duchon basis requires at least one covariate dimension");
}
if periodic_per_axis.len() != d {
crate::bail_invalid_basis!(
"periodic_per_axis must have length d={d}, got {}",
periodic_per_axis.len()
);
}
if periods.len() != d {
crate::bail_invalid_basis!("periods must have length d={d}, got {}", periods.len());
}
for (j, (&per, &period)) in periodic_per_axis.iter().zip(periods.iter()).enumerate() {
if per && !(period.is_finite() && period > 0.0) {
crate::bail_invalid_basis!(
"axis {j} is periodic but period={period} is not finite & positive"
);
}
}
if centers.ncols() != d {
crate::bail_invalid_basis!(
"centers ncols={} does not match data ncols={d}",
centers.ncols()
);
}
if spec.length_scale.is_some() {
crate::bail_invalid_basis!(
"mixed-periodicity Duchon basis currently only supports the pure polyharmonic spectrum (length_scale=None)"
);
}
if spec.power != 0.0 {
crate::bail_invalid_basis!(
"mixed-periodicity Duchon basis currently requires power = 0 (pure polyharmonic); got power={}",
spec.power
);
}
let user_m = duchon_p_from_nullspace_order(spec.nullspace_order);
let effective_nullspace_order = DuchonNullspaceOrder::Zero;
let p_order = duchon_p_from_nullspace_order(effective_nullspace_order);
let s_order_int = 0usize;
validate_duchon_kernel_orders(None, p_order, s_order_int as f64, d)?;
let z = kernel_constraint_nullspace(
centers.view(),
effective_nullspace_order,
&mut workspace.cache,
)?;
let kernel_cols = z.ncols();
let m_kernel = pure_duchon_block_order(user_m, s_order_int as f64);
let ppc = PolyharmonicBlockCoeff::new(m_kernel, d);
let centers_owned = centers.clone();
let k_centers = centers_owned.nrows();
let n_data = data.nrows();
let mut raw_kernel = Array2::<f64>::zeros((n_data, k_centers));
raw_kernel
.axis_chunks_iter_mut(ndarray::Axis(0), 1024)
.into_par_iter()
.enumerate()
.for_each(|(chunk_idx, mut block)| {
let row_offset = chunk_idx * 1024;
for (local_i, mut out_row) in block.outer_iter_mut().enumerate() {
let i = row_offset + local_i;
let x_row = data.row(i);
for j in 0..k_centers {
let c_row = centers_owned.row(j);
let r =
duchon_mixed_periodicity_distance(x_row, c_row, periodic_per_axis, periods);
out_row[j] = ppc.eval(r);
}
}
});
let design_kernel = fast_ab(&raw_kernel, &z);
let mut basis = Array2::<f64>::zeros((n_data, kernel_cols + 1));
basis
.slice_mut(s![.., 0..kernel_cols])
.assign(&design_kernel);
basis.column_mut(kernel_cols).fill(1.0);
let mut center_kernel = Array2::<f64>::zeros((k_centers, k_centers));
fill_symmetric_from_row_kernel(&mut center_kernel, |i, j| {
let r = duchon_mixed_periodicity_distance(
centers_owned.row(i),
centers_owned.row(j),
periodic_per_axis,
periods,
);
Ok(ppc.eval(r))
})?;
let omega = fast_ab(&fast_atb(&z, ¢er_kernel), &z);
let mut penalty = Array2::<f64>::zeros((basis.ncols(), basis.ncols()));
penalty
.slice_mut(s![0..kernel_cols, 0..kernel_cols])
.assign(&omega);
let base_design = DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(basis));
let identifiability_transform = spatial_identifiability_transform_from_design_matrix(
data,
&base_design,
&spec.identifiability,
"mixed-periodicity Duchon",
)?;
let (design, primary) = if let Some(transform) = identifiability_transform.as_ref() {
let design =
wrap_dense_design_with_transform(base_design, transform, "mixed-periodicity Duchon")?;
let transformed = fast_ab(&fast_atb(transform, &penalty), transform);
(design, transformed)
} else {
(base_design, penalty)
};
let candidates = vec![normalize_penalty_candidate(
primary,
1,
PenaltySource::Primary,
)];
let (penalties, nullspace_dims, penaltyinfo, null_eigenvectors, ops) =
filter_active_penalty_candidates_with_ops(candidates)?;
Ok(BasisBuildResult {
design,
penalties,
nullspace_dims,
penaltyinfo,
ops,
null_eigenvectors,
joint_null_rotation: None,
metadata: BasisMetadata::Duchon {
centers: centers_owned,
length_scale: None,
periodic: Some(
periodic_per_axis
.iter()
.zip(periods.iter())
.map(|(&is_periodic, &period)| if is_periodic { Some(period) } else { None })
.collect(),
),
power: spec.power,
nullspace_order: effective_nullspace_order,
identifiability_transform,
input_scales: None,
aniso_log_scales: None,
operator_collocation_points: None,
},
kronecker_factored: None,
})
}
pub fn build_duchon_basis_mixed_periodicity_auto(
data: ArrayView2<'_, f64>,
spec: &DuchonBasisSpec,
periodic_per_axis: &[bool],
periods: Option<&[f64]>,
) -> Result<BasisBuildResult, BasisError> {
let mut workspace = BasisWorkspace::default();
let centers = select_centers_by_strategy(data, &spec.center_strategy)?;
assert_spatial_centers_below_large_scale_cap(data.ncols(), centers.view())?;
let d = data.ncols();
if periodic_per_axis.len() != d {
crate::bail_invalid_basis!(
"periodic_per_axis must have length d={d}, got {}",
periodic_per_axis.len()
);
}
let resolved_periods: Vec<f64> = match periods {
Some(p) => {
if p.len() != d {
crate::bail_invalid_basis!("periods must have length d={d}, got {}", p.len());
}
p.to_vec()
}
None => {
let mut out = vec![1.0_f64; d];
for j in 0..d {
if periodic_per_axis[j] {
let col = centers.column(j);
let left = col.iter().fold(f64::INFINITY, |a, &b| a.min(b));
let right = col.iter().fold(f64::NEG_INFINITY, |a, &b| a.max(b));
if !left.is_finite() || !right.is_finite() || left >= right {
return Err(BasisError::InvalidRange(left, right));
}
out[j] = right - left;
}
}
out
}
};
if d == 1 && periodic_per_axis[0] {
let mut periodic_spec = spec.clone();
periodic_spec.periodic = Some(vec![Some(resolved_periods[0])]);
return build_periodic_duchon_basis_1d(data, &periodic_spec, centers, &mut workspace);
}
build_duchon_basis_mixed_periodicity(
data,
spec,
centers,
periodic_per_axis,
&resolved_periods,
&mut workspace,
)
}
pub fn duchon_cubic_default(dim: usize) -> (DuchonNullspaceOrder, f64) {
(DuchonNullspaceOrder::Linear, (dim as f64 - 1.0) / 2.0)
}
fn duchon_native_penalty_candidates(
centers: ArrayView2<'_, f64>,
length_scale: Option<f64>,
power: f64,
nullspace_order: DuchonNullspaceOrder,
aniso_log_scales: Option<&[f64]>,
kernel_transform: &Array2<f64>,
outer_identifiability: Option<&Array2<f64>>,
poly_cols: usize,
) -> Result<Vec<PenaltyCandidate>, BasisError> {
let dim = centers.ncols();
if dim == 0 {
crate::bail_invalid_basis!(
"duchon_native_penalty_candidates: centers must have at least one column"
);
}
let k = centers.nrows();
let z = kernel_transform;
let n_kernel = z.ncols();
let p_order = duchon_p_from_nullspace_order(nullspace_order);
let s_int = duchon_power_to_usize(power);
let pure = length_scale.is_none();
let pure_poly_coeff = if pure {
Some(PolyharmonicBlockCoeff::new(
pure_duchon_block_order(p_order, power),
dim,
))
} else {
None
};
let coeffs =
length_scale.map(|ls| duchon_partial_fraction_coeffs(p_order, s_int, 1.0 / ls.max(1e-300)));
let kernel_amp = duchon_kernel_amplification(
centers,
length_scale,
p_order,
s_int,
dim,
aniso_log_scales,
coeffs.as_ref(),
pure_poly_coeff.as_ref(),
);
let axis_scales = aniso_log_scales.map(aniso_axis_scales);
let mut center_kernel = Array2::<f64>::zeros((k, k));
fill_symmetric_from_row_kernel(&mut center_kernel, |i, j| {
let r = if let Some(scales) = axis_scales.as_deref() {
aniso_distance_rows_with_scales(centers, i, centers, j, scales)
} else {
euclidean_distance_rows(centers, i, centers, j)
};
if let Some(ppc) = pure_poly_coeff.as_ref() {
Ok(ppc.eval(r))
} else {
duchon_matern_kernel_general_from_distance(
r,
length_scale,
p_order,
s_int,
dim,
coeffs.as_ref(),
)
}
})?;
let amp2 = kernel_amp * kernel_amp;
let omega = {
let zt_k = fast_atb(z, ¢er_kernel);
fast_ab(&zt_k, z).mapv(|v| v * amp2)
};
let n_pre = n_kernel + poly_cols;
let mut primary_pre = Array2::<f64>::zeros((n_pre, n_pre));
primary_pre
.slice_mut(s![..n_kernel, ..n_kernel])
.assign(&omega);
let primary = symmetrize(&project_penalty_matrix(&primary_pre, outer_identifiability));
let shrink = if poly_cols > 1 {
let mut shrink_pre = Array2::<f64>::zeros((n_pre, n_pre));
for col in (n_kernel + 1)..n_pre {
shrink_pre[[col, col]] = 1.0;
}
let shrink = symmetrize(&project_penalty_matrix(&shrink_pre, outer_identifiability));
Some(shrink)
} else {
None
};
let mut out = Vec::new();
out.push(normalize_penalty_candidate(
primary,
0,
PenaltySource::Primary,
));
if let Some(shrink) = shrink {
out.push(normalize_penalty_candidate(
shrink,
0,
PenaltySource::DoublePenaltyNullspace,
));
}
Ok(out)
}
const DUCHON_COLLOCATION_OVERSAMPLE: usize = 3;
fn duchon_operator_penalty_candidates(
collocation_points: ArrayView2<'_, f64>,
centers: ArrayView2<'_, f64>,
operator_penalties: &DuchonOperatorPenaltySpec,
length_scale: Option<f64>,
power: f64,
nullspace_order: DuchonNullspaceOrder,
per_axis_relevance: bool,
identifiability_transform: Option<&Array2<f64>>,
workspace: &mut BasisWorkspace,
) -> Result<Vec<PenaltyCandidate>, BasisError> {
let want_mass = matches!(operator_penalties.mass, OperatorPenaltySpec::Active { .. });
let mut want_tension = matches!(
operator_penalties.tension,
OperatorPenaltySpec::Active { .. }
);
let mut want_stiffness = matches!(
operator_penalties.stiffness,
OperatorPenaltySpec::Active { .. }
);
let effective_order = duchon_effective_nullspace_order(centers, nullspace_order);
let p_order = duchon_p_from_nullspace_order(effective_order);
let dim = centers.ncols();
let two_pps = 2.0 * (p_order as f64 + power);
want_tension = want_tension && two_pps > dim as f64 + 1.0;
want_stiffness = want_stiffness && two_pps > dim as f64 + 2.0;
if !want_mass && !want_tension && !want_stiffness {
return Ok(Vec::new());
}
let mut effective_spec = operator_penalties.clone();
if !want_tension {
effective_spec.tension = OperatorPenaltySpec::Disabled;
}
if !want_stiffness {
effective_spec.stiffness = OperatorPenaltySpec::Disabled;
}
let max_op = duchon_max_active_operator_derivative_order(&effective_spec);
let ops = build_duchon_collocation_operator_matriceswithworkspace(
centers,
collocation_points,
None,
length_scale,
power,
nullspace_order,
None,
identifiability_transform.map(|t| t.view()),
max_op,
workspace,
)?;
let kernel_nullspace = ops.kernel_nullspace_transform.as_ref();
let poly_cols = ops.polynomial_block_cols;
let split_tension = per_axis_relevance && want_tension;
let factory_spec = if split_tension {
let mut spec = effective_spec.clone();
spec.tension = OperatorPenaltySpec::Disabled;
spec
} else {
effective_spec
};
let mut candidates = if let Some(length_scale) = length_scale {
operator_penalty_candidates_closed_form(
centers,
&ops.d0,
&ops.d1,
&ops.d2,
&factory_spec,
p_order,
duchon_power_to_usize(power),
length_scale,
None,
kernel_nullspace,
poly_cols,
identifiability_transform,
)
} else {
operator_penalty_candidates_closed_form_pure(
centers,
&ops.d0,
&ops.d1,
&ops.d2,
&factory_spec,
p_order,
power,
None,
kernel_nullspace,
poly_cols,
identifiability_transform,
)
};
if split_tension {
for axis in 0..dim {
let d1_axis = ops.d1.slice(s![axis..; dim, ..]).to_owned();
candidates.push(normalize_penalty_candidate(
symmetrize(&fast_ata(&d1_axis)),
0,
PenaltySource::OperatorRelevance { axis },
));
}
}
Ok(candidates)
}