#pragma once

// Include after importing mean_field. These maps change the discrete shape
// space; they are not an invertible preconditioner of the original equations.
#include <algorithm>
#include <cmath>
#include <concepts>
#include <cstdint>
#include <limits>
#include <stdexcept>
#include <string>
#include <type_traits>
#include <vector>

#include <mfem.hpp>
#include <mpi.h>

namespace experiment::plateau_diagnostics {
    struct ReducedSurfaceCoordinateOptions final {
        double relativeRankTolerance{1.e-12};
        double orthogonalityTolerance{5.e-12};
    };

    struct ReducedSurfaceCoordinateReport final {
        double stateOrthogonalityError{0.0};
        double residualOrthogonalityError{0.0};
        double stateReconstructionError{0.0};
        double residualReconstructionError{0.0};
        double stateMinimumIndependentFraction{1.0};
        double residualMinimumIndependentFraction{1.0};
    };

    namespace reduced_surface_detail {
        inline void RequireFinite(const mfem::Vector &value, const char *description) {
            for (int i = 0; i < value.Size(); ++i)
                if (!std::isfinite(value[i])) throw std::invalid_argument(description);
        }

        inline double StableNorm(const mfem::Vector &value) {
            double norm = 0.0;
            for (int i = 0; i < value.Size(); ++i) norm = std::hypot(norm, value[i]);
            return norm;
        }

        inline double DotColumn(const mfem::DenseMatrix &matrix, int column, const mfem::Vector &value) {
            long double result = 0.0L;
            for (int i = 0; i < matrix.Height(); ++i)
                result += static_cast<long double>(matrix(i, column)) * value[i];
            return static_cast<double>(result);
        }

        struct QRReport final {
            double orthogonalityError{0.0};
            double reconstructionError{0.0};
            double minimumIndependentFraction{1.0};
        };

        // Reorthogonalization preserves the advertised coordinate identities
        // without a new dense linear-algebra dependency. Rejecting a dependent
        // mode keeps the scientific degree/mode ordering explicit in artifacts.
        inline QRReport ThinQR(
            const mfem::DenseMatrix &input,
            mfem::DenseMatrix &basis,
            mfem::DenseMatrix &transform,
            ReducedSurfaceCoordinateOptions options
        ) {
            const int rows = input.Height(), columns = input.Width();
            basis.SetSize(rows, columns);
            transform.SetSize(columns, columns);
            basis = 0.0;
            transform = 0.0;
            QRReport report;
            mfem::Vector column(rows);
            for (int j = 0; j < columns; ++j) {
                for (int i = 0; i < rows; ++i) column[i] = input(i, j);
                RequireFinite(column, "The reduced surface basis contains nonfinite entries.");
                const double originalNorm = StableNorm(column);
                if (!std::isfinite(originalNorm) || originalNorm <= 0.0)
                    throw std::invalid_argument("The reduced surface basis has a zero or overflowing column.");
                for (int pass = 0; pass < 2; ++pass) {
                    for (int previous = 0; previous < j; ++previous) {
                        const double projection = DotColumn(basis, previous, column);
                        if (!std::isfinite(projection))
                            throw std::invalid_argument("The reduced surface QR projection is nonfinite.");
                        transform(previous, j) += projection;
                        for (int i = 0; i < rows; ++i) column[i] -= projection * basis(i, previous);
                    }
                }
                const double independentNorm = StableNorm(column);
                const double fraction = independentNorm / originalNorm;
                if (!std::isfinite(fraction) || fraction <= options.relativeRankTolerance)
                    throw std::invalid_argument("The reduced surface basis is numerically rank deficient.");
                report.minimumIndependentFraction = std::min(report.minimumIndependentFraction, fraction);
                transform(j, j) = independentNorm;
                for (int i = 0; i < rows; ++i) basis(i, j) = column[i] / independentNorm;
            }
            double inputNorm = 0.0, errorNorm = 0.0;
            for (int j = 0; j < columns; ++j) {
                for (int l = 0; l < columns; ++l) {
                    long double product = 0.0L;
                    for (int i = 0; i < rows; ++i)
                        product += static_cast<long double>(basis(i, j)) * basis(i, l);
                    report.orthogonalityError = std::max(
                        report.orthogonalityError, std::abs(static_cast<double>(product) - (j == l ? 1.0 : 0.0))
                    );
                }
                for (int i = 0; i < rows; ++i) {
                    long double reconstructed = 0.0L;
                    for (int l = 0; l <= j; ++l)
                        reconstructed += static_cast<long double>(basis(i, l)) * transform(l, j);
                    inputNorm = std::hypot(inputNorm, input(i, j));
                    errorNorm = std::hypot(errorNorm, static_cast<double>(reconstructed) - input(i, j));
                }
            }
            report.reconstructionError = errorNorm / inputNorm;
            if (!std::isfinite(report.reconstructionError) ||
                !std::isfinite(report.orthogonalityError) ||
                report.orthogonalityError > options.orthogonalityTolerance ||
                report.reconstructionError > options.orthogonalityTolerance)
                throw std::invalid_argument("The reduced surface QR failed its coordinate consistency check.");
            return report;
        }
    } // namespace reduced_surface_detail

