#ifndef NEURAL_NET_MQH
#define NEURAL_NET_MQH
#include "Statistics.mqh"

class CNeuralNet {
private:
   int      m_inputs;
   int      m_hidden;
   int      m_outputs;

   matrix   m_W1;
   vector   m_b1;
   matrix   m_W2;
   vector   m_b2;

   matrix   m_mW1, m_vW1;
   vector   m_mb1, m_vb1;
   matrix   m_mW2, m_vW2;
   vector   m_mb2, m_vb2;
    int      m_t;
    double   m_adamBeta1;
    double   m_adamBeta2;
    double   m_adamEps;
    double   m_beta1T;
    double   m_beta2T;

    double   m_l2;
   int      m_epochsTrained;
   double   m_lastLoss;
   bool     m_initialized;

   // HeInit per ReLU (NeuroBook §1.3)
   double HeScale(int fanIn) const {
      return MathSqrt(6.0 / MathMax(1, fanIn));
   }

   // Dropout (NeuroBook §6.2)
   double   m_dropoutRate;
   double   m_dropoutMask[];

   // BatchNorm semplificata su hidden (NeuroBook §6.3)
   // Running stats (EMA) + learnable affine
   vector   m_bnGamma;         // scala learnable
   vector   m_bnBeta;          // shift learnable
   vector   m_bnRunningMean;   // media mobile
   vector   m_bnRunningVar;    // varianza mobile
   double   m_bnMomentum;
   double   m_bnEps;

   // Debug: loss history
   double   m_lossHistory[];
   string   m_lossCsvFn;

   void ApplyDropout(vector &h) {
      for(int i = 0; i < h.Size(); i++) {
         m_dropoutMask[i] = ((double)MathRand() / 32767.0) > m_dropoutRate ? 1.0 : 0.0;
         h[i] *= m_dropoutMask[i];
      }
   }

   void ApplyBN(vector &h) {
      // Normalizza h IN PLACE: per-neuron stats (NeuroBook §6.3)
      for(int i = 0; i < h.Size(); i++) {
         double denom = MathSqrt(m_bnRunningVar[i] + m_bnEps);
         double standardized = (h[i] - m_bnRunningMean[i]) / denom;
         h[i] = m_bnGamma[i] * standardized + m_bnBeta[i];
      }
   }

   // Normalizza E AGGIORNA running stats (training, batch_size=1).
   // Usa running stats per normalizzare (stessa modalità di inference) per
   // evitare degenerazione con batch=1 (var=0). (NeuroBook §6.3, online BN)
   void TrainBN(const vector &h_in, vector &h_out) {
      h_out.Resize(h_in.Size());
      for(int i = 0; i < h_in.Size(); i++) {
         // Update running stats (per-neuron EMA)
         m_bnRunningMean[i] = m_bnMomentum * m_bnRunningMean[i] + (1.0 - m_bnMomentum) * h_in[i];
         double dx = h_in[i] - m_bnRunningMean[i];
         m_bnRunningVar[i]  = m_bnMomentum * m_bnRunningVar[i]  + (1.0 - m_bnMomentum) * dx * dx;
         // Normalize using running stats (same as ApplyBN)
         double denom = MathSqrt(m_bnRunningVar[i] + m_bnEps);
         double standardized = (h_in[i] - m_bnRunningMean[i]) / denom;
         h_out[i] = m_bnGamma[i] * standardized + m_bnBeta[i];
      }
   }

   // Backward attraverso BN: dh_out[i] = dh_norm[i] * gamma[i] / sqrt(runningVar[i] + eps)
   void BNBackward(vector &dh_out, const vector &h_in, const vector &dh_norm) {
      dh_out.Resize(h_in.Size());
      for(int i = 0; i < h_in.Size(); i++) {
         double denom = MathSqrt(m_bnRunningVar[i] + m_bnEps);
         dh_out[i] = dh_norm[i] * m_bnGamma[i] / denom;
      }
   }

   double ReLU(double x) const {
      return (x > 0.0) ? x : 0.0;
   }

   double ReLUDeriv(double x) const {
      return (x > 0.0) ? 1.0 : 0.0;
   }

