///|
pub fn finite_difference_derivative(
  f : (Double) -> Double raise VleError,
  point : Double,
  step : Double,
) -> Double raise VleError {
  if step <= 0.0 {
    raise VleError::InvalidParameter("finite-difference step must be positive")
  }
  (f(point + step) - f(point - step)) / (2.0 * step)
}

///|
pub fn forward_difference_derivative(
  f : (Double) -> Double raise VleError,
  point : Double,
  step : Double,
) -> Double raise VleError {
  if step <= 0.0 {
    raise VleError::InvalidParameter("finite-difference step must be positive")
  }
  (f(point + step) - f(point)) / step
}

///|
pub fn saturation_pressure_sensitivity(
  component : Component,
  temperature_k : Double,
  step_k : Double,
) -> SensitivityPoint raise VleError {
  assert_temperature(temperature_k)
  let pressure = component.saturation_pressure_bar(temperature_k)
  let derivative = finite_difference_derivative(
    fn(t : Double) -> Double raise VleError {
      component.saturation_pressure_bar(t)
    },
    temperature_k,
    step_k,
  )
  SensitivityPoint::{
    temperature_k,
    pressure_bar: pressure,
    derivative_bar_per_k: derivative,
    relative_sensitivity: derivative * temperature_k / clamp_positive(pressure),
  }
}

///|
pub fn temperature_sensitivity_profile(
  component : Component,
  temperatures_k : Array[Double],
  step_k : Double,
) -> Array[SensitivityPoint] raise VleError {
  if temperatures_k.length() == 0 {
    raise VleError::EmptyMixture
  }
  [
    for temperature in temperatures_k => {
      saturation_pressure_sensitivity(component, temperature, step_k)
    }
  ]
}

///|
pub fn pressure_sensitivity_ratio(
  component : Component,
  low_temperature_k : Double,
  high_temperature_k : Double,
  step_k : Double,
) -> Double raise VleError {
  let low = saturation_pressure_sensitivity(
    component, low_temperature_k, step_k,
  )
  let high = saturation_pressure_sensitivity(
    component, high_temperature_k, step_k,
  )
  high.derivative_bar_per_k / clamp_positive(low.derivative_bar_per_k)
}

///|
pub fn relative_sensitivity_from_perturbation(
  base_value : Double,
  perturbed_value : Double,
  relative_perturbation : Double,
) -> Double raise VleError {
  if base_value == 0.0 || relative_perturbation == 0.0 {
    raise VleError::InvalidParameter(
      "sensitivity base and perturbation cannot be zero",
    )
  }
  (perturbed_value - base_value) / (base_value * relative_perturbation)
}

///|
pub fn pressure_uncertainty_from_temperature(
  sensitivity : SensitivityPoint,
  temperature_uncertainty_k : Double,
) -> Double raise VleError {
  if temperature_uncertainty_k < 0.0 {
    raise VleError::InvalidParameter(
      "temperature uncertainty cannot be negative",
    )
  }
  abs_double(sensitivity.derivative_bar_per_k) * temperature_uncertainty_k
}

///|
pub fn combine_independent_uncertainties(
  uncertainties : Array[Double],
) -> Double raise VleError {
  if uncertainties.length() == 0 {
    raise VleError::EmptyMixture
  }
  for uncertainty in uncertainties {
    if uncertainty < 0.0 {
      raise VleError::InvalidParameter("uncertainty cannot be negative")
    }
  }
  @math.pow(
    for i = 0, value = 0.0; i < uncertainties.length(); {
      continue i + 1, value + uncertainties[i] * uncertainties[i]
    } nobreak {
      value
    },
    0.5,
  )
}