    class ReducedSurfaceCoordinates final {
    public:
        ReducedSurfaceCoordinates(
            const mfem::DenseMatrix &physicalBasis,
            const mfem::Vector &stateFactors,
            const mfem::Vector &residualFactors,
            int shapeValueOffset,
            int shapeResidualOffset,
            MPI_Comm communicator,
            ReducedSurfaceCoordinateOptions options = {}
        )
            : m_physicalBasis(physicalBasis), m_stateFactors(stateFactors), m_residualFactors(residualFactors),
              m_shapeValueOffset(shapeValueOffset), m_shapeResidualOffset(shapeResidualOffset) {
            int ranks = 0;
            if (MPI_Comm_size(communicator, &ranks) != MPI_SUCCESS || ranks != 1)
                throw std::invalid_argument("The reduced surface experiment currently requires one MPI rank.");
            if (!std::isfinite(options.relativeRankTolerance) || options.relativeRankTolerance <= 0.0 ||
                options.relativeRankTolerance >= 1.0 || !std::isfinite(options.orthogonalityTolerance) ||
                options.orthogonalityTolerance <= 0.0 || options.orthogonalityTolerance >= 1.0)
                throw std::invalid_argument("Invalid reduced surface coordinate tolerances.");
            const int size = stateFactors.Size(), shape = physicalBasis.Height(), modes = physicalBasis.Width();
            if (size != residualFactors.Size() || shape <= 0 || modes <= 0 || modes > shape || shape > size ||
                shapeValueOffset < 0 || shapeValueOffset > size - shape ||
                shapeResidualOffset < 0 || shapeResidualOffset > size - shape)
                throw std::invalid_argument("Invalid reduced surface dimensions or distinct block offsets.");
            for (int i = 0; i < size; ++i)
                if (!std::isfinite(stateFactors[i]) || stateFactors[i] <= 0.0 ||
                    !std::isfinite(residualFactors[i]) || residualFactors[i] <= 0.0)
                    throw std::invalid_argument("Reduced surface normalization factors must be finite and positive.");

            mfem::DenseMatrix stateImage(shape, modes), residualDualImage(shape, modes);
            for (int j = 0; j < modes; ++j) {
                for (int i = 0; i < shape; ++i) {
                    stateImage(i, j) = stateFactors[shapeValueOffset + i] * physicalBasis(i, j);
                    residualDualImage(i, j) = physicalBasis(i, j) / residualFactors[shapeResidualOffset + i];
                }
            }
            // State and residual scales act on different layouts and are not
            // interchangeable. Q_r^T Fhat = T_r^{-T} B^T Fphysical preserves
            // precisely the physical Galerkin equations with an invertible
            // reduced residual scaling, even when Q_x differs from Q_r.
            const auto state = reduced_surface_detail::ThinQR(stateImage, m_stateBasis, m_stateTransform, options);
            const auto residual = reduced_surface_detail::ThinQR(
                residualDualImage, m_residualBasis, m_residualTransform, options
            );
            m_report = {
                state.orthogonalityError, residual.orthogonalityError,
                state.reconstructionError, residual.reconstructionError,
                state.minimumIndependentFraction, residual.minimumIndependentFraction
            };
        }

