Files
tensor/README.md
T
petrbalvin af4ee19703
Release / gates (push) Successful in 4m38s
Test / test (push) Successful in 5m16s
Release / release (push) Successful in 35s
feat: initial release
Assisted-by: GLM 5.3 Flash
2026-09-03 10:00:00 +02:00

25 KiB
Raw Permalink Blame History

Tensor

Tensor is a scientific computing library for Go: n-dimensional arrays, dense and sparse linear algebra, differential equation solvers, quadrature, statistics, signal transforms, optimisation, deterministic scientific charts and a reverse-mode differentiable core, built for the natural sciences: cosmology, astronomy, quantum, particle and nuclear physics, condensed matter and materials science, chemistry, biology and genetics. It is pure Go with no cgo, no GPU stack and no third-party dependency, and it is deterministic by contract: parallel kernels reduce in a fixed order, so the same program with the same seed produces bit-identical output run to run.

What Tensor optimises for

Tensor is not built to win a speed contest. It keeps no benchmark tables against other libraries, NumPy included: they are not the competition, and racing them would settle nothing. The only performance numbers in this repository compare Tensor against its own previous revisions, under the discipline in docs/BENCHMARKING.md, so that a change is judged by what it costs and never by an impression. Speed is an engineering duty, not the goal.

The goal is a tool a scientist can trust with their numbers:

  • Determinism, everywhere, without exceptions. The same program with the same seed produces bit-identical output, run to run, on any core count and under any worker setting. Parallel kernels reduce in an order derived from the data alone, never from the machine, and the generator is pinned and stable across Go releases. A number that moves between runs is not a result.
  • Portability, without variants. Pure Go on the standard library: no cgo, no third-party dependency. One portable build, no build flags, no CPU feature probes; there is no fast edition and no slow edition of the truth to keep in agreement.
  • Precision, measured rather than claimed. Narrow element types accumulate in float64, and where two algorithms round differently the choice is settled against high-precision referents, big.Float arithmetic at hundreds of bits. Every solver carries its exact-reference or residual check in the test suite.

This is also why the GPU is not the path for Tensor. GPU execution is non-deterministic by construction: the order of a reduction follows the hardware topology and the device scheduler instead of the data, and the arithmetic carries no exactness guarantee to hold, the graphics APIs themselves permitting results several ULP off. A library contracted to exact reproducibility does not hand its numbers to hardware that cannot sign for them.

