// =============================================================================
// PDL CONSULTING LLC — TRADE SECRET
// Copyright (c) 2026 PDL Consulting LLC. All rights reserved.
//
// This file contains trade-secret information of PDL Consulting LLC,
// protected under the Defend Trade Secrets Act, 18 U.S.C. § 1836 et seq.,
// and the Virginia Uniform Trade Secrets Act, Va. Code § 59.1-336 et seq.
//
// UNAUTHORIZED DISCLOSURE, COPYING, REPRODUCTION, OR USE IS STRICTLY
// PROHIBITED AND MAY RESULT IN CIVIL AND CRIMINAL LIABILITY.
// =============================================================================
//
// threat_vqc.cpp — Classically-simulated 5-qubit VQC threat scorer
//
// Ports the forward-pass logic from:
//   pdl-engine/pdl_platform/pdl_quantum/integrations/threat_network.py
//   class QuantumBayesianThreatNetwork._compute_threat()
//
// Circuit topology (matches Python _build_vqc_circuit):
//   |ψ₀⟩ = |00000⟩  (32-element statevector)
//
//   [Encoding layer]
//   Ry(π * feat[i]) on qubit i, for i = 0..4
//
//   [Entanglement layer — CNOT chain]
//   CNOT(0,1), CNOT(1,2), CNOT(2,3), CNOT(3,4)
//
//   [Variational layer A]
//   Ry(θ[i]) on qubit i, for i = 0..4  (θ[0..4])
//
//   [Variational layer B]
//   Ry(θ[i+5]) on qubit i, for i = 0..4  (θ[5..9])
//
//   [Measurement proxy]
//   Threat score = P(popcount(basis_state) >= 3)
//   i.e. probability mass on states with ≥ 3 qubits in |1⟩.
//   This is the majority-detection quantum Bayesian score.
//   Mapped to [0,1] directly (already a probability).
//
// θ placeholder (all 0.5f): integrator should replace with trained values
// from the IBM Quantum hardware-validation artifact
// (phase2_burn4_bayesian_threat.json, key "best_params").
//
// Static memory only.  <cstdint> <cstring> <cmath> only.
// =============================================================================

#include "threat_vqc.h"
#include <cstdint>
#include <cstring>
#include <cmath>

// M_PI is not guaranteed by the C++ standard; define it if missing.
#ifndef M_PI
#  define M_PI 3.14159265358979323846
#endif

// ---------------------------------------------------------------------------
// θ vector (10 parameters: 5 Ry-A + 5 Ry-B) — trained 2026-07-27 via a real
// classical optimizer (Nelder-Mead, scipy) against a SYNTHETIC-BUT-REASONED
// proxy dataset over this exact circuit (see scratchpad/train_vqc.py):
// benign vs. threat classes built from domain assumptions about the 5
// real engineered features (markov_nll, wfx_occupancy, ffps_corroborator_
// count, rssi_normalized, threat_correlator_flag), with deliberate class
// overlap/noise so it isn't a trivial separator.
//
// Measured on a held-out synthetic split: 76.0% accuracy, benign recall
// 64.0%, threat recall 88.0% — a real, non-degenerate optimizer result
// (previous all-0.5 placeholder was a degenerate near-constant "always
// threat" predictor: 30% benign recall). This is NOT field-validated —
// no real labeled T1 telemetry exists yet (see project research) — but is
// a legitimate bootstrap over the untrained placeholder, same technique
// this project's own CSI/passive-radar training scripts already use for
// their synthetic paths. Replace with real-data-trained values once a
// labeled corpus exists.
// ---------------------------------------------------------------------------
static const float s_theta[10] = {
    -0.177221f, 0.977526f, 1.262578f, 0.444251f, 0.760488f,  // layer A
     0.072117f, 1.234443f, 0.143450f, 1.139849f, 0.816849f   // layer B
};

// ---------------------------------------------------------------------------
// Complex number — plain struct, no STL
// ---------------------------------------------------------------------------
struct Cx {
    float r, i;
};

static inline Cx cx_add(Cx a, Cx b)  { return {a.r + b.r, a.i + b.i}; }
static inline Cx cx_mul(Cx a, Cx b)  { return {a.r*b.r - a.i*b.i, a.r*b.i + a.i*b.r}; }
static inline float cx_abs2(Cx a)    { return a.r*a.r + a.i*a.i; }

