// Object recognition entirely in KNP with AltAI neurons: retina → learned feature detectors → readout.
//
// altai_objects features <patches.bin> <count> <K> <out_weights.bin> <backend.so>
//     Learns K feature detectors on 12×12 RGB crops (uint8, 3×12×12 each). Network: pixels → retina (fixed
//     ON/OFF centre–surround neurons, 6 maps 6×6) → K competing AltAI neurons with trainable synapses. Each
//     detector inhibits the others and its own spike gives it dopamine, so only winners learn; with the delta
//     backend (KNP_DELTA_RULE, KNP_DELTA_DECAY, KNP_DELTA_LAG) the winner moves its weights towards the input.
//     Writes float32 weights [K][216] (retina channel, y, x order).
//
// altai_objects run <images.bin> <labels.txt> <train> <eval> <features.bin> <K> <classes> <backend.so> <prefix>
//     Full network: pixels (3×32×32) → retina (6×32×32) → detectors (K per position, stride 2, weights from
//     features.bin minus their mean, integers) → readout (one neuron per class, delta-rule training with the
//     label dopamine +a every active step and own-spike dopamine −b, KNP_DELTA=a,b). After <train> images the
//     readout weights are read back from the backend and fixed; then <eval> images are run with a plain AltAI
//     readout. Writes <prefix>_readout.bin (float32 [classes][neurons]) and <prefix>_eval.csv
//     (index,label,spike counts per class).
//
// Environment: KNP_RETINA_THR (default 250), KNP_DET_THR (detector threshold, default 2000),
// KNP_OUT_THR (readout threshold, default 4000), KNP_INHIB (detector lateral inhibition, default 5000),
// KNP_SUMTHR (features: threshold increase per unit of positive weight sum, default 0).

#include <knp/framework/model_executor.h>
#include <knp/framework/network.h>
#include <knp/neuron-traits/altai_lif.h>
#include <knp/synapse-traits/delta.h>
#include <knp/synapse-traits/stdp_synaptic_resource_rule.h>

#include <chrono>
#include <cmath>
#include <fstream>
#include <iostream>
#include <map>
#include <numeric>
#include <random>
#include <string>
#include <vector>


using DeltaSynapse = knp::synapse_traits::DeltaSynapse;
using DeltaProjection = knp::core::Projection<DeltaSynapse>;
using ResourceSynapse = knp::synapse_traits::SynapticResourceSTDPDeltaSynapse;
using ResourceProjection = knp::core::Projection<ResourceSynapse>;
using AltAINeuron = knp::neuron_traits::AltAILIF;
using StdpNeuron = knp::neuron_traits::SynapticResourceSTDPAltAILIFNeuron;
using AltAIPopulation = knp::core::Population<AltAINeuron>;
using StdpPopulation = knp::core::Population<StdpNeuron>;

constexpr int steps_per_image = 15, active_steps = 10, patch = 6, dog_r = 3, maps = 6;
constexpr float w_min = -253.122F, w_max = 92.3957F;


float env_f(const char *name, float def)
{
    const char *v = std::getenv(name);
    return v ? std::stof(v) : def;
}


std::vector<float> env_list(const char *name, std::vector<float> def)
{
    const char *v = std::getenv(name);
    if (!v) return def;
    std::vector<float> out;
    std::string s(v);
    size_t pos = 0;
    while (pos <= s.size())
    {
        size_t next = s.find(',', pos);
        if (next == std::string::npos) next = s.size();
        out.push_back(std::stof(s.substr(pos, next - pos)));
        pos = next + 1;
    }
    return out;
}


