// fp32 throughout - this project targets consumer GPUs (e.g. AMD Polaris/RX 580) where double
// precision runs at a small fraction of float throughput. No kernel below needs fp64 range or
// precision: weight/delta clamps stay within +-100, gradients/activations are float32-safe, and
// the Adam epsilon guard (vt > 0 ? vt : l*10) isn't a tiny-epsilon comparison sensitive to
// fp32 vs fp64 underflow. Removing cl_khr_fp64 means any accidental double slipping back in fails
// kernel compilation loudly instead of silently running at 1/16th speed.
// Tightened from 1.0e6 - see AI\Network.mqh's matching MAX_WEIGHT comment for the full rationale
// (that ceiling was so loose it never engaged before training had already gone unstable). 100.0
// matches the equivalent clamp in Dmitriy Gizlyk's reference NeuroNet.cl engine.
#define MAX_WEIGHT 100.0f
// Floor on the magnitude of a saturated tanh/sigmoid unit's activation derivative - see AI\Network.mqh's
// matching MIN_ACTIVATION_DERIVATIVE comment for the full rationale (a unit pinned near its activation
// extremes must still receive some corrective gradient, however small, rather than exactly zero).
// 2026-07-27: Increased from 1.0e-4f to 1.0e-3f to strengthen the escape signal through saturated hidden
// neurons. On the fp32 OpenCL backend, a saturated sigmoid output layer (the "neutral collapse from era 0"
// bug) propagated only this floor value backward through every hidden layer, which at 1e-4 was too weak
// to ever pull weights out of saturation. The output-layer saturation root cause is now fixed (NONE/logit
// activation instead of SIGMOID — see ExpertSignalAIBase.mqh's BuildFreshTopology), but this increase
// provides a 10× stronger safety net against any remaining fp32-specific saturation in hidden layers.
#define MIN_ACTIVATION_DERIVATIVE 1.0e-3f
// Decoupled (AdamW-style) weight decay applied inside every Adam kernel below - see WarriorCPU.cpp's
// matching WEIGHT_DECAY comment for the full rationale (unbounded slow weight growth across hundreds
// of oversampled training eras causing periodic collapse-then-recover cycles). 0.01 is the standard
// AdamW default.
#define WEIGHT_DECAY 0.001f
// Per-step update clip - see WarriorCPU.cpp's matching MAX_WEIGHT_DELTA comment for the full
// rationale (weight decay alone didn't stop the collapse cycles - they turned out to be sudden Adam
// overshoot events, most likely from 5x back-to-back oversampling replay building artificially
// correlated momentum). Applied to the raw delta BEFORE it's added to the weight, unlike MAX_WEIGHT.
#define MAX_WEIGHT_DELTA 0.1f
//+------------------------------------------------------------------+
//|                                                                  |
//+------------------------------------------------------------------+
__kernel void FeedForward(__global float *matrix_w,
                          __global float *matrix_i,
                          __global float *matrix_o,
                          int inputs, int activation)
  {
   int i = get_global_id(0);
   float sum = 0.0f;
   float4 inp, weight;
   int shift = (inputs + 1) * i;
   for(int k = 0; k <= inputs; k = k + 4)
     {
      switch(inputs - k)
        {
         case 0:
            inp = (float4)(1, 0, 0, 0);
            weight = (float4)(matrix_w[shift + k], 0, 0, 0);
            break;
         case 1:
            inp = (float4)(matrix_i[k], 1, 0, 0);
            weight = (float4)(matrix_w[shift + k], matrix_w[shift + k + 1], 0, 0);
            break;
         case 2:
            inp = (float4)(matrix_i[k], matrix_i[k + 1], 1, 0);
            weight = (float4)(matrix_w[shift + k], matrix_w[shift + k + 1], matrix_w[shift + k + 2], 0);
            break;
         case 3:
            inp = (float4)(matrix_i[k], matrix_i[k + 1], matrix_i[k + 2], 1);
            weight = (float4)(matrix_w[shift + k], matrix_w[shift + k + 1], matrix_w[shift + k + 2], matrix_w[shift + k + 3]);
            break;
         default:
            inp = (float4)(matrix_i[k], matrix_i[k + 1], matrix_i[k + 2], matrix_i[k + 3]);
            weight = (float4)(matrix_w[shift + k], matrix_w[shift + k + 1], matrix_w[shift + k + 2], matrix_w[shift + k + 3]);
            break;
        }
      sum += dot(inp, weight);
     }
   switch(activation)
     {
      case 0:
         sum = tanh(sum);
         break;
      case 1:
         sum = 1 / (1 + exp(-clamp(sum, -50.0f, 50.0f)));
         break;
      case 2: // PReLU, param=0.01 - matches FeedForwardConv's case 2 below
         if(sum < 0)
            sum *= 0.01f;
         break;
     }
   matrix_o[i] = sum;
   barrier(CLK_GLOBAL_MEM_FENCE);
  }
//+------------------------------------------------------------------+
//|                                                                  |
//+------------------------------------------------------------------+
__kernel void CaclOutputGradient(__global float *matrix_t,
                                 __global float *matrix_o,
                                 __global float *matrix_ig,
                                 int activation)
  {
   int i = get_global_id(0);
   float temp = 0;
   float out = matrix_o[i];
   switch(activation)
     {
      case 0:
         // Deliberately NOT multiplied by the tanh derivative (1-out^2): that
         // factor vanishes as out approaches +-1, which is exactly where a
         // TANH output neuron needs to converge for a +-1 target (e.g. the
         // buy/sell extremes of a single-neuron regression head), stalling
         // training right when it matters most. Using the raw (target-out)
         // error here is the same fix as pairing sigmoid with cross-entropy -
         // a bounded activation whose output error isn't re-damped by its
         // own saturating derivative.
         temp = clamp(matrix_t[i], -1.0f, 1.0f) - out;
         break;
      case 1:
         // Also deliberately NOT multiplied by the sigmoid derivative out*(1-out) - same
         // reasoning as case 0 above. This is the classification output layer (3 neurons,
         // one-hot 0/1 targets - see BuildFreshTopology()); (target-out) unscaled is exactly
         // the sigmoid+binary-cross-entropy gradient (the loss derivative and the sigmoid
         // derivative algebraically cancel), whereas re-multiplying by out*(1-out) here is
         // the MSE-with-sigmoid formula, which damps toward zero as out approaches 0 or 1 -
         // exactly where a 0/1-target neuron needs to converge. With this damping (plus the
         // 0.00000001 underflow guard at the saturated extremes), the classification output
         // layer's effective learning rate was throttled to a small fraction of the
         // regression path's, observed in practice as all 3 output neurons converging to
         // (and never escaping) an identical value - softmax of three equal logits is always
         // exactly 1/3 each, with ties resolved to index 0 ("Buy" wins by default).
         temp = clamp(matrix_t[i], 0.0f, 1.0f) - out;
         break;
      case 2:
         // PReLU, param=0.01 - unlike tanh/sigmoid, PReLU is unbounded and never saturates,
         // so none of the clamp-to-range/derivative-skip workarounds above are needed here:
         // plain (target-out) scaled by the derivative at this output is exact.
         temp = (matrix_t[i] - out) * (out >= 0 ? 1.0f : 0.01f);
         break;
      default:
         // NONE (raw logits - see the 3-class softmax classification output layer in
         // BuildFreshTopology()/NativeActivationCode()'s NONE case, which maps here): identity
         // activation, derivative 1 everywhere, so the plain (target-out) error passes through
         // unscaled - same formula CNeuron::calcOutputGradients() already uses unconditionally
         // on the un-accelerated CPU fallback path. Leaving this case unhandled left `temp` at
         // its 0.0 initializer, i.e. every softmax-classification output neuron got a
         // permanently zero gradient on this backend - which, since calcHiddenGradients()
         // propagates backward FROM the output layer's gradient, silently froze the entire
         // network (not just the output layer) at its random initial weights whenever an
         // OpenCL-accelerated tier was active.
         temp = matrix_t[i] - out;
         break;
     }
   matrix_ig[i] = temp;
  }
