"""Reproduce the five-by-five Hilbert mixed-precision refinement example."""
import json
import numpy as np
import scipy
from scipy import linalg


def run():
    h5 = linalg.hilbert(5)
    x_ref = np.ones(5, dtype=np.float64)
    b = h5 @ x_ref
    factors = linalg.lu_factor(h5.astype(np.float32))
    x = linalg.lu_solve(factors, b.astype(np.float32)).astype(np.float64)
    rows = []
    for k in range(5):
        residual = b - h5 @ x
        rows.append({
            "iteration": k,
            "relative_residual_inf": float(linalg.norm(residual, np.inf) / linalg.norm(b, np.inf)),
            "relative_reference_error_inf": float(linalg.norm(x - x_ref, np.inf) / linalg.norm(x_ref, np.inf)),
        })
        if k < 4:
            correction = linalg.lu_solve(factors, residual.astype(np.float32))
            x += correction.astype(np.float64)
    return {
        "numpy": np.__version__,
        "scipy": scipy.__version__,
        "condition_2": float(np.linalg.cond(h5)),
        "rows": rows,
    }


if __name__ == "__main__":
    print(json.dumps(run(), indent=2))