   void Softmax(vector &v) const {
      double maxVal = v[0];
      for(int i = 1; i < v.Size(); i++)
         if(v[i] > maxVal) maxVal = v[i];

      double sum = 0.0;
      for(int i = 0; i < v.Size(); i++) {
         v[i] = MathExp(v[i] - maxVal);
         sum += v[i];
      }
      double eps = DATA_EPS(sum);
      if(sum > eps) {
         for(int i = 0; i < v.Size(); i++)
            v[i] /= sum;
      } else {
         double eq = 1.0 / v.Size();
         for(int i = 0; i < v.Size(); i++)
            v[i] = eq;
      }
   }

   void BiasInit(vector &v) {
      for(int i = 0; i < v.Size(); i++)
         v[i] = 0.0;
   }

   void ZeroMatrix(matrix &m) {
      for(int r = 0; r < (int)m.Rows(); r++)
         for(int c = 0; c < (int)m.Cols(); c++)
            m[r][c] = 0.0;
   }

   void ZeroVector(vector &v) {
      for(int i = 0; i < v.Size(); i++)
         v[i] = 0.0;
   }

    void AdamUpdate(matrix &param, matrix &m, matrix &v, matrix &grad, double lr) {
       for(int r = 0; r < (int)param.Rows(); r++) {
          for(int c = 0; c < (int)param.Cols(); c++) {
             double g = grad[r][c] + m_l2 * param[r][c];

             m[r][c] = m_adamBeta1 * m[r][c] + (1.0 - m_adamBeta1) * g;
             v[r][c] = m_adamBeta2 * v[r][c] + (1.0 - m_adamBeta2) * g * g;

             double mHat = m[r][c] / (1.0 - m_beta1T);
             double vHat = v[r][c] / (1.0 - m_beta2T);

             param[r][c] -= lr * mHat / (MathSqrt(vHat) + m_adamEps);
          }
       }
    }

    void AdamUpdate(vector &param, vector &m, vector &v, vector &grad, double lr) {
       for(int i = 0; i < param.Size(); i++) {
          double g = grad[i] + m_l2 * param[i];

          m[i] = m_adamBeta1 * m[i] + (1.0 - m_adamBeta1) * g;
          v[i] = m_adamBeta2 * v[i] + (1.0 - m_adamBeta2) * g * g;

          double mHat = m[i] / (1.0 - m_beta1T);
          double vHat = v[i] / (1.0 - m_beta2T);

          param[i] -= lr * mHat / (MathSqrt(vHat) + m_adamEps);
       }
    }

public:
    CNeuralNet()
       : m_inputs(0), m_hidden(0), m_outputs(0),
         m_t(0), m_beta1T(1.0), m_beta2T(1.0),
         m_l2(0), m_epochsTrained(0), m_lastLoss(0.0),
         m_initialized(false),
         m_dropoutRate(0), m_bnMomentum(0), m_bnEps(0) {}