//+------------------------------------------------------------------+
//|                                                                  |
//+------------------------------------------------------------------+
__kernel void CaclHiddenGradient(__global float *matrix_w,
                                 __global float *matrix_g,
                                 __global float *matrix_o,
                                 __global float *matrix_ig,
                                 int outputs, int activation)
  {
   int i = get_global_id(0);
   // This layer's neuron count. Dispatch is exactly Neurons() (no +1: biases receive no backprop
   // gradient, and the old +1 work-item only ever produced an out-of-bounds read of matrix_o).
   int inputs = get_global_size(0);
   float sum = 0;
   float out = matrix_o[i];
   // matrix_w is THIS layer's outgoing weight matrix, laid out exactly as FeedForward consumes it:
   // one row per NEXT-layer neuron k, row stride (inputs + 1), so the weight FROM this neuron i
   // INTO next-layer neuron k is matrix_w[k * (inputs + 1) + i] - a column read, dL/dout_i =
   // sum_k g[k] * W[k][i]. Until 2026-08-11 this kernel read matrix_w[(outputs + 1) * i + k]: the
   // TRANSPOSE for square layers and, for the non-square boundaries this EA actually builds
   // (tapered stacks, the 3-neuron head), a mis-strided walk that ran past the buffer's end. Every
   // gradient crossing a dense boundary on its way down - the entire learning signal reaching the
   // BN/conv/LSTM front ends - passed through that fixed wrong matrix, i.e. feedback-alignment
   // dynamics rather than backprop (which is why nets still "learned something" and it survived).
   // The book reference (references\...\NeuroNet_DNG\NeuroNet.cl CalcHiddenGradient) has the
   // corrected form; this repo's kernel descended from the earlier article version. Mirrored in
   // DirectML\WarriorCPU.cpp (CPU_CalcHiddenGradient) and WarriorDML.cpp (kHlslHiddenGradient) -
   // the three must stay in lockstep.
   for(int k = 0; k < outputs; k++)
      sum += matrix_g[k] * matrix_w[k * (inputs + 1) + i];
   switch(activation)
     {
      case 0:
         sum = clamp(sum + out, -1.0f, 1.0f) - out;
         sum = sum * fmax(MIN_ACTIVATION_DERIVATIVE, 1 - pow(out, 2.0f));
         break;
      case 1:
         sum = clamp(sum + out, 0.0f, 1.0f) - out;
         sum = sum * fmax(MIN_ACTIVATION_DERIVATIVE, out * (1 - out));
         break;
      case 2:
         // PReLU, param=0.01 - unbounded, non-saturating, so the plain chain rule applies
         // directly with no clamp/implied-target reformulation needed (contrast with the
         // tanh/sigmoid cases above, which exist specifically to dodge saturation).
         sum = sum * (out >= 0 ? 1.0f : 0.01f);
         break;
     }
   matrix_ig[i] = sum;
  }
//+------------------------------------------------------------------+
//|                                                                  |
//+------------------------------------------------------------------+
__kernel void UpdateWeightsMomentum(__global float *matrix_w,
                                    __global float *matrix_g,
                                    __global float *matrix_i,
                                    __global float *matrix_dw,
                                    int inputs, float learning_rates, float momentum, int optimizer)
  {
   int i = get_global_id(0);
   int j = get_global_id(1);
   int wi = i * (inputs + 1) + j;
   // `optimizer` is retained (unused) only to keep this kernel's argument list binary-compatible with
   // the already-deployed WarriorDML.dll/WarriorCPU.dll exports. It used to select an index-parity sign
   // flip on the gradient ("DFA"), which put half of every weight tensor into permanent gradient ASCENT
   // - see ENUM_OPTIMIZATION's comment in AI\Network.mqh. Always plain descent now; callers pass 0.
   float delta = learning_rates * (matrix_g[i] * (j < inputs ? matrix_i[j] : 1)) + momentum * matrix_dw[wi];
   matrix_dw[wi] = delta;
   // Unlike UpdateWeightsConvMomentum/UpdateWeightsConvAdam/LSTM_UpdateWeightsAdam below, this dense
   // kernel (and UpdateWeightsAdam right after it) had never gotten the MAX_WEIGHT clamp - a gradient
   // spike could drive a weight to +-Infinity with nothing to stop it, and the next Adam step turns
   // Infinity into NaN (Inf/Inf), which then propagates through every FeedForward sum touching that
   // weight and never recovers. Dense layers are exactly where the classification output head lives,
   // so this is the most consequential of the three missing clamps, not the least.
   matrix_w[wi] = clamp(matrix_w[wi] + delta, -MAX_WEIGHT, MAX_WEIGHT);
  };
