hessboost 0.2.2

//! Count and positive-continuous regression objectives with a log link:
//! Poisson, Gamma, and Tweedie. All predict `exp(margin)`.

use super::{
    GradPair, Loss, OutputDomain, check_base_score_domain, check_label_domain, log_link,
    newton_intercepts, weighted_label_mean,
};
use crate::K_RT_EPS_F32;
use crate::data::MetaInfo;
use crate::error::Result;
use crate::metric::EvalMetric;
use crate::objective::Tweedie;

/// Half the Poisson unit deviance at log-mean `margin`: `y ln(y/μ) − (y − μ)`,
/// whose margin derivative is the Poisson gradient `μ − y`.
fn poisson_deviance(margin: f32, label: f32) -> f64 {
    let (m, y) = (f64::from(margin), f64::from(label));
    let mu = m.exp();
    if y > 0.0 {
        y * (y.ln() - m) - (y - mu)
    } else {
        mu
    }
}

/// The log-link intercept of XGBoost's `FitInterceptGlmLike`: the
/// (weighted) label mean through the link.
fn log_label_mean(info: &MetaInfo) -> Vec<f32> {
    let mut margin = [weighted_label_mean(info.label_values(), info.weights)];
    log_link(&mut margin);
    margin.to_vec()
}

/// Emit the `pred_transform`/`probs_to_margins`/`validate_base_score`
/// hooks shared by the log-link objectives (all predict `exp(margin)`). The
/// link is XGBoost's `ProbToMargin`, `ln(v)` in `f32`.
macro_rules! log_link_objective {
    () => {
        fn pred_transform(&self, preds: &mut [f32]) {
            crate::simd::exp_inplace(preds);
        }

        fn probs_to_margins(&self, scores: &mut [f32]) {
            log_link(scores);
        }

        fn validate_base_score(&self, base_score: f64) -> Result<()> {
            check_base_score_domain(base_score, OutputDomain::Positive)
        }
    };
}

/// Poisson regression (`count:poisson`). Gradient is `exp(m) − y`. The Hessian is
/// stabilized by `max_delta_step` (default 0.7 in XGBoost) via
/// `exp(m + max_delta_step)`.
#[derive(Debug, Clone, Copy)]
pub struct Poisson {
    max_delta_step: f32,
}

impl Poisson {
    /// Create with the given Hessian-stabilizing max delta step.
    pub fn new(max_delta_step: f32) -> Self {
        Poisson { max_delta_step }
    }
}

impl Default for Poisson {
    fn default() -> Self {
        Poisson {
            max_delta_step: 0.7,
        }
    }
}

impl Loss for Poisson {
    fn name(&self) -> &'static str {
        "count:poisson"
    }

    fn gradient(
        &self,
        preds: &[f32],
        labels: &[f32],
        weights: Option<&[f32]>,
        out: &mut [GradPair],
    ) {
        let max_delta_step = self.max_delta_step;
        super::rowwise_gradient(
            labels.len(),
            1,
            preds,
            labels,
            weights,
            out,
            |p, l, w, o| {
                crate::simd::poisson_gradient(p, l, w, max_delta_step, o);
            },
        );
    }

    log_link_objective!();

    fn base_margins_info(&self, info: &MetaInfo) -> Vec<f32> {
        log_label_mean(info)
    }

    fn pointwise_loss(&self) -> Option<super::PointwiseLoss<'_>> {
        Some(Box::new(poisson_deviance))
    }

    fn validate_info(&self, info: &MetaInfo) -> Result<()> {
        check_label_domain(info, |y| y < 0.0)
    }

    fn default_metric(&self) -> EvalMetric {
        EvalMetric::PoissonNLogLik
    }
}

/// Gamma regression (`reg:gamma`), a log-link objective for positive
/// targets: XGBoost's `RegLossObj<GammaDeviance>`. With `p = exp(m)` in
/// `f32`, gradient `1 − y/p` and Hessian `y/p`, both times the row weight,
/// which `scale_pos_weight` multiplies for a label of exactly `1`. The
/// intercept is the (weighted) label mean through the link, or with
/// `scale_pos_weight != 1` XGBoost's Newton step (`FitIntercept`) on the
/// reweighted gradients.
#[derive(Debug, Clone, Copy)]
pub(crate) struct Gamma {
    scale_pos_weight: f32,
}

impl Gamma {
    /// The loss with positive-label weight `scale_pos_weight`.
    pub(crate) fn new(scale_pos_weight: f32) -> Self {
        Gamma { scale_pos_weight }
    }
}

