// Signal/gradient flow diagnostic for the sequence-LSTM kernels, at PRODUCTION shapes.
//
// Why this exists (separate from lstm_seq_gradcheck.cpp): that check proves the BPTT
// math is CORRECT. Correct is not the same as USABLE. A textbook-correct LSTM whose
// forget gate sits at sigmoid(0)=0.5 multiplies the cell state by 0.5 every step, so
// over a 20-bar window it retains 1e-6 of bar 0 in the forward pass and propagates
// 1e-6 of the gradient back to it. It trains to a constant and the training log shows
// a flat IS error - exactly what LSTM/HYBRID did on 2026-07-30.
//
// This measures both directions at the real shapes and sweeps the forget-gate bias:
//   fwd  - how much h_T moves when ONLY bar t is perturbed  (memory reach)
//   bwd  - ||dL/dx_t||                                       (credit assignment reach)
//
// Build (VS x64 dev prompt, this folder):
//   cl /nologo /EHsc /O2 /std:c++17 lstm_seq_flowcheck.cpp WarriorCPU.lib
#include <cstdio>
#include <cmath>
#include <vector>
#include <random>
#include "WarriorCPU.h"

static const int H = 64;   // derived hidden size at H1 defaults (see ComputeLstmHiddenSize)
static const int IW = 21;  // per-bar feature count for the plain-LSTM topology
static const int T = 20;   // m_historyBars

int main()
{
   CpuHandle h = CPU_Init(2);
   const int perGate = H * (H + IW + 1);
   const int wTotal = 4 * perGate;
   const int bias = H + IW;   // bias column within a row

   int wB = CPU_BufferCreate(h, wTotal);
   int inB = CPU_BufferCreate(h, T * IW);
   int cgB = CPU_BufferCreate(h, T * 4 * H);
   int ccB = CPU_BufferCreate(h, T * H);
   int chB = CPU_BufferCreate(h, T * H);
   int outB = CPU_BufferCreate(h, H);
   int ogB = CPU_BufferCreate(h, H);
   int wgB = CPU_BufferCreate(h, wTotal);
   int igB = CPU_BufferCreate(h, T * IW);

   std::mt19937 rng(1234);
   // LeCun-uniform, matching CNeuronLSTMOCL::SetInputs exactly.
   const double scale = 1.0 / std::sqrt((double)(H + IW) + 1.0);
   std::uniform_real_distribution<double> wd(-scale, scale);
   std::vector<double> w0((size_t)wTotal);
   for(int i = 0; i < wTotal; i++)
      w0[i] = wd(rng);

   std::normal_distribution<double> xd(0.0, 1.0);
   std::vector<double> x((size_t)(T * IW));
   for(int i = 0; i < T * IW; i++)
      x[i] = xd(rng);

   std::vector<double> dL((size_t)H, 0.0);
   for(int k = 0; k < H; k++)
      dL[k] = ((k % 7) - 3) / 3.0;   // asymmetric upstream gradient

   // Must track LSTM_FORGET_BIAS_INIT in AI\Network.mqh. The sweep exists to show WHY that value is
   // what it is; the verdict at the end is what makes this a regression test rather than a printout.
   const double SHIPPED = 2.0;
   const double MIN_REACH = 0.05;   // oldest bar must retain >=5% of the newest bar's influence
   double shippedFwd = 0.0, shippedBwd = 0.0;

   const double fb[] = { 0.0, 1.0, SHIPPED };
   std::printf("shapes: H=%d IW=%d T=%d   weights=%d\n\n", H, IW, T, wTotal);

   for(double b : fb)
   {
      std::vector<double> w = w0;
      for(int hid = 0; hid < H; hid++)
         w[(size_t)(0 * perGate + hid * (H + IW + 1) + bias)] = b;   // gate 0 == forget

      CPU_BufferWrite(h, wB, w.data(), wTotal);
      CPU_BufferWrite(h, inB, x.data(), T * IW);
      CPU_LSTMSeqForward(h, wB, inB, cgB, ccB, chB, outB, H, IW, T);
      std::vector<double> base((size_t)H);
      CPU_BufferRead(h, outB, base.data(), H);

      // --- forward reach: perturb ONLY bar t, measure ||dh_T||
      double fwd0 = 0.0, fwdLast = 0.0;
      for(int t : {0, T - 1})
      {
         std::vector<double> xp = x;
         for(int j = 0; j < IW; j++)
            xp[(size_t)(t * IW + j)] += 0.01;
         CPU_BufferWrite(h, inB, xp.data(), T * IW);
         CPU_LSTMSeqForward(h, wB, inB, cgB, ccB, chB, outB, H, IW, T);
         std::vector<double> pert((size_t)H);
         CPU_BufferRead(h, outB, pert.data(), H);
         double s = 0.0;
         for(int k = 0; k < H; k++)
            s += (pert[k] - base[k]) * (pert[k] - base[k]);
         (t == 0 ? fwd0 : fwdLast) = std::sqrt(s);
      }

      // --- backward reach: ||dL/dx_t|| at the first and last step
      CPU_BufferWrite(h, inB, x.data(), T * IW);
      CPU_LSTMSeqForward(h, wB, inB, cgB, ccB, chB, outB, H, IW, T);
      CPU_BufferWrite(h, ogB, dL.data(), H);
      CPU_LSTMSeqBackward(h, wB, inB, cgB, ccB, chB, ogB, wgB, igB, H, IW, T);
      std::vector<double> dX((size_t)(T * IW));
      CPU_BufferRead(h, igB, dX.data(), T * IW);
      double n0 = 0.0, nL = 0.0;
      for(int j = 0; j < IW; j++)
      {
         n0 += dX[(size_t)j] * dX[(size_t)j];
         nL += dX[(size_t)((T - 1) * IW + j)] * dX[(size_t)((T - 1) * IW + j)];
      }
      n0 = std::sqrt(n0);
      nL = std::sqrt(nL);

      std::printf("forget-gate bias %+.1f (sigmoid=%.3f)\n", b, 1.0 / (1.0 + std::exp(-b)));
      std::printf("   fwd  ||dh_T|| from bar 0  = %.3e   from bar %d = %.3e   ratio %.2e\n",
                  fwd0, T - 1, fwdLast, fwdLast > 0 ? fwd0 / fwdLast : 0.0);
      std::printf("   bwd  ||dL/dx_0||          = %.3e   ||dL/dx_%d|| = %.3e   ratio %.2e\n\n",
                  n0, T - 1, nL, nL > 0 ? n0 / nL : 0.0);

      if(b == SHIPPED)
      {
         shippedFwd = fwdLast > 0 ? fwd0 / fwdLast : 0.0;
         shippedBwd = nL > 0 ? n0 / nL : 0.0;
      }
   }

   CPU_Shutdown(h);

   bool ok = (shippedFwd >= MIN_REACH && shippedBwd >= MIN_REACH);
   std::printf("shipped forget-gate bias %+.1f: fwd reach %.3f, bwd reach %.3f (need >= %.2f) -> %s\n",
               SHIPPED, shippedFwd, shippedBwd, MIN_REACH, ok ? "PASS" : "FAIL");
   if(!ok)
      std::printf("   The recurrence cannot see across its own window. Expect a flat IS error and\n"
                  "   OOS recall Neutral:100%%. Check LSTM_FORGET_BIAS_INIT and the cell-state update.\n");
   return ok ? 0 : 1;
}