//+------------------------------------------------------------------+
//|                                                                  |
//+------------------------------------------------------------------+
__kernel void UpdateWeightsAdam(__global float *matrix_w,
                                __global const float *matrix_g,
                                __global const float *matrix_i,
                                __global float *matrix_m,
                                __global float *matrix_v,
                                const int inputs, const float l, const float b1, const float b2)
  {
   const int i = get_global_id(0);
   const int j = get_global_id(1);
   // Each work-item owns up to 4 consecutive slots of weight row i. The row is (inputs + 1) wide -
   // the last column is the bias, whose implicit input is 1 - and dispatch dim 1 is sized on
   // ceil((inputs + 1) / 4), so the bias column is always reachable. Until 2026-08-11 this kernel
   // had two defects the book reference had already fixed: the input for slot group j was read at
   // matrix_i[j] instead of matrix_i[j * 4] (every group past the first paired its weights with the
   // WRONG input - a corrupted outer product), and the tail switch tested (inputs - j * 4) with
   // dim 1 sized on ceil(inputs / 4), which made the bias column unreachable whenever
   // inputs % 4 == 0 (32/64-wide layers: dense biases simply never trained on OpenCL). The batched
   // path (AccumulateWeightGrad + ApplyAccumAdam) never had either bug; this per-sample kernel is
   // what SetBatchSize(1) paths run - including online continual learning on client machines,
   // where OpenCL is the only tier. The per-lane guard keeps the tail group off the next row.
   for(int n = 0; n < 4; n++)
     {
      int col = j * 4 + n;
      if(col > inputs)
         break;
      int wi = i * (inputs + 1) + col;
      float inp = (col < inputs ? matrix_i[col] : 1.0f);
      float g = matrix_g[i] * inp;
      float mt = b1 * matrix_m[wi] + (1 - b1) * g;
      // v is STORED already square-rooted (a standard deviation) so it can be the denominator
      // directly; it must therefore be squared back before re-entering the recursion. Feeding the
      // stored sqrt in as if it were the variance - which every Adam kernel here did until
      // 2026-08-09 - mixes a std dev with a variance, and the resulting recursion has a fixed
      // point at v ~= b2 for ANY |g| below 1. The denominator then stops tracking the gradient
      // scale, delta becomes proportional to |g|, and Adam silently degrades to plain SGD. That is
      // not academic: it is why the conv/LSTM front end (which sits behind a batch-norm and so
      // receives gradients divided by sqrt(var) ~ 500) froze at dW/W 0.00% while the dense stack
      // trained. Measured with DirectML\batch_accum_check.cpp against the shipped kernel: 3285x
      // less displacement at |g|=1e-5 than at |g|=1, where a scale-invariant optimizer gives the
      // same distance for both. See NeuronBatchNorm.mqh, which got this right first.
      float vt = sqrt(b2 * matrix_v[wi] * matrix_v[wi] + (1 - b2) * g * g);
      float w = matrix_w[wi];
      float delta = clamp(l * mt / (vt > 0 ? vt : l * 10) - l * WEIGHT_DECAY * w, -MAX_WEIGHT_DELTA, MAX_WEIGHT_DELTA);
      // See UpdateWeightsMomentum's comment just above for why the matrix_w write is clamped to
      // MAX_WEIGHT - this dense Adam kernel is the one that updates the classification output
      // layer. No sign-agreement gate on the matrix_w write (removed 2026-07, all four backends in
      // sync) - see AI\NeuronCPU.mqh's CNeuron::updateInputWeights comment for why: it rectified
      // the one-hot softmax-CCE gradient stream into a permanent downward ratchet (the all-Neutral
      // collapse).
      matrix_w[wi] = clamp(w + delta, -MAX_WEIGHT, MAX_WEIGHT);
      matrix_m[wi] = mt;
      matrix_v[wi] = vt;
     }
  };
//+------------------------------------------------------------------+
//| Max-pooling layer (CNeuronPoolOCL) - no weights, just a sliding   |
//| max over the previous layer's output. Ported from the NeuroNet_DNG|
//| reference's CNeuronProofOCL kernels (these two had no bug).       |
//+------------------------------------------------------------------+
__kernel void FeedForwardProof(__global float *matrix_i,
                               __global float *matrix_o,
                               int inputs, int window, int step)
  {
   int i = get_global_id(0);
   int pos = i * step;
   // Mirrors the bounds guard in this kernel's pure-MQL5 host counterpart
   // (CNeuronPoolOCL::feedForwardCPU in AI\NeuronOCLConvPool.mqh) - outputs is sized so pos shouldn't
   // reach inputs in practice, but without this guard an out-of-bounds pos would read past matrix_i.
   if(pos >= inputs)
      return;
   float result = matrix_i[pos];
   for(int k = 1; k < window; k++)
     {
      int shift = k + pos;
      if(shift >= inputs)
         break;
      result = max(result, matrix_i[shift]);
     }
   matrix_o[i] = result;
  }
//+------------------------------------------------------------------+
//|                                                                  |
//+------------------------------------------------------------------+
__kernel void CalcInputGradientProof(__global float *matrix_i,
                                     __global float *matrix_g,
                                     __global float *matrix_o,
                                     __global float *matrix_ig,
                                     int outputs, int window, int step)
  {
   int i = get_global_id(0);
   float prev_gradient = 0.0f;
   float value = matrix_i[i];
   int start = i - window + step;
   start = (start - start % step) / step;
   int stop = (i - i % step) / step + 1;
   for(int out = max(0, start); out < min(outputs, stop); out++)
      if(value == matrix_o[out])
         prev_gradient += matrix_g[out];
   matrix_ig[i] = prev_gradient;
  }
//+------------------------------------------------------------------+
//| Convolution layer (CNeuronConvOCL) - ported from the NeuroNet_DNG |
//| reference library and adapted to this project's buffer convention.|
//| window_in+1 weights (incl. bias) are shared across all sliding    |
//| positions; window_out lets one position produce several output    |
//| channels (filters) at once.                                       |
//+------------------------------------------------------------------+
__kernel void FeedForwardConv(__global float *matrix_w,
                              __global float *matrix_i,
                              __global float *matrix_o,
                              int inputs, int step, int window_in, int window_out,
                              int activation)
  {
   int i = get_global_id(0);
   int shift_out = window_out * i;
   int shift_in = step * i;
   for(int out = 0; out < window_out; out++)
     {
      int shift = (window_in + 1) * out;
      int stop = (window_in <= (inputs - shift_in) ? window_in : (inputs - shift_in));
      float sum = 0.0f;
      for(int k = 0; k < stop; k++)
         sum += matrix_i[shift_in + k] * matrix_w[shift + k];
      sum += matrix_w[shift + window_in]; // bias
      switch(activation)
        {
         case 0:
            sum = tanh(sum);
            break;
         case 1:
            sum = 1 / (1 + exp(-clamp(sum, -50.0f, 50.0f)));
            break;
         case 2: // PReLU, param=0.01
            if(sum < 0)
               sum *= 0.01f;
            break;
        }
      matrix_o[out + shift_out] = sum;
     }
  }
//+------------------------------------------------------------------+
//|                                                                  |
//+------------------------------------------------------------------+
__kernel void CalcHiddenGradientConv(__global float *matrix_w,
                                     __global float *matrix_g,
                                     __global float *matrix_o,
                                     __global float *matrix_ig,
                                     int outputs, int step, int window_in, int window_out,
                                     int activation)
  {
   int i = get_global_id(0);
   float out = matrix_o[i];
   int start = i - window_in + step;
   start = max((start - start % step) / step, 0);
   int stop = (i - i % step) / step + 1;
   if(stop > (outputs / window_out))
      stop = outputs / window_out;
   float sum = 0.0f;
   for(int h = 0; h < window_out; h++)
     {
      for(int k = start; k < stop; k++)
        {
         int shift_w = (stop - k - 1) * step + i % step + h * (window_in + 1);
         int shift_g = k * window_out + h;
         if(shift_g >= outputs || shift_w >= (window_in + 1) * window_out)
            break;
         sum += matrix_g[shift_g] * matrix_w[shift_w];
        }
     }
   switch(activation)
     {
      case 0:
         sum = clamp(sum + out, -1.0f, 1.0f) - out;
         sum = sum * fmax(MIN_ACTIVATION_DERIVATIVE, 1 - pow(out, 2.0f));
         break;
      case 1:
         sum = clamp(sum + out, 0.0f, 1.0f) - out;
         sum = sum * fmax(MIN_ACTIVATION_DERIVATIVE, out * (1 - out));
         break;
      case 2: // PReLU, param=0.01
         if(out < 0)
            sum *= 0.01f;
         break;
     }
   matrix_ig[i] = sum;
  }