impl Default for Gamma {
    fn default() -> Self {
        Gamma::new(1.0)
    }
}

impl Loss for Gamma {
    fn name(&self) -> &'static str {
        "reg:gamma"
    }

    fn gradient(
        &self,
        preds: &[f32],
        labels: &[f32],
        weights: Option<&[f32]>,
        out: &mut [GradPair],
    ) {
        let scale_pos_weight = self.scale_pos_weight;
        super::rowwise_gradient(
            labels.len(),
            1,
            preds,
            labels,
            weights,
            out,
            |p, l, w, o| {
                crate::simd::gamma_gradient(p, l, w, scale_pos_weight, o);
            },
        );
    }

    log_link_objective!();

    fn base_margins_info(&self, info: &MetaInfo) -> Vec<f32> {
        // XGBoost `RegLossObj::InitEstimation`, as for `reg:squarederror`.
        if (self.scale_pos_weight - 1.0).abs() > K_RT_EPS_F32 {
            return newton_intercepts(self, info);
        }
        log_label_mean(info)
    }

    fn pointwise_loss(&self) -> Option<super::PointwiseLoss<'_>> {
        // Half the Gamma unit deviance: `y/μ − ln(y/μ) − 1`, with a label of
        // 1 reweighted by `scale_pos_weight` exactly as in the gradient.
        let scale_pos_weight = f64::from(self.scale_pos_weight);
        Some(Box::new(move |margin, label| {
            let (m, y) = (f64::from(margin), f64::from(label));
            let weight = if label == 1.0 { scale_pos_weight } else { 1.0 };
            weight * (y * (-m).exp() + m - y.ln() - 1.0)
        }))
    }

    fn validate_info(&self, info: &MetaInfo) -> Result<()> {
        check_label_domain(info, |y| y <= 0.0)
    }

    fn default_metric(&self) -> EvalMetric {
        EvalMetric::GammaNLogLik
    }
}

/// Tweedie regression (`reg:tweedie`) with variance power `rho ∈ [1, 2)`
/// (XGBoost `tweedie_variance_power`; 1 is Poisson, 2 would be Gamma).
#[derive(Debug, Clone, Copy)]
pub(crate) struct TweedieLoss {
    param: Tweedie,
    rho: f32,
}

impl TweedieLoss {
    /// The loss at variance power `param`.
    pub(crate) fn new(param: Tweedie) -> Self {
        TweedieLoss {
            param,
            rho: param.variance_power() as f32,
        }
    }
}

impl Default for TweedieLoss {
    fn default() -> Self {
        TweedieLoss::new(Tweedie::default())
    }
}