    void Init(int inputs, int hidden, int outputs, double l2 = -1, double dropoutRate = -1) {
       m_inputs  = inputs;
       m_hidden  = hidden;
       m_outputs = outputs;

       // Iperparametri dall'architettura:
       //   hidden min. = 2 (un neurone non può fare computazione utile)
       //   β₂/β₁ ratio = 100 (Adam originale: 0.999/0.9)
       //   BN window   = 5× hidden (momentum più lento per layer più grandi)
       int totalParams = inputs * hidden + hidden + hidden * outputs + outputs;
       int minHidden = MathMax(hidden, 2);
       m_l2 = 1.0 / MathMax(totalParams, 1);              // L2 = 1/param (1 = nessun parametro → L2=1)
       m_dropoutRate = 1.0 / MathSqrt((double)minHidden); // Dropout: 1/sqrt(hidden)
       m_adamBeta1 = 1.0 - 1.0 / (double)minHidden;       // β₁ = 1 - 1/hidden
       m_adamBeta2 = 1.0 - 1.0 / (double)MathMax(minHidden * 100, 200);  // β₂ = 1 - 1/(hidden×100)
       m_adamEps = DATA_EPS(1.0);                          // Adam ε: machine epsilon
       m_bnMomentum = 1.0 - 1.0 / (double)MathMax(minHidden * 5, 10);  // BN momentum = 1 - 1/(hidden×5)
       m_bnEps = DATA_EPS(1.0);                            // BN ε: machine epsilon

       // Se chiamata esterna vuole override, usa quelli
       if(l2 >= 0)         m_l2 = l2;
       if(dropoutRate >= 0) m_dropoutRate = dropoutRate;

      m_W1.Init(inputs, hidden);
      m_b1.Init(hidden);
      m_W2.Init(hidden, outputs);
      m_b2.Init(outputs);

      // He Init per ReLU (NeuroBook §1.3)
      {
         double scale = HeScale(inputs);
         MathSrand(GetTickCount());
         for(int r = 0; r < inputs; r++)
            for(int c = 0; c < hidden; c++)
               m_W1[r][c] = ((double)MathRand() / 32767.0 * 2.0 - 1.0) * scale;
      }
      BiasInit(m_b1);
      {
         double scale = HeScale(hidden);
         for(int r = 0; r < hidden; r++)
            for(int c = 0; c < outputs; c++)
               m_W2[r][c] = ((double)MathRand() / 32767.0 * 2.0 - 1.0) * scale;
      }
      BiasInit(m_b2);

      m_mW1.Init(inputs, hidden);  ZeroMatrix(m_mW1);
      m_vW1.Init(inputs, hidden);  ZeroMatrix(m_vW1);
      m_mb1.Init(hidden);          ZeroVector(m_mb1);
      m_vb1.Init(hidden);          ZeroVector(m_vb1);
      m_mW2.Init(hidden, outputs); ZeroMatrix(m_mW2);
      m_vW2.Init(hidden, outputs); ZeroMatrix(m_vW2);
      m_mb2.Init(outputs);         ZeroVector(m_mb2);
      m_vb2.Init(outputs);         ZeroVector(m_vb2);

      // Dropout mask
      ArrayResize(m_dropoutMask, hidden);

      // BatchNorm (NeuroBook §6.3)
      m_bnGamma.Init(hidden);
      m_bnBeta.Init(hidden);
      m_bnRunningMean.Init(hidden);
      m_bnRunningVar.Init(hidden);
      for(int i = 0; i < hidden; i++) {
         m_bnGamma[i] = 1.0;
         m_bnBeta[i] = 0.0;
         m_bnRunningMean[i] = 0.0;
         m_bnRunningVar[i] = 1.0;
      }

      m_t        = 0;
      m_beta1T   = 1.0;
      m_beta2T   = 1.0;
      m_epochsTrained = 0;
      m_lastLoss = 0.0;
      m_initialized = true;
      ArrayResize(m_lossHistory, 0);
      m_lossCsvFn = "";
   }

    void ForwardPass(vector &inp, vector &output, vector &h1Cache) {
       h1Cache.Resize(m_hidden);
       for(int i = 0; i < m_hidden; i++) {
          double sum = m_b1[i];
          for(int j = 0; j < m_inputs; j++)
             sum += inp[j] * m_W1[j][i];
          h1Cache[i] = ReLU(sum);
       }
       // BatchNorm inferenza: normalizza con running stats (NeuroBook §6.3)
       ApplyBN(h1Cache);
       // Dropout scaling in inferenza: scale = 1 - dropoutRate (NeuroBook §6.2)
       if(m_dropoutRate > 0) {
          for(int i = 0; i < h1Cache.Size(); i++)
             h1Cache[i] *= (1.0 - m_dropoutRate);
       }

       output.Resize(m_outputs);
       for(int i = 0; i < m_outputs; i++) {
          double sum = m_b2[i];
          for(int j = 0; j < m_hidden; j++)
             sum += h1Cache[j] * m_W2[j][i];
          output[i] = sum;
       }
       Softmax(output);
    }

    void Forward(vector &inp, vector &output) {
       vector h1Cache;
       ForwardPass(inp, output, h1Cache);
    }