        [[nodiscard]] int FullSize() const noexcept { return m_stateFactors.Size(); }
        [[nodiscard]] int ReducedSize() const noexcept { return FullSize() - ShapeSize() + ModeCount(); }
        [[nodiscard]] int ShapeSize() const noexcept { return m_physicalBasis.Height(); }
        [[nodiscard]] int ModeCount() const noexcept { return m_physicalBasis.Width(); }
        [[nodiscard]] int ShapeValueOffset() const noexcept { return m_shapeValueOffset; }
        [[nodiscard]] int ShapeResidualOffset() const noexcept { return m_shapeResidualOffset; }
        [[nodiscard]] const mfem::DenseMatrix &PhysicalBasis() const noexcept { return m_physicalBasis; }
        [[nodiscard]] const mfem::DenseMatrix &StateBasis() const noexcept { return m_stateBasis; }
        [[nodiscard]] const mfem::DenseMatrix &ResidualBasis() const noexcept { return m_residualBasis; }
        [[nodiscard]] const mfem::DenseMatrix &StateTransform() const noexcept { return m_stateTransform; }
        [[nodiscard]] const mfem::DenseMatrix &ResidualTransform() const noexcept { return m_residualTransform; }
        [[nodiscard]] const ReducedSurfaceCoordinateReport &Report() const noexcept { return m_report; }

        [[nodiscard]] bool MatchesNormalization(
            const mfem::Vector &stateFactors, const mfem::Vector &residualFactors
        ) const noexcept {
            if (stateFactors.Size() != FullSize() || residualFactors.Size() != FullSize()) return false;
            for (int i = 0; i < FullSize(); ++i)
                if (stateFactors[i] != m_stateFactors[i] || residualFactors[i] != m_residualFactors[i]) return false;
            return true;
        }

        void InjectState(const mfem::Vector &reduced, mfem::Vector &full) const {
            Inject(reduced, full, m_shapeValueOffset, m_stateBasis);
        }
        void RestrictState(const mfem::Vector &full, mfem::Vector &reduced) const {
            Restrict(full, reduced, m_shapeValueOffset, m_stateBasis);
        }
        void RestrictResidual(const mfem::Vector &full, mfem::Vector &reduced) const {
            Restrict(full, reduced, m_shapeResidualOffset, m_residualBasis);
        }
        void InjectResidual(const mfem::Vector &reduced, mfem::Vector &full) const {
            Inject(reduced, full, m_shapeResidualOffset, m_residualBasis);
        }

    private:
        void Inject(
            const mfem::Vector &reduced, mfem::Vector &full, int offset, const mfem::DenseMatrix &basis
        ) const {
            if (reduced.Size() != ReducedSize() || (&reduced == &full))
                throw std::invalid_argument("Reduced surface injection requires distinct correctly sized vectors.");
            reduced_surface_detail::RequireFinite(reduced, "Cannot inject a nonfinite reduced surface vector.");
            full.SetSize(FullSize());
            for (int i = 0; i < offset; ++i) full[i] = reduced[i];
            for (int i = 0; i < ShapeSize(); ++i) {
                long double value = 0.0L;
                for (int j = 0; j < ModeCount(); ++j)
                    value += static_cast<long double>(basis(i, j)) * reduced[offset + j];
                full[offset + i] = static_cast<double>(value);
            }
            for (int i = offset + ShapeSize(); i < FullSize(); ++i)
                full[i] = reduced[i - ShapeSize() + ModeCount()];
            reduced_surface_detail::RequireFinite(full, "Reduced surface injection overflowed.");
        }

