// Copyright (c) 2026 Petr Balvín <opensource@petrbalvin.org> (https://petrbalvin.org)
// SPDX-License-Identifier: MIT

package tensor_test

// The godoc examples: every flagship workflow as a runnable, checked
// snippet. pkg.go.dev renders these beside the API, and `go test`
// executes them, so the documentation cannot rot.

import (
	"fmt"
	"log"
	"math"

	tensor "sourcedock.dev/petrbalvin/tensor"
	grad "sourcedock.dev/petrbalvin/tensor/grad"
)

// A rank-1 array from literals, element-wise arithmetic, a reduction.
func ExampleAdd() {
	a, err := tensor.FromFloats([]float64{1, 2, 3, 4}, 4)
	if err != nil {
		log.Fatal(err)
	}
	b, _ := tensor.FromFloats([]float64{10, 20, 30, 40}, 4)
	sum, _ := tensor.Add(a, b)
	mean, _ := tensor.Mean(sum)
	fmt.Println(sum, mean)
	// Output: float (4) [11, 22, 33, 44] 27.5
}

// The 2-D matrix product.
func ExampleMatMul2D() {
	a, _ := tensor.FromFloats([]float64{1, 2, 3, 4}, 2, 2)
	b, _ := tensor.FromFloats([]float64{5, 6, 7, 8}, 2, 2)
	prod, err := tensor.MatMul2D(a, b)
	if err != nil {
		log.Fatal(err)
	}
	fmt.Println(prod)
	// Output: float (2, 2) [19, 22, 43, 50]
}

// Einstein summation, batched over an ellipsis axis.
func ExampleEinsum() {
	a, _ := tensor.FromFloats([]float64{1, 2, 3, 4, 5, 6, 7, 8}, 2, 2, 2)
	b, _ := tensor.FromFloats([]float64{1, 0, 0, 1, 1, 0, 0, 1}, 2, 2, 2)
	// Batched matrix product: each batch times its identity.
	got, err := tensor.Einsum("...ij,...jk->...ik", a, b)
	if err != nil {
		log.Fatal(err)
	}
	fmt.Println(got.Shape(), got.FloatAt(0), got.FloatAt(3))
	// Output: [2 2 2] 1 4
}

// Reverse-mode autograd: the gradient of Σ (x·w)² arrives exact.
func ExampleTensor_Backward() {
	x, _ := grad.FromFloat64s([]float64{1, 2, 3}, true, 3)
	w, _ := grad.FromFloat64s([]float64{0.5, -1, 2}, false, 3)
	prod, _ := x.Mul(w)
	sq, _ := prod.Pow(2)
	loss, _ := sq.Sum()
	if err := loss.Backward(); err != nil {
		log.Fatal(err)
	}
	// dL/dx = 2·x·w², non-negative because w enters squared.
	g := x.Grad()
	fmt.Printf("%.4g %.4g %.4g\n", g.FloatAt(0), g.FloatAt(1), g.FloatAt(2))
	// Output: 0.5 4 24
}

// Ordinary least squares with the full classical inference.
func ExampleLinearRegression() {
	// y = 2 + 3x on x = 0..5, the intercept column first.
	design, _ := tensor.FromFloats([]float64{
		1, 0, 1, 1, 1, 2, 1, 3, 1, 4, 1, 5,
	}, 6, 2)
	y, _ := tensor.FromFloats([]float64{2.1, 4.9, 8.2, 11.1, 13.8, 17.2}, 6)
	res, err := tensor.LinearRegression(design, y)
	if err != nil {
		log.Fatal(err)
	}
	fmt.Printf("slope %.3f ± %.3f, R2 %.4f\n",
		res.Coefficients[1], res.StandardErrors[1], res.RSquared)
	// Output: slope 3.003 ± 0.044, R2 0.9991
}

// Solving an initial value problem: the exponential decay y' = −y.
func ExampleIntegrateODE() {
	f := func(t float64, y *tensor.Array) (*tensor.Array, error) {
		return tensor.MulF(y, -1), nil
	}
	y0, _ := tensor.FromFloats([]float64{1}, 1)
	end, err := tensor.IntegrateODE(f, 0, 1, y0, tensor.ODEOptions{})
	if err != nil {
		log.Fatal(err)
	}
	fmt.Printf("y(1) = %.6f\n", end.FloatAt(0))
	// Output: y(1) = 0.367880
}

