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

package signal

import (
	"math"
	"testing"

	"sourcedock.dev/petrbalvin/tensor/internal/core"
)

// The spectral Poisson solves read their source through poissonFloats,
// which must materialise a rebased view in element order. A view's
// payload is longer than its element count and starts past offset zero,
// so a raw payload read would solve a different problem: this pin
// compares a 2-D view answer against the same values in a fresh array,
// bit for bit, on both solvers.
func TestSolvePoissonViewsMatchFresh(t *testing.T) {
	const rows, cols = 9, 9
	win := make([]float64, rows*cols)
	for i := range win {
		win[i] = math.Sin(float64(i)*0.37)*0.5 + float64(i%7)*0.1 - 0.3
	}
	// Neumann's compatibility gate measures the source with the
	// trapezoidal rule: interior weight 1, edge weight 1/2, corner
	// weight 1/4. Subtract the trapezoidal mean so the window's measure
	// is exactly zero.
	weights := make([]float64, rows*cols)
	wSum := 0.0
	for r := range rows {
		for c := range cols {
			w := 1.0
			if r == 0 || r == rows-1 {
				w /= 2
			}
			if c == 0 || c == cols-1 {
				w /= 2
			}
			weights[r*cols+c] = w
			wSum += w
		}
	}
	mean := 0.0
	for i := range win {
		mean += weights[i] * win[i]
	}
	mean /= wSum
	for i := range win {
		win[i] -= mean
	}
	padded := make([]float64, (rows+2)*(cols+2))
	for r := range rows {
		copy(padded[(r+1)*(cols+2)+1:(r+1)*(cols+2)+1+cols], win[r*cols:(r+1)*cols])
	}
	back, err := core.FromFloats(padded, rows+2, cols+2)
	if err != nil {
		t.Fatal(err)
	}
	v1, err := core.Slice(back, 0, 1, rows+1)
	if err != nil {
		t.Fatal(err)
	}
	view, err := core.Slice(v1, 1, 1, cols+1)
	if err != nil {
		t.Fatal(err)
	}
	fresh, err := core.FromFloats(win, rows, cols)
	if err != nil {
		t.Fatal(err)
	}
	solvers := map[string]func(*core.Array) (*core.Array, error){
		"SolvePoissonDirichlet": func(a *core.Array) (*core.Array, error) {
			return SolvePoissonDirichlet(a, 1, 1)
		},
		"SolvePoissonNeumann": func(a *core.Array) (*core.Array, error) {
			return SolvePoissonNeumann(a, 1, 1)
		},
	}
	for name, solve := range solvers {
		gv, err := solve(view)
		if err != nil {
			t.Fatalf("%s(view): %v", name, err)
		}
		gf, err := solve(fresh)
		if err != nil {
			t.Fatalf("%s(fresh): %v", name, err)
		}
		if gv.Len() != gf.Len() {
			t.Fatalf("%s: view length %d, fresh %d", name, gv.Len(), gf.Len())
		}
		for i := range gv.Len() {
			if math.Float64bits(gv.FloatAt(i)) != math.Float64bits(gf.FloatAt(i)) {
				t.Fatalf("%s element %d: view %#x, fresh %#x",
					name, i, math.Float64bits(gv.FloatAt(i)), math.Float64bits(gf.FloatAt(i)))
			}
		}
	}
}