impl Loss for TweedieLoss {
    fn name(&self) -> &'static str {
        "reg:tweedie"
    }

    fn gradient(
        &self,
        preds: &[f32],
        labels: &[f32],
        weights: Option<&[f32]>,
        out: &mut [GradPair],
    ) {
        let rho = self.rho;
        super::rowwise_gradient(
            labels.len(),
            1,
            preds,
            labels,
            weights,
            out,
            |p, l, w, o| {
                crate::simd::tweedie_gradient(p, l, w, rho, o);
            },
        );
    }

    log_link_objective!();

    fn base_margins_info(&self, info: &MetaInfo) -> Vec<f32> {
        log_label_mean(info)
    }

    fn pointwise_loss(&self) -> Option<super::PointwiseLoss<'_>> {
        // Half the Tweedie unit deviance for `1 < ρ < 2` (Poisson at `ρ = 1`):
        // `y^(2−ρ)/((1−ρ)(2−ρ)) − y μ^(1−ρ)/(1−ρ) + μ^(2−ρ)/(2−ρ)`.
        let rho = f64::from(self.rho);
        if (rho - 1.0).abs() < 1e-9 {
            return Some(Box::new(poisson_deviance));
        }
        let (a, b) = (1.0 - rho, 2.0 - rho);
        Some(Box::new(move |margin, label| {
            let (m, y) = (f64::from(margin), f64::from(label));
            y.powf(b) / (a * b) - y * (a * m).exp() / a + (b * m).exp() / b
        }))
    }

    fn validate_info(&self, info: &MetaInfo) -> Result<()> {
        check_label_domain(info, |y| y < 0.0)
    }

    fn default_metric(&self) -> EvalMetric {
        // XGBoost `TweedieRegression::Configure` names the metric with the
        // power the loss trains with (its `float`), and the metric reads
        // that name back: the `f32` power, as its shortest decimal in `f64`.
        let trained = format!("{}", self.rho).parse::<f64>().map(Tweedie::new);
        EvalMetric::TweedieNLogLik(match trained {
            Ok(Ok(power)) => power,
            _ => self.param,
        })
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::objective::{base_margins, gradient_pairs};
    use approx::assert_relative_eq;

    /// Tweedie's default metric evaluates at the power the loss trains with
    /// (its `f32`), as XGBoost names and reads it back, so a power that
    /// rounds to 1 in `f32` trains and evaluates at 1.
    #[test]
    fn tweedie_default_metric_uses_the_trained_power() {
        let loss = TweedieLoss::new(Tweedie::new(1.000_000_04).unwrap());
        let at_one = EvalMetric::TweedieNLogLik(Tweedie::new(1.0).unwrap());
        assert_eq!(loss.default_metric(), at_one);
        let (preds, labels) = ([1.0f32, 2.0], [1.0f32, 3.0]);
        let score = |metric: EvalMetric| metric.build(1).unwrap().eval(&preds, &labels, None);
        assert_eq!(
            score(loss.default_metric()).to_bits(),
            score(at_one).to_bits()
        );
        let default = TweedieLoss::new(Tweedie::default());
        assert_eq!(default.default_metric().name(), "tweedie-nloglik@1.5");
    }

    #[test]
    fn poisson_gradient_at_log_mean_is_zero_sum() {
        // At margin = log(y) the gradient exp(m)-y = 0.
        let obj = Poisson::default();
        let labels = [2.0f32, 5.0];
        let preds = [2.0f32.ln(), 5.0f32.ln()];
        let out = gradient_pairs(&obj, &preds, &labels, None);
        assert_relative_eq!(out[0].grad, 0.0, epsilon = 1e-5);
        assert_relative_eq!(out[1].grad, 0.0, epsilon = 1e-5);
        assert!(out[0].hess > 0.0);
    }

    #[test]
    fn gamma_gradient_zero_at_log_y() {
        let obj = Gamma::default();
        let labels = [3.0f32];
        let preds = [3.0f32.ln()];
        let out = gradient_pairs(&obj, &preds, &labels, None);
        // 1 - y*exp(-log y) = 1 - 1 = 0
        assert_relative_eq!(out[0].grad, 0.0, epsilon = 1e-5);
    }

    /// `scale_pos_weight` multiplies the weight of the rows labeled exactly
    /// 1 only, and the intercept becomes the Newton step from margin 0,
    /// `Σ w'(y − 1) / Σ w'y`, through `exp` and back.
    #[test]
    fn gamma_scale_pos_weight_reweights_rows_labeled_one() {
        let obj = Gamma::new(2.0);
        let out = gradient_pairs(&obj, &[0.0, 0.0], &[1.0, 4.0], Some(&[3.0, 3.0]));
        assert_eq!(out[0], GradPair::new(0.0, 6.0)); // (1 - 1) * 6, 1 * 6
        assert_eq!(out[1], GradPair::new(-9.0, 12.0)); // (1 - 4) * 3, 4 * 3
        // Labels 1 (weight 2) and 3: (0 + 2) / (2 + 3) = 0.4.
        let mut expected = [0.4f32];
        obj.pred_transform(&mut expected);
        obj.probs_to_margins(&mut expected);
        assert_eq!(base_margins(&obj, &[1.0, 3.0], None), expected.to_vec());
    }

    #[test]
    fn tweedie_transform_is_exp() {
        let obj = TweedieLoss::default();
        let mut p = [0.0f32, 1.0];
        obj.pred_transform(&mut p);
        assert_relative_eq!(p[0], 1.0, epsilon = 1e-6);
        assert_relative_eq!(p[1], 1.0f32.exp(), epsilon = 1e-6);
    }

    /// The log-link intercept is `ln(mean)` in f32 (XGBoost
    /// `FitInterceptGlmLike` + `ProbToMargin`); with weights it is the
    /// weighted mean, and a zero mean maps to `-inf` like XGBoost (which the
    /// trainer rejects rather than floors).
    #[test]
    fn log_link_intercept_is_ln_of_mean() {
        let obj = Poisson::default();
        assert_eq!(base_margins(&obj, &[2.0, 6.0], None), vec![4f32.ln()]);
        let w = [3.0f32, 1.0];
        assert_eq!(base_margins(&obj, &[2.0, 6.0], Some(&w)), vec![3f32.ln()]);
        assert_eq!(
            base_margins(&Gamma::default(), &[0.0, 0.0], None),
            vec![f32::NEG_INFINITY]
        );
    }
}