// 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).
#define MIN_ACTIVATION_DERIVATIVE 1.0e-4f
// 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);
   float sum = 0;
   float out = matrix_o[i];
   float4 grad, weight;
   int shift = (outputs + 1) * i;
   for(int k = 0; k < outputs; k += 4)
     {
      switch(outputs - k)
        {
         case 0:
            grad = (float4)(matrix_g[k], 0, 0, 0);
            weight = (float4)(matrix_w[shift + k], 0, 0, 0);
            break;
         case 1:
            grad = (float4)(matrix_g[k], matrix_g[k + 1], 0, 0);
            weight = (float4)(matrix_w[shift + k], matrix_w[shift + k + 1], 0, 0);
            break;
         case 2:
            grad = (float4)(matrix_g[k], matrix_g[k + 1], matrix_g[k + 2], 0);
            weight = (float4)(matrix_w[shift + k], matrix_w[shift + k + 1], matrix_w[shift + k + 2], 0);
            break;
         default:
            grad = (float4)(matrix_g[k], matrix_g[k + 1], matrix_g[k + 2], matrix_g[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(grad, weight);
     }
   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 i = get_global_id(0);
   int j = get_global_id(1);
   int wi = i * (inputs + 1) + j;
   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);
   const int wi = i * (inputs + 1) + j * 4;
   float4 m, v, weight, inp;
   switch(inputs - j * 4)
     {
      case 0:
         inp = (float4)(1, 0, 0, 0);
         weight = (float4)(matrix_w[wi], 0, 0, 0);
         m = (float4)(matrix_m[wi], 0, 0, 0);
         v = (float4)(matrix_v[wi], 0, 0, 0);
         break;
      case 1:
         inp = (float4)(matrix_i[j], 1, 0, 0);
         weight = (float4)(matrix_w[wi], matrix_w[wi + 1], 0, 0);
         m = (float4)(matrix_m[wi], matrix_m[wi + 1], 0, 0);
         v = (float4)(matrix_v[wi], matrix_v[wi + 1], 0, 0);
         break;
      case 2:
         inp = (float4)(matrix_i[j], matrix_i[j + 1], 1, 0);
         weight = (float4)(matrix_w[wi], matrix_w[wi + 1], matrix_w[wi + 2], 0);
         m = (float4)(matrix_m[wi], matrix_m[wi + 1], matrix_m[wi + 2], 0);
         v = (float4)(matrix_v[wi], matrix_v[wi + 1], matrix_v[wi + 2], 0);
         break;
      case 3:
         inp = (float4)(matrix_i[j], matrix_i[j + 1], matrix_i[j + 2], 1);
         weight = (float4)(matrix_w[wi], matrix_w[wi + 1], matrix_w[wi + 2], matrix_w[wi + 3]);
         m = (float4)(matrix_m[wi], matrix_m[wi + 1], matrix_m[wi + 2], matrix_m[wi + 3]);
         v = (float4)(matrix_v[wi], matrix_v[wi + 1], matrix_v[wi + 2], matrix_v[wi + 3]);
         break;
      default:
         inp = (float4)(matrix_i[j], matrix_i[j + 1], matrix_i[j + 2], matrix_i[j + 3]);
         weight = (float4)(matrix_w[wi], matrix_w[wi + 1], matrix_w[wi + 2], matrix_w[wi + 3]);
         m = (float4)(matrix_m[wi], matrix_m[wi + 1], matrix_m[wi + 2], matrix_m[wi + 3]);
         v = (float4)(matrix_v[wi], matrix_v[wi + 1], matrix_v[wi + 2], matrix_v[wi + 3]);
         break;
     }
   float4 g = matrix_g[i] * inp;
   float4 mt = b1 * m + (1 - b1) * g;
   float4 vt = sqrt(b2 * v + (1 - b2) * pow(g, 2.0f));
   float4 delta = clamp(l * mt / (vt > 0 ? vt : l * 10) - l * WEIGHT_DECAY * weight, -MAX_WEIGHT_DELTA, MAX_WEIGHT_DELTA);
   // See UpdateWeightsMomentum's comment just above for why every matrix_w write below 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 writes (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).
   switch(inputs - j * 4)
     {
      case 2:
         matrix_w[wi + 2] = clamp(matrix_w[wi + 2] + delta.s2, -MAX_WEIGHT, MAX_WEIGHT);
         matrix_m[wi + 2] = mt.s2;
         matrix_v[wi + 2] = vt.s2;
      case 1:
         matrix_w[wi + 1] = clamp(matrix_w[wi + 1] + delta.s1, -MAX_WEIGHT, MAX_WEIGHT);
         matrix_m[wi + 1] = mt.s1;
         matrix_v[wi + 1] = vt.s1;
      case 0:
         matrix_w[wi] = clamp(matrix_w[wi] + delta.s0, -MAX_WEIGHT, MAX_WEIGHT);
         matrix_m[wi] = mt.s0;
         matrix_v[wi] = vt.s0;
         break;
      default:
         matrix_w[wi] = clamp(matrix_w[wi] + delta.s0, -MAX_WEIGHT, MAX_WEIGHT);
         matrix_m[wi] = mt.s0;
         matrix_v[wi] = vt.s0;
         matrix_w[wi + 1] = clamp(matrix_w[wi + 1] + delta.s1, -MAX_WEIGHT, MAX_WEIGHT);
         matrix_m[wi + 1] = mt.s1;
         matrix_v[wi + 1] = vt.s1;
         matrix_w[wi + 2] = clamp(matrix_w[wi + 2] + delta.s2, -MAX_WEIGHT, MAX_WEIGHT);
         matrix_m[wi + 2] = mt.s2;
         matrix_v[wi + 2] = vt.s2;
         matrix_w[wi + 3] = clamp(matrix_w[wi + 3] + delta.s3, -MAX_WEIGHT, MAX_WEIGHT);
         matrix_m[wi + 3] = mt.s3;
         matrix_v[wi + 3] = vt.s3;
         break;
     }
  };
//+------------------------------------------------------------------+
//| 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;
   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 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]);
     }
   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;
      float vt = sqrt(b2 * 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;
     }
  };
//+------------------------------------------------------------------+
//| 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;
  }
//+------------------------------------------------------------------+
//|                                                                  |
//+------------------------------------------------------------------+
__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;
   float vt = sqrt(b2 * 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 wi = get_global_id(0);
   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);
  };