    double TrainSample(vector &inp, vector &target, double lr, double weight = 1.0) {
      if(!m_initialized) return -1.0;

      // ── Forward ──
      // Layer 1: W1*x + b1 → z1 → ReLU → h1_raw → Dropout → h1_drop → BN → h1_norm
      vector z1(m_hidden);
      vector h1_raw(m_hidden);
for(int i = 0; i < m_hidden; i++) {
          z1[i] = m_b1[i];
          for(int j = 0; j < m_inputs; j++)
             z1[i] += inp[j] * m_W1[j][i];
          h1_raw[i] = ReLU(z1[i]);
       }

      // Dropout (NeuroBook §6.2): salva in m_dropoutMask, applica a h1_drop
      vector h1_drop(m_hidden);
      if(m_dropoutRate > 0.0) {
         for(int i = 0; i < m_hidden; i++) {
            m_dropoutMask[i] = ((double)MathRand() / 32767.0) > m_dropoutRate ? 1.0 : 0.0;
            h1_drop[i] = h1_raw[i] * m_dropoutMask[i];
         }
      } else {
         for(int i = 0; i < m_hidden; i++)
            h1_drop[i] = h1_raw[i];
      }

      // BatchNorm (NeuroBook §6.3): normalizza h1_drop → h1_norm, aggiorna running stats
      vector h1_cache = h1_drop;  // copia per backward
      vector h1_norm(m_hidden);
      TrainBN(h1_drop, h1_norm);

      // Layer 2: W2*h1_norm + b2 → z2 → Softmax → output
      vector z2(m_outputs);
      for(int i = 0; i < m_outputs; i++) {
         z2[i] = m_b2[i];
         for(int j = 0; j < m_hidden; j++)
            z2[i] += h1_norm[j] * m_W2[j][i];
      }
      vector output(m_outputs);
      for(int i = 0; i < m_outputs; i++) output[i] = z2[i];
      Softmax(output);

      // ── Loss: CCE pesata + L2 ──
      double loss = 0.0;
      for(int i = 0; i < m_outputs; i++) {
         double p = MathMax(DATA_EPS(output[i]), output[i]);
         loss -= weight * target[i] * MathLog(p);
      }
      loss += 0.5 * m_l2 * (L2Norm(m_W1) + L2Norm(m_W2));
      m_lastLoss = loss;

      m_t++;
      m_beta1T *= m_adamBeta1;
      m_beta2T *= m_adamBeta2;

      // ── Backward ──
      // dL/dz2 = weight * (output - target) (CCE+Softmax combinata, pesata)
      vector dL_dz2(m_outputs);
      for(int i = 0; i < m_outputs; i++)
         dL_dz2[i] = weight * (output[i] - target[i]);

      // dL/dW2 = h1_norm ⊗ dL/dz2
      matrix dL_dW2(m_hidden, m_outputs);
      for(int i = 0; i < m_hidden; i++)
         for(int j = 0; j < m_outputs; j++)
            dL_dW2[i][j] = h1_norm[i] * dL_dz2[j];

      // dL/db2 = dL/dz2
      vector dL_db2(m_outputs);
      for(int i = 0; i < m_outputs; i++)
         dL_db2[i] = dL_dz2[i];

      // dL/dh1_norm = W2^T * dL/dz2
      vector dL_dh1_norm(m_hidden);
      for(int i = 0; i < m_hidden; i++) {
         double sum = 0.0;
         for(int j = 0; j < m_outputs; j++)
            sum += dL_dz2[j] * m_W2[i][j];
         dL_dh1_norm[i] = sum;
      }

      // BatchNorm backward: dL/dh1_drop = gamma/sqrt(var) * dL/dh1_norm
      // + update gamma, beta
      vector dL_dh1_drop(m_hidden);
      for(int i = 0; i < m_hidden; i++) {
          double denom = MathSqrt(m_bnRunningVar[i] + m_bnEps);
          double h_std = (h1_cache[i] - m_bnRunningMean[i]) / denom;
          double gOld = m_bnGamma[i];  // salva gamma pre-update per backward pass-through
          m_bnGamma[i] -= lr * dL_dh1_norm[i] * h_std;
          m_bnBeta[i]  -= lr * dL_dh1_norm[i];
          dL_dh1_drop[i] = dL_dh1_norm[i] * gOld / denom;
      }

      // Dropout backward: applica stessa mask
      if(m_dropoutRate > 0.0) {
         for(int i = 0; i < m_hidden; i++)
            dL_dh1_drop[i] *= m_dropoutMask[i];
      }

      // dL/dz1 = dL/dh1_drop * ReLU'(z1)
      vector dL_dz1(m_hidden);
      for(int i = 0; i < m_hidden; i++)
         dL_dz1[i] = dL_dh1_drop[i] * ReLUDeriv(z1[i]);

      // dL/dW1 = inp ⊗ dL/dz1
      matrix dL_dW1(m_inputs, m_hidden);
      for(int i = 0; i < m_inputs; i++)
         for(int j = 0; j < m_hidden; j++)
            dL_dW1[i][j] = inp[i] * dL_dz1[j];

      // dL/db1 = dL/dz1
      vector dL_db1(m_hidden);
      for(int i = 0; i < m_hidden; i++)
         dL_db1[i] = dL_dz1[i];

      // Adam update
      AdamUpdate(m_W1, m_mW1, m_vW1, dL_dW1, lr);
      AdamUpdate(m_b1, m_mb1, m_vb1, dL_db1, lr);
      AdamUpdate(m_W2, m_mW2, m_vW2, dL_dW2, lr);
      AdamUpdate(m_b2, m_mb2, m_vb2, dL_db2, lr);

      return loss;
   }