//+------------------------------------------------------------------+
//|                                                                  |
//+------------------------------------------------------------------+
__kernel void UpdateWeightsConvMomentum(__global float *matrix_w,
                                        __global float *matrix_g,
                                        __global float *matrix_i,
                                        __global float *matrix_dw,
                                        int inputs, float learning_rates, float momentum,
                                        int window_in, int window_out, int step, int optimizer)
  {
   int i = get_global_id(0);
   int shift = i % (window_in + 1);
   int shift_out = (i - shift) / (window_in + 1);
   int total = (inputs - window_in) % step;
   total = (inputs - window_in - total) / step + (total > 0 ? 1 : 0);
   float grad = 0.0f;
   for(int t = 0; t < total; t++)
     {
      if(shift != window_in && (shift + t * step) >= inputs)
         break;
      grad += matrix_g[t * window_out + shift_out] * (shift == window_in ? 1 : matrix_i[shift + t * step]);
     }
   // `optimizer` retained unused - see UpdateWeightsMomentum's comment above.
   float delta = learning_rates * grad + momentum * matrix_dw[i];
   matrix_dw[i] = delta;
   matrix_w[i] = clamp(matrix_w[i] + delta, -MAX_WEIGHT, MAX_WEIGHT);
  };
//+------------------------------------------------------------------+
//|                                                                  |
//+------------------------------------------------------------------+
__kernel void UpdateWeightsConvAdam(__global float *matrix_w,
                                    __global const float *matrix_g,
                                    __global const float *matrix_i,
                                    __global float *matrix_m,
                                    __global float *matrix_v,
                                    const int inputs, const float l, const float b1, const float b2,
                                    int window_in, int window_out, int step)
  {
   int i = get_global_id(0); // weight offset within a single output channel's (window_in+1) block
   if(i > window_in)
      return;
   int total = (inputs - (window_in - step)) % step;
   total = (inputs - (window_in - step) - total) / step + (total > 0 ? 1 : 0);
   for(int out = 0; out < window_out; out++)
     {
      int shift_w = i + out * (window_in + 1);
      float grad = 0.0f;
      for(int t = 0; t < total; t++)
        {
         if(i != window_in && (i + t * step) >= inputs)
            break;
         grad += matrix_g[t * window_out + out] * (i == window_in ? 1 : matrix_i[i + t * step]);
        }
      float mt = b1 * matrix_m[shift_w] + (1 - b1) * grad;
      // v squared back before the recursion - see UpdateWeightsAdam's comment above for why.
      float vt = sqrt(b2 * matrix_v[shift_w] * matrix_v[shift_w] + (1 - b2) * pow(grad, 2.0f));
      float delta = clamp(l * mt / (vt > 0 ? vt : l * 10) - l * WEIGHT_DECAY * matrix_w[shift_w], -MAX_WEIGHT_DELTA, MAX_WEIGHT_DELTA);
      // Sign-agreement gate removed - see UpdateWeightsAdam's comment above for why.
      matrix_w[shift_w] = clamp(matrix_w[shift_w] + delta, -MAX_WEIGHT, MAX_WEIGHT);
      matrix_m[shift_w] = mt;
      matrix_v[shift_w] = vt;
     }
  };
//+------------------------------------------------------------------+
//| MINI-BATCH GRADIENT ACCUMULATION (2026-08-09 audit, F4).          |
//|                                                                   |
//| Training here is otherwise pure online SGD - one weight update    |
//| per sample - which makes the per-update gradient noise maximal    |
//| and the end-of-era weight state a high-variance random draw. The  |
//| fix is to sum each weight's gradient over TRAIN_BATCH_SIZE        |
//| samples and take ONE optimizer step on the mean.                  |
//|                                                                   |
//| Only the ACCUMULATION lives here. The optimizer step itself is    |
//| elementwise O(weights), so - exactly like the batch-norm layer    |
//| and the softmax+CCE output gradient before it - it is computed    |
//| once per BATCH in MQL5 against the host mirrors of these buffers  |
//| (CNeuronBaseOCL::ApplyAccumulatedGradients). That keeps ONE       |
//| optimizer implementation for all four compute tiers instead of    |
//| four that can drift apart, which this engine has been bitten by   |
//| repeatedly.                                                       |
//|                                                                   |
//| The per-weight gradient expressions below are copied verbatim     |
//| from UpdateWeightsAdam / UpdateWeightsConvAdam above - they MUST  |
//| stay identical, because batch size 1 has to reproduce the         |
//| unbatched path exactly.                                           |
//+------------------------------------------------------------------+
__kernel void AccumulateWeightGrad(__global float *matrix_acc,
                                   __global const float *matrix_g,
                                   __global const float *matrix_i,
                                   const int inputs)
  {
   const int i = get_global_id(0);          // destination neuron
   const int j = get_global_id(1);          // source slot, j == inputs is the bias
   const int wi = i * (inputs + 1) + j;
   // Deliberately scalar, not the float4 form UpdateWeightsAdam uses: this kernel replaces that
   // dispatch on batched runs rather than adding to it, so the dispatch count is unchanged, and a
   // plain expression is far easier to hold against the host-side apply it has to agree with.
   matrix_acc[wi] += matrix_g[i] * (j < inputs ? matrix_i[j] : 1.0f);
  };
//+------------------------------------------------------------------+
//| Convolution counterpart. Flat over (window_in+1)*window_out, i.e. |
//| UpdateWeightsConvMomentum's indexing, but with the POSITION COUNT |
//| taken from UpdateWeightsConvAdam.                                 |
//|                                                                   |
//| Those two disagree, and it is worth being explicit about which is |
//| reproduced here. Positions of a window_in-wide window stepping by |
//| `step` over `inputs` values is floor((inputs-window_in)/step)+1.  |
//| The Adam kernel's total - ceil((inputs-window_in+step)/step) -    |
//| equals that; the Momentum kernel's ceil((inputs-window_in)/step)  |
//| is one SHORT whenever step divides the span, so the SGD conv path |
//| silently drops the last sliding position's gradient. That is a    |
//| pre-existing defect in the SGD path (ADAM is the default and the  |
//| only optimizer these models ship with) and is NOT fixed here -    |
//| fixing it would change SGD's unbatched behaviour, which is not    |
//| this change's business. The correct count is used because this    |
//| accumulator feeds the shared host-side apply for BOTH optimizers. |
//+------------------------------------------------------------------+
__kernel void AccumulateWeightGradConv(__global float *matrix_acc,
                                       __global const float *matrix_g,
                                       __global const float *matrix_i,
                                       const int inputs, int window_in, int window_out, int step)
  {
   int i = get_global_id(0);
   int shift = i % (window_in + 1);
   int shift_out = (i - shift) / (window_in + 1);
   int total = (inputs - (window_in - step)) % step;
   total = (inputs - (window_in - step) - total) / step + (total > 0 ? 1 : 0);
   float grad = 0.0f;
   for(int t = 0; t < total; t++)
     {
      if(shift != window_in && (shift + t * step) >= inputs)
         break;
      grad += matrix_g[t * window_out + shift_out] * (shift == window_in ? 1 : matrix_i[shift + t * step]);
     }
   matrix_acc[i] += grad;
  };
//+------------------------------------------------------------------+
//| dst += src, elementwise. Used by the LSTM's mini-batch path: that |
//| layer's WeightsGradient already holds one sample's complete dW    |
//| (summed over the BPTT timesteps), but the seq-backward kernels    |
//| MEMSET it at the top of every call, so the running batch total    |
//| has to live in a second buffer. Deliberately generic rather than  |
//| LSTM-specific - it is just a vector add.                          |
//+------------------------------------------------------------------+
__kernel void AccumulateBufferInto(__global float *dst,
                                   __global const float *src)
  {
   const int i = get_global_id(0);
   dst[i] += src[i];
  };
