// SPDX-License-Identifier: MPL-2.0
#pragma once

#include "vibeqc/tensor/structure.hpp"
#include <chrono>
#include <cmath>
#include <cstring>
#if defined(_OPENMP)
#include <omp.h>
#endif
#include <Eigen/Core>
#include <array>

namespace vibeqc::tensor {

enum class TileState { evaluated, exact_zero, unavailable };
struct FactorTile {
    const double* data = nullptr;
    std::size_t bytes = 0;
    TileState state = TileState::unavailable;
};

// Compare layout-aware/explicit packing as in Matthews,
// doi:10.1137/16M108968X, sections 4-7/9.1. This path deliberately retains
// the existing small Eigen kernels; it does not implement GETT or TBLIS.
// First provider capability: immutable, replayable, resident FP64 column-major
// panels. A snapshot's pointers remain valid until the synchronous scope ends.
// Transient/single-pass providers need a different admitted schedule (#797).
class ResidentFactors {
    Budget::Reservation inventory_;
    std::weak_ptr<BorrowScope::State> lifetime_;
public:
    const Budget budget;
    const Domain auxiliary, left, right;
    const Space source_revision;
    ResidentFactors(const Budget& b, BorrowScope& scope, Domain p, Domain m, Domain n,
                    std::size_t live_bytes, bool replayable = true, std::size_t workspace_bytes = 0, Space source = {})
        : inventory_(b.reserve(add(sizeof(ResidentFactors), live_bytes))), lifetime_(scope.token()),
          budget(b), auxiliary(std::move(p)), left(std::move(m)), right(std::move(n)),
          source_revision(source.valid() ? source : Space::create()) {
        if (!scope.same_budget(b)) throw std::invalid_argument("tensor provider ledger mismatch");
        auxiliary.validate(b); left.validate(b); right.validate(b);
        if (!replayable || workspace_bytes != 0)
            throw std::invalid_argument("RI-J requires resident replayable tiles with zero provider workspace");
        signed_size(mul(left.extent, right.extent)); signed_size(auxiliary.extent);
    }
    virtual ~ResidentFactors() = default;
    void check_lifetime() const {
        auto token = lifetime_.lock();
        if (!token || !token->active) throw std::runtime_error("expired tensor provider scope");
    }
    virtual FactorTile tile(std::size_t p) const = 0;
    virtual std::uint64_t value_revision() const { return 0; }
    virtual Space pattern_revision() const { return source_revision; }
};

class PerPFactors final : public ResidentFactors {
    const std::vector<Eigen::MatrixXd>& factors_;
    static std::size_t inventory(const std::vector<Eigen::MatrixXd>& factors, std::size_t external) {
        auto bytes = add(external, add(sizeof(PerPFactors) - sizeof(ResidentFactors),
                                      mul(factors.capacity(), sizeof(Eigen::MatrixXd))));
        for (const auto& matrix : factors) bytes = add(bytes, mul(std::size_t(matrix.size()), sizeof(double)));
        return bytes;
    }
public:
    PerPFactors(const Budget& b, BorrowScope& scope, const std::vector<Eigen::MatrixXd>& factors,
                Domain p, Domain m, Domain n, std::size_t external_live_bytes = 0, Space source = {})
        : ResidentFactors(b, scope, p, m, n, inventory(factors, external_live_bytes),true,0,source), factors_(factors) {
        if (factors.size() != p.extent) throw std::invalid_argument("RI-J auxiliary extent mismatch");
        for (const auto& matrix : factors)
            if (std::size_t(matrix.rows()) != m.extent || std::size_t(matrix.cols()) != n.extent)
                throw std::invalid_argument("RI-J per-P shape mismatch");
    }
    FactorTile tile(std::size_t p) const override {
        const auto& matrix = factors_.at(p);
        return {matrix.data(), mul(std::size_t(matrix.size()), sizeof(double)), TileState::evaluated};
    }
};

// Fixed-rank schemas at this registered kernel boundary still enter the same
// runtime-rank index checker as general affine contractions.

enum class RiJSelection { existing, measured, panel };
struct RiJOptions {
    NumericalMode numerical_mode=NumericalMode::preserve_order;
    RiJSelection selection=RiJSelection::existing;
    // Resource/search ceilings, not machine-specific performance thresholds.
    std::size_t maximum_packing_bytes=8*1024*1024;
    std::size_t maximum_tuning_terms=1000000000;
    int tuning_passes=3; // plus one warmup for each of two candidates
    Space hardware=local_hardware(),placement=local_placement();
};

namespace ri_j_detail {
// Audited allocation-free Eigen boundary for a single RHS. Explicit BuiltIn
// bypasses GeneralMatrixVector_BLAS even when EIGEN_USE_BLAS is enabled.
// Caller validates extents and owns/admits pack, gamma, D and private output.
inline void panel_products(std::size_t elements,std::size_t count,const double* pack,
                           const double* density,double* gamma,double* output) {
#if EIGEN_VERSION_AT_LEAST(3,4,0)
    using I=Eigen::Index;
    using RM=Eigen::internal::const_blas_data_mapper<double,I,Eigen::RowMajor>;
    using CM=Eigen::internal::const_blas_data_mapper<double,I,Eigen::ColMajor>;
    std::fill_n(gamma,count,0.0);
    Eigen::internal::general_matrix_vector_product<I,double,RM,Eigen::RowMajor,false,double,CM,false,Eigen::BuiltIn>::run(
        signed_size(count),signed_size(elements),RM(pack,signed_size(elements)),CM(density,1),gamma,1,1.0);
    Eigen::internal::general_matrix_vector_product<I,double,CM,Eigen::ColMajor,false,double,RM,false,Eigen::BuiltIn>::run(
        signed_size(elements),signed_size(count),CM(pack,signed_size(elements)),RM(gamma,1),output,1,1.0);
#else
    (void)elements;(void)count;(void)pack;(void)density;(void)gamma;(void)output;
    throw std::logic_error("Eigen built-in panel backend unavailable");
#endif
}
// The Eigen 5 (internal version 3.5) packet-mapper ABI was inspected and
// allocation-tested. Older versions keep the canonical packed-data mapper.
#if EIGEN_WORLD_VERSION == 3 && EIGEN_MAJOR_VERSION == 5 && EIGEN_MINOR_VERSION == 0
template<int Order> struct FactorMapper {
    const FactorTile* tiles;
    Eigen::Index leading;
    Eigen::Index stride() const { return leading; }
    double operator()(Eigen::Index i,Eigen::Index j) const {
        const auto& tile=tiles[Order==Eigen::ColMajor ? j : i];
        return tile.state==TileState::exact_zero ? 0.0 : tile.data[Order==Eigen::ColMajor ? i : j];
    }
    template<class Packet,int Alignment> Packet load(Eigen::Index i,Eigen::Index j) const {
        const auto& tile=tiles[Order==Eigen::ColMajor ? j : i];
        return tile.state==TileState::exact_zero ? Eigen::internal::pset1<Packet>(0.0) :
            Eigen::internal::ploadu<Packet>(tile.data+(Order==Eigen::ColMajor ? i : j));
    }
};
inline void mapped_products(std::size_t elements,std::size_t count,const FactorTile* tiles,
                            const double* density,double* gamma,double* output) {
    using I=Eigen::Index;
    using RM=Eigen::internal::const_blas_data_mapper<double,I,Eigen::RowMajor>;
    using CM=Eigen::internal::const_blas_data_mapper<double,I,Eigen::ColMajor>;
    std::fill_n(gamma,count,0.0);
    Eigen::internal::general_matrix_vector_product<I,double,FactorMapper<Eigen::RowMajor>,Eigen::RowMajor,false,double,CM,false,Eigen::BuiltIn>::run(
        signed_size(count),signed_size(elements),FactorMapper<Eigen::RowMajor>{tiles,signed_size(elements)},CM(density,1),gamma,1,1.0);
    Eigen::internal::general_matrix_vector_product<I,double,FactorMapper<Eigen::ColMajor>,Eigen::ColMajor,false,double,RM,false,Eigen::BuiltIn>::run(
        signed_size(elements),signed_size(count),FactorMapper<Eigen::ColMajor>{tiles,signed_size(elements)},RM(gamma,1),output,1,1.0);
}
#endif
} // namespace ri_j_detail

class RiJPlan {
    Budget::Reservation control_;
    std::shared_ptr<const ResidentFactors> factors_;
    View density_;
    IndexProgram gamma_ir_, j_ir_;
    std::uint64_t density_version_;
    int workers_;
    RiJOptions options_;
    StructuralCache* cache_=nullptr; // synchronous caller-owned schedule cache
    StructuralKey key_;
    Schedule schedule_;
    std::uint64_t factor_version_;
    std::size_t panel_columns_=0;
    bool tune_=false,cache_hit_=false,mapped_=false;
    static constexpr std::size_t backend_id() {
        return 1000000 + EIGEN_WORLD_VERSION*10000 + EIGEN_MAJOR_VERSION*100 + EIGEN_MINOR_VERSION;
    }
    static Schema<3> factor_schema(const ResidentFactors& f) {
        return {{Axis{f.auxiliary, 0}, Axis{f.left, 0}, Axis{f.right, 0}}};
    }
    static Schema<1> gamma_schema(const ResidentFactors& f) { return {{Axis{f.auxiliary, 1}}}; }
    static Schema<2> output_schema(const ResidentFactors& f) {
        return {{Axis{f.left, 1}, Axis{f.right, signed_size(f.left.extent)}}};
    }
public:
    struct Result {
        Budget::Reservation reservation;
        Eigen::MatrixXd value;
        int workers_used = 1;
        ScheduleKernel kernel_used=ScheduleKernel::per_p;
        std::size_t panel_columns=0,tile_visits=0;
        bool structural_cache_hit=false;
        double packing_seconds=0,per_p_seconds=0,panel_seconds=0;
        const char* selection_reason="existing numerical/thread policy";
        Result(const Budget& b, std::size_t m, std::size_t n)
            : reservation(b.reserve(add(sizeof(Result), mul(mul(m, n), sizeof(double))))),
              value(signed_size(m), signed_size(n)) { value.setZero(); }
    };
    RiJPlan(std::shared_ptr<const ResidentFactors> factors, View density, int max_workers = 1,
            RiJOptions options = {}, StructuralCache* cache = nullptr)
        : control_(checked_handle(factors)->budget.reserve(sizeof(RiJPlan))), factors_(std::move(factors)), density_(std::move(density)),
          gamma_ir_(factors_->budget, factor_schema(*factors_), {0, 1, 2}, density_, {1, 2}, gamma_schema(*factors_), {0}),
          j_ir_(factors_->budget, factor_schema(*factors_), {0, 1, 2}, gamma_schema(*factors_), {0}, output_schema(*factors_), {1, 2}),
          density_version_(density_.allocation()->version), workers_(1),options_(options),cache_(cache),
          factor_version_(factors_->value_revision()) {
        factors_->check_lifetime();
        if (!factors_->budget.same(density_.allocation()->budget)) throw std::invalid_argument("RI-J ledger mismatch");
        if (density_.allocation()->type != DType::fp64 || density_.access() == Access::write || density_.conjugate())
            throw std::invalid_argument("RI-J requires readable FP64 density");
        if (max_workers < 1) throw std::invalid_argument("RI-J worker count must be positive");
        auto terms = mul(factors_->auxiliary.extent, mul(factors_->left.extent, factors_->right.extent));
        // Preserve the existing DF threshold/team ownership; never nest teams.
#if defined(_OPENMP)
        if (options_.selection==RiJSelection::existing && terms >= 10000000 && !omp_in_parallel())
            workers_ = std::max(1, static_cast<int>(std::min(
                std::min(factors_->auxiliary.extent, factors_->right.extent),
                static_cast<std::size_t>(std::min(max_workers, omp_get_max_threads())))));
#endif
        if(options_.tuning_passes<1 || options_.tuning_passes>3)
            throw std::invalid_argument("RI-J tuning passes must be between one and three");
        if(options_.selection==RiJSelection::panel && options_.numerical_mode!=NumericalMode::reassociated_fp64)
            throw std::invalid_argument("panel RI-J requires explicit reassociated FP64 mode");
        schedule_.workers=workers_;
        key_.operation=1;key_.backend=backend_id();key_.source=factors_->source_revision;key_.pattern=factors_->pattern_revision();
        key_.hardware=options_.hardware;key_.placement=options_.placement;
        key_.mode=options_.numerical_mode;key_.requested_workers=max_workers;key_.actual_workers=workers_;
        key_.thread_policy=options_.selection==RiJSelection::existing ? 0 : 1;
        key_.planner_policy=static_cast<std::size_t>(options_.selection);
        key_.pack_limit=options_.maximum_packing_bytes;key_.tuning_work_limit=options_.maximum_tuning_terms;
        key_.tuning_passes=options_.tuning_passes;
        key_.budget=factors_->budget.limit();key_.execution_allowance=factors_->budget.available();
#if defined(_OPENMP)
        key_.nested=omp_in_parallel();
#endif
        key_.append_schema(factor_schema(*factors_),DType::fp64,Access::read);
        key_.append(density_);
        // The audited Eigen BuiltIn GEMV boundary bypasses opaque BLAS and
        // dispatches one RHS without any dynamic backend packing workspace.
#if EIGEN_VERSION_AT_LEAST(3,4,0)
        if(options_.selection!=RiJSelection::existing && options_.numerical_mode==NumericalMode::reassociated_fp64 &&
            factors_->auxiliary.extent && density_.size()) {
            const auto extra_output=options_.selection==RiJSelection::measured ?
                add(sizeof(Result),mul(density_.size(),sizeof(double))) : 0;
            const auto base=add(base_execution_bytes(),extra_output);
            const auto column=mul(add(density_.size(),1),sizeof(double));
            if(key_.execution_allowance>=base)
                panel_columns_=std::min(factors_->auxiliary.extent,
                    std::min(options_.maximum_packing_bytes,key_.execution_allowance-base)/column);
        }
#endif
#if EIGEN_WORLD_VERSION == 3 && EIGEN_MAJOR_VERSION == 5 && EIGEN_MINOR_VERSION == 0
        // Stable separate per-P buffers have unit-stride packet runs. The
        // measured zero-copy mapper avoids the packing cost entirely; no
        // representation/ownership fiction makes these one affine allocation.
        if(options_.selection==RiJSelection::measured && options_.numerical_mode==NumericalMode::reassociated_fp64 &&
            factors_->auxiliary.extent && density_.size() && key_.execution_allowance>=
                add(base_execution_bytes(),add(sizeof(Result),mul(density_.size(),sizeof(double))))) {
            mapped_=true;panel_columns_=factors_->auxiliary.extent;
        }
#endif
        if(options_.selection==RiJSelection::panel) {
            if(!panel_columns_) throw std::runtime_error("RI-J panel workspace unavailable");
            schedule_.kernel=ScheduleKernel::panel_product;schedule_.panel_columns=panel_columns_;
        } else if(options_.selection==RiJSelection::measured && panel_columns_ &&
                  terms<=options_.maximum_tuning_terms/std::size_t(4*(options_.tuning_passes+1))) {
            tune_=true;
        }
        if(cache_ && cache_->lookup(key_,schedule_)) {
            cache_hit_=true;tune_=false;
            // A schedule never waives current shape/capability/admission checks.
            if(schedule_.workers!=workers_ ||
                (schedule_.kernel!=ScheduleKernel::per_p && schedule_.kernel!=ScheduleKernel::panel_product && schedule_.kernel!=ScheduleKernel::mapped_panel) ||
                (schedule_.kernel==ScheduleKernel::mapped_panel && !mapped_))
                throw std::logic_error("cached RI-J schedule capability mismatch");
            if((schedule_.kernel==ScheduleKernel::panel_product || schedule_.kernel==ScheduleKernel::mapped_panel) &&
                (!panel_columns_ || schedule_.panel_columns>panel_columns_))
                throw std::logic_error("cached RI-J panel exceeds current capability");
        }
        if(!tune_ && schedule_.kernel!=ScheduleKernel::panel_product && schedule_.kernel!=ScheduleKernel::mapped_panel)panel_columns_=0;
        if(schedule_.kernel==ScheduleKernel::panel_product || schedule_.kernel==ScheduleKernel::mapped_panel)panel_columns_=schedule_.panel_columns;
    }
    const StructuralKey& structural_key() const { return key_; }
    bool structural_cache_hit() const { return cache_hit_; }
    std::size_t panel_columns() const { return panel_columns_; }
    const char* kernel() const { return panel_columns_ ? "cpu.eigen.ri_j.bounded_selection.v1" : "cpu.eigen.ri_j.per_p.v1"; }
    const IndexProgram& gamma_program() const { return gamma_ir_; }
    const IndexProgram& j_program() const { return j_ir_; }
    const char* equations() const { return "gamma[p]=B[p,m,n]*D[m,n]; J[m,n]=B[p,m,n]*gamma[p]"; }
    const char* reduction_policy() const { return panel_columns_ ? "explicit reassociated FP64 panel products; per-P fallback" : "Eigen Frobenius; ascending auxiliary output accumulation"; }
    std::size_t packing_bytes() const {
        return dense_density() ? 0 : mul(density_.size(), sizeof(double));
    }
    std::size_t base_execution_bytes() const {
        auto p = factors_->auxiliary.extent;
        return add(executor_control_bytes(), add(add(mul(p, sizeof(FactorTile)), mul(p, sizeof(double))),
                   add(packing_bytes(), add(sizeof(Result), mul(density_.size(), sizeof(double))))));
    }
    std::size_t factor_packing_bytes() const { return mapped_ ? 0 : mul(panel_columns_,mul(density_.size(),sizeof(double))); }
    std::size_t execution_bytes() const {
        return add(base_execution_bytes(),add(factor_packing_bytes(),add(mul(mapped_ ? 0 : panel_columns_,sizeof(double)),
            tune_ ? add(sizeof(Result),mul(density_.size(),sizeof(double))) : 0)));
    }
    static constexpr std::size_t executor_control_bytes() {
        return sizeof(Array<FactorTile>) + 4*sizeof(Array<double>) + 1024;
    }
    bool dense_density() const {
        return density_.axis(0).stride == 1 && density_.axis(1).stride == signed_size(factors_->left.extent);
    }
    Result execute() const {
        factors_->check_lifetime();
        Lease density_lease(density_.allocation(), false);
        if (density_.allocation()->version != density_version_ || factors_->value_revision()!=factor_version_ || factors_->pattern_revision()!=key_.pattern)
            throw std::runtime_error("stale RI-J density/factor binding");
        // Preflight all execution storage before provider visits. Individual
        // tickets below remain live; concurrent admission can still refuse safely
        // before numerical execution. Inputs and plan stay reserved throughout.
        { auto preflight = factors_->budget.reserve(execution_bytes()); }
        const auto p_count = factors_->auxiliary.extent, rows = factors_->left.extent, cols = factors_->right.extent;
        const auto elements = mul(rows, cols), bytes = mul(elements, sizeof(double));
        auto executor_control = factors_->budget.reserve(executor_control_bytes());
        Array<FactorTile> tiles(factors_->budget, p_count);
        for (std::size_t p = 0; p < p_count; ++p) {
            tiles[p] = factors_->tile(p);
            const auto& tile = tiles[p];
            add(reinterpret_cast<std::uintptr_t>(tile.data), tile.bytes);
            if (tile.state == TileState::exact_zero) {
                if (tile.data || tile.bytes) throw std::invalid_argument("exact-zero tile must have no payload");
            } else if (tile.state != TileState::evaluated || tile.bytes != bytes ||
                       (bytes && (!tile.data || reinterpret_cast<std::uintptr_t>(tile.data) % alignof(double)))) {
                throw std::runtime_error("unavailable or invalid RI-J provider tile");
            }
        }
        Array<double> gamma(factors_->budget, p_count);
        Array<double> packed(factors_->budget, packing_bytes() / sizeof(double));
        Array<double> factor_pack(factors_->budget,factor_packing_bytes()/sizeof(double));
        Array<double> panel_gamma(factors_->budget,mapped_ ? 0 : panel_columns_);
        Result result(factors_->budget, rows, cols); // private until success
        result.structural_cache_hit=cache_hit_;result.tile_visits=p_count;
        if (!elements) return result;
        const auto* base = static_cast<const double*>(density_.allocation()->base);
        const double* d = base + density_.offset();
        if (!dense_density()) {
            for (std::size_t n = 0; n < cols; ++n)
                for (std::size_t m = 0; m < rows; ++m)
                    packed[n * rows + m] = base[signed_add(density_.offset(), signed_add(
                        stride_span(m, density_.axis(0).stride), stride_span(n, density_.axis(1).stride)))];
            d = packed.data();
        }
        using Matrix = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>;
        using Map = Eigen::Map<const Matrix>;
        Map D(d, signed_size(rows), signed_size(cols));
        auto factor = [&](std::size_t p) { return Map(tiles[p].data, signed_size(rows), signed_size(cols)); };
        // These coefficient-wise Eigen kernels allocate no backend packing or
        // GEMM workspace. No global BLAS thread setting is modified.
        int effective_workers = workers_;
#if defined(_OPENMP)
        if (omp_in_parallel()) effective_workers = 1;
#endif
        auto per_p=[&](Result& result) {
            result.value.setZero();
            if (effective_workers == 1) {
                for (std::size_t p = 0; p < p_count; ++p)
                    gamma[p] = tiles[p].state == TileState::exact_zero ? 0.0 : (factor(p).array() * D.array()).sum();
                for (std::size_t p = 0; p < p_count; ++p)
                    if (tiles[p].state != TileState::exact_zero) result.value.noalias() += gamma[p] * factor(p);
            } else {
                int actual_workers = 1;
            #pragma omp parallel num_threads(effective_workers)
                {
                #pragma omp single
                    {
#if defined(_OPENMP)
                        actual_workers = omp_get_num_threads();
#endif
                    }
                #pragma omp for schedule(static)
                    for (Eigen::Index p = 0; p < signed_size(p_count); ++p)
                        gamma[std::size_t(p)] = tiles[std::size_t(p)].state == TileState::exact_zero ? 0.0 :
                            (factor(std::size_t(p)).array() * D.array()).sum();
                #pragma omp for schedule(static)
                    for (Eigen::Index n = 0; n < signed_size(cols); ++n)
                        for (std::size_t p = 0; p < p_count; ++p)
                            if (tiles[p].state != TileState::exact_zero)
                                result.value.col(n).noalias() += gamma[p] * factor(p).col(n);
                }
                result.workers_used = actual_workers;
            }
        };
        auto panel=[&](Result& target) {
            target.value.setZero();target.packing_seconds=0;
#if EIGEN_WORLD_VERSION == 3 && EIGEN_MAJOR_VERSION == 5 && EIGEN_MINOR_VERSION == 0
            if(mapped_) {
                ri_j_detail::mapped_products(elements,p_count,tiles.data(),d,gamma.data(),target.value.data());
                target.kernel_used=ScheduleKernel::mapped_panel;target.panel_columns=p_count;return;
            }
#endif
#if EIGEN_VERSION_AT_LEAST(3,4,0)
            for(std::size_t first=0;first<p_count;first+=std::min(panel_columns_,p_count-first)) {
                const auto count=std::min(panel_columns_,p_count-first);
                const auto start=std::chrono::steady_clock::now();
                for(std::size_t p=0;p<count;++p) {
                    auto* dest=factor_pack.data()+p*elements;
                    if(tiles[first+p].state==TileState::exact_zero)std::fill_n(dest,elements,0.0);
                    else std::copy_n(tiles[first+p].data,elements,dest);
                }
                target.packing_seconds+=std::chrono::duration<double>(std::chrono::steady_clock::now()-start).count();
                ri_j_detail::panel_products(elements,count,factor_pack.data(),d,panel_gamma.data(),target.value.data());
            }
#endif
            target.kernel_used=ScheduleKernel::panel_product;target.panel_columns=panel_columns_;
        };
        auto measured=schedule_;
        bool do_tune=tune_;
        if(do_tune && cache_ && cache_->lookup(key_,measured)) {
            do_tune=false;result.structural_cache_hit=true;
        }
        if(do_tune) {
            Result candidate(factors_->budget,rows,cols);
            std::array<double,3> old_times{},new_times{};
            auto timed=[](auto operation) {
                const auto start=std::chrono::steady_clock::now();operation();
                return std::chrono::duration<double>(std::chrono::steady_clock::now()-start).count();
            };
            per_p(result);panel(candidate); // bounded warmup; never stored numeric data
            for(int i=0;i<options_.tuning_passes;++i) {
                if(i%2) { new_times[i]=timed([&]{panel(candidate);});old_times[i]=timed([&]{per_p(result);}); }
                else { old_times[i]=timed([&]{per_p(result);});new_times[i]=timed([&]{panel(candidate);}); }
            }
            std::sort(old_times.begin(),old_times.begin()+options_.tuning_passes);
            std::sort(new_times.begin(),new_times.begin()+options_.tuning_passes);
            measured.per_p_seconds=old_times[options_.tuning_passes/2];measured.panel_seconds=new_times[options_.tuning_passes/2];
            double maximum=0,error=0;bool finite=true;
            for(std::size_t i=0;i<elements;++i) {
                finite=finite && std::isfinite(result.value.data()[i]) && std::isfinite(candidate.value.data()[i]);
                maximum=std::max(maximum,std::abs(result.value.data()[i]));
                error=std::max(error,std::abs(result.value.data()[i]-candidate.value.data()[i]));
            }
            const double bound=128*std::numeric_limits<double>::epsilon()*double(add(elements,p_count))*std::max(1.0,maximum);
            // Five percent is a portable measurement-noise guard, not a shape
            // crossover. No hardware threshold or numeric pack is cached.
            if(finite && error<=bound && measured.panel_seconds<0.95*measured.per_p_seconds) {
                measured.kernel=mapped_ ? ScheduleKernel::mapped_panel : ScheduleKernel::panel_product;measured.panel_columns=panel_columns_;
                result.value.swap(candidate.value);result.kernel_used=measured.kernel;result.panel_columns=panel_columns_;
                result.packing_seconds=candidate.packing_seconds;
                result.selection_reason="bounded measurement selected panel product";
            } else result.selection_reason="per-P wins measurement or numerical comparison";
            if(cache_)cache_->store(key_,measured);
        } else if(measured.kernel==ScheduleKernel::panel_product || measured.kernel==ScheduleKernel::mapped_panel) {
            panel(result);result.selection_reason=result.structural_cache_hit ? "cached bounded panel schedule" : "explicit bounded panel schedule";
        } else {
            per_p(result);result.selection_reason=cache_hit_ ? "cached per-P schedule" : "ordered mode, capability or tuning-work fallback";
            if(cache_ && !cache_hit_)cache_->store(key_,measured);
        }
        if(result.kernel_used!=ScheduleKernel::per_p && !result.value.allFinite()) {
            per_p(result);result.kernel_used=ScheduleKernel::per_p;result.panel_columns=0;
            result.selection_reason="non-finite panel result uses existing IEEE path";
        }
        result.per_p_seconds=measured.per_p_seconds;result.panel_seconds=measured.panel_seconds;
        return result;
    }
};

} // namespace vibeqc::tensor