///|
/// Feature matrix with reproducible column metadata.
pub struct FeatureSet {
  values : Array[Array[Double]]
  names : Array[String]
  row_count : Int
  column_count : Int
  fingerprint : UInt64
}

///|
/// Integer-category encoding map.
pub struct CategoryEncoding {
  categories : Array[Int]
  values : Array[Double]
  unknown_value : Double
  ordered : Bool
}

///|
/// Feature screening result.
pub struct FeatureScreening {
  selected : Array[Int]
  dropped : Array[Int]
  variance : Array[Double]
  outcome_correlation : Array[Double]
  reasons : Array[String]
}

///|
/// Train-validation split with explicit row indices.
pub struct FeatureSplit {
  train_indices : Array[Int]
  validation_indices : Array[Int]
  train_values : Array[Array[Double]]
  validation_values : Array[Array[Double]]
}

///|
/// Monotonic transformation specification.
pub struct TransformSpec {
  center : Double
  scale : Double
  power : Double
  lower : Double
  upper : Double
  missing_value : Double
}

///|
fn fe_unique_int(values : Array[Int]) -> Array[Int] {
  let result : Array[Int] = Array::new()
  for value in values {
    if !result.contains(value) {
      result.push(value)
    }
  }
  result
}

///|
fn fe_sorted_int(values : Array[Int]) -> Array[Int] {
  let result = values.copy()
  for i in 1.. 0 && result[cursor - 1] > value {
      result[cursor] = result[cursor - 1]
      cursor -= 1
    }
    result[cursor] = value
  }
  result
}

///|
fn fe_column(matrix : Array[Array[Double]], index : Int) -> Array[Double] {
  let result = Array::new(capacity=matrix.length())
  for row in matrix {
    if index < row.length() {
      result.push(row[index])
    }
  }
  result
}

