<?php

declare(strict_types=1);

namespace Dochia;

/**
 * Dochia v1 — reference implementation (pure PHP, no deps).
 *
 * Pure PHP reference implementation; mirrors python/dochia.py. The golden
 * fixtures in tests/fixtures/golden.json bind the two implementations.
 */
final class MicroModel
{
    public function __construct(public readonly ModelParams $params)
    {
    }

    /**
     * Personalise a user's cycle-length distribution from observed lengths.
     *
     * @param  list<int|float>  $lengths
     * @return array<string, mixed>
     */
    public function personalState(array $lengths): array
    {
        $lengths = $this->filterLengths($lengths);
        $n = count($lengths);

        if ($n === 0) {
            $withinVar = $this->params->withinSdPrior ** 2;
            $betweenVar = $this->params->popBetweenSd ** 2;
            $posteriorVar = ($withinVar * $betweenVar) / max(1e-9, $withinVar);
            return [
                'mean' => $this->params->popMean,
                'within_var' => $withinVar,
                'cycle_sd' => sqrt($withinVar + $posteriorVar),
                'n' => 0,
                'tail_prob' => $this->params->tailProbBase,
                'trend' => 0.0,
                'recent_sd' => 0.0,
                'recent_range' => 0,
                'gap' => false,
                'abstain' => false,
                'confidence' => 'new',
            ];
        }

        $mean = array_sum($lengths) / $n;
        $sampleVar = 0.0;
        if ($n > 1) {
            $sumSq = 0.0;
            foreach ($lengths as $y) {
                $sumSq += ($y - $mean) ** 2;
            }
            $sampleVar = $sumSq / ($n - 1);
        }
        $median = self::median($lengths);
        $absDev = [];
        foreach ($lengths as $y) {
            $absDev[] = abs($y - $median);
        }
        $robustSd = 1.4826 * self::median($absDev);
        $robustVar = $robustSd * $robustSd;
        $sampleVar = $this->params->varianceBlend * $sampleVar + (1.0 - $this->params->varianceBlend) * $robustVar;

        $priorVar = $this->params->withinSdPrior ** 2;
        $priorWeight = $this->params->priorWeight;
        // Trust the user's data more when observed variance clearly exceeds the
        // prior; this widens intervals for genuinely variable users instead of
        // shrinking them back to the population prior.
        if ($sampleVar > $priorVar) {
            $priorWeight = max(0.5, $this->params->priorWeight * $priorVar / $sampleVar);
        }
        $withinVar = ($n * $sampleVar + $priorWeight * $priorVar) / ($n + $priorWeight);
        $withinVar = min(max($withinVar, $this->params->minWithinVar), $this->params->maxWithinVar);

        $betweenVar = $this->params->popBetweenSd ** 2;
        $weight = ($n * $betweenVar) / max(1e-9, $n * $betweenVar + $withinVar);
        $shrunkMean = $weight * $mean + (1.0 - $weight) * $this->params->popMean;
        $shrunkMean = min(max($shrunkMean, $this->params->minMean), $this->params->maxMean);

        $recent = array_slice($lengths, -6);
        $trend = 0.0;
        if ($n >= 4 && count($recent) >= 3) {
            $weights = [];
            $w = 1.0;
            foreach (array_reverse($recent) as $ignored) {
                $weights[] = $w;
                $w *= (1.0 - $this->params->trendEwma);
            }
            $weights = array_reverse($weights);
            $totalW = array_sum($weights);
            $acc = 0.0;
            foreach ($recent as $i => $y) {
                $acc += ($y - $shrunkMean) * $weights[$i];
            }
            $trend = $acc / $totalW;
            $shift = abs($trend) < $this->params->trendThreshold
                ? 0.0
                : max(-$this->params->trendMaxShift, min($this->params->trendMaxShift, $trend));
            $shrunkMean = min(max($shrunkMean + $shift, $this->params->minMean), $this->params->maxMean);
        }

        $recentMean = array_sum($recent) / count($recent);
        $recentSumSq = 0.0;
        foreach ($recent as $y) {
            $recentSumSq += ($y - $recentMean) ** 2;
        }
        $recentSd = sqrt($recentSumSq / max(1, count($recent) - 1));
        $recentRange = (int) (max($recent) - min($recent));
        $gap = max($recent) > $this->params->abstainGapDays;

        $pi = $this->params->tailProbBase
            + $this->params->tailProbSdGain * max(0.0, $recentSd - 2.0)
            + $this->params->tailProbRangeGain * max(0.0, $recentRange - 7);
        if ($gap) {
            $pi += $this->params->tailProbGapBonus;
        }
        $pi = min(max($pi, $this->params->tailProbBase), $this->params->tailProbMax);

        $abstain = $gap || $recentSd > $this->params->abstainRecentSd || $n < $this->params->abstainMinCycles;
        $confidence = $abstain
            ? 'abstain'
            : ($n >= 6 ? 'established' : ($n >= 3 ? 'building' : ($n >= 1 ? 'learning' : 'new')));

        return [
            'mean' => $shrunkMean,
            'within_var' => $withinVar,
            'cycle_sd' => sqrt($withinVar + ($withinVar * $betweenVar) / max(1e-9, $n * $betweenVar + $withinVar)),
            'n' => $n,
            'tail_prob' => $pi,
            'trend' => $trend,
            'recent_sd' => $recentSd,
            'recent_range' => $recentRange,
            'gap' => $gap,
            'abstain' => $abstain,
            'confidence' => $confidence,
        ];
    }