// Difference of Gaussians, sigma 0.5 and 1.5, 7×7, scaled by 1000 (as in altai-objects/retina_proto.py).
std::vector<float> dog_kernel()
{
    auto gauss = [](float sigma)
    {
        std::vector<float> g((2 * dog_r + 1) * (2 * dog_r + 1));
        float sum = 0;
        for (int y = -dog_r; y <= dog_r; ++y)
            for (int x = -dog_r; x <= dog_r; ++x)
                sum += g[(y + dog_r) * (2 * dog_r + 1) + x + dog_r] = std::exp(-(x * x + y * y) / (2 * sigma * sigma));
        for (auto &v : g) v /= sum;
        return g;
    };
    auto a = gauss(0.5F), b = gauss(1.5F);
    std::vector<float> d(a.size());
    for (size_t i = 0; i < a.size(); ++i) d[i] = std::round((a[i] - b[i]) * 1000);
    return d;
}


float resource_from_weight(float w)
{
    const double diff = w_max - w_min, over = w - w_min;
    return static_cast<float>(diff * over / (diff - over));
}


knp::neuron_traits::neuron_parameters<AltAINeuron> plain_neuron(float threshold)
{
    knp::neuron_traits::neuron_parameters<AltAINeuron> p;
    p.activation_threshold_ = static_cast<uint16_t>(std::min(threshold, 65535.F));
    p.negative_activation_threshold_ = 0;
    p.potential_reset_value_ = 0;
    p.potential_leak_ = 0;
    return p;
}


knp::neuron_traits::neuron_parameters<StdpNeuron> learning_neuron(float threshold)
{
    knp::neuron_traits::neuron_parameters<StdpNeuron> p;
    p.activation_threshold_ = static_cast<uint16_t>(std::min(threshold, 65535.F));
    p.negative_activation_threshold_ = 0;
    p.potential_reset_value_ = 0;
    p.potential_leak_ = 0;
    p.dopamine_plasticity_time_ = 10;
    p.isi_max_ = 10;
    p.d_h_ = 0;
    p.stability_change_parameter_ = 0;
    p.synapse_sum_threshold_coefficient_ = 0;
    return p;
}


// Retina projection: input channel (pixels c×h×w) → retina population (maps×oh×ow): ON map c uses +DoG,
// OFF map c uses −DoG around pixel (oy + off, ox + off) of channel c; taps outside the image are skipped.
DeltaProjection make_retina(const knp::core::UID &post, int h, int w, int oh, int ow, int off)
{
    const auto dog = dog_kernel();
    const int taps = (2 * dog_r + 1) * (2 * dog_r + 1);
    std::vector<DeltaProjection::Synapse> synapses;
    for (int m = 0; m < maps; ++m)
        for (int oy = 0; oy < oh; ++oy)
            for (int ox = 0; ox < ow; ++ox)
                for (int t = 0; t < taps; ++t)
                {
                    const int y = oy + off + t / (2 * dog_r + 1) - dog_r, x = ox + off + t % (2 * dog_r + 1) - dog_r;
                    if (y < 0 || y >= h || x < 0 || x >= w || dog[t] == 0) continue;
                    const int c = m % 3;
                    const float weight = m < 3 ? dog[t] : -dog[t];
                    synapses.push_back(DeltaProjection::Synapse{
                        {weight, 1, knp::synapse_traits::OutputType::EXCITATORY},
                        static_cast<size_t>(c * h * w + y * w + x), static_cast<size_t>((m * oh + oy) * ow + ox)});
                }
    return DeltaProjection(
        knp::core::UID{false}, post, [&](size_t i) { return std::optional(synapses[i]); }, synapses.size());
}


// Spike raster of an image (uint8 values), incrementing converter from rest, active steps only.
std::vector<knp::core::messaging::SpikeData> raster(const uint8_t *px, size_t n)
{
    std::vector<knp::core::messaging::SpikeData> out(active_steps);
    std::vector<float> st(n, 0.F);
    for (int t = 0; t < active_steps; ++t)
        for (size_t i = 0; i < n; ++i)
        {
            st[i] += px[i] / 255.F;
            if (st[i] >= 1.F)
            {
                st[i] -= 1.F;
                out[t].push_back(static_cast<uint32_t>(i));
            }
        }
    return out;
}