///|
pub fn bubble_pressure_temperature_sensitivity(
  components : Array[Component],
  liquid : Array[Double],
  temperature_k : Double,
  step_k : Double,
  model? : ActivityModel = ActivityModel::Ideal,
) -> SensitivityPoint raise VleError {
  let pressure = bubble_pressure_bar(components, liquid, temperature_k, model~).pressure_bar
  let derivative = finite_difference_derivative(
    fn(t : Double) -> Double raise VleError {
      bubble_pressure_bar(components, liquid, t, model~).pressure_bar
    },
    temperature_k,
    step_k,
  )
  SensitivityPoint::{
    temperature_k,
    pressure_bar: pressure,
    derivative_bar_per_k: derivative,
    relative_sensitivity: derivative * temperature_k / clamp_positive(pressure),
  }
}

///|
pub fn bubble_pressure_composition_sensitivity(
  components : Array[Component],
  liquid : Array[Double],
  temperature_k : Double,
  fraction_step : Double,
) -> Array[Double] raise VleError {
  let x = normalize(liquid)
  if fraction_step <= 0.0 || fraction_step >= 0.5 {
    raise VleError::InvalidParameter("composition step must be in (0, 0.5)")
  }
  if components.length() != 2 || x.length() != 2 {
    raise VleError::LengthMismatch(expected=2, actual=components.length())
  }
  let low = bubble_pressure_bar(
      components,
      [
        clamp_positive(x[0] - fraction_step),
        1.0 - clamp_positive(x[0] - fraction_step),
      ],
      temperature_k,
    ).pressure_bar
  let high = bubble_pressure_bar(
      components,
      [
        clamp_positive(x[0] + fraction_step),
        1.0 - clamp_positive(x[0] + fraction_step),
      ],
      temperature_k,
    ).pressure_bar
  [low, high, (high - low) / (2.0 * fraction_step)]
}

///|
pub fn sensitivity_profile_mean(
  profile : Array[SensitivityPoint],
) -> Double raise VleError {
  if profile.length() == 0 {
    raise VleError::EmptyMixture
  }
  for i = 0, value = 0.0; i < profile.length(); {
    continue i + 1, value + profile[i].relative_sensitivity
  } nobreak {
    value / profile.length().to_double()
  }
}

///|
pub fn sensitivity_profile_maximum(
  profile : Array[SensitivityPoint],
) -> SensitivityPoint raise VleError {
  if profile.length() == 0 {
    raise VleError::EmptyMixture
  }
  for i = 1, index = 0; i < profile.length(); {
    continue i + 1,
      if profile[i].relative_sensitivity > profile[index].relative_sensitivity {
        i
      } else {
        index
      }
  } nobreak {
    profile[index]
  }
}

///|
pub fn sensitivity_is_monotone_positive(
  profile : Array[SensitivityPoint],
) -> Bool {
  for point in profile {
    if point.derivative_bar_per_k <= 0.0 {
      return false
    }
  }
  true
}

///|
pub fn normalized_sensitivity(
  derivative : Double,
  value : Double,
  scale : Double,
) -> Double raise VleError {
  if value == 0.0 || scale <= 0.0 {
    raise VleError::InvalidParameter("normalization values are invalid")
  }
  derivative * scale / value
}

///|
pub fn perturb_value(
  value : Double,
  relative_change : Double,
) -> Double raise VleError {
  if relative_change <= -1.0 {
    raise VleError::InvalidParameter(
      "relative change would make value nonpositive",
    )
  }
  value * (1.0 + relative_change)
}

///|
pub fn symmetric_perturbations(
  value : Double,
  relative_change : Double,
) -> (Double, Double) raise VleError {
  if relative_change <= 0.0 || relative_change >= 1.0 {
    raise VleError::InvalidParameter("symmetric perturbation must be in (0, 1)")
  }
  (value * (1.0 - relative_change), value * (1.0 + relative_change))
}

///|
pub fn finite_difference_condition_number(
  base_value : Double,
  low_value : Double,
  high_value : Double,
  relative_step : Double,
) -> Double raise VleError {
  if base_value == 0.0 || relative_step <= 0.0 {
    raise VleError::InvalidParameter("condition-number inputs are invalid")
  }
  abs_double(high_value - low_value) /
  (2.0 * abs_double(base_value) * relative_step)
}