   double L2Norm(matrix &m) {
      double sum = 0.0;
      for(int r = 0; r < (int)m.Rows(); r++)
         for(int c = 0; c < (int)m.Cols(); c++)
            sum += m[r][c] * m[r][c];
      return sum;
   }

   // Overload senza pesi (compatibilità)
   double Train(matrix &features, matrix &targets, int epochs, double lr) {
      vector empty;
      return Train(features, targets, epochs, lr, empty);
   }

   double Train(matrix &features, matrix &targets, int epochs, double lr, vector &weights) {
      if(!m_initialized || features.Rows() == 0) return -1.0;

      int n = (int)features.Rows();
      double avgLoss = 0.0;
      ArrayResize(m_lossHistory, epochs);
      int logEvery = MathMax(1, epochs / 10);

      for(int epoch = 0; epoch < epochs; epoch++) {
         int indices[];
         ArrayResize(indices, n);
         for(int i = 0; i < n; i++) indices[i] = i;

         for(int i = n - 1; i > 0; i--) {
            int j = MathRand() % (i + 1);
            int tmp = indices[i];
            indices[i] = indices[j];
            indices[j] = tmp;
         }

         double epochLoss = 0.0;
         for(int s = 0; s < n; s++) {
            int idx = indices[s];
            vector inp = features.Row(idx);
            vector tgt = targets.Row(idx);
            double w = (weights.Size() > idx) ? weights[idx] : 1.0;
            epochLoss += TrainSample(inp, tgt, lr, w);
         }
         epochLoss /= (double)n;
         m_lossHistory[epoch] = epochLoss;

         if(epoch == epochs - 1) avgLoss = epochLoss;

         if(epoch == 0 || epoch == epochs - 1 || (epoch+1) % logEvery == 0) {
            Print("  Epoch ", epoch+1, "/", epochs,
                  " | loss: ", StringFormat("%.6f", epochLoss),
                  " | lr: ", StringFormat("%.5f", lr));
         }
      }

      m_epochsTrained += epochs;
      m_lastLoss = avgLoss;
      return avgLoss;
   }

   // Salva loss history in CSV nella cartella Common
   bool SaveLossCsv(string filename = "") {
      if(ArraySize(m_lossHistory) == 0) return false;
      if(filename == "") filename = "NN_LossHistory.csv";
      int fh = FileOpen(filename, FILE_TXT|FILE_WRITE|FILE_COMMON);
      if(fh == INVALID_HANDLE) return false;

      FileWriteString(fh, "epoch,loss\r\n");
      for(int i = 0; i < ArraySize(m_lossHistory); i++) {
         FileWriteString(fh, (string)(i+1) + "," + StringFormat("%.8f", m_lossHistory[i]) + "\r\n");
      }
      FileClose(fh);
      Print("Loss history saved to ", filename, " (", ArraySize(m_lossHistory), " epochs)");
      return true;
   }

   string LossHistorySummary() const {
      if(ArraySize(m_lossHistory) == 0) return "no history";
      double first = m_lossHistory[0];
      double last  = m_lossHistory[ArraySize(m_lossHistory)-1];
      double best  = first;
      int bestEpoch = 0;
      for(int i = 0; i < ArraySize(m_lossHistory); i++) {
         if(m_lossHistory[i] < best) { best = m_lossHistory[i]; bestEpoch = i; }
      }
      return StringFormat("loss: %.6f → %.6f (best: %.6f @ epoch %d, %d epochs)",
         first, last, best, bestEpoch+1, ArraySize(m_lossHistory));
   }

int Predict(vector &inp) {
       vector output;
       Forward(inp, output);
      int bestIdx = 0;
      double bestVal = output[0];
      for(int i = 1; i < m_outputs; i++) {
         if(output[i] > bestVal) {
            bestVal = output[i];
            bestIdx = i;
         }
      }
      return bestIdx;
   }

double GetCombinedZ(vector &inp) {
       vector output;
       Forward(inp, output);
      return output[0] - output[2];
   }

