// SPDX-License-Identifier: Unlicense OR MIT

// Scroll fling physics.
// Ported from gioui.org/internal/fling (animation.go + extrapolation.go).

package gio

import (
	"math"
	"time"
)

// FlingAnimation models exponential deceleration of a fling gesture.
type FlingAnimation struct {
	x  float32
	t0 time.Time
	v0 float32
}

const (
	minFlingVelocityDp  Dp  = 50
	maxFlingVelocityDp  Dp  = 8000
	thresholdVelocityPx     = 1
	flingDecay              = -4.2 // pixels/second^2 drag constant
)

// FlingStart initiates a fling from the given velocity (pixels/second).
// Returns true if velocity is above the minimum threshold.
func (f *FlingAnimation) FlingStart(m Metric, now time.Time, velocity float32) bool {
	min := float32(m.Dp(minFlingVelocityDp))
	v := velocity
	if -min <= v && v <= min {
		return false
	}
	max := float32(m.Dp(maxFlingVelocityDp))
	if v > max {
		v = max
	} else if v < -max {
		v = -max
	}
	f.t0 = now
	f.v0 = v
	f.x = 0
	return true
}

// FlingActive reports whether a fling is in progress.
func (f *FlingAnimation) FlingActive() bool { return f.v0 != 0 }

// FlingStop stops any active fling.
func (f *FlingAnimation) FlingStop() { f.v0 = 0 }

// FlingTick advances the fling and returns the distance to scroll.
func (f *FlingAnimation) FlingTick(now time.Time) int {
	if !f.FlingActive() {
		return 0
	}
	k := float32(flingDecay)
	t := now.Sub(f.t0)
	ekt := float32(math.Exp(float64(k) * t.Seconds()))
	x := f.v0*ekt/k - f.v0/k
	dist := x - f.x
	idist := int(dist)
	f.x += float32(idist)
	v := f.v0 * ekt
	if -thresholdVelocityPx < v && v < thresholdVelocityPx {
		f.v0 = 0
	}
	return idist
}

// --- Velocity extrapolation ---

// FlingExtrapolation estimates fling velocity from a series of position samples.
// Uses Android's least-squares polynomial fit algorithm.
type FlingExtrapolation struct {
	idx       int
	samples   []flingSample
	lastValue float32
	cache     [flingHistorySize]flingSample
	values    [flingHistorySize]float32
	times     [flingHistorySize]float32
}

type flingSample struct {
	t time.Duration
	v float32
}

// FlingEstimate is the estimated velocity and distance from extrapolation.
type FlingEstimate struct {
	Velocity float32
	Distance float32
}

const (
	flingDegree      = 2
	flingHistorySize = 20
	flingMaxAge      = 100 * time.Millisecond
	flingMaxGap      = 40 * time.Millisecond
)

type flingCoeffs [flingDegree + 1]float32

// Sample records a position observation at time t.
func (e *FlingExtrapolation) Sample(t time.Duration, val float32) {
	e.lastValue = val
	if e.samples == nil {
		e.samples = e.cache[:0]
	}
	s := flingSample{t: t, v: val}
	if e.idx == len(e.samples) && e.idx < cap(e.samples) {
		e.samples = append(e.samples, s)
	} else {
		e.samples[e.idx] = s
	}
	e.idx++
	if e.idx == cap(e.samples) {
		e.idx = 0
	}
}

// SampleDelta records a delta sample relative to the last value.
func (e *FlingExtrapolation) SampleDelta(t time.Duration, delta float32) {
	e.Sample(t, delta+e.lastValue)
}

// Estimate computes the velocity and distance from recorded samples.
func (e *FlingExtrapolation) Estimate() FlingEstimate {
	if len(e.samples) == 0 {
		return FlingEstimate{}
	}
	values := e.values[:0]
	times := e.times[:0]
	first := e.flingGet(0)
	prev := first.t
	for i := range e.samples {
		p := e.flingGet(-i)
		age := first.t - p.t
		if age >= flingMaxAge || prev-p.t >= flingMaxGap {
			break
		}
		prev = p.t
		values = append(values, p.v)
		times = append(times, float32(p.t-first.t)/float32(time.Second))
	}
	if len(values) == 1 {
		return FlingEstimate{}
	}
	degree := flingDegree
	if len(values) <= degree {
		degree = len(values) - 1
	}
	coeffs, ok := flingPolyFit(times, values, degree)
	if !ok {
		return FlingEstimate{}
	}
	// Velocity is the derivative at t=0.
	velocity := coeffs[1]
	// Distance is the polynomial value at t=0 minus current.
	distance := values[0] - first.v
	return FlingEstimate{Velocity: velocity, Distance: distance}
}

func (e *FlingExtrapolation) flingGet(i int) flingSample {
	n := len(e.samples)
	idx := (e.idx + n + i - 1) % n
	return e.samples[idx]
}

// flingPolyFit fits a polynomial of the given degree to (times, values).
// Returns coefficients [c0, c1, c2, ...] and whether the fit succeeded.
func flingPolyFit(times, values []float32, degree int) (flingCoeffs, bool) {
	// Build the Vandermonde matrix A where A[i][j] = times[i]^j.
	n := len(values)
	cols := degree + 1
	A := []float32{:n * cols}
	for i, t := range times {
		for j := 0; j < cols; j++ {
			A[i*cols+j] = float32(math.Pow(float64(t), float64(j)))
		}
	}
	// Solve the normal equations A^T A c = A^T b using Cholesky-like method.
	// For simplicity, use Gaussian elimination on the normal matrix.
	ATA := []float32{:cols * cols}
	ATb := []float32{:cols}
	for i := 0; i < cols; i++ {
		for k := 0; k < n; k++ {
			ATb[i] += A[k*cols+i] * values[k]
			for j := 0; j < cols; j++ {
				ATA[i*cols+j] += A[k*cols+i] * A[k*cols+j]
			}
		}
	}
	// Gaussian elimination with partial pivoting.
	aug := []float32{:cols * (cols + 1)}
	for i := 0; i < cols; i++ {
		for j := 0; j < cols; j++ {
			aug[i*(cols+1)+j] = ATA[i*cols+j]
		}
		aug[i*(cols+1)+cols] = ATb[i]
	}
	for col := 0; col < cols; col++ {
		// Find pivot.
		pivot := col
		for row := col + 1; row < cols; row++ {
			if math.Abs(float64(aug[row*(cols+1)+col])) > math.Abs(float64(aug[pivot*(cols+1)+col])) {
				pivot = row
			}
		}
		// Swap rows.
		for j := 0; j <= cols; j++ {
			aug[col*(cols+1)+j], aug[pivot*(cols+1)+j] = aug[pivot*(cols+1)+j], aug[col*(cols+1)+j]
		}
		if aug[col*(cols+1)+col] == 0 {
			return flingCoeffs{}, false
		}
		// Eliminate.
		for row := 0; row < cols; row++ {
			if row == col {
				continue
			}
			factor := aug[row*(cols+1)+col] / aug[col*(cols+1)+col]
			for j := col; j <= cols; j++ {
				aug[row*(cols+1)+j] -= factor * aug[col*(cols+1)+j]
			}
		}
	}
	var coeffs flingCoeffs
	for i := 0; i < cols; i++ {
		coeffs[i] = aug[i*(cols+1)+cols] / aug[i*(cols+1)+i]
	}
	return coeffs, true
}