        void Restrict(
            const mfem::Vector &full, mfem::Vector &reduced, int offset, const mfem::DenseMatrix &basis
        ) const {
            if (full.Size() != FullSize() || (&full == &reduced))
                throw std::invalid_argument("Reduced surface restriction requires distinct correctly sized vectors.");
            reduced_surface_detail::RequireFinite(full, "Cannot restrict a nonfinite full surface vector.");
            reduced.SetSize(ReducedSize());
            for (int i = 0; i < offset; ++i) reduced[i] = full[i];
            for (int j = 0; j < ModeCount(); ++j) {
                long double value = 0.0L;
                for (int i = 0; i < ShapeSize(); ++i)
                    value += static_cast<long double>(basis(i, j)) * full[offset + i];
                reduced[offset + j] = static_cast<double>(value);
            }
            for (int i = offset + ShapeSize(); i < FullSize(); ++i)
                reduced[i - ShapeSize() + ModeCount()] = full[i];
            reduced_surface_detail::RequireFinite(reduced, "Reduced surface restriction overflowed.");
        }

        mfem::DenseMatrix m_physicalBasis, m_stateBasis, m_residualBasis, m_stateTransform, m_residualTransform;
        mfem::Vector m_stateFactors, m_residualFactors;
        int m_shapeValueOffset, m_shapeResidualOffset;
        ReducedSurfaceCoordinateReport m_report;
    };

    template <std::derived_from<mfem::Operator> FullOperator>
    class ReducedSurfaceOperator final : public mfem::Operator {
    public:
        ReducedSurfaceOperator(const FullOperator &operation, const ReducedSurfaceCoordinates &coordinates)
            : mfem::Operator(coordinates.ReducedSize()), m_operation(&operation), m_coordinates(&coordinates),
              m_fullDirection(coordinates.FullSize()), m_fullAction(coordinates.FullSize()) {
            Refresh();
        }

        // Refresh does not prepare the physical problem. It adopts an already
        // prepared generation while requiring the frozen coordinate scales.
        void Refresh() {
            m_current = false;
            if (m_operation->Height() != m_coordinates->FullSize() ||
                m_operation->Width() != m_coordinates->FullSize())
                throw std::invalid_argument("Reduced surface operator needs the matching square full Jacobian.");
            if constexpr (requires { m_operation->IsPrepared(); })
                if (!m_operation->IsPrepared()) throw std::invalid_argument("The full Jacobian is not prepared.");
            if constexpr (requires { m_operation->GetNormalization(); }) {
                const auto &normalization = m_operation->GetNormalization();
                if (!m_coordinates->MatchesNormalization(normalization.StateFactors(), normalization.ResidualFactors()))
                    throw std::invalid_argument("The frozen reduced surface normalization has changed.");
            }
            if constexpr (requires { m_operation->GetProblem().GetPreparationGeneration(); })
                m_generation = m_operation->GetProblem().GetPreparationGeneration();
            m_current = true;
        }

        [[nodiscard]] bool IsCurrent() const {
            if (!m_current || m_operation->Height() != m_coordinates->FullSize() ||
                m_operation->Width() != m_coordinates->FullSize()) return false;
            if constexpr (requires { m_operation->IsPrepared(); })
                if (!m_operation->IsPrepared()) return false;
            if constexpr (requires { m_operation->GetProblem().GetPreparationGeneration(); })
                if (m_generation != m_operation->GetProblem().GetPreparationGeneration()) return false;
            if constexpr (requires { m_operation->GetNormalization(); }) {
                const auto &normalization = m_operation->GetNormalization();
                if (!m_coordinates->MatchesNormalization(normalization.StateFactors(), normalization.ResidualFactors()))
                    return false;
            }
            return true;
        }