//+------------------------------------------------------------------+
//| MINI-BATCH APPLY - one optimizer step on the accumulated mean.    |
//|                                                                   |
//| These replace a HOST-SIDE apply that read the whole weight matrix,|
//| the accumulator and both Adam moments back over the bus, stepped  |
//| them in MQL5, and wrote four buffers out again: eight full-matrix |
//| transfers per batch per weight block, each a blocking sync. At    |
//| batch 8 that was ~1 weight-matrix-worth of PCIe traffic PER       |
//| SAMPLE, which is why an RX 580 lost to a CPU thread pool on this  |
//| workload. Market builds cannot import a DLL, so OpenCL is the     |
//| tier clients actually run and this path is the product's speed.   |
//|                                                                   |
//| Flat and elementwise: the accumulator already holds the summed    |
//| per-weight gradient (AccumulateWeightGrad above), so there is no  |
//| outer product to reconstruct and no float4 dispatch to mirror.    |
//|                                                                   |
//| `scale` is 1/batchCount - the MEAN. The arithmetic is a           |
//| line-for-line transcription of UpdateWeightsAdam above, INCLUDING |
//| the squared-back second moment, both clamps and the decoupled     |
//| decay, so batch size 1 still reproduces the unbatched kernel. The |
//| host implementation in NeuronOCLBase.mqh remains the reference    |
//| for the DLL and pure-MQL5 tiers and must be edited in step.       |
//|                                                                   |
//| The accumulator is zeroed HERE, saving a separate clear dispatch  |
//| and guaranteeing it cannot be left dirty for the next batch by an |
//| early return between the step and the clear.                      |
//+------------------------------------------------------------------+
__kernel void ApplyAccumAdam(__global float *matrix_w,
                             __global float *matrix_acc,
                             __global float *matrix_m,
                             __global float *matrix_v,
                             const float scale, const float l,
                             const float b1, const float b2)
  {
   const int i = get_global_id(0);
   const float grad = matrix_acc[i] * scale;
   const float wv = matrix_w[i];
   const float mt = b1 * matrix_m[i] + (1 - b1) * grad;
   // Squared back before re-entering the recursion - see UpdateWeightsAdam for the full account of
   // what feeding the stored square root in as the variance did to this engine.
   const float vt = sqrt(b2 * matrix_v[i] * matrix_v[i] + (1 - b2) * grad * grad);
   const float delta = clamp(l * mt / (vt > 0 ? vt : l * 10) - l * WEIGHT_DECAY * wv,
                             -MAX_WEIGHT_DELTA, MAX_WEIGHT_DELTA);
   matrix_w[i] = clamp(wv + delta, -MAX_WEIGHT, MAX_WEIGHT);
   matrix_m[i] = mt;
   matrix_v[i] = vt;
   matrix_acc[i] = 0.0f;
  };
//+------------------------------------------------------------------+
//| SGD/momentum counterpart. Note the delta is NOT clamped to        |
//| MAX_WEIGHT_DELTA here - matching UpdateWeightsMomentum and the    |
//| host apply, both of which clamp only the resulting weight. Kept   |
//| that way deliberately rather than "improved": the batched path    |
//| has to reproduce the unbatched one, not diverge from it.          |
//+------------------------------------------------------------------+
__kernel void ApplyAccumMomentum(__global float *matrix_w,
                                 __global float *matrix_acc,
                                 __global float *matrix_dw,
                                 const float scale, const float lr, const float momentum)
  {
   const int i = get_global_id(0);
   const float grad = matrix_acc[i] * scale;
   const float delta = lr * grad + momentum * matrix_dw[i];
   matrix_dw[i] = delta;
   matrix_w[i] = clamp(matrix_w[i] + delta, -MAX_WEIGHT, MAX_WEIGHT);
   matrix_acc[i] = 0.0f;
  };
//+------------------------------------------------------------------+
//| BATCH NORM (CNeuronBatchNormOCL), device-side since 2026-08-09.   |
//|                                                                   |
//| This layer computed host-side on every tier, which on OpenCL      |
//| meant a blocking read of the previous layer's output and a        |
//| blocking write of its own - per sample, per layer, forward AND    |
//| backward. These kernels are LINE-FOR-LINE transcriptions of the   |
//| host implementations in AI\NeuronBatchNorm.mqh (NormalizeHost,    |
//| calcInputGradients, StepGammaBeta), including every NaN guard and |
//| clamp - that layer latched NaN into its persisted statistics once |
//| (2026-08-02) and the guards are what keep a single garbage value  |
//| from poisoning the .nnw forever. The host copies remain the       |
//| reference AND the runtime for the DLL / pure-MQL5 tiers; edit the |
//| two together. The MQL5 side also self-checks each kernel against  |
//| its host twin on first use and latches back to the host path on   |
//| any disagreement, so a transcription error costs a warning and    |
//| some speed, never a corrupted model.                              |
//|                                                                   |
//| Slot layout below MUST match NeuronBatchNorm.mqh's BN_OPT_*.      |
//+------------------------------------------------------------------+
#define BN_OPT_STRIDE   9
#define BN_OPT_MEAN     0
#define BN_OPT_VAR      1
#define BN_OPT_NX       2
#define BN_OPT_GAMMA    3
#define BN_OPT_BETA     4
#define BN_OPT_MG       5
#define BN_OPT_MB       6
#define BN_OPT_VG       7
#define BN_OPT_VB       8
#define BN_EPSILON      1.0e-10f
#define BN_MIN_STD      1.0e-4f
#define BN_MAX_INPUT    1.0e6f
//+------------------------------------------------------------------+
//| Forward transform for one unit. `w` is the host-computed          |
//| effective EMA window (the 1,2,3..iBatchSize ramp lives host-side  |
//| with iSamplesSeen); `frozen` mirrors bStatsFrozen - statistics    |
//| are USED but not UPDATED, classic inference semantics.            |
//+------------------------------------------------------------------+
__kernel void BatchNormForward(__global const float *matrix_i,
                               __global float *matrix_o,
                               __global float *options,
                               const float w, const int frozen)
  {
   const int i = get_global_id(0);
   const int shift = i * BN_OPT_STRIDE;
   float x = matrix_i[i];
   // Guards in the same order as NormalizeHost: non-finite -> 0, then the range clamp that keeps the
   // statistics arithmetic bounded (a merely-huge finite x overflows (x-mean)^2 long before infinity).
   if(isnan(x) || isinf(x))
      x = 0.0f;
   x = clamp(x, -BN_MAX_INPUT, BN_MAX_INPUT);
   float mean = options[shift + BN_OPT_MEAN];
   float variance = options[shift + BN_OPT_VAR];
   // Self-heal statistics poisoned before the guards existed - Load faithfully restores NaN.
   if(isnan(mean) || isinf(mean))
      mean = x;
   if(isnan(variance) || isinf(variance) || variance < 0.0f)
      variance = 0.0f;
   if(frozen == 0)
     {
      // EMA update, new mean used by the variance term - exact order of the host code.
      mean = (mean * (w - 1.0f) + x) / w;
      variance = (variance * (w - 1.0f) + (x - mean) * (x - mean)) / w;
     }
   const float delta = x - mean;
   const float sd = fmax(sqrt(variance + BN_EPSILON), BN_MIN_STD);
   const float nx = delta / sd;
   float gamma = options[shift + BN_OPT_GAMMA];
   float beta = options[shift + BN_OPT_BETA];
   // Local substitution only, exactly like the host forward: the PERSISTED repair happens in the
   // gamma/beta step, which is the copy that makes it stick.
   if(isnan(gamma) || isinf(gamma))
      gamma = 1.0f;
   if(isnan(beta) || isinf(beta))
      beta = 0.0f;
   if(frozen == 0)
     {
      options[shift + BN_OPT_MEAN] = mean;
      options[shift + BN_OPT_VAR] = variance;
     }
   // nx cached even when frozen - costs nothing, keeps the buffer consistent with the output.
   options[shift + BN_OPT_NX] = nx;
   matrix_o[i] = gamma * nx + beta;
  };