std::vector<uint8_t> read_bin(const std::string &path)
{
    std::ifstream f(path, std::ios::binary);
    return std::vector<uint8_t>((std::istreambuf_iterator<char>(f)), std::istreambuf_iterator<char>());
}


// Read the trainable projection back from the backend: float weights [post][pre-index-in-population].
template <class Backend>
std::vector<float> trained_weights(Backend &backend, const knp::core::UID &proj_uid, size_t pre_n, size_t post_n)
{
    std::vector<float> w(pre_n * post_n, 0.F);
    auto data = backend->get_network_data();
    for (auto &it = *data.projection_range.first; it != *data.projection_range.second; ++it)
    {
        auto proj = *it;
        std::visit(
            [&](auto &p)
            {
                using P = std::decay_t<decltype(p)>;
                if constexpr (std::is_same_v<P, ResourceProjection>)
                    if (p.get_uid() == proj_uid)
                        for (const auto &s : p)
                            w[std::get<2>(s) * pre_n + std::get<1>(s)] = std::get<knp::core::synapse_data>(s).weight_;
            },
            proj);
    }
    return w;
}


int features_mode(int argc, const char *const argv[])
{
    const std::string patches_path = argv[2], out_path = argv[5], backend_path = argv[6];
    const size_t count = std::stoul(argv[3]), k = std::stoul(argv[4]);
    constexpr int crop = 12, in_n = 3 * crop * crop, ret_n = maps * patch * patch;
    const auto px = read_bin(patches_path);
    if (px.size() < count * in_n) throw std::runtime_error("not enough patches");

    AltAIPopulation retina([&](size_t) { return plain_neuron(env_f("KNP_RETINA_THR", 250)); }, ret_n);
    StdpPopulation detectors(
        [&](size_t)
        {
            auto p = learning_neuron(env_f("KNP_DET_THR", 2000));
            // AltAI homeostasis: the threshold grows with the sum of positive weights (KNP_SUMTHR × sum), which
            // works like weight normalisation, so a detector wins by pattern match rather than by weight mass.
            p.synapse_sum_threshold_coefficient_ = env_f("KNP_SUMTHR", 0);
            // Lasting inhibition: the potential may go down to -KNP_DET_NEG, and a leak of KNP_DET_LEAK per step
            // brings it back towards 0 (in the 5 silent steps between crops).
            p.negative_activation_threshold_ = static_cast<uint16_t>(env_f("KNP_DET_NEG", 0));
            p.potential_leak_ = static_cast<int16_t>(env_f("KNP_DET_LEAK", 0));
            return p;
        },
        k);
    auto retina_proj = make_retina(retina.get_uid(), crop, crop, patch, patch, dog_r);

    std::mt19937 rng(static_cast<unsigned>(env_f("KNP_SEED", 1)));
    std::uniform_real_distribution<float> init(0.F, env_f("KNP_INIT_W", 30));
    std::vector<ResourceProjection::Synapse> learn_syn;
    for (size_t j = 0; j < k; ++j)
        for (size_t i = 0; i < ret_n; ++i)
        {
            knp::synapse_traits::synapse_parameters<ResourceSynapse> s;
            s.rule_.w_min_ = w_min;
            s.rule_.w_max_ = w_max;
            s.rule_.dopamine_plasticity_period_ = 10;
            s.weight_ = init(rng);
            s.rule_.synaptic_resource_ = resource_from_weight(s.weight_);
            s.delay_ = 1;
            s.output_type_ = knp::synapse_traits::OutputType::EXCITATORY;
            learn_syn.push_back(ResourceProjection::Synapse{s, i, j});
        }
    ResourceProjection learn_proj(
        retina.get_uid(), detectors.get_uid(), [&](size_t i) { return std::optional(learn_syn[i]); },
        learn_syn.size());
    learn_proj.unlock_weights();

    // Lateral inhibition and own-spike reward.
    const float inhib = env_f("KNP_INHIB", 5000), reward = env_f("KNP_WTA_REWARD", 1);
    DeltaProjection lateral(
        detectors.get_uid(), detectors.get_uid(),
        [&](size_t i) -> std::optional<DeltaProjection::Synapse>
        {
            const size_t pre = i / k, post = i % k;
            if (pre == post)
                return DeltaProjection::Synapse{{reward, 1, knp::synapse_traits::OutputType::DOPAMINE}, pre, post};
            return DeltaProjection::Synapse{{-inhib, 1, knp::synapse_traits::OutputType::EXCITATORY}, pre, post};
        },
        k * k);

    // Self-inhibition after a spike (KNP_SELF_INHIB), so a winner fires about once per crop.
    const float self_inhib = env_f("KNP_SELF_INHIB", 0);
    DeltaProjection self_proj(
        detectors.get_uid(), detectors.get_uid(),
        [&](size_t i) {
            return std::optional(DeltaProjection::Synapse{
                {-self_inhib, 1, knp::synapse_traits::OutputType::EXCITATORY}, i, i});
        },
        self_inhib > 0 ? k : 0);
    // Imprinting (KNP_IMPRINT=1): crop number c teaches detector c mod K through a dopamine channel on every active
    // step; with KNP_DELTA_WDECAY the detector's weights settle at the retina activity of that crop.
    const bool imprint = env_f("KNP_IMPRINT", 0) > 0;
    DeltaProjection teacher(
        knp::core::UID{false}, detectors.get_uid(),
        [&](size_t i) {
            return std::optional(DeltaProjection::Synapse{
                {env_f("KNP_TEACH", 1), 1, knp::synapse_traits::OutputType::DOPAMINE}, i, i});
        },
        k);
    const auto teacher_uid = teacher.get_uid();
    const auto retina_uid = retina_proj.get_uid(), learn_uid = learn_proj.get_uid(), det_uid = detectors.get_uid();
    knp::framework::Network network;
    network.add_population(std::move(retina));
    network.add_population(std::move(detectors));
    network.add_projection(std::move(retina_proj));
    network.add_projection(std::move(learn_proj));
    if (imprint)
        network.add_projection(std::move(teacher));
    else
        network.add_projection(std::move(lateral));
    if (self_inhib > 0 && !imprint) network.add_projection(std::move(self_proj));
    knp::core::UID in_channel, out_channel, teach_channel;
    knp::framework::Model model(std::move(network));
    model.add_input_channel(in_channel, retina_uid);
    model.add_output_channel(out_channel, det_uid);
    if (imprint) model.add_input_channel(teach_channel, teacher_uid);

    std::vector<knp::core::messaging::SpikeData> current;
    auto input_gen = [&](knp::core::Step step) -> knp::core::messaging::SpikeData
    {
        const size_t img = step / steps_per_image, local = step % steps_per_image;
        if (local == 0) current = raster(px.data() + img * in_n, in_n);
        if (local >= active_steps) return {};
        return current[local];
    };
    knp::framework::BackendLoader loader;
    auto backend = loader.load(backend_path);
    auto teach_gen = [&](knp::core::Step step) -> knp::core::messaging::SpikeData
    {
        if (step % steps_per_image >= active_steps) return {};
        return {static_cast<uint32_t>((step / steps_per_image) % k)};
    };
    knp::framework::ModelLoader::InputChannelMap channels{{in_channel, input_gen}};
    if (imprint) channels.insert({teach_channel, teach_gen});
    knp::framework::ModelExecutor executor(model, backend, std::move(channels));
    auto &out = executor.get_loader().get_output_channel(out_channel);
    const auto t0 = std::chrono::steady_clock::now();
    executor.start([&](size_t step) { return step < count * steps_per_image; });
    std::vector<size_t> wins(k, 0);
    for (const auto &m : out.update())
        for (auto i : m.neuron_indexes_) ++wins[i];
    auto w = trained_weights(backend, learn_uid, ret_n, k);
    std::ofstream(out_path, std::ios::binary).write(reinterpret_cast<const char *>(w.data()), w.size() * 4);
    size_t used = 0;
    for (auto v : wins) used += v > 0;
    std::cout << "features: " << count << " patches, " << k << " detectors, "
              << std::chrono::duration<double>(std::chrono::steady_clock::now() - t0).count() << " s; detectors that "
              << "fired: " << used << "; spikes in last batch: " << std::accumulate(wins.begin(), wins.end(), 0UL)
              << std::endl;
    return 0;
}


