//+------------------------------------------------------------------+
//| dense_backprop_check.cpp                                          |
//|                                                                   |
//| Offline math check for the 2026-08-11 dense hidden-gradient        |
//| transpose fix. Same shape as batch_accum_check beside it: link     |
//| straight against WarriorCPU.dll, drive the real exports, compare   |
//| against an independent reference computed here.                    |
//|                                                                   |
//| WHAT THIS PROVES, precisely:                                       |
//|   1. CPU_CalcHiddenGradient is CONSISTENT WITH CPU_FeedForward:    |
//|      with no activation on either side, the backward kernel's      |
//|      ig[i] equals the central finite difference of                  |
//|      J(x) = Sum_k g_k * forward(x)_k through the REAL forward      |
//|      kernel, on square AND non-square layer pairs. This is the     |
//|      exact property the pre-fix kernel violated: it read the       |
//|      weight matrix transposed (square case) or mis-strided         |
//|      (non-square case), so its output was NOT the gradient of any  |
//|      function the forward pass computes.                           |
//|   2. The activation-derivative branches (tanh / sigmoid / PReLU)   |
//|      match a direct transcription of the kernel's documented        |
//|      clamp-reformulated formulas applied to the correctly-indexed  |
//|      propagated sum.                                               |
//|                                                                   |
//| WHAT IT DOES NOT PROVE - see [[feedback_verify_in_situ_not_       |
//| offline]]: that a layer TRAINS in the assembled network, and       |
//| nothing about the OpenCL/DirectML mirrors (those share the fixed   |
//| indexing by transcription; the in-situ dW/W report on a real era   |
//| is the check that covers whichever tier actually runs).            |
//|                                                                   |
//| Build (from a plain cmd, after build_cpu.bat):                     |
//|   cl /nologo /EHsc /O2 /std:c++17 dense_backprop_check.cpp WarriorCPU.lib
//+------------------------------------------------------------------+
#include <cstdio>
#include <cmath>
#include <vector>
#include <random>
#include <algorithm>

// WarriorCPU.h is deliberately NOT included - see batch_accum_check.cpp for why.
// Signatures must match WarriorCPU.h exactly.
typedef long long CpuHandle;
extern "C"
{
   __declspec(dllimport) CpuHandle __stdcall CPU_Init(int threads);
   __declspec(dllimport) void __stdcall CPU_Shutdown(CpuHandle ctx);
   __declspec(dllimport) int __stdcall CPU_BufferCreate(CpuHandle ctx, int elementCount);
   __declspec(dllimport) int __stdcall CPU_BufferWrite(CpuHandle ctx, int handle, const double *data, int count);
   __declspec(dllimport) int __stdcall CPU_BufferRead(CpuHandle ctx, int handle, double *data, int count);
   __declspec(dllimport) void __stdcall CPU_BufferFree(CpuHandle ctx, int handle);
   __declspec(dllimport) int __stdcall CPU_FeedForward(CpuHandle ctx, int wHandle, int iHandle, int oHandle,
         int inputs, int activation);
   __declspec(dllimport) int __stdcall CPU_CalcHiddenGradient(CpuHandle ctx, int wHandle, int gHandle, int oHandle,
         int igHandle, int outputs, int activation, int count);
}