//+------------------------------------------------------------------+
//| Backward for one unit: dL/dx = gamma/sd * dL/dy, then the         |
//| PREVIOUS layer's activation derivative applied with the same      |
//| clamp-to-range "implied target" treatment CaclHiddenGradient      |
//| uses, so from below a batch-norm layer is indistinguishable from  |
//| any other. `activation` is the NATIVE code (0=TANH 1=SIGMOID      |
//| 2=PRELU, anything else passes through) - NativeActivationCode()   |
//| at the dispatch site, never a raw enum cast.                      |
//+------------------------------------------------------------------+
__kernel void BatchNormHiddenGrad(__global const float *matrix_g,
                                  __global const float *prev_o,
                                  __global float *prev_g,
                                  __global const float *options,
                                  const int activation)
  {
   const int i = get_global_id(0);
   const int shift = i * BN_OPT_STRIDE;
   const float sd = fmax(sqrt(options[shift + BN_OPT_VAR] + BN_EPSILON), BN_MIN_STD);
   float g = matrix_g[i] * options[shift + BN_OPT_GAMMA] / sd;
   const float out = prev_o[i];
   switch(activation)
     {
      case 0:
         g = clamp(g + out, -1.0f, 1.0f) - out;
         g = g * fmax(MIN_ACTIVATION_DERIVATIVE, 1.0f - out * out);
         break;
      case 1:
         g = clamp(g + out, 0.0f, 1.0f) - out;
         g = g * fmax(MIN_ACTIVATION_DERIVATIVE, out * (1.0f - out));
         break;
      case 2:
         g = g * (out >= 0.0f ? 1.0f : 0.01f);
         break;
      default:
         break;
     }
   prev_g[i] = g;
  };
//+------------------------------------------------------------------+
//| Mini-batch accumulate for gamma/beta: dL/dgamma = dL/dy * nx,     |
//| dL/dbeta = dL/dy, summed per unit into acc[2i]/acc[2i+1]. Formed  |
//| per sample because BN_OPT_NX is overwritten by every forward.     |
//+------------------------------------------------------------------+
__kernel void BatchNormAccumGammaBeta(__global const float *matrix_g,
                                      __global const float *options,
                                      __global float *acc)
  {
   const int i = get_global_id(0);
   const float g = matrix_g[i];
   acc[2 * i]     += g * options[i * BN_OPT_STRIDE + BN_OPT_NX];
   acc[2 * i + 1] += g;
  };
//+------------------------------------------------------------------+
//| One gamma/beta optimizer step on the accumulated mean, then the   |
//| accumulator is zeroed. Transcribed from StepGammaBeta including   |
//| its exact write ordering: moments are stored BEFORE the           |
//| finiteness check, so a skipped unit still advances its optimizer  |
//| state, and the persisted gamma/beta self-heal happens HERE (this  |
//| is the copy that sticks). NO WEIGHT DECAY - decaying gamma toward |
//| zero is the collapse this layer exists to prevent. The per-sample |
//| path is this same pair of kernels at scale=1.                     |
//| `optimizer`: 0 = SGD/momentum (lr=eta, momentum=alpha), else Adam |
//| (lt = bias-corrected step, b1, b2; v stored square-rooted, same   |
//| convention as every Adam kernel in this file).                    |
//+------------------------------------------------------------------+
__kernel void BatchNormApplyGammaBeta(__global float *options,
                                      __global float *acc,
                                      const float scale, const float lt,
                                      const float b1, const float b2,
                                      const float lr, const float momentum,
                                      const int optimizer)
  {
   const int i = get_global_id(0);
   const int shift = i * BN_OPT_STRIDE;
   const float gGamma = acc[2 * i] * scale;
   const float gBeta = acc[2 * i + 1] * scale;
   acc[2 * i] = 0.0f;
   acc[2 * i + 1] = 0.0f;
   float gamma = options[shift + BN_OPT_GAMMA];
   float beta = options[shift + BN_OPT_BETA];
   if(isnan(gamma) || isinf(gamma))
      gamma = 1.0f;
   if(isnan(beta) || isinf(beta))
      beta = 0.0f;
   float dGamma, dBeta;
   if(optimizer == 0)
     {
      dGamma = lr * gGamma + momentum * options[shift + BN_OPT_MG];
      dBeta  = lr * gBeta  + momentum * options[shift + BN_OPT_MB];
      options[shift + BN_OPT_MG] = dGamma;
      options[shift + BN_OPT_MB] = dBeta;
     }
   else
     {
      const float mg = b1 * options[shift + BN_OPT_MG] + (1.0f - b1) * gGamma;
      const float mb = b1 * options[shift + BN_OPT_MB] + (1.0f - b1) * gBeta;
      const float vgOld = options[shift + BN_OPT_VG];
      const float vbOld = options[shift + BN_OPT_VB];
      const float vg = sqrt(b2 * vgOld * vgOld + (1.0f - b2) * gGamma * gGamma);
      const float vb = sqrt(b2 * vbOld * vbOld + (1.0f - b2) * gBeta * gBeta);
      dGamma = lt * mg / (vg > 0.0f ? vg : lt * 10.0f);
      dBeta  = lt * mb / (vb > 0.0f ? vb : lt * 10.0f);
      options[shift + BN_OPT_MG] = mg;
      options[shift + BN_OPT_MB] = mb;
      options[shift + BN_OPT_VG] = vg;
      options[shift + BN_OPT_VB] = vb;
     }
   dGamma = clamp(dGamma, -MAX_WEIGHT_DELTA, MAX_WEIGHT_DELTA);
   dBeta  = clamp(dBeta, -MAX_WEIGHT_DELTA, MAX_WEIGHT_DELTA);
   // A non-finite delta means SKIP this unit's parameter write, never poison it - the moments above
   // are already stored, matching the host's ordering exactly.
   if(isnan(dGamma) || isinf(dGamma) || isnan(dBeta) || isinf(dBeta))
      return;
   options[shift + BN_OPT_GAMMA] = clamp(gamma + dGamma, -MAX_WEIGHT, MAX_WEIGHT);
   options[shift + BN_OPT_BETA]  = clamp(beta + dBeta, -MAX_WEIGHT, MAX_WEIGHT);
  };