Features

  • Arrays: immutable, row-major, n-dimensional arrays of int64, IEEE 754 half precision, float32, float64 and complex128 with a strict promotion ladder and loud shape errors; slicing views, gathers, scatters, sorting, Einsum, interpolation and the special functions from the gamma family to the elliptic integrals.
  • Linear algebra (linalg): dense factorisations and eigenproblems in real, complex and general arithmetic, matrix functions, regularised and truncated solves, the rank-revealing QR, sparse CSR solvers with ILU preconditioning, sparse direct Cholesky and LU with fill-reducing orderings and the rank-one update and downdate, LSQR and LSMR least squares, Lanczos and Arnoldi eigensolvers, the Krylov action of a matrix exponential, polynomial fitting and natural cubic splines.
  • Differential equations (integrate): adaptive Dormand-Prince, stiff systems from variable-step BDF2 up to variable-order BDF 1 to 5, the L-stable Rosenbrock-Wanner ROS4, index-1 differential-algebraic systems in mass-matrix form, symplectic integrators from velocity Verlet through Yoshida's fourth order to the implicit midpoint, event detection, boundary values by shooting and by adaptive Lobatto collocation, Gauss-Legendre quadrature, globally adaptive cubature, heat and wave evolution in one and two space dimensions, flux-limited advection, and finite-element Poisson solves on triangular and tetrahedral meshes.
  • Signal and transforms (signal): FFTs of any length, multi-dimensional and real-input transforms, cosine and sine transforms, the NUFFT, Welch, spectrogram and Lomb-Scargle spectra, the Hilbert envelope, decimation and resampling, Butterworth, Chebyshev, inverse Chebyshev and elliptic filter design, the window catalogue, zero-phase filtfilt, median and rank filters, Haar and Daubechies wavelets, continuous wavelets, Kalman filtering, AR and ARMA estimation, convolutions, pooling and stencils, and spectral Poisson solves.
  • Statistics (stats): CDFs, quantiles and draws for the normal, exponential, gamma, chi-square, Student t, Poisson and binomial laws; density, CDF and quantile for Weibull, lognormal and Pareto; the negative binomial PMF, CDF and quantile; the Dirichlet density, mean, mode and draws; the noncentral chi-square, F and t families through their density, CDF and quantile; histograms and rolling windows, robust descriptives, Welch's t-test, Kolmogorov-Smirnov, Mann-Whitney U, one-way ANOVA, bootstrap intervals, rank correlations and multiple-testing corrections, multivariate normals, kernel density, Gaussian processes, clustering, PCA, and linear, weighted, logistic, Poisson, lasso, elastic-net, Huber and quantile regression with the classical inference beside Theil-Sen's robust pair.
  • Optimisation (optim): Levenberg-Marquardt with an optional analytic Jacobian, L-BFGS with box bounds, linearly constrained minimisation through the augmented Lagrangian with nonlinear equality and inequality rows, Nelder-Mead, differential evolution, CMA-ES and simulated annealing, the revised simplex and active-set quadratic programming, and Brent, Newton, Broyden quasi-Newton and damped-Newton root finding.
  • Differentiable core (grad): a reverse-mode graph over the arithmetic surface, the matrix products and the transforms, complex Wirtinger differentiation, Hessians, Newton-CG on the graph, Hamiltonian Monte Carlo and adjoint ODE sensitivities.
  • Data (io): CSV, FITS images and tables, HDF5 datasets read and written (both superblock generations, chunked storage, deflate and shuffle filters, attributes and string data), the NetCDF classic model, and memory-mapped arrays that let a data cube far larger than RAM open instantly.
  • Worked examples: thirteen runnable workflows in examples/, each solving its problem end to end: ODE parameter fitting by adjoint sensitivities, PSF deconvolution, Hamiltonian Monte Carlo sampling, spectral analysis, wavelet denoising, the exact pendulum period through the elliptic integral, a Helmholtz system on the complex sparse solvers, quasi-Monte Carlo integration, heat and wave evolution, regression inference, a FITS star field, a NetCDF climate round trip, and an FFT tour.

Experimental distributed computing

This is an experiment. The spmd package carries a revolutionary technological concept: one program running on many ranks, in one process or across machines, where the order of a reduction is a function of the data alone, so a sharded reduction carries the single-array reduction's exact bits at any world size. The concept is not fully verified and remains the subject of research. The single-machine library is the settled product; the distributed surface is new, its behaviour and performance on real networks are not yet measured, and it is expected to change as the research moves.

What the experiment carries today, and what already holds: the movement collectives (Broadcast, Scatter, Gather, AllGather) and two families of reductions whose answers the test suite, the oracle digests and the loopback cluster tests in CI pin against the single-array reductions they must equal. What the suite cannot reach, the bandwidth and behaviour of real networks above all, is exactly what the research is for. Read the contract in docs/API.md before you build on it, and treat its edges as open questions rather than finished answers.

The shortest form of the experiment, four ranks over one million values, with the sharded sum equal to the single-array sum bit for bit:

package main