   bool Save(string filename) {
      int fh = FileOpen(filename, FILE_WRITE | FILE_BIN | FILE_COMMON);
      if(fh == INVALID_HANDLE) return false;
      bool ok = Save(fh);
      FileClose(fh);
      return ok;
   }

   bool Save(int fh) {
      if(fh == INVALID_HANDLE) return false;

      FileWriteInteger(fh, 3);
      FileWriteInteger(fh, m_inputs);
      FileWriteInteger(fh, m_hidden);
      FileWriteInteger(fh, m_outputs);
      FileWriteInteger(fh, m_epochsTrained);
      FileWriteDouble(fh, m_lastLoss);
      FileWriteDouble(fh, m_l2);
      FileWriteDouble(fh, m_dropoutRate);
      FileWriteDouble(fh, m_bnMomentum);
      FileWriteDouble(fh, m_bnEps);

      for(int r = 0; r < m_inputs; r++)
         for(int c = 0; c < m_hidden; c++)
            FileWriteDouble(fh, m_W1[r][c]);

      for(int i = 0; i < m_hidden; i++)
         FileWriteDouble(fh, m_b1[i]);

      for(int r = 0; r < m_hidden; r++)
         for(int c = 0; c < m_outputs; c++)
            FileWriteDouble(fh, m_W2[r][c]);

      for(int i = 0; i < m_outputs; i++)
         FileWriteDouble(fh, m_b2[i]);

      for(int r = 0; r < m_inputs; r++)
         for(int c = 0; c < m_hidden; c++) {
            FileWriteDouble(fh, m_mW1[r][c]);
            FileWriteDouble(fh, m_vW1[r][c]);
         }

      for(int i = 0; i < m_hidden; i++) {
         FileWriteDouble(fh, m_mb1[i]);
         FileWriteDouble(fh, m_vb1[i]);
      }

      for(int r = 0; r < m_hidden; r++)
         for(int c = 0; c < m_outputs; c++) {
            FileWriteDouble(fh, m_mW2[r][c]);
            FileWriteDouble(fh, m_vW2[r][c]);
         }

      for(int i = 0; i < m_outputs; i++) {
         FileWriteDouble(fh, m_mb2[i]);
         FileWriteDouble(fh, m_vb2[i]);
      }

      FileWriteDouble(fh, m_beta1T);
      FileWriteDouble(fh, m_beta2T);
      FileWriteInteger(fh, m_t);

      // BatchNorm (v3)
      for(int i = 0; i < m_hidden; i++) {
         FileWriteDouble(fh, m_bnGamma[i]);
         FileWriteDouble(fh, m_bnBeta[i]);
         FileWriteDouble(fh, m_bnRunningMean[i]);
         FileWriteDouble(fh, m_bnRunningVar[i]);
      }

      return true;
   }

   bool Load(string filename) {
      if(!FileIsExist(filename, FILE_COMMON)) return false;

      int fh = FileOpen(filename, FILE_READ | FILE_BIN | FILE_COMMON);
      if(fh == INVALID_HANDLE) return false;
      bool ok = Load(fh);
      FileClose(fh);
      return ok;
   }