//+------------------------------------------------------------------+
//| LSTM (CNeuronLSTMOCL) - derived from scratch from the standard    |
//| LSTM equations (NOT ported from any reference implementation).    |
//| Single-timestep-truncated BPTT: gradient does not flow back into  |
//| h_prev/c_prev from a previous step - only this step's gates and   |
//| weights receive gradient. Adam-only optimizer.                    |
//| Weight layout: 4 gates (order: forget, input, output, candidate), |
//| each a block of H*(I+H+1) floats - per hidden unit: H weights on  |
//| h_prev, I weights on the current input, 1 bias.                   |
//+------------------------------------------------------------------+
__kernel void LSTM_Gates(__global float *matrix_w,
                         __global float *hidden_prev,
                         __global float *inputs,
                         __global float *concatenated,
                         int hidden_size, int input_size)
  {
   int id = get_global_id(0);
   int gate = get_global_id(1);
   int per_gate = hidden_size * (hidden_size + input_size + 1);
   int shift = gate * per_gate + id * (hidden_size + input_size + 1);
   float sum = 0.0f;
   for(int k = 0; k < hidden_size; k++)
      sum += hidden_prev[k] * matrix_w[shift + k];
   for(int k = 0; k < input_size; k++)
      sum += inputs[k] * matrix_w[shift + hidden_size + k];
   sum += matrix_w[shift + hidden_size + input_size];
   float val = (gate < 3) ? (1.0f / (1.0f + exp(-sum))) : tanh(sum);
   concatenated[gate * hidden_size + id] = val;
  }
//+------------------------------------------------------------------+
//|                                                                  |
//+------------------------------------------------------------------+
__kernel void LSTM_State(__global float *concatenated,
                         __global float *memory,
                         __global float *hidden_prev,
                         __global float *hidden_cache,
                         __global float *output,
                         int hidden_size)
  {
   int id = get_global_id(0);
   float f = concatenated[id];
   float ii = concatenated[hidden_size + id];
   float o = concatenated[2 * hidden_size + id];
   float g = concatenated[3 * hidden_size + id];
   float c_prev = memory[id];
   memory[hidden_size + id] = c_prev;
   float c_t = f * c_prev + ii * g;
   memory[id] = c_t;
   hidden_cache[id] = hidden_prev[id];
   output[id] = o * tanh(c_t);
  }
//+------------------------------------------------------------------+
//|                                                                  |
//+------------------------------------------------------------------+
__kernel void LSTM_GateGradient(__global float *gradient,
                                __global float *memory,
                                __global float *concatenated,
                                __global float *concatenated_gradient,
                                int hidden_size)
  {
   int id = get_global_id(0);
   float c_t = memory[id];
   float c_prev = memory[hidden_size + id];
   float f = concatenated[id];
   float ii = concatenated[hidden_size + id];
   float o = concatenated[2 * hidden_size + id];
   float g = concatenated[3 * hidden_size + id];
   float t = tanh(c_t);
   float dh = gradient[id];
   // Derivative floor on every gate below - see MIN_ACTIVATION_DERIVATIVE's declaration comment.
   // LSTM gates are sigmoid/tanh by design and routinely sit near their saturated extremes (a fully
   // "open" or "closed" gate), which is exactly where an unfloored derivative goes to zero and
   // blocks gradient from flowing through the recurrent path at all.
   float dc = dh * o * fmax(MIN_ACTIVATION_DERIVATIVE, 1 - t * t);
   concatenated_gradient[2 * hidden_size + id] = dh * t * fmax(MIN_ACTIVATION_DERIVATIVE, o * (1 - o));       // output gate
   concatenated_gradient[id] = dc * c_prev * fmax(MIN_ACTIVATION_DERIVATIVE, f * (1 - f));                    // forget gate
   concatenated_gradient[hidden_size + id] = dc * g * fmax(MIN_ACTIVATION_DERIVATIVE, ii * (1 - ii));         // input gate
   concatenated_gradient[3 * hidden_size + id] = dc * ii * fmax(MIN_ACTIVATION_DERIVATIVE, 1 - g * g);        // candidate
  }
//+------------------------------------------------------------------+
//|                                                                  |
//+------------------------------------------------------------------+
__kernel void LSTM_WeightsGradient(__global float *concatenated_gradient,
                                   __global float *hidden_cache,
                                   __global float *inputs,
                                   __global float *weights_gradient,
                                   int hidden_size, int input_size)
  {
   int wi = get_global_id(0);
   int per_gate = hidden_size * (hidden_size + input_size + 1);
   int gate = wi / per_gate;
   int rem = wi % per_gate;
   int id = rem / (hidden_size + input_size + 1);
   int k = rem % (hidden_size + input_size + 1);
   float inp = (k < hidden_size ? hidden_cache[k] : (k < hidden_size + input_size ? inputs[k - hidden_size] : 1.0f));
   weights_gradient[wi] = concatenated_gradient[gate * hidden_size + id] * inp;
  }
//+------------------------------------------------------------------+
//|                                                                  |
//+------------------------------------------------------------------+
__kernel void LSTM_InputsGradient(__global float *concatenated_gradient,
                                  __global float *matrix_w,
                                  __global float *inputs_gradient,
                                  int hidden_size, int input_size)
  {
   int j = get_global_id(0);
   int per_gate = hidden_size * (hidden_size + input_size + 1);
   float sum = 0.0f;
   for(int gate = 0; gate < 4; gate++)
      for(int id = 0; id < hidden_size; id++)
         sum += concatenated_gradient[gate * hidden_size + id] * matrix_w[gate * per_gate + id * (hidden_size + input_size + 1) + hidden_size + j];
   inputs_gradient[j] = sum;
  }
//+------------------------------------------------------------------+
//|                                                                  |
//+------------------------------------------------------------------+
//+------------------------------------------------------------------+
//| Sequence LSTM - ONE LAUNCH PER TIMESTEP.                          |
//|                                                                   |
//| The recurrence is sequential: step t needs h/c from step t-1, and |
//| every hidden unit's gates read ALL of h_{t-1}. OpenCL can only    |
//| barrier WITHIN a work-group, and nothing here controls how the    |
//| runtime partitions the global size, so a single-launch loop with  |
//| barrier() would be correct only by luck. Driving the loop from    |
//| the host instead makes each launch an implicit global barrier -   |
//| slower per step, correct on every device.                         |
//|                                                                   |
//| Layouts match CPU_LSTMSeqForward in DirectML\WarriorCPU.cpp (the  |
//| authoritative contract): gate order [f,i,o,g]; weight row `hid`   |
//| of a gate is [H recurrent | Iw input | 1 bias]; caches are        |
//| timestep-major. h_{-1} and c_{-1} are zero.                       |
//+------------------------------------------------------------------+
__kernel void LSTM_SeqStepForward(__global float *matrix_w,
                                  __global float *inputs,
                                  __global float *cache_gates,
                                  __global float *cache_cell,
                                  __global float *cache_hidden,
                                  __global float *output,
                                  int hidden_size,
                                  int step_inputs,
                                  int steps,
                                  int t)
  {
   int id = get_global_id(0);
   int H = hidden_size;
   int Iw = step_inputs;
   int row = H + Iw + 1;
   int per_gate = H * row;
   __global float *x = inputs + t * Iw;
   float g[4];
   for(int gate = 0; gate < 4; gate++)
     {
      int shift = gate * per_gate + id * row;
      float sum = 0;
      if(t > 0)
        {
         __global float *h_prev = cache_hidden + (t - 1) * H;
         for(int k = 0; k < H; k++)
            sum += h_prev[k] * matrix_w[shift + k];
        }
      for(int k = 0; k < Iw; k++)
         sum += x[k] * matrix_w[shift + H + k];
      sum += matrix_w[shift + H + Iw];
      g[gate] = (gate < 3) ? (1 / (1 + exp(-sum))) : tanh(sum);
      cache_gates[t * 4 * H + gate * H + id] = g[gate];
     }
   float c_prev = (t > 0) ? cache_cell[(t - 1) * H + id] : 0;
   float c_t = g[0] * c_prev + g[1] * g[3];
   cache_cell[t * H + id] = c_t;
   float h_t = g[2] * tanh(c_t);
   cache_hidden[t * H + id] = h_t;
   //--- the layer's visible output is the LAST hidden state
   if(t == steps - 1)
      output[id] = h_t;
  }