        void Mult(const mfem::Vector &direction, mfem::Vector &action) const override {
            if (!IsCurrent()) throw std::invalid_argument("The reduced surface Jacobian is stale.");
            m_coordinates->InjectState(direction, m_fullDirection);
            m_operation->Mult(m_fullDirection, m_fullAction);
            m_coordinates->RestrictResidual(m_fullAction, action);
        }

        [[nodiscard]] const FullOperator &GetFullOperator() const noexcept { return *m_operation; }

    private:
        const FullOperator *m_operation;
        const ReducedSurfaceCoordinates *m_coordinates;
        mutable mfem::Vector m_fullDirection, m_fullAction;
        std::uint64_t m_generation{0};
        bool m_current{false};
    };

    template <class FullInverse>
        requires std::derived_from<FullInverse, mfem::Solver> &&
                 mean_field::solver::LinearPreconditionerApplicationContractAvailable<FullInverse>
    class ReducedSurfacePreconditioner final : public mfem::Solver {
    public:
        // Projecting an inverse is an approximation, not an inverse identity.
        // In particular a flexible underlying application remains flexible.
        static constexpr auto applicationContract =
            mean_field::solver::linearPreconditionerApplicationContract<FullInverse>;

        ReducedSurfacePreconditioner(FullInverse &inverse, const ReducedSurfaceCoordinates &coordinates)
            : mfem::Solver(coordinates.ReducedSize()), m_inverse(&inverse), m_coordinates(&coordinates),
              m_fullResidual(coordinates.FullSize()), m_fullCorrection(coordinates.FullSize()) {
            if (!IsCurrent()) throw std::invalid_argument("The full reduced-surface inverse is not current.");
            iterative_mode = false;
        }

        void SetOperator(const mfem::Operator &operation) override {
            if (operation.Height() != Height() || operation.Width() != Width())
                throw std::invalid_argument("The projected inverse received an incompatible reduced operator.");
            // The physical inverse is already bound to the authoritative full
            // operator. Passing this reduced operator into it would be invalid.
            if (!IsCurrent()) throw std::invalid_argument("The projected inverse is stale.");
        }

        [[nodiscard]] bool IsCurrent() const {
            if (m_inverse->Height() != m_coordinates->FullSize() ||
                m_inverse->Width() != m_coordinates->FullSize() || m_inverse->iterative_mode) return false;
            if constexpr (requires { m_inverse->IsCurrent(); })
                if (!m_inverse->IsCurrent()) return false;
            return true;
        }

        void Mult(const mfem::Vector &residual, mfem::Vector &correction) const override {
            if (!IsCurrent()) throw std::invalid_argument("The projected inverse is stale.");
            m_coordinates->InjectResidual(residual, m_fullResidual);
            m_fullCorrection = 0.0;
            m_inverse->Mult(m_fullResidual, m_fullCorrection);
            m_coordinates->RestrictState(m_fullCorrection, correction);
        }

    private:
        FullInverse *m_inverse;
        const ReducedSurfaceCoordinates *m_coordinates;
        mutable mfem::Vector m_fullResidual, m_fullCorrection;
    };

    struct ReducedSurfaceAlgebraSelfCheckReport final {
        int checks{0};
        double maximumError{0.0};
        bool rejectedDependentBasis{false};
        bool rejectedNonfiniteBasis{false};
        bool rejectedStaleOperator{false};
    };

    namespace reduced_surface_detail {
        class SelfCheckOperator final : public mfem::Operator {
        public:
            explicit SelfCheckOperator(int size) : mfem::Operator(size), matrix(size) {
                matrix = 0.0;
                for (int i = 0; i < size; ++i) {
                    matrix(i, i) = 2.0 + i;
                    for (int j = 0; j < size; ++j)
                        if (i != j) matrix(i, j) = 0.03 * (1 + i - 2 * j);
                }
            }
            void Mult(const mfem::Vector &x, mfem::Vector &y) const override {
                ++applications;
                y.SetSize(Height());
                matrix.Mult(x, y);
            }
            [[nodiscard]] bool IsPrepared() const noexcept { return true; }
            [[nodiscard]] const SelfCheckOperator &GetProblem() const noexcept { return *this; }
            [[nodiscard]] std::uint64_t GetPreparationGeneration() const noexcept { return generation; }
            mfem::DenseMatrix matrix;
            std::uint64_t generation{1};
            mutable int applications{0};
        };