namespace
{
// Must match WarriorCPU.cpp (NOT AI\Network.cl's 1e-3f - the CPU DLL floors at 1e-4).
const double MIN_ACTIVATION_DERIVATIVE = 1.0e-4;

int g_failures = 0;

void CheckMax(const char *what, double maxDiff, double tol)
{
   if(!(maxDiff <= tol))
   {
      std::printf("  FAIL %-52s max |diff| %.3g > %.3g\n", what, maxDiff, tol);
      ++g_failures;
   }
   else
      std::printf("  ok   %-52s max |diff| %.3g\n", what, maxDiff);
}

// J(x) = Sum_k g_k * forward(x)_k through the real forward kernel, activation NONE (-1).
double ScalarJ(CpuHandle ctx, int wH, int iH, int oH, int thisN, int nextN,
               const std::vector<double> &x, const std::vector<double> &g)
{
   CPU_BufferWrite(ctx, iH, x.data(), thisN);
   if(!CPU_FeedForward(ctx, wH, iH, oH, thisN, -1))
   {
      std::printf("  FAIL CPU_FeedForward returned 0\n");
      ++g_failures;
      return 0.0;
   }
   std::vector<double> out(nextN);
   CPU_BufferRead(ctx, oH, out.data(), nextN);
   double j = 0.0;
   for(int k = 0; k < nextN; k++)
      j += g[k] * out[k];
   return j;
}

void RunPair(CpuHandle ctx, int thisN, int nextN, std::mt19937 &rng)
{
   std::uniform_real_distribution<double> uni(-1.0, 1.0);
   std::vector<double> w((size_t)(thisN + 1) * nextN), x(thisN), g(nextN);
   for(auto &v : w) v = uni(rng);
   for(auto &v : x) v = uni(rng);
   for(auto &v : g) v = uni(rng);

   int wH  = CPU_BufferCreate(ctx, (int)w.size());
   int iH  = CPU_BufferCreate(ctx, thisN);
   int oH  = CPU_BufferCreate(ctx, nextN);
   int gH  = CPU_BufferCreate(ctx, nextN);
   int outH = CPU_BufferCreate(ctx, thisN);   // this layer's own output vector (= x here)
   int igH = CPU_BufferCreate(ctx, thisN);
   CPU_BufferWrite(ctx, wH, w.data(), (int)w.size());
   CPU_BufferWrite(ctx, gH, g.data(), nextN);
   CPU_BufferWrite(ctx, outH, x.data(), thisN);

   //--- 1) forward/backward consistency by central finite differences, activation NONE
   if(!CPU_CalcHiddenGradient(ctx, wH, gH, outH, igH, nextN, -1, thisN))
   {
      std::printf("  FAIL CPU_CalcHiddenGradient returned 0 (%dx%d)\n", thisN, nextN);
      ++g_failures;
      return;
   }
   std::vector<double> ig(thisN);
   CPU_BufferRead(ctx, igH, ig.data(), thisN);

   const double h = 1.0e-6;
   double maxDiff = 0.0;
   for(int i = 0; i < thisN; i++)
   {
      std::vector<double> xp = x, xm = x;
      xp[i] += h;
      xm[i] -= h;
      double fd = (ScalarJ(ctx, wH, iH, oH, thisN, nextN, xp, g)
                   - ScalarJ(ctx, wH, iH, oH, thisN, nextN, xm, g)) / (2.0 * h);
      maxDiff = std::max(maxDiff, std::fabs(fd - ig[i]));
   }
   char label[96];
   std::snprintf(label, sizeof(label), "fd-vs-kernel %dx%d (activation NONE)", thisN, nextN);
   CheckMax(label, maxDiff, 1e-7);

   //--- 2) activation branches vs a direct transcription on the correctly-indexed sum
   for(int act = 0; act <= 2; act++)
   {
      //--- for tanh/sigmoid the stored output must be inside the activation's range or the
      //--- clamp-reformulation departs from the plain derivative by design; keep |o| < 1.
      std::vector<double> o(thisN);
      for(auto &v : o) v = 0.9 * uni(rng);
      CPU_BufferWrite(ctx, outH, o.data(), thisN);
      if(!CPU_CalcHiddenGradient(ctx, wH, gH, outH, igH, nextN, act, thisN))
      {
         std::printf("  FAIL CPU_CalcHiddenGradient returned 0 (act %d)\n", act);
         ++g_failures;
         continue;
      }
      CPU_BufferRead(ctx, igH, ig.data(), thisN);
      double maxD = 0.0;
      for(int i = 0; i < thisN; i++)
      {
         double sum = 0.0;
         for(int k = 0; k < nextN; k++)
            sum += g[k] * w[(size_t)k * (thisN + 1) + i];
         double ov = o[i];
         if(act == 0)
         {
            sum = std::max(-1.0, std::min(1.0, sum + ov)) - ov;
            sum *= std::max(MIN_ACTIVATION_DERIVATIVE, 1.0 - ov * ov);
         }
         else if(act == 1)
         {
            sum = std::max(0.0, std::min(1.0, sum + ov)) - ov;
            sum *= std::max(MIN_ACTIVATION_DERIVATIVE, ov * (1.0 - ov));
         }
         else if(act == 2 && ov < 0.0)
            sum *= 0.01;
         maxD = std::max(maxD, std::fabs(sum - ig[i]));
      }
      std::snprintf(label, sizeof(label), "activation branch %d, %dx%d", act, thisN, nextN);
      CheckMax(label, maxD, 1e-12);
   }

   CPU_BufferFree(ctx, wH);
   CPU_BufferFree(ctx, iH);
   CPU_BufferFree(ctx, oH);
   CPU_BufferFree(ctx, gH);
   CPU_BufferFree(ctx, outH);
   CPU_BufferFree(ctx, igH);
}
} // namespace

int main()
{
   CpuHandle ctx = CPU_Init(2);
   if(!ctx)
   {
      std::printf("CPU_Init failed\n");
      return 1;
   }
   std::mt19937 rng(20260811);
   std::printf("dense_backprop_check:\n");
   //--- square (the transpose case), the real head boundary (wide -> 3), a taper boundary, and a
   //--- width that is a multiple of 4 (the OpenCL Adam bias-skip shape, kept here for coverage).
   RunPair(ctx, 8, 8, rng);
   RunPair(ctx, 64, 3, rng);
   RunPair(ctx, 33, 64, rng);
   RunPair(ctx, 5, 3, rng);
   CPU_Shutdown(ctx);
   if(g_failures)
   {
      std::printf("%d FAILURE(S)\n", g_failures);
      return 1;
   }
   std::printf("all checks passed\n");
   return 0;
}