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

// Command qmc compares quasi-random integration against plain Monte
// Carlo on the same two-dimensional integral. Sobol points are a
// digital lattice: every block of 2^m points stratifies each
// coordinate exactly, so the error decays far faster than the
// 1/sqrt(n) of random sampling, and Halton sits in between.
//
// Usage: go run ./examples/qmc
package main

import (
	"fmt"
	"log"
	"math"

	"sourcedock.dev/petrbalvin/tensor"
)

// f is the integrand: smooth, with its curvature spread over the unit
// square. The exact value is (1-e^-1)*sqrt(pi)/2*erf(1), the product
// of the x integral and the error function integral over y.
func f(x, y float64) float64 { return math.Exp(-x - y*y) }

const exact = 0.4720828881800443 // (1-e^-1)*sqrt(pi)/2*erf(1)

// estimate integrates f over [0,1]^2 from an (n,2) point set.
func estimate(pts *tensor.Array, n int) float64 {
	s := 0.0
	for i := range n {
		x, err := tensor.FloatAt(pts, i, 0)
		if err != nil {
			log.Fatal(err)
		}
		y, err := tensor.FloatAt(pts, i, 1)
		if err != nil {
			log.Fatal(err)
		}
		s += f(x, y)
	}
	return s / float64(n)
}

func main() {
	fmt.Printf("integral of exp(-x - y^2) over the unit square, exact %.10f\n\n", exact)
	fmt.Println("  points     Monte Carlo      Halton        Sobol")
	for _, n := range []int{64, 256, 1024, 4096, 16384} {
		// Monte Carlo: uniform draws from the seeded generator.
		g := tensor.NewGenerator(int64(n))
		mc, err := tensor.Floats(g, 2*n)
		if err != nil {
			log.Fatal(err)
		}
		mcPts, err := tensor.Reshape(mc, n, 2)
		if err != nil {
			log.Fatal(err)
		}
		// Halton and Sobol from the first point on; Sobol skips its
		// origin point exactly as Halton does.
		hal, err := tensor.HaltonPoints(n, 2, 0)
		if err != nil {
			log.Fatal(err)
		}
		sob, err := tensor.SobolPoints(n, 2, 0)
		if err != nil {
			log.Fatal(err)
		}
		eMC := math.Abs(estimate(mcPts, n) - exact)
		eHal := math.Abs(estimate(hal, n) - exact)
		eSob := math.Abs(estimate(sob, n) - exact)
		fmt.Printf("  %6d     %.3e    %.3e    %.3e\n", n, eMC, eHal, eSob)
	}
	fmt.Println()
	fmt.Println("the quasi-random errors collapse with n; the Monte Carlo")
	fmt.Println("error only shrinks as 1/sqrt(n) and stays noisy on top")
}