<?php

declare(strict_types=1);

namespace Dochia;

/**
 * Predictive distribution with body scale and separate tail stretching.
 *
 * The body (between tail anchors) is a lognormal mixture; the lower and upper
 * tails are stretched so that the 80% interval and the 95% interval can be
 * calibrated independently. This is a coherent piecewise distribution.
 */
final class CalibratedPredictive
{
    public function __construct(
        public readonly LognormalMixture $base,
        public readonly float $anchorLow,
        public readonly float $anchorHigh,
        public readonly float $scaleLow,
        public readonly float $scaleHigh,
        public readonly float $qLow,
        public readonly float $qHigh,
    ) {
    }

    public static function fromBase(LognormalMixture $base, ModelParams $params): self
    {
        return new self(
            base: $base,
            anchorLow: $params->tailAnchorLow,
            anchorHigh: $params->tailAnchorHigh,
            scaleLow: max($params->tailScaleLow, 1e-6),
            scaleHigh: max($params->tailScaleHigh, 1e-6),
            qLow: $base->quantile($params->tailAnchorLow),
            qHigh: $base->quantile($params->tailAnchorHigh),
        );
    }

    public function quantile(float $p): float
    {
        if ($p <= 0.0) {
            return 0.0;
        }
        if ($p >= 1.0) {
            return 1000.0;
        }
        if ($p <= $this->anchorLow) {
            return $this->qLow + ($this->base->quantile($p) - $this->qLow) * $this->scaleLow;
        }
        if ($p >= $this->anchorHigh) {
            return $this->qHigh + ($this->base->quantile($p) - $this->qHigh) * $this->scaleHigh;
        }
        return $this->base->quantile($p);
    }

    public function cdf(float $x): float
    {
        if ($x <= $this->qLow) {
            $y = $this->qLow + ($x - $this->qLow) / $this->scaleLow;
            return min(max($this->base->cdf($y), 0.0), 1.0);
        }
        if ($x >= $this->qHigh) {
            $y = $this->qHigh + ($x - $this->qHigh) / $this->scaleHigh;
            return min(max($this->base->cdf($y), 0.0), 1.0);
        }
        return $this->base->cdf($x);
    }

    public function pdf(float $x): float
    {
        if ($x <= 0.0) {
            return 0.0;
        }
        if ($x < $this->qLow) {
            return $this->base->pdf($this->qLow + ($x - $this->qLow) / $this->scaleLow) / $this->scaleLow;
        }
        if ($x > $this->qHigh) {
            return $this->base->pdf($this->qHigh + ($x - $this->qHigh) / $this->scaleHigh) / $this->scaleHigh;
        }
        return $this->base->pdf($x);
    }

    public function mean(): float
    {
        return $this->base->mean;
    }

    public function sd(): float
    {
        return $this->base->sd;
    }
}