// The globally adaptive cubature over a box: the 2-D Gaussian.
func ExampleIntegrateND() {
	got, err := tensor.IntegrateND(func(x []float64) float64 {
		return math.Exp(-x[0]*x[0] - x[1]*x[1])
	}, []float64{-3, -3}, []float64{3, 3}, tensor.CubatureOptions{Tolerance: 1e-11})
	if err != nil {
		log.Fatal(err)
	}
	fmt.Printf("%.6f\n", got)
	// Output: 3.141454
}

// The Crank-Nicolson heat equation holds its eigenmode shape while it
// decays.
func ExampleIntegrateHeat1D() {
	const (
		n     = 49
		kappa = 0.1
		dx    = 1.0 / 50
	)
	u0 := make([]float64, n)
	for i := range n {
		u0[i] = math.Sin(math.Pi * float64(i+1) * dx)
	}
	u0Arr, _ := tensor.FromFloats(u0, n)
	states, err := tensor.IntegrateHeat1D(u0Arr, kappa, dx, 1.0, 0.001, 2, 0, 0)
	if err != nil {
		log.Fatal(err)
	}
	last := (states.Shape()[0] - 1) * n
	// The mode's centre started at 1 and decays by exp(−κπ²t).
	decay := math.Exp(-kappa * math.Pi * math.Pi)
	fmt.Printf("decay %.4f, centre %.4f\n", decay, states.FloatAt(last+n/2))
	// Output: decay 0.3727, centre 0.3728
}

// The Lanczos eigensolver on a sparse symmetric matrix.
func ExampleSpEigen() {
	dense, _ := tensor.FromFloats([]float64{2, 1, 1, 2}, 2, 2)
	sp, _ := tensor.SparseFrom(dense)
	vals, _, err := tensor.SpEigen(sp, 2, tensor.NewGenerator(7))
	if err != nil {
		log.Fatal(err)
	}
	fmt.Printf("%.6f %.6f\n", vals.FloatAt(0), vals.FloatAt(1))
	// Output: 3.000000 1.000000
}

// The complete elliptic integral of the first kind by the AGM.
func ExampleEllipticK() {
	m, _ := tensor.FromFloats([]float64{0, 0.5}, 2)
	k, err := tensor.EllipticK(m)
	if err != nil {
		log.Fatal(err)
	}
	fmt.Printf("%.10f %.10f\n", k.FloatAt(0), k.FloatAt(1))
	// Output: 1.5707963268 1.8540746773
}

// Haar wavelets turn an off-grid step into one detail coefficient,
// and the inverse restores the signal exactly.
func ExampleDWT() {
	vals := make([]float64, 16)
	for i := 5; i < 16; i++ {
		vals[i] = 1
	}
	x, _ := tensor.FromFloats(vals, 16)
	c, err := tensor.DWT(x, 1)
	if err != nil {
		log.Fatal(err)
	}
	restored, _ := tensor.IDWT(c, 1)
	ok := math.Abs(restored.FloatAt(7)-1) < 1e-12 && restored.FloatAt(4) == 0
	fmt.Printf("detail[2] = %.4f, restored: %v\n", c.FloatAt(10), ok)
	// Output: detail[2] = -0.7071, restored: true
}

// The autocorrelation of a periodic signal is itself periodic.
func ExampleAutocorrelate() {
	vals := make([]float64, 64)
	for i := range vals {
		vals[i] = math.Cos(2 * math.Pi * float64(i) / 16)
	}
	x, _ := tensor.FromFloats(vals, 64)
	acf, err := tensor.Autocorrelate(x, 16)
	if err != nil {
		log.Fatal(err)
	}
	// The biased normalisation tapers the tail: 48/64 of the pairs
	// overlap at lag 16.
	fmt.Printf("acf(0) = %.4f, acf(16) = %.4f\n", acf.FloatAt(0), acf.FloatAt(16))
	// Output: acf(0) = 1.0000, acf(16) = 0.7500
}

// Differential evolution crosses Rastrigin's minefield of local
// minima to the global one.
func ExampleMinimiseDifferentialEvolution() {
	lower, _ := tensor.FromFloats([]float64{-5.12, -5.12}, 2)
	upper, _ := tensor.FromFloats([]float64{5.12, 5.12}, 2)
	_, fv, err := tensor.MinimiseDifferentialEvolution(func(a *tensor.Array) (float64, error) {
		s := 0.0
		for i := range 2 {
			z := a.FloatAt(i)
			s += z*z - 10*math.Cos(2*math.Pi*z) + 10
		}
		return s, nil
	}, lower, upper, tensor.DifferentialEvolutionOptions{Generations: 800})
	if err != nil {
		log.Fatal(err)
	}
	fmt.Printf("%.1e\n", fv)
	// Output: 0.0e+00
}