   bool Load(int fh) {
      if(fh == INVALID_HANDLE) return false;

      int version = FileReadInteger(fh);
      int inputs  = FileReadInteger(fh);
      int hidden  = FileReadInteger(fh);
      int outputs = FileReadInteger(fh);

      Init(inputs, hidden, outputs);

      m_epochsTrained = FileReadInteger(fh);
      m_lastLoss = FileReadDouble(fh);
      if(version >= 2) m_l2 = FileReadDouble(fh);
      if(version >= 3) {
         m_dropoutRate = FileReadDouble(fh);
         m_bnMomentum  = FileReadDouble(fh);
         m_bnEps       = FileReadDouble(fh);
      }

      for(int r = 0; r < m_inputs; r++)
         for(int c = 0; c < m_hidden; c++)
            m_W1[r][c] = FileReadDouble(fh);

      for(int i = 0; i < m_hidden; i++)
         m_b1[i] = FileReadDouble(fh);

      for(int r = 0; r < m_hidden; r++)
         for(int c = 0; c < m_outputs; c++)
            m_W2[r][c] = FileReadDouble(fh);

      for(int i = 0; i < m_outputs; i++)
         m_b2[i] = FileReadDouble(fh);

      if(version >= 1) {
         for(int r = 0; r < m_inputs; r++)
            for(int c = 0; c < m_hidden; c++) {
               m_mW1[r][c] = FileReadDouble(fh);
               m_vW1[r][c] = FileReadDouble(fh);
            }

         for(int i = 0; i < m_hidden; i++) {
            m_mb1[i] = FileReadDouble(fh);
            m_vb1[i] = FileReadDouble(fh);
         }

         for(int r = 0; r < m_hidden; r++)
            for(int c = 0; c < m_outputs; c++) {
               m_mW2[r][c] = FileReadDouble(fh);
               m_vW2[r][c] = FileReadDouble(fh);
            }

         for(int i = 0; i < m_outputs; i++) {
            m_mb2[i] = FileReadDouble(fh);
            m_vb2[i] = FileReadDouble(fh);
         }

           if(version >= 1) {
              m_beta1T = FileReadDouble(fh);
              m_beta2T = FileReadDouble(fh);
              m_t      = FileReadInteger(fh);
           }

           // BatchNorm (v3)
           if(version >= 3) {
              m_bnGamma.Init(m_hidden);
              m_bnBeta.Init(m_hidden);
              m_bnRunningMean.Init(m_hidden);
              m_bnRunningVar.Init(m_hidden);
              for(int i = 0; i < m_hidden; i++) {
                 m_bnGamma[i]       = FileReadDouble(fh);
                 m_bnBeta[i]        = FileReadDouble(fh);
                 m_bnRunningMean[i] = FileReadDouble(fh);
                 m_bnRunningVar[i]  = FileReadDouble(fh);
              }
           }
        }

        return true;
     }

   bool IsInitialized() const { return m_initialized; }
   int  EpochsTrained()  const { return m_epochsTrained; }
   double LastLoss()     const { return m_lastLoss; }
   int  Inputs()         const { return m_inputs; }
   int  Hidden()         const { return m_hidden; }
   int  Outputs()        const { return m_outputs; }

   string Info() const {
      if(!m_initialized) return "NN: uninitialized";
      return StringFormat("NN: %d→%d→%d (epoche=%d, loss=%.6f, dropout=%.2f, BN=%s)",
         m_inputs, m_hidden, m_outputs, m_epochsTrained, m_lastLoss,
         m_dropoutRate, m_initialized ? "on" : "off");
   }

   string WeightsSummary() const {
      if(!m_initialized) return "NN: uninitialized";
      double w1Min = 1e99, w1Max = -1e99, w1Sum = 0;
      double w2Min = 1e99, w2Max = -1e99, w2Sum = 0;
      int w1Count = 0, w2Count = 0;
      for(int r = 0; r < m_inputs; r++) {
         for(int c = 0; c < m_hidden; c++) {
            double v = m_W1[r][c];
            if(v < w1Min) w1Min = v; if(v > w1Max) w1Max = v;
            w1Sum += v; w1Count++;
         }
      }
      for(int r = 0; r < m_hidden; r++) {
         for(int c = 0; c < m_outputs; c++) {
            double v = m_W2[r][c];
            if(v < w2Min) w2Min = v; if(v > w2Max) w2Max = v;
            w2Sum += v; w2Count++;
         }
      }
      return StringFormat("  W1 [%.4f, %.4f] μ=%.4f | W2 [%.4f, %.4f] μ=%.4f",
         w1Min, w1Max, (w1Count>0?w1Sum/w1Count:0),
         w2Min, w2Max, (w2Count>0?w2Sum/w2Count:0));
   }
};

struct NNTrainSample {
   double features[NN_FEATURES];
   double target[NN_TARGETS];
   double weight;    // peso del sample: trades con grosso impatto pesano di più
};

class NNTrainBuffer {
private:
   NNTrainSample m_samples[];
   int m_count;
public:
   NNTrainBuffer() : m_count(0) {}