int run_mode(int argc, const char *const argv[])
{
    const std::string images_path = argv[2], labels_path = argv[3], feat_path = argv[6], backend_path = argv[9],
                      prefix = argv[10];
    const size_t train = std::stoul(argv[4]), eval = std::stoul(argv[5]), k = std::stoul(argv[7]),
                 classes = std::stoul(argv[8]);
    constexpr int side = 32, in_n = 3 * side * side, ret_side = side, ret_n = maps * side * side, stride = 2;
    constexpr int per_side = (side - patch) / stride + 1, positions = per_side * per_side, fan_in = maps * patch * patch;
    const size_t det_n = k * positions;
    const auto px = read_bin(images_path);
    std::vector<size_t> labels;
    {
        std::ifstream f(labels_path);
        for (size_t v; f >> v;) labels.push_back(v);
    }
    // Detector weights: learned prototypes minus their mean, scaled to integers in ±KNP_DET_WMAX.
    std::vector<float> proto(k * fan_in);
    {
        std::ifstream f(feat_path, std::ios::binary);
        f.read(reinterpret_cast<char *>(proto.data()), proto.size() * 4);
    }
    std::vector<float> mean(fan_in, 0.F);
    for (size_t j = 0; j < k; ++j)
        for (int i = 0; i < fan_in; ++i) mean[i] += proto[j * fan_in + i] / k;
    float max_abs = 1e-6F;
    for (size_t j = 0; j < k; ++j)
        for (int i = 0; i < fan_in; ++i) max_abs = std::max(max_abs, std::fabs(proto[j * fan_in + i] - mean[i]));
    const float det_wmax = env_f("KNP_DET_WMAX", 127);
    std::vector<float> det_w(k * fan_in);
    for (size_t j = 0; j < k; ++j)
        for (int i = 0; i < fan_in; ++i)
            det_w[j * fan_in + i] = std::round((proto[j * fan_in + i] - mean[i]) / max_abs * det_wmax);

    const float out_thr = env_f("KNP_OUT_THR", 4000);
    auto build = [&](bool learning, const std::vector<float> &readout_w, knp::core::UID &image_proj_uid,
                     knp::core::UID &label_proj_uid, knp::core::UID &readout_proj_uid, knp::core::UID &out_uid)
    {
        AltAIPopulation retina([&](size_t) { return plain_neuron(env_f("KNP_RETINA_THR", 250)); }, ret_n);
        AltAIPopulation detectors([&](size_t) { return plain_neuron(env_f("KNP_DET_THR", 2000)); }, det_n);
        auto retina_proj = make_retina(retina.get_uid(), side, side, ret_side, ret_side, 0);
        // Detector index = proto * positions + position; synapse index = detector * fan_in + element.
        DeltaProjection det_proj(
            retina.get_uid(), detectors.get_uid(),
            [&](size_t index) -> std::optional<DeltaProjection::Synapse>
            {
                const size_t neuron = index / fan_in, e = index % fan_in;
                const size_t j = neuron / positions, pos = neuron % positions;
                const int py = static_cast<int>(pos / per_side) * stride, pxx = static_cast<int>(pos % per_side) * stride;
                const int m = static_cast<int>(e / (patch * patch)), dy = static_cast<int>(e % (patch * patch)) / patch,
                          dx = static_cast<int>(e % patch);
                const float weight = det_w[j * fan_in + e];
                if (weight == 0) return std::nullopt;
                return DeltaProjection::Synapse{
                    {weight, 1, knp::synapse_traits::OutputType::EXCITATORY},
                    static_cast<size_t>((m * ret_side + py + dy) * ret_side + pxx + dx), neuron};
            },
            det_n * fan_in);
        image_proj_uid = retina_proj.get_uid();
        knp::framework::Network network;
        if (learning)
        {
            StdpPopulation readout([&](size_t) { return learning_neuron(out_thr); }, classes);
            ResourceProjection rproj(
                detectors.get_uid(), readout.get_uid(),
                [&](size_t i) -> std::optional<ResourceProjection::Synapse>
                {
                    knp::synapse_traits::synapse_parameters<ResourceSynapse> s;
                    s.rule_.w_min_ = w_min;
                    s.rule_.w_max_ = w_max;
                    s.rule_.dopamine_plasticity_period_ = 10;
                    s.weight_ = 0;
                    s.rule_.synaptic_resource_ = resource_from_weight(0);
                    s.delay_ = 1;
                    s.output_type_ = knp::synapse_traits::OutputType::EXCITATORY;
                    return ResourceProjection::Synapse{s, i % det_n, i / det_n};
                },
                det_n * classes);
            rproj.unlock_weights();
            const auto delta = env_list("KNP_DELTA", {1, 4});
            DeltaProjection label_proj(
                knp::core::UID{false}, readout.get_uid(),
                [&](size_t i) {
                    return std::optional(DeltaProjection::Synapse{
                        {delta[0], 1, knp::synapse_traits::OutputType::DOPAMINE}, i, i});
                },
                classes);
            DeltaProjection self_proj(
                readout.get_uid(), readout.get_uid(),
                [&](size_t i) {
                    return std::optional(DeltaProjection::Synapse{
                        {-delta[1], 1, knp::synapse_traits::OutputType::DOPAMINE}, i, i});
                },
                classes);
            label_proj_uid = label_proj.get_uid();
            readout_proj_uid = rproj.get_uid();
            out_uid = readout.get_uid();
            network.add_population(std::move(readout));
            network.add_projection(std::move(rproj));
            network.add_projection(std::move(label_proj));
            network.add_projection(std::move(self_proj));
        }
        else
        {
            AltAIPopulation readout([&](size_t) { return plain_neuron(out_thr); }, classes);
            DeltaProjection rproj(
                detectors.get_uid(), readout.get_uid(),
                [&](size_t i) -> std::optional<DeltaProjection::Synapse>
                {
                    const float weight = readout_w[(i / det_n) * det_n + i % det_n];
                    if (weight == 0) return std::nullopt;
                    return DeltaProjection::Synapse{
                        {weight, 1, knp::synapse_traits::OutputType::EXCITATORY}, i % det_n, i / det_n};
                },
                det_n * classes);
            out_uid = readout.get_uid();
            network.add_population(std::move(readout));
            network.add_projection(std::move(rproj));
        }
        network.add_population(std::move(retina));
        network.add_population(std::move(detectors));
        network.add_projection(std::move(retina_proj));
        network.add_projection(std::move(det_proj));
        return network;
    };

    std::vector<knp::core::messaging::SpikeData> current;
    auto image_gen = [&](size_t first)
    {
        return [&, first](knp::core::Step step) -> knp::core::messaging::SpikeData
        {
            const size_t img = first + step / steps_per_image, local = step % steps_per_image;
            if (local == 0) current = raster(px.data() + img * in_n, in_n);
            if (local >= active_steps) return {};
            return current[local];
        };
    };
    knp::framework::BackendLoader loader;
    std::vector<float> readout_w(det_n * classes, 0.F);
    if (train > 0)
    {
        knp::core::UID image_proj, label_proj, readout_proj, out_uid;
        knp::framework::Model model(build(true, readout_w, image_proj, label_proj, readout_proj, out_uid));
        knp::core::UID in_channel, label_channel;
        model.add_input_channel(in_channel, image_proj);
        model.add_input_channel(label_channel, label_proj);
        auto label_gen = [&](knp::core::Step step) -> knp::core::messaging::SpikeData
        {
            if (step % steps_per_image >= active_steps) return {};
            return {static_cast<uint32_t>(labels[step / steps_per_image])};
        };
        auto backend = loader.load(backend_path);
        knp::framework::ModelExecutor executor(model, backend, {{in_channel, image_gen(0)}, {label_channel, label_gen}});
        const auto t0 = std::chrono::steady_clock::now();
        executor.start(
            [&](size_t step)
            {
                if (step % (1000 * steps_per_image) == 0 && step)
                    std::cout << "trained " << step / steps_per_image << " images, "
                              << std::chrono::duration<double>(std::chrono::steady_clock::now() - t0).count() << " s"
                              << std::endl;
                return step < train * steps_per_image;
            });
        readout_w = trained_weights(backend, readout_proj, det_n, classes);
        std::ofstream(prefix + "_readout.bin", std::ios::binary)
            .write(reinterpret_cast<const char *>(readout_w.data()), readout_w.size() * 4);
        std::cout << "training done in " << std::chrono::duration<double>(std::chrono::steady_clock::now() - t0).count()
                  << " s" << std::endl;
    }
    else
    {
        std::ifstream f(prefix + "_readout.bin", std::ios::binary);
        f.read(reinterpret_cast<char *>(readout_w.data()), readout_w.size() * 4);
    }
    if (eval == 0) return 0;
    knp::core::UID image_proj, label_proj, readout_proj, out_uid;
    knp::framework::Model model(build(false, readout_w, image_proj, label_proj, readout_proj, out_uid));
    knp::core::UID in_channel, out_channel;
    model.add_input_channel(in_channel, image_proj);
    model.add_output_channel(out_channel, out_uid);
    // Evaluation starts right after the training images unless KNP_EVAL_FIRST is given.
    const size_t first = std::getenv("KNP_EVAL_FIRST") ? static_cast<size_t>(env_f("KNP_EVAL_FIRST", 0)) : train;
    auto backend = loader.load(backend_path);
    knp::framework::ModelExecutor executor(model, backend, {{in_channel, image_gen(first)}});
    auto &out = executor.get_loader().get_output_channel(out_channel);
    const auto t0 = std::chrono::steady_clock::now();
    executor.start([&](size_t step) { return step < eval * steps_per_image; });
    std::vector<std::vector<int>> counts(eval, std::vector<int>(classes, 0));
    for (const auto &m : out.update())
    {
        // A spike sent on step s belongs to image s / steps_per_image (output arrives within the same image window).
        const size_t img = m.header_.send_time_ / steps_per_image;
        if (img < eval)
            for (auto i : m.neuron_indexes_) ++counts[img][i];
    }
    std::ofstream csv(prefix + "_eval.csv");
    csv << "index,label";
    for (size_t c = 0; c < classes; ++c) csv << ",c" << c;
    csv << "\n";
    size_t correct = 0, silent = 0;
    for (size_t i = 0; i < eval; ++i)
    {
        const auto &c = counts[i];
        const int best = static_cast<int>(std::max_element(c.begin(), c.end()) - c.begin());
        const bool none = *std::max_element(c.begin(), c.end()) == 0;
        silent += none;
        correct += !none && static_cast<size_t>(best) == labels[first + i];
        csv << first + i << ',' << labels[first + i];
        for (auto v : c) csv << ',' << v;
        csv << "\n";
    }
    std::cout << "eval: " << correct << " / " << eval << " correct, " << silent << " silent, "
              << std::chrono::duration<double>(std::chrono::steady_clock::now() - t0).count() << " s" << std::endl;
    return 0;
}


int main(int argc, const char *const argv[])
{
    const std::string mode = argc > 1 ? argv[1] : "";
    if (mode == "features" && argc == 7) return features_mode(argc, argv);
    if (mode == "run" && argc == 11) return run_mode(argc, argv);
    std::cerr << "usage: see the header of altai_objects.cpp" << std::endl;
    return 1;
}