//+------------------------------------------------------------------+
//| BPTT step 1/3: gate gradients for step t, and dc for step t-1.    |
//| dh comes from the layer above at t == steps-1 and from the        |
//| recurrence (dh_buf, written by LSTM_SeqStepInputGrad) otherwise.  |
//| dc_buf must be zeroed by the host before the first step.          |
//+------------------------------------------------------------------+
__kernel void LSTM_SeqStepGateGrad(__global float *out_gradient,
                                   __global float *dh_buf,
                                   __global float *dc_buf,
                                   __global float *cache_gates,
                                   __global float *cache_cell,
                                   __global float *gate_grad,
                                   int hidden_size,
                                   int steps,
                                   int t)
  {
   int id = get_global_id(0);
   int H = hidden_size;
   __global float *gt = cache_gates + t * 4 * H;
   float f = gt[id];
   float ii = gt[H + id];
   float o = gt[2 * H + id];
   float g = gt[3 * H + id];
   float c_t = cache_cell[t * H + id];
   float c_prev = (t > 0) ? cache_cell[(t - 1) * H + id] : 0;
   float tc = tanh(c_t);
   float dh = (t == steps - 1) ? out_gradient[id] : dh_buf[id];
   //--- the term the per-step kernels could never express: dc arriving from step t+1
   float dcTot = dh * o * fmax(MIN_ACTIVATION_DERIVATIVE, 1 - tc * tc) + dc_buf[id];
   gate_grad[2 * H + id] = dh * tc * fmax(MIN_ACTIVATION_DERIVATIVE, o * (1 - o));
   gate_grad[id] = dcTot * c_prev * fmax(MIN_ACTIVATION_DERIVATIVE, f * (1 - f));
   gate_grad[H + id] = dcTot * g * fmax(MIN_ACTIVATION_DERIVATIVE, ii * (1 - ii));
   gate_grad[3 * H + id] = dcTot * ii * fmax(MIN_ACTIVATION_DERIVATIVE, 1 - g * g);
   dc_buf[id] = dcTot * f;
  }
//+------------------------------------------------------------------+
//| BPTT step 2/3: ACCUMULATE dW for step t. Accumulates rather than  |
//| assigns because the gate weights are shared by every timestep -   |
//| the host zeroes weights_gradient once before the backward loop.   |
//+------------------------------------------------------------------+
__kernel void LSTM_SeqStepWeightGrad(__global float *gate_grad,
                                     __global float *cache_hidden,
                                     __global float *inputs,
                                     __global float *weights_gradient,
                                     int hidden_size,
                                     int step_inputs,
                                     int t)
  {
   int wi = get_global_id(0);
   int H = hidden_size;
   int Iw = step_inputs;
   int row = H + Iw + 1;
   int per_gate = H * row;
   int gate = wi / per_gate;
   int rem = wi % per_gate;
   int hid = rem / row;
   int k = rem % row;
   float inp;
   if(k < H)
      inp = (t > 0) ? cache_hidden[(t - 1) * H + k] : 0;
   else
      if(k < H + Iw)
         inp = inputs[t * Iw + (k - H)];
      else
         inp = 1;
   weights_gradient[wi] += gate_grad[gate * H + hid] * inp;
  }
//+------------------------------------------------------------------+
//| BPTT step 3/3: dx_t (this step's slice of the input gradient) and |
//| dh_{t-1} (the recurrent path). One flat range over both so a      |
//| narrow per-step input still saturates the device.                 |
//+------------------------------------------------------------------+
__kernel void LSTM_SeqStepInputGrad(__global float *gate_grad,
                                    __global float *matrix_w,
                                    __global float *inputs_gradient,
                                    __global float *dh_buf,
                                    int hidden_size,
                                    int step_inputs,
                                    int t)
  {
   int j = get_global_id(0);
   int H = hidden_size;
   int Iw = step_inputs;
   int row = H + Iw + 1;
   int per_gate = H * row;
   int col = (j < Iw) ? (H + j) : (j - Iw);
   float sum = 0;
   for(int gate = 0; gate < 4; gate++)
      for(int hid = 0; hid < H; hid++)
         sum += gate_grad[gate * H + hid] * matrix_w[gate * per_gate + hid * row + col];
   if(j < Iw)
      inputs_gradient[t * Iw + j] = sum;
   else
      dh_buf[j - Iw] = sum;
  }
//+------------------------------------------------------------------+
//|                                                                  |
//+------------------------------------------------------------------+
__kernel void LSTM_UpdateWeightsAdam(__global float *matrix_w,
                                     __global float *weights_gradient,
                                     __global float *matrix_m,
                                     __global float *matrix_v,
                                     float l, float b1, float b2)
  {
   int wi = get_global_id(0);
   float g = weights_gradient[wi];
   float mt = b1 * matrix_m[wi] + (1 - b1) * g;
   // v squared back before the recursion - see UpdateWeightsAdam's comment above for why.
   float vt = sqrt(b2 * matrix_v[wi] * matrix_v[wi] + (1 - b2) * g * g);
   float delta = clamp(l * mt / (vt > 0 ? vt : l * 10) - l * WEIGHT_DECAY * matrix_w[wi], -MAX_WEIGHT_DELTA, MAX_WEIGHT_DELTA);
   // Sign-agreement gate removed - see UpdateWeightsAdam's comment above for why.
   matrix_w[wi] = clamp(matrix_w[wi] + delta, -MAX_WEIGHT, MAX_WEIGHT);
   matrix_m[wi] = mt;
   matrix_v[wi] = vt;
  };
//+------------------------------------------------------------------+
//| SGD+momentum counterpart to LSTM_UpdateWeightsAdam above - same   |
//| flat (already-elementwise) weights_gradient/matrix_dw layout, just|
//| the classic heavy-ball update instead of Adam's per-parameter step|
//+------------------------------------------------------------------+
__kernel void LSTM_UpdateWeightsMomentum(__global float *matrix_w,
                                         __global float *weights_gradient,
                                         __global float *matrix_dw,
                                         float learning_rates, float momentum, int optimizer)
  {
   int wi = get_global_id(0);
   // `optimizer` retained unused - see UpdateWeightsMomentum's comment above.
   float delta = learning_rates * weights_gradient[wi] + momentum * matrix_dw[wi];
   matrix_dw[wi] = delta;
   matrix_w[wi] = clamp(matrix_w[wi] + delta, -MAX_WEIGHT, MAX_WEIGHT);
  };