    /** @param array<string, mixed> $state */
    public function predictive(array $state): CalibratedPredictive
    {
        $sd = ((float) $state['cycle_sd']) * $this->params->bodyScale;
        $base = LognormalMixture::fromMoments(
            mean: (float) $state['mean'],
            sd: $sd,
            pi: (float) $state['tail_prob'],
            tailShift: $this->params->tailShiftDays,
            tailSdMultiplier: $this->params->tailSdMultiplier,
            momentMatch: $this->params->momentMatch,
        );
        return CalibratedPredictive::fromBase($base, $this->params);
    }

    /**
     * @param  list<int|float>  $lengths
     * @return array<string, mixed>
     */
    public function forecast(array $lengths, int $currentDay = 1): array
    {
        $state = $this->personalState($lengths);
        $dist = $this->predictive($state);
        $hazard = $this->onsetProbabilities($dist, $currentDay);
        $ovulation = $this->ovulationDistribution($dist);

        $quantiles = [];
        foreach ([0.025, 0.05, 0.10, 0.25, 0.50, 0.75, 0.90, 0.95, 0.975] as $q) {
            $quantiles[(string) $q] = $dist->quantile($q);
        }

        $explanations = [];
        if ($state['abstain']) {
            $explanations[] = 'Insufficient or highly variable history for a confident forecast.';
        }
        if ($state['gap']) {
            $explanations[] = 'A long gap was detected; the forecast is intentionally wide.';
        }
        if ($state['recent_sd'] > 4.0) {
            $explanations[] = 'Recent cycle lengths are unusually variable; the interval is wide.';
        }
        if ($state['n'] < 3) {
            $explanations[] = 'Only a few cycles recorded; the population prior dominates.';
        }
        if ($state['trend'] != 0.0) {
            $explanations[] = 'A recent trend was applied to the expected length.';
        }

        return [
            'params_version' => $this->params->version,
            'mean_length' => (float) $state['mean'],
            'sd_length' => ((float) $state['cycle_sd']) * $this->params->bodyScale,
            'cycles_used' => (int) $state['n'],
            'confidence' => (string) $state['confidence'],
            'abstain' => (bool) $state['abstain'],
            'abstain_reason' => $this->abstainReason($state),
            'tail_prob' => (float) $state['tail_prob'],
            'trend' => (float) $state['trend'],
            'quantiles' => $quantiles,
            'hazard' => $hazard,
            'ovulation' => [
                'mean_day' => $ovulation['mean_day'],
                'sd_day' => $ovulation['sd_day'],
                'fertile_probabilities' => $ovulation['fertile_probabilities'],
                'peak_days' => $ovulation['peak_days'],
                'window_start' => $ovulation['window_start'],
                'window_end' => $ovulation['window_end'],
            ],
            'prompt' => $this->prompt($state, $hazard),
            'explanations' => $explanations,
        ];
    }

    /** @return array<int, float> day => hazard */
    public function hazardCurve(CalibratedPredictive $dist): array
    {
        $curve = [];
        $prevSurvival = 1.0;
        for ($day = 1; $day <= $this->params->hazardMaxDay; $day++) {
            if ($prevSurvival <= 1e-12) {
                $curve[$day] = 1.0;
                continue;
            }
            $f = $dist->cdf((float) $day) - $dist->cdf((float) ($day - 1));
            $h = max(0.0, $f) / $prevSurvival;
            $curve[$day] = min(max($h, 0.0), 1.0);
            $prevSurvival -= max(0.0, $f);
            if ($prevSurvival < 0.0) {
                $prevSurvival = 0.0;
            }
        }
        return $curve;
    }

    /** @return array<string, mixed> */
    public function onsetProbabilities(CalibratedPredictive $dist, int $currentDay): array
    {
        $currentDay = max(1, $currentDay);
        $survivalNow = max(1e-12, 1.0 - $dist->cdf((float) ($currentDay - 1)));
        $within = function (int $k) use ($dist, $currentDay, $survivalNow): float {
            $survivalThen = max(0.0, 1.0 - $dist->cdf((float) ($currentDay - 1 + $k)));
            return min(max(1.0 - $survivalThen / $survivalNow, 0.0), 1.0);
        };
        $curve = $this->hazardCurve($dist);
        $pToday = $curve[$currentDay] ?? 0.0;
        $window = [];
        for ($d = $currentDay; $d <= min($this->params->hazardMaxDay, $currentDay + 13); $d++) {
            $window[$d] = $curve[$d] ?? 0.0;
        }
        return [
            'current_day' => $currentDay,
            'p_today' => $pToday,
            'p_within_3' => $within(3),
            'p_within_7' => $within(7),
            'p_within_14' => $within(14),
            'late_probability' => 1.0 - $dist->cdf((float) ($currentDay - 1)),
            'curve' => $window,
        ];
    }