///|
/// Creates a feature set and computes a stable fingerprint.
pub fn feature_set(
  values : Array[Array[Double]],
  names? : Array[String] = [],
) -> FeatureSet {
  let row_count = values.length()
  let column_count = if row_count == 0 { 0 } else { values[0].length() }
  let actual_names : Array[String] = Array::new(capacity=column_count)
  if names.length() == column_count {
    for name in names {
      actual_names.push(name)
    }
  } else {
    for i in 0.. FeatureSet {
  if matrix.length() == 0 {
    return feature_set([])
  }
  let width = matrix[0].length()
  let result : Array[Array[Double]] = Array::new(capacity=matrix.length())
  for row in matrix {
    let expanded = row.copy()
    for value in row {
      expanded.push(if is_finite(value) { 0.0 } else { 1.0 })
    }
    result.push(expanded)
  }
  let names : Array[String] = Array::new(capacity=width * 2)
  for i in 0.. FeatureSet {
  let actual_levels = if levels.length() > 0 {
    fe_sorted_int(fe_unique_int(levels))
  } else {
    fe_sorted_int(fe_unique_int(categories))
  }
  let result : Array[Array[Double]] = Array::new(capacity=categories.length())
  for category in categories {
    let row = Array::make(actual_levels.length(), 0.0)
    for j in 0.. Array[Double] {
  let levels = fe_unique_int(categories)
  let counts = Array::make(levels.length(), 0)
  for category in categories {
    for j in 0.. CategoryEncoding {
  let n = categories.length().min(values.length())
  let levels = fe_unique_int(categories[:n].to_owned())
  let encoded = Array::make(levels.length(), unknown_value)
  let counts = Array::make(levels.length(), 0)
  for i in 0.. 0 {
      encoded[j] /= counts[j].to_double()
    }
  }
  { categories: levels, values: encoded, unknown_value, ordered }
}

///|
/// Applies a category encoding map to new category values.
pub fn apply_category_encoding(
  encoding : CategoryEncoding,
  categories : Array[Int],
) -> Array[Double] {
  let result = Array::new(capacity=categories.length())
  for category in categories {
    let mut encoded = encoding.unknown_value
    for j in 0.. Array[Double] {
  let n = categories.length().min(outcome.length())
  let global = mean_or(outcome[:n].to_owned(), 0.0)
  let levels = fe_unique_int(categories[:n].to_owned())
  let sums = Array::make(levels.length(), 0.0)
  let counts = Array::make(levels.length(), 0)
  for i in 0.. Array[Int] {
  let width = if bins > 0 { bins } else { 1 }
  let observed = values.filter(fn(value) { is_finite(value) })
  let sorted = observed.copy()
  for i in 1.. 0 && sorted[cursor - 1] > value {
      sorted[cursor] = sorted[cursor - 1]
      cursor -= 1
    }
    sorted[cursor] = value
  }
  let cutpoints = Array::new(capacity=width - 1)
  for bin in 1..= cutpoint {
          bucket += 1
        }
      }
      result.push(bucket.min(width - 1))
    }
  }
  result
}

///|
/// Computes a piecewise-linear spline basis with knots.
pub fn linear_spline_basis(
  values : Array[Double],
  knots : Array[Double],
) -> FeatureSet {
  let result : Array[Array[Double]] = Array::new(capacity=values.length())
  for value in values {
    let row = Array::new(capacity=knots.length() + 2)
    row.push(value)
    row.push(value * value)
    for knot in knots {
      row.push(if value > knot { value - knot } else { 0.0 })
    }
    result.push(row)
  }
  let names : Array[String] = ["linear", "quadratic"]
  for i in 0.. FeatureSet {
  let order = if degree > 0 { degree } else { 1 }
  let result : Array[Array[Double]] = Array::new(capacity=values.length())
  for value in values {
    let row = Array::new(capacity=order)
    let mut power = 1.0
    for _ in 1..<=order {
      power *= value
      row.push(power)
    }
    result.push(row)
  }
  let names : Array[String] = Array::new(capacity=order)
  for i in 0.. FeatureSet {
  if matrix.length() == 0 {
    return feature_set([])
  }
  let width = matrix[0].length()
  let result : Array[Array[Double]] = Array::new(capacity=matrix.length())
  for row in matrix {
    let expanded : Array[Double] = Array::new(
      capacity=if include_original {
        width + width * (width - 1) / 2
      } else {
        width * (width - 1) / 2
      },
    )
    if include_original {
      for value in row {
        expanded.push(value)
      }
    }
    for i in 0.. Array[Double] {
  let result = Array::new(capacity=values.length())
  for value in values {
    if value >= 0.0 {
      result.push(@math.ln(1.0 + value))
    } else {
      result.push(-@math.ln(1.0 + value.abs()))
    }
  }
  result
}

///|
/// Applies a signed square-root transform.
pub fn signed_sqrt_transform(values : Array[Double]) -> Array[Double] {
  let result = Array::new(capacity=values.length())
  for value in values {
    result.push(if value >= 0.0 { value.sqrt() } else { -value.abs().sqrt() })
  }
  result
}

///|
/// Fits a robust center-scale-power transformation specification.
pub fn transform_spec(
  values : Array[Double],
  power? : Double = 1.0,
  lower? : Double = -1.0e12,
  upper? : Double = 1.0e12,
  missing_value? : Double = 0.0,
) -> TransformSpec {
  let observed = values.filter(fn(value) { is_finite(value) })
  let center = mean_or(observed, 0.0)
  let scale = std_dev(observed).max(1.0e-12)
  { center, scale, power, lower, upper, missing_value }
}

///|
/// Applies a configured center-scale-power transformation.
pub fn apply_transform_spec(
  spec : TransformSpec,
  values : Array[Double],
) -> Array[Double] {
  let result = Array::new(capacity=values.length())
  for value in values {
    if !is_finite(value) {
      result.push(spec.missing_value)
    } else {
      let bounded = clamp(value, spec.lower, spec.upper)
      let standardized = (bounded - spec.center) / spec.scale
      result.push(
        if spec.power == 1.0 {
          standardized
        } else if standardized >= 0.0 {
          @math.pow(standardized, spec.power)
        } else {
          -@math.pow(standardized.abs(), spec.power)
        },
      )
    }
  }
  result
}

///|
/// Clips every matrix column to supplied quantile limits.
pub fn clip_matrix(
  matrix : Array[Array[Double]],
  lower_probability? : Double = 0.01,
  upper_probability? : Double = 0.99,
) -> Array[Array[Double]] {
  if matrix.length() == 0 {
    return []
  }
  let width = matrix[0].length()
  let lower = Array::new(capacity=width)
  let upper = Array::new(capacity=width)
  for j in 0.. Array[Array[Double]] {
  if matrix.length() == 0 {
    return []
  }
  let width = matrix[0].length()
  let result : Array[Array[Double]] = Array::new(capacity=width)
  for j in 0.. Array[Array[Double]] {
  if matrix.length() == 0 {
    return []
  }
  let scales = robust_feature_scales(matrix)
  let result : Array[Array[Double]] = Array::new(capacity=matrix.length())
  for row in matrix {
    let scaled = Array::new(capacity=row.length())
    for j in 0.. FeatureScreening {
  let width = if matrix.length() == 0 { 0 } else { matrix[0].length() }
  let variances = feature_variances(matrix)
  let correlations = Array::new(capacity=width)
  let selected : Array[Int] = Array::new()
  let dropped : Array[Int] = Array::new()
  let reasons : Array[String] = Array::new()
  for j in 0.. minimum_variance &&
      correlation_value.abs() < maximum_correlation {
      selected.push(j)
    } else {
      dropped.push(j)
      reasons.push(
        if variances[j] <= minimum_variance {
          "low variance"
        } else {
          "near-perfect outcome correlation"
        },
      )
    }
  }
  {
    selected,
    dropped,
    variance: variances,
    outcome_correlation: correlations,
    reasons,
  }
}

///|
/// Computes a variance inflation proxy from pairwise correlations.
pub fn variance_inflation_proxy(matrix : Array[Array[Double]]) -> Array[Double] {
  if matrix.length() == 0 {
    return []
  }
  let width = matrix[0].length()
  let result = Array::make(width, 1.0)
  for i in 0.. Array[Int] {
  let vif = variance_inflation_proxy(matrix)
  let result : Array[Int] = Array::new()
  for i in 0.. FeatureSplit {
  let n = matrix.length()
  let validation_size = (n.to_double() * clamp(validation_fraction, 0.0, 0.9))
    .round()
    .to_int()
  let order = shuffled_indices(n, seed)
  let validation_indices = order[:validation_size].to_owned()
  let train_indices = order[validation_size:].to_owned()
  let train_values = Array::new(capacity=train_indices.length())
  let validation_values = Array::new(capacity=validation_indices.length())
  for index in train_indices {
    train_values.push(matrix[index])
  }
  for index in validation_indices {
    validation_values.push(matrix[index])
  }
  { train_indices, validation_indices, train_values, validation_values }
}

///|
/// Creates deterministic feature names for a matrix width.
pub fn generated_feature_names(
  width : Int,
  prefix? : String = "x",
) -> Array[String] {
  let result : Array[String] = Array::new(
    capacity=if width > 0 { width } else { 0 },
  )
  for i in 0.. FeatureSet {
  let result : Array[Array[Double]] = Array::new(capacity=matrix.length())
  for row in matrix {
    let selected = Array::new(capacity=indexes.length())
    for index in indexes {
      if index >= 0 && index < row.length() {
        selected.push(row[index])
      }
    }
    result.push(selected)
  }
  let names = generated_feature_names(indexes.length())
  feature_set(result, names~)
}

///|
/// Concatenates two feature matrices column-wise.
pub fn concatenate_features(
  first : Array[Array[Double]],
  second : Array[Array[Double]],
) -> FeatureSet {
  let n = first.length().min(second.length())
  let result : Array[Array[Double]] = Array::new(capacity=n)
  for i in 0.. Array[Array[Double]] {
  let result : Array[Array[Double]] = Array::new(capacity=matrix.length())
  for row in matrix {
    let centered = Array::new(capacity=row.length())
    for j in 0.. Array[Array[Double]] {
  let result : Array[Array[Double]] = Array::new(capacity=matrix.length())
  for row in matrix {
    let scaled = Array::new(capacity=row.length())
    for j in 0.. Array[Double] {
  let width = if buckets > 0 { buckets } else { 1 }
  let result = Array::make(width, 0.0)
  for value in values {
    let index = value.abs() % width
    result[index] += 1.0
  }
  result
}

///|
/// Computes row norms for regularization diagnostics.
pub fn feature_row_norms(matrix : Array[Array[Double]]) -> Array[Double] {
  let result = Array::new(capacity=matrix.length())
  for row in matrix {
    let mut norm = 0.0
    for value in row {
      norm += value * value
    }
    result.push(norm.sqrt())
  }
  result
}

///|
/// Computes column L2 norms.
pub fn feature_column_norms(matrix : Array[Array[Double]]) -> Array[Double] {
  if matrix.length() == 0 {
    return []
  }
  let width = matrix[0].length()
  let result = Array::make(width, 0.0)
  for row in matrix {
    for j in 0.. Array[Array[Double]] {
  let norms = feature_column_norms(matrix)
  let result : Array[Array[Double]] = Array::new(capacity=matrix.length())
  for row in matrix {
    let normalized = Array::new(capacity=row.length())
    for j in 0.. Double {
  let mut nonzero = 0
  let mut total = 0
  for row in matrix {
    for value in row {
      total += 1
      if value.abs() > tolerance {
        nonzero += 1
      }
    }
  }
  if total == 0 {
    0.0
  } else {
    nonzero.to_double() / total.to_double()
  }
}

///|
/// Returns a compact feature engineering summary vector.
pub fn feature_summary(features : FeatureSet) -> Array[Double] {
  [
    features.row_count.to_double(),
    features.column_count.to_double(),
    feature_density(features.values),
    features.fingerprint.to_double(),
  ]
}