        class SelfCheckInverse final : public mfem::Solver {
        public:
            static constexpr auto applicationContract =
                mean_field::preconditioning::ApplicationContract::stationary_linear;
            explicit SelfCheckInverse(int size) : mfem::Solver(size) { iterative_mode = false; }
            void SetOperator(const mfem::Operator &) override {}
            [[nodiscard]] bool IsCurrent() const noexcept { return true; }
            void Mult(const mfem::Vector &x, mfem::Vector &y) const override {
                ++applications;
                y.SetSize(Height());
                for (int i = 0; i < Height(); ++i) y[i] = x[i] / (2.0 + i);
            }
            mutable int applications{0};
        };
    } // namespace reduced_surface_detail

    // Tiny deterministic algebra only: this exercises different value/row
    // offsets and non-proportional scales without constructing any FEM state.
    [[nodiscard]] inline ReducedSurfaceAlgebraSelfCheckReport RunReducedSurfaceAlgebraSelfCheck(MPI_Comm comm) {
        ReducedSurfaceAlgebraSelfCheckReport report;
        const auto check = [&report](double error) {
            ++report.checks;
            if (!std::isfinite(error) || error > 2.e-12)
                throw std::runtime_error("Reduced surface algebra self-check failed.");
            report.maximumError = std::max(report.maximumError, std::abs(error));
        };
        const auto vectorError = [&check](const mfem::Vector &a, const mfem::Vector &b) {
            if (a.Size() != b.Size()) throw std::runtime_error("Reduced surface self-check size mismatch.");
            double error = 0.0;
            for (int i = 0; i < a.Size(); ++i) error = std::max(error, std::abs(a[i] - b[i]));
            check(error);
        };
        mfem::Vector stateFactors(7), residualFactors(7), full(7), reduced(6);
        for (int i = 0; i < 7; ++i) {
            stateFactors[i] = 0.7 + 0.3 * i;
            residualFactors[i] = 1.9 - 0.13 * i;
            full[i] = 0.2 - 0.17 * i;
        }
        for (int i = 0; i < 6; ++i) reduced[i] = 0.11 * (i + 1);
        mfem::DenseMatrix physical(3, 2);
        physical(0, 0) = 1.0; physical(1, 0) = 0.5; physical(2, 0) = 0.0;
        physical(0, 1) = 0.2; physical(1, 1) = -0.3; physical(2, 1) = 1.0;
        ReducedSurfaceCoordinates coordinates(physical, stateFactors, residualFactors, 1, 4, comm);
        check(coordinates.Report().stateOrthogonalityError);
        check(coordinates.Report().residualOrthogonalityError);
        mfem::Vector injected, restricted, roundtrip, other;
        coordinates.InjectState(reduced, injected);
        coordinates.RestrictState(injected, roundtrip);
        vectorError(reduced, roundtrip);
        coordinates.RestrictState(full, restricted);
        check(std::abs(injected * full - reduced * restricted));
        coordinates.InjectResidual(reduced, injected);
        coordinates.RestrictResidual(injected, roundtrip);
        vectorError(reduced, roundtrip);
        coordinates.RestrictResidual(full, restricted);
        check(std::abs(injected * full - reduced * restricted));

        // Verify the physical dual identity independently of InjectResidual:
        // T_r^T rhat_shape = B^T (D_r^{-1} Fhat)_shape.
        for (int j = 0; j < 2; ++j) {
            double lhs = 0.0, rhs = 0.0;
            for (int l = 0; l < 2; ++l) lhs += coordinates.ResidualTransform()(l, j) * restricted[4 + l];
            for (int i = 0; i < 3; ++i) rhs += physical(i, j) * full[4 + i] / residualFactors[4 + i];
            check(std::abs(lhs - rhs));
        }
        // The two QR bases really differ in this fixture; transpose reuse
        // would therefore fail the physical dual identity above.
        double mapDifference = 0.0;
        for (int j = 0; j < 2; ++j)
            for (int i = 0; i < 3; ++i)
                mapDifference += std::abs(coordinates.StateBasis()(i, j) - coordinates.ResidualBasis()(i, j));
        if (mapDifference < 1.e-3) throw std::runtime_error("Reduced surface dual-map fixture is ineffective.");
        ++report.checks;

        reduced_surface_detail::SelfCheckOperator fullOperator(7);
        ReducedSurfaceOperator operation(fullOperator, coordinates);
        operation.Mult(reduced, roundtrip);
        coordinates.InjectState(reduced, injected);
        fullOperator.Mult(injected, other);
        coordinates.RestrictResidual(other, restricted);
        vectorError(roundtrip, restricted);
        reduced_surface_detail::SelfCheckInverse fullInverse(7);
        ReducedSurfacePreconditioner inverse(fullInverse, coordinates);
        inverse.SetOperator(operation);
        const int oldApplications = fullOperator.applications;
        inverse.Mult(reduced, roundtrip);
        coordinates.InjectResidual(reduced, injected);
        fullInverse.Mult(injected, other);
        coordinates.RestrictState(other, restricted);
        vectorError(roundtrip, restricted);
        check(fullOperator.applications - oldApplications);

        ++fullOperator.generation;
        try { operation.Mult(reduced, roundtrip); }
        catch (const std::invalid_argument &) { report.rejectedStaleOperator = true; }
        if (!report.rejectedStaleOperator || operation.IsCurrent())
            throw std::runtime_error("Reduced surface stale-generation guard failed.");
        ++report.checks;
        operation.Refresh();
        operation.Mult(reduced, roundtrip);

        mfem::DenseMatrix identity(3);
        identity = 0.0;
        for (int i = 0; i < 3; ++i) identity(i, i) = 1.0;
        ReducedSurfaceCoordinates fullCoordinates(identity, stateFactors, residualFactors, 1, 4, comm);
        fullCoordinates.InjectState(full, injected);
        vectorError(full, injected);
        fullCoordinates.RestrictResidual(full, restricted);
        vectorError(full, restricted);
        ReducedSurfaceOperator fullRecovery(fullOperator, fullCoordinates);
        fullRecovery.Mult(full, roundtrip);
        fullOperator.Mult(full, other);
        vectorError(roundtrip, other);
        ReducedSurfacePreconditioner inverseRecovery(fullInverse, fullCoordinates);
        inverseRecovery.Mult(full, roundtrip);
        fullInverse.Mult(full, other);
        vectorError(roundtrip, other);

        mfem::DenseMatrix invalid(physical);
        for (int i = 0; i < 3; ++i) invalid(i, 1) = invalid(i, 0);
        try { ReducedSurfaceCoordinates rejected(invalid, stateFactors, residualFactors, 1, 4, comm); }
        catch (const std::invalid_argument &) { report.rejectedDependentBasis = true; }
        invalid = physical;
        invalid(0, 0) = std::numeric_limits<double>::quiet_NaN();
        try { ReducedSurfaceCoordinates rejected(invalid, stateFactors, residualFactors, 1, 4, comm); }
        catch (const std::invalid_argument &) { report.rejectedNonfiniteBasis = true; }
        if (!report.rejectedDependentBasis || !report.rejectedNonfiniteBasis)
            throw std::runtime_error("Reduced surface basis validity guards failed.");
        report.checks += 2;
        return report;
    }
} // namespace experiment::plateau_diagnostics