   void Add(double &features[], double &target[], double weight = 1.0) {
      int idx = m_count;
      ArrayResize(m_samples, idx + 1);
      for(int i = 0; i < NN_FEATURES; i++) m_samples[idx].features[i] = features[i];
      for(int i = 0; i < NN_TARGETS; i++) m_samples[idx].target[i]   = target[i];
      m_samples[idx].weight = weight;
      m_count++;
   }

   int Count() const { return m_count; }

   void ToMatrices(matrix &features, matrix &targets) {
      if(m_count == 0) return;
      features.Init(m_count, NN_FEATURES);
      targets.Init(m_count, NN_TARGETS);
      for(int i = 0; i < m_count; i++) {
         for(int j = 0; j < NN_FEATURES; j++)
            features[i][j] = m_samples[i].features[j];
         for(int j = 0; j < NN_TARGETS; j++)
            targets[i][j] = m_samples[i].target[j];
      }
   }

   vector GetWeights() {
      vector w(m_count);
      for(int i = 0; i < m_count; i++)
         w[i] = m_samples[i].weight;
      return w;
   }

   void Clear() {
      ArrayResize(m_samples, 0);
      m_count = 0;
   }

    void Trim(int maxSamples) {
       if(m_count <= maxSamples) return;
       int remove = m_count - maxSamples;
       for(int i = 0; i < maxSamples; i++)
          m_samples[i] = m_samples[i + remove];
       ArrayResize(m_samples, maxSamples);
       m_count = maxSamples;
    }

    // Dump training samples to CSV nella cartella Common per debug
    void SaveCsv(string filename = "NN_TrainingSamples.csv") {
       int fh = FileOpen(filename, FILE_TXT|FILE_WRITE|FILE_COMMON);
       if(fh == INVALID_HANDLE) return;
        FileWriteString(fh, "sample");
        for(int f = 0; f < NN_FEATURES; f++) FileWriteString(fh, ",feat_" + (string)f);
        for(int t = 0; t < NN_TARGETS; t++) FileWriteString(fh, ",target_" + (string)t);
        FileWriteString(fh, ",label\r\n");
        for(int i = 0; i < m_count; i++) {
           string line = (string)i;
           for(int f = 0; f < NN_FEATURES; f++) line += "," + StringFormat("%+.6f", m_samples[i].features[f]);
           for(int t = 0; t < NN_TARGETS; t++) line += "," + StringFormat("%.0f", m_samples[i].target[t]);
          // Label leggibile
          if(m_samples[i].target[0] == 1) line += ",BUY";
          else if(m_samples[i].target[2] == 1) line += ",SELL";
          else line += ",FLAT";
          line += "\r\n";
          FileWriteString(fh, line);
       }
       FileClose(fh);
       Print("Training samples saved to ", filename, " (", m_count, " samples)");
    }

    // Stampa statistiche riassuntive del dataset
    void PrintStats() {
       int buyCount = 0, sellCount = 0, flatCount = 0;
       for(int i = 0; i < m_count; i++) {
          if(m_samples[i].target[0] == 1) buyCount++;
          else if(m_samples[i].target[2] == 1) sellCount++;
          else flatCount++;
       }
       Print("  Dataset: ", m_count, " samples | Buy: ", buyCount,
             " Sell: ", sellCount, " Flat: ", flatCount);
    }

    // Media feature per classe (utile per vedere se le feature discriminano)
    void PrintFeatureStats() {
       if(m_count < 3) { Print("  Feature stats: too few samples"); return; }
       double meanFeat[NN_FEATURES] = {};
       double absMean[NN_FEATURES] = {};
       for(int i = 0; i < m_count; i++) {
          for(int f = 0; f < NN_FEATURES; f++) {
             meanFeat[f] += m_samples[i].features[f];
             absMean[f] += MathAbs(m_samples[i].features[f]);
          }
       }
       Print("  Feature means (signal strength per class):");
       for(int f = 0; f < NN_FEATURES; f++) {
          meanFeat[f] /= m_count;
          absMean[f] /= m_count;
       }
       string labels[NN_FEATURES] = {"Hurst","ADX","MA","Momentum","Consensus","Hunter","Agreement","TrendStr"};
       for(int f = 0; f < NN_FEATURES; f++) {
          Print("    [", f, "] ", labels[f], ": μ=", StringFormat("%+.4f", meanFeat[f]),
                " |avg|=", StringFormat("%.4f", absMean[f]));
       }
    }
 };
#endif