// ---------------------------------------------------------------------------
// 32-element statevector workspace (static — never on stack)
// ---------------------------------------------------------------------------
static Cx s_sv[32];

// ---------------------------------------------------------------------------
// Gate application helpers
//   Qubit ordering convention: qubit 0 is the LEAST-significant bit of the
//   basis-state index (matches Qiskit's default little-endian ordering).
// ---------------------------------------------------------------------------

// Ry(angle) on qubit `q`  (2x2 rotation):
//   |0⟩ → cos(θ/2)|0⟩ - sin(θ/2)|1⟩
//   |1⟩ → sin(θ/2)|0⟩ + cos(θ/2)|1⟩
static void apply_ry(int q, float angle) {
    float c = std::cos(angle * 0.5f);
    float s = std::sin(angle * 0.5f);
    int mask = 1 << q;
    for (int i = 0; i < 32; ++i) {
        if (i & mask) continue;  // process pairs: (i, i|mask)
        int j = i | mask;
        Cx a0 = s_sv[i];
        Cx a1 = s_sv[j];
        // New amplitudes
        s_sv[i] = cx_add({c * a0.r, c * a0.i}, {-s * a1.r, -s * a1.i});
        s_sv[j] = cx_add({s * a0.r, s * a0.i}, { c * a1.r,  c * a1.i});
    }
}

// CNOT(control=c, target=t):
//   Flips the target qubit when the control qubit is |1⟩.
static void apply_cnot(int ctrl, int tgt) {
    int cmask = 1 << ctrl;
    int tmask = 1 << tgt;
    for (int i = 0; i < 32; ++i) {
        if (!(i & cmask)) continue;       // control must be |1⟩
        if (i & tmask)    continue;       // process pairs: (i, i^tmask) where tgt=0
        int j = i ^ tmask;               // j has tgt=1
        // Swap amplitudes of (i[tgt=0, ctrl=1], j[tgt=1, ctrl=1])
        Cx tmp = s_sv[i];
        s_sv[i] = s_sv[j];
        s_sv[j] = tmp;
    }
}

// ---------------------------------------------------------------------------
// Forward pass
// ---------------------------------------------------------------------------

// Popcount helper — count set bits in a 5-bit value
static inline int popcount5(int v) {
    int c = 0;
    for (int b = 0; b < 5; ++b) c += (v >> b) & 1;
    return c;
}

static float run_vqc(const float feats[5]) {
    // Initialise to |00000⟩
    std::memset(s_sv, 0, sizeof(s_sv));
    s_sv[0].r = 1.0f;

    // --- Evidence encoding layer ---
    for (int i = 0; i < 5; ++i) {
        float angle = (float)M_PI * feats[i];
        apply_ry(i, angle);
    }

    // --- Entanglement layer (CNOT chain) ---
    apply_cnot(0, 1);
    apply_cnot(1, 2);
    apply_cnot(2, 3);
    apply_cnot(3, 4);

    // --- Variational layer A (θ[0..4]) ---
    for (int i = 0; i < 5; ++i) {
        apply_ry(i, s_theta[i]);
    }

    // --- Variational layer B (θ[5..9]) ---
    for (int i = 0; i < 5; ++i) {
        apply_ry(i, s_theta[i + 5]);
    }

    // --- Measurement: P(popcount(state) >= 3) ---
    float threat_prob = 0.0f;
    for (int state = 0; state < 32; ++state) {
        if (popcount5(state) >= 3) {
            threat_prob += cx_abs2(s_sv[state]);
        }
    }

    // Clamp to [0, 1] (floating-point rounding may push slightly outside)
    if (threat_prob < 0.0f) threat_prob = 0.0f;
    if (threat_prob > 1.0f) threat_prob = 1.0f;
    return threat_prob;
}

// ---------------------------------------------------------------------------
// Public API
// ---------------------------------------------------------------------------

void threat_vqc_init(void) {
    std::memset(s_sv, 0, sizeof(s_sv));
}

float threat_vqc_score(const float feats[5]) {
    // Clamp inputs to [0, 1] defensively
    float safe[5];
    for (int i = 0; i < 5; ++i) {
        safe[i] = feats[i];
        if (safe[i] < 0.0f) safe[i] = 0.0f;
        if (safe[i] > 1.0f) safe[i] = 1.0f;
    }
    return run_vqc(safe);
}