import (
	"fmt"

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

func main() {
	err := spmd.Launch(4, func(w *spmd.World) error {
		const n = 1000000
		vals := make([]float64, n)
		for i := range vals {
			vals[i] = float64((i*7919)%2001-1000) / 7.0
		}
		whole, _ := tensor.FromFloats(vals, n)
		span, err := spmd.Partition(n, w.Size(), w.Rank())
		if err != nil {
			return err
		}
		local, err := tensor.Slice(whole, 0, span.Lo, span.Hi)
		if err != nil {
			return err
		}
		got, err := w.AllReduceShards(local, span, spmd.Sum)
		if err != nil {
			return err
		}
		if w.Rank() == 0 {
			want := tensor.Sum(whole)
			fmt.Printf("size %d: sharded %v equals single-array %v: %v\n",
				w.Size(), got.Float(), want.Float(), got.Float() == want.Float())
		}
		return w.Barrier()
	})
	if err != nil {
		panic(err)
	}
}

Launch swaps for Listen and Join when the ranks are processes on different machines, and nothing else in the program moves.

Install

As a library:

go get sourcedock.dev/petrbalvin/tensor

Requires Go 1.27.1 or newer, the exact version go.mod declares.

One import covers everything: the root package re-exports the exported surface of every domain package, so tensor.SVD and linalg.SVD name the same function. A domain package may also be imported on its own for a narrow dependency graph; the general array constructors live in the root package (linalg.ArrayFromFloatsSafe is the one exported outside it), and the arrays are the same type either way, because tensor.Array is an alias for the core array, not a wrapper around it. Nothing is lost by mixing the two styles.

Quick start

A stiff relaxation with a slow forcing, integrated to the analytic answer:

package main

import (
	"fmt"
	"math"

	"sourcedock.dev/petrbalvin/tensor"
)

func main() {
	// y' = -1e5·(y - cos t): a transient of width 1e-5 under a slow
	// forcing. The adaptive BDF2 steps over the transient and then
	// follows the forcing; an explicit scheme is pinned to the
	// stability limit h < 2e-5 for the whole run.
	f := func(t float64, y *tensor.Array) (*tensor.Array, error) {
		return tensor.FromFloats([]float64{-1e5 * (y.FloatAt(0) - math.Cos(t))}, 1)
	}
	y0, _ := tensor.FromFloats([]float64{0}, 1)
	end, err := tensor.IntegrateBDF2(f, 0, 1, y0, tensor.ODEOptions{MaxSteps: 2000})
	if err != nil {
		panic(err)
	}
	exact := (1e10*math.Cos(1) + 1e5*math.Sin(1)) / (1e10 + 1)
	fmt.Printf("y(1) = %.12f, error %.2e\n", end.FloatAt(0),
		math.Abs(end.FloatAt(0)-exact))
}

The stiff solver lands on the exact value (k²·cos 1 + k·sin 1)/(k² + 1) inside a 2000-step budget where the explicit pair, pinned to its stability limit, cannot follow the run at all. It prints y(1) = 0.540310721007, error 4.83e-10.

Usage

Tensor's behaviour is a contract, not a convention:

  • Immutable arrays. No operation mutates its inputs; results are fresh arrays, and views never alias a buffer a later step could rewrite. A returned array is yours alone, which is also what makes them safe to share between goroutines.
  • Errors, not lies. Singular systems, exhausted step budgets, malformed files and impossible shapes come back as errors prefixed tensor: that say what happened. Nothing is silently truncated, clamped or filled.
  • Determinism. The generator is xoshiro256++ seeded through splitmix64, stable across Go releases, because the standard library does not promise stable output and reproducibility is the point of a seed. Parallel kernels keep their reduction order fixed, so a result does not move with the core count; SetNumCPU(n) pins the worker count for containers and small machines.
  • Scientific scope. Tensor is for the natural sciences and for nothing else. It carries nothing for artificial intelligence, economics or finance: no neural-network machinery, no training loops, no market or portfolio helpers. The convolution and pooling functions are signal-processing stencils (PSF deconvolution, image filtering), not model layers. Differentiation exists because fitting parameters to data and sensitivity analysis are scientific tools.

What follows is one short program per package, each complete and runnable as written. The full surface, option by option, is in docs/API.md, and the longer workflows are the thirteen programs in examples/.

Arrays

package main

import (
	"fmt"

	"sourcedock.dev/petrbalvin/tensor"
)

func main() {
	y, err := tensor.FromFloats([]float64{1, 2, 3, 4, 5, 6}, 2, 3)
	if err != nil {
		panic(err)
	}
	// A slice is a read-only view: no copy, and no way for a later
	// step to write through it.
	cols, _ := tensor.Slice(y, 1, 1, 3)
	// Reductions name the axis, and the axis disappears from the shape.
	rowSum, _ := tensor.SumAxis(y, 1)
	// The dtype ladder promotes on request, never implicitly downward.
	halves, _ := tensor.Astype(y, tensor.Float32)

	fmt.Println(y.Shape(), rowSum, halves.Dtype())
	fmt.Println(cols)
	fmt.Println(y)
}

Shape reports the extents, Dtype the element type, Len the element count. Element-wise operations (Add, Mul, Exp, Sqrt, the comparisons, Where) and the shape moves (Reshape, Transpose, Concat, Stack, Pad) all take and return whole arrays, so a formula reads as one expression per line.

Linear algebra

package main

import (
	"fmt"

	"sourcedock.dev/petrbalvin/tensor"
)

func main() {
	a, _ := tensor.FromFloats([]float64{4, 1, 1, 3}, 2, 2)
	b, _ := tensor.FromFloats([]float64{1, 2}, 2)

	// One LU with partial pivoting; a singular matrix is an error.
	x, err := tensor.Solve(a, b)
	if err != nil {
		panic(err)
	}
	// The symmetric eigenproblem returns ascending values and the
	// orthonormal eigenvectors as columns.
	values, vectors, _ := tensor.Eigen(a)
	// The Cholesky factor for reuse across right-hand sides.
	l, _ := tensor.Cholesky(a)

	fmt.Println(x, values)
	fmt.Println(l, vectors.Shape())
}

Sparse systems go through the same array type: SparseFrom builds the COO form, CSRFromCOO and CSCFromCOO the compressed views, NewSparseCholesky and NewSparseLU the direct factorisations under a fill-reducing ordering, NewSparseILU the preconditioner, and SpSolve, SpSolveBiCGSTAB, SpLSQR and SpEigen the iterative solvers.

Differential equations

package main

import (
	"fmt"
	"math"

	"sourcedock.dev/petrbalvin/tensor"
)

func main() {
	// The harmonic oscillator as a first-order system: y = (position,
	// velocity), so y' = (velocity, -position).
	f := func(t float64, y *tensor.Array) (*tensor.Array, error) {
		return tensor.FromFloats([]float64{y.FloatAt(1), -y.FloatAt(0)}, 2)
	}
	y0, _ := tensor.FromFloats([]float64{0, 1}, 2)

	// Event detection is a watch on the trajectory, so the crossing
	// time comes from the interpolant rather than from the step grid.
	hits, end, err := tensor.IntegrateODEEvents(f, 0, 4, y0,
		[]tensor.ODEWatch{{
			Function: func(t float64, y *tensor.Array) (float64, error) {
				return y.FloatAt(0), nil
			},
			Direction: -1, // falling crossings only
		}}, tensor.ODEOptions{})
	if err != nil {
		panic(err)
	}
	fmt.Printf("first minimum at t = %.6f, want pi = %.6f\n", hits[0].Time, math.Pi)
	fmt.Printf("y(4) = %.6f\n", end.FloatAt(0))
}

IntegrateODE is the adaptive Dormand-Prince pair, IntegrateBDF2 and IntegrateBDFVar the stiff routes, IntegrateROS4 the L-stable one, IntegrateDAE the mass-matrix form. Quadrature is IntegrateFunction and IntegrateND, the boundary value problems are IntegrateBoundary and SolveBoundaryCollocation, and the turnkey time steppers are the IntegrateHeat1D/2D, IntegrateWave1D/2D and advection families.

Transforms, filters and spectra

package main

import (
	"fmt"
	"math"

	"sourcedock.dev/petrbalvin/tensor"
)

func main() {
	const (
		fs = 1000.0
		n  = 1000
	)
	samples := make([]float64, n)
	for i := range n {
		t := float64(i) / fs
		samples[i] = math.Sin(2*math.Pi*50*t) + 0.25*math.Sin(2*math.Pi*120*t)
	}
	x, _ := tensor.FromFloats(samples, n)

	// Welch's averaged periodogram: windowed segments, one-sided.
	freqs, psd, err := tensor.WelchPSD(x, fs, 256, 128, "hann")
	if err != nil {
		panic(err)
	}
	// Bin 0 is the mean; the peak above it is the tone.
	rest, _ := tensor.Slice(psd, 0, 1, psd.Len())
	idx, _ := tensor.ArgMax(rest)
	fmt.Printf("peak at %.1f Hz, bin spacing %.1f Hz\n",
		freqs.FloatAt(idx+1), freqs.FloatAt(1))

	// A Butterworth design plus filtfilt: zero phase, so no lag.
	b, a, _ := tensor.ButterworthLowPass(4, fs, 60)
	clean, _ := tensor.Filtfilt(b, a, x)
	peakIn, _ := tensor.Max(x)
	peakOut, _ := tensor.Max(clean)
	fmt.Printf("input peak %.3f, filtered peak %.3f\n", peakIn.Float(), peakOut.Float())
}

The Fourier family covers FFT/IFFT of any length, FFT2, FFT3, FFTN, the real-input RFFT/IRFFT, the cosine and sine transforms, the Haar and Daubechies wavelets, the analytic CWT, and the STFT, spectrogram and Lomb-Scargle estimators beside Welch.

Statistics

package main

import (
	"fmt"

	"sourcedock.dev/petrbalvin/tensor"
)

func main() {
	// y = 2 + 3x with noise, 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)

	fit, err := tensor.LinearRegression(design, y)
	if err != nil {
		panic(err)
	}
	fmt.Printf("slope %.3f ± %.3f, t = %.2f, p = %.2g, R2 = %.4f\n",
		fit.Coefficients[1], fit.StandardErrors[1],
		fit.TStatistics[1], fit.PValues[1], fit.RSquared)

	// Distributions answer one call each, with the tail the caller asks for.
	fmt.Printf("P(Z <= 1.96) = %.4f, t(10) 97.5%% = %.4f\n",
		tensor.NormalCDF(1.96), mustQuantile(tensor.StudentTQuantile(0.975, 10)))
}

// Quantiles can fail on a parameter outside their domain, so the
// example carries the error rather than dropping it.
func mustQuantile(q float64, err error) float64 {
	if err != nil {
		panic(err)
	}
	return q
}

The models are LinearRegression, WeightedLinearRegression, LogisticRegression, PoissonRegression, Lasso, ElasticNet, HuberRegression and QuantileRegression, each returning its coefficient table with the standard errors and p-values beside it; TheilSenRegression returns the robust intercept and slope as a pair. The tests are WelchTTest, KolmogorovSmirnovTest, MannWhitneyU, ANOVAOneWay and ChiSquareGoodnessOfFit; the multivariate tools are PCA, KMeans, GaussianMixture and GaussianProcessRegression.

Optimisation

package main

import (
	"fmt"
	"math"

	"sourcedock.dev/petrbalvin/tensor"
)

func main() {
	// Fit A·exp(−k·t) to six noisy observations by least squares.
	ts := []float64{0, 1, 2, 3, 4, 5}
	obs := []float64{2.0, 1.22, 0.74, 0.45, 0.27, 0.17}
	residual := func(p *tensor.Array) (*tensor.Array, error) {
		out := make([]float64, len(ts))
		for i, t := range ts {
			out[i] = p.FloatAt(0)*math.Exp(-p.FloatAt(1)*t) - obs[i]
		}
		return tensor.FromFloats(out, len(out))
	}
	p0, _ := tensor.FromFloats([]float64{1, 0.5}, 2)

	p, ss, err := tensor.LevenbergMarquardt(residual, p0, tensor.LMOptions{})
	if err != nil {
		panic(err)
	}
	fmt.Printf("A = %.4f, k = %.4f, sum of squares %.3e\n",
		p.FloatAt(0), p.FloatAt(1), ss)
}

MinimiseLBFGS with optional box walls, Minimise (Nelder-Mead), MinimiseConstrained for the linear rows, MinimiseNonlinearConstrained for constraint functions, MinimiseDifferentialEvolution, MinimiseCMAES and MinimiseSimulatedAnnealing for the multimodal landscapes, MinimiseLinear and MinimiseQP for the programs, and FindRoot, FindRootNewton and FindRootSystem for the roots. A solver that runs out of budget refuses rather than reporting a converged answer, and AllowBudgetExit opts into the best point instead; MinimiseDifferentialEvolution is the exception, returning its best generation without an error.

Data in and out

package main

import (
	"fmt"
	"os"
	"path/filepath"

	"sourcedock.dev/petrbalvin/tensor"
)

func main() {
	dir, err := os.MkdirTemp("", "tensor-readme")
	if err != nil {
		panic(err)
	}
	defer os.RemoveAll(dir)

	scan, _ := tensor.FromFloats([]float64{1, 2, 3, 4}, 2, 2)
	path := filepath.Join(dir, "field.h5")
	err = tensor.SaveHDF5(path, []tensor.HDF5Dataset{{
		Path:   "/scan/temperature",
		Values: scan,
		Attrs:  map[string]string{"units": "K"},
	}}, nil, tensor.HDF5WriteOptions{Gzip: 6, Shuffle: true})
	if err != nil {
		panic(err)
	}

	sets, err := tensor.LoadHDF5(path)
	if err != nil {
		panic(err)
	}
	fmt.Println(sets[0].Path, sets[0].Shape, sets[0].Attrs["units"], sets[0].Values)
}

LoadCSV/SaveCSV handle the tabular case, LoadFITS/SaveFITS images and LoadFITSTable/SaveFITSTable the table extensions, LoadNetCDF/SaveNetCDF the classic model, and MapFloats, MapFloat32s and MapInts open a native-endian file as a read-only array without reading it; SaveNativeFloats writes the float64 form MapFloats reads.

Charts

package main

import (
	"fmt"
	"math"

	"sourcedock.dev/petrbalvin/tensor"
)

func main() {
	wavelengths := make([]float64, 200)
	intensity := make([]float64, 200)
	for i := range wavelengths {
		wavelengths[i] = 400 + 2*float64(i)
		intensity[i] = 100 + 40*math.Exp(-math.Pow(wavelengths[i]-589, 2)/25)
	}
	xs, _ := tensor.FromFloats(wavelengths, 200)
	ys, _ := tensor.FromFloats(intensity, 200)
	series, _ := tensor.Line("sodium D line", xs, ys)
	chart := tensor.Chart{
		Title:  "Absorption spectrum",
		XLabel: "wavelength [nm]", YLabel: "intensity",
		Series: []tensor.Series{series},
	}
	if err := chart.WriteSVG("spectrum.svg"); err != nil {
		panic(err)
	}
	fmt.Println("spectrum.svg written")
}

Deterministic SVG line charts: linear axes, five ticks each, one legend line per series, and a byte-identical file on every run, so a figure in a paper is compared exactly like any other computed number. The package is small by intent; it draws the figures, it does not stage a cinema.

Automatic differentiation

package main

import (
	"fmt"

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

func main() {
	x, _ := grad.FromFloat64s([]float64{1, 2, 3}, true, 3)
	w, _ := grad.FromFloat64s([]float64{0.5, -1, 2}, true, 3)

	prod, _ := x.Mul(w)
	sq, _ := prod.Pow(2)
	loss, _ := sq.Sum()
	if err := loss.Backward(); err != nil {
		panic(err)
	}
	// dL/dx = 2·x·w² and dL/dw = 2·x²·w, exact to rounding.
	fmt.Println(x.Grad(), w.Grad())

	// Second-order questions come from the same graph: H·v in two
	// gradient evaluations, the tool Newton-CG scales on. The Hessian
	// of Σz² is 2·I, so H·x is 2·x.
	quadratic := func(z *grad.Tensor) (*grad.Tensor, error) {
		sq, err := z.Pow(2)
		if err != nil {
			return nil, err
		}
		return sq.Sum()
	}
	hv, err := grad.HessianVectorProduct(quadratic, x, x, grad.HessianOptions{})
	if err != nil {
		panic(err)
	}
	fmt.Println(hv)
}

MinimiseNewtonCG minimises a graph function with truncated-CG Newton steps, SampleHMC runs Hamiltonian Monte Carlo on any differentiable unnormalised density, and AdjointODE differentiates an ODE solution at the cost of one extra solve. Complex graphs follow the Wirtinger convention, with Real, Imag, Conj and Abs2 bridging into a real loss.

at the cost of one extra solve. Complex graphs follow the Wirtinger convention, with Real, Imag, Conj and Abs2 bridging into a real loss.

Where to look next

Development

just build   # compile everything, examples included
just test    # the suite with the coverage floor
just gates   # the definition of done: build, fmt-check, vet, test, race

CI (Gitea Actions) enforces the same gates on every push to development, race excepted: the race detector is dispatched by hand. See docs/DEVELOPMENT.md for the full workflow and CONTRIBUTING.md for how to contribute.

Documentation

Licence

MIT. See LICENSE for the text.

Copyright © 2026 Petr Balvín