    /** @return array<string, mixed> */
    public function ovulationDistribution(CalibratedPredictive $dist): array
    {
        $sdO = sqrt(max($dist->sd() ** 2 - $this->params->lutealSd ** 2, $this->params->ovulationSdFloor ** 2));
        $meanO = $dist->mean() - $this->params->lutealMean;

        $step = $this->params->ovulationGridStep;
        $maxY = 365.0;
        $edges = [];
        $y = 0.0;
        while ($y <= $maxY) {
            $edges[] = $y;
            $y += $step;
        }
        $cdfEdges = [];
        foreach ($edges as $edge) {
            $cdfEdges[] = $dist->cdf($edge);
        }

        $cdfO = function (float $d) use ($edges, $cdfEdges): float {
            $total = 0.0;
            $count = count($edges);
            for ($i = 0; $i < $count - 1; $i++) {
                $mass = $cdfEdges[$i + 1] - $cdfEdges[$i];
                if ($mass <= 0.0) {
                    continue;
                }
                $yMid = 0.5 * ($edges[$i] + $edges[$i + 1]);
                $z = ($d - $yMid + $this->params->lutealMean) / $this->params->lutealSd;
                $total += $mass * Erf::normCdf($z);
            }
            return $total;
        };

        $probs = [];
        $prev = 0.0;
        for ($day = 1; $day <= $this->params->hazardMaxDay; $day++) {
            $cur = $cdfO((float) $day);
            $probs[$day] = min(max($cur - $prev, 0.0), 1.0);
            $prev = $cur;
        }
        $probs[$this->params->hazardMaxDay] += max(0.0, 1.0 - $prev);

        $fertile = [];
        for ($day = 1; $day <= $this->params->hazardMaxDay; $day++) {
            $window = 0.0;
            $last = min($this->params->hazardMaxDay, $day + $this->params->fertileWindowDays - 1);
            for ($o = $day; $o <= $last; $o++) {
                $window += $probs[$o];
            }
            $fertile[$day] = min(max($window, 0.0), 1.0);
        }

        $maxProb = max($fertile);
        $peakDays = [];
        foreach ($fertile as $day => $p) {
            if ($p >= max(0.25, 0.5 * $maxProb)) {
                $peakDays[] = $day;
            }
        }
        $windowDays = [];
        foreach ($fertile as $day => $p) {
            if ($p >= 0.05) {
                $windowDays[] = $day;
            }
        }

        return [
            'mean_day' => $meanO,
            'sd_day' => $sdO,
            'probabilities' => $probs,
            'fertile_probabilities' => $fertile,
            'peak_days' => $peakDays,
            'window_start' => $windowDays === [] ? null : min($windowDays),
            'window_end' => $windowDays === [] ? null : max($windowDays),
        ];
    }

    /** @param array<string, mixed> $state */
    private function abstainReason(array $state): ?string
    {
        if (! $state['abstain']) {
            return null;
        }
        if ($state['gap']) {
            return 'long_gap';
        }
        if ($state['recent_sd'] > $this->params->abstainRecentSd) {
            return 'high_variability';
        }
        return 'insufficient_history';
    }

    /**
     * @param  list<int|float>  $lengths
     * @return list<float>
     */
    private function filterLengths(array $lengths): array
    {
        $out = [];
        foreach ($lengths as $length) {
            $y = (float) $length;
            if ($y >= $this->params->minLength && $y <= $this->params->maxLength) {
                $out[] = $y;
            }
        }
        return $out;
    }

    /** @param list<int|float> $values */
    private static function median(array $values): float
    {
        $n = count($values);
        if ($n === 0) {
            return 0.0;
        }
        sort($values);
        return $n % 2 === 1 ? (float) $values[intdiv($n, 2)] : 0.5 * ((float) $values[$n / 2 - 1] + (float) $values[$n / 2]);
    }

    /** @param array<string, mixed> $state @param array<string, mixed> $hazard */
    private function prompt(array $state, array $hazard): string
    {
        if ($state['abstain']) {
            return 'Not enough signal for a reliable forecast.';
        }
        $p3 = (float) $hazard['p_within_3'];
        if ($p3 < 0.05) {
            return 'Your period is unlikely in the next three days.';
        }
        if ($p3 < 0.25) {
            return 'There is a small chance your period starts in the next three days.';
        }
        if ($p3 < 0.5) {
            return 'Your period is fairly likely within three days.';
        }
        if ($p3 < 0.8) {
            return 'Your period is likely within three days.';
        }
        return 'Your period is very likely within three days.';
    }
}