// Slicing a contiguous selection returns a read-only view sharing the
// storage; an interior range is copied.
func ExampleSlice() {
	m, _ := tensor.FromFloats([]float64{
		1, 2, 3, 4,
		5, 6, 7, 8,
		9, 10, 11, 12,
	}, 3, 4)
	// Whole rows: a view.
	rows, err := tensor.Slice(m, 0, 1, 3)
	if err != nil {
		log.Fatal(err)
	}
	// An interior column range: a copy.
	block, _ := tensor.Slice(rows, 1, 1, 3)
	fmt.Println(rows.Shape(), rows.FloatAt(0), block.Shape(), block.FloatAt(0))
	// Output: [2 4] 5 [2 2] 6
}

// A size-1 dimension replicates over the target, and anything that
// would not is an error rather than a silent broadcast.
func ExampleBroadcastTo() {
	col, _ := tensor.FromFloats([]float64{1, 2, 3}, 3, 1)
	wide, err := tensor.BroadcastTo(col, 3, 4)
	if err != nil {
		log.Fatal(err)
	}
	_, bad := tensor.BroadcastTo(col, 2, 4)
	fmt.Println(wide)
	fmt.Println("refused:", bad != nil)
	// Output: float (3, 4) [1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3]
	// refused: true
}

// The gamma function, including the negative half line.
func ExampleGamma() {
	x, _ := tensor.FromFloats([]float64{-0.5, 0.5, 5}, 3)
	g, err := tensor.Gamma(x)
	if err != nil {
		log.Fatal(err)
	}
	fmt.Printf("%.7f %.7f %.1f\n", g.FloatAt(0), g.FloatAt(1), g.FloatAt(2))
	// Output: -3.5449077 1.7724539 24.0
}

// Linear interpolation clamps outside the sampled range; the monotone
// cubic passes through the same knots with a bounded slope.
func ExampleInterpolate() {
	xs, _ := tensor.FromFloats([]float64{0, 1, 3}, 3)
	ys, _ := tensor.FromFloats([]float64{0, 2, 2.5}, 3)
	query, _ := tensor.FromFloats([]float64{0.5, 2, 5}, 3)
	lin, err := tensor.Interpolate(xs, ys, query)
	if err != nil {
		log.Fatal(err)
	}
	mono, _ := tensor.InterpolateMonotone(xs, ys, query)
	// Both stay inside the bracketing samples, and the tail clamps to
	// the last knot.
	fmt.Printf("linear %.4f %.4f %.4f\n", lin.FloatAt(0), lin.FloatAt(1), lin.FloatAt(2))
	fmt.Printf("pchip  %.4f %.4f %.4f\n", mono.FloatAt(0), mono.FloatAt(1), mono.FloatAt(2))
	// Output: linear 1.0000 2.2500 2.5000
	// pchip  1.2621 2.3716 2.5000
}

// The Halton sequence stratifies progressively in every base; the
// index-zero origin is dropped, as it carries no information.
func ExampleHaltonPoints() {
	pts, err := tensor.HaltonPoints(3, 2, 0)
	if err != nil {
		log.Fatal(err)
	}
	fmt.Println(pts.Shape())
	for i := range 3 {
		fmt.Printf("(%.4f, %.4f) ", pts.FloatAt(2*i), pts.FloatAt(2*i+1))
	}
	fmt.Println()
	// Output: [3 2]
	// (0.5000, 0.3333) (0.2500, 0.6667) (0.7500, 0.1111)
}

// A fixed seed replays the same stream, which is what makes a random
// workflow testable.
func ExampleNewGenerator() {
	u, err := tensor.Floats(tensor.NewGenerator(7), 3)
	if err != nil {
		log.Fatal(err)
	}
	again, _ := tensor.Floats(tensor.NewGenerator(7), 3)
	fmt.Printf("%.6f %.6f %.6f, equal: %v\n",
		u.FloatAt(0), u.FloatAt(1), u.FloatAt(2), tensor.Equal(u, again))
	// Output: 0.055360 0.172116 0.717576, equal: true
}