// AutoARIMA_liteC.c
// Input convention:
//   Close[0] = most recent price
//   Close[1] = previous price
//   Close[2] = older price
//
// Notes:
// - No #include needed in Zorro lite-C.
// - Uses var/var* instead of double/double*.
// - Uses malloc/free/memset instead of calloc.
// - No ternary operator, no const qualifiers.

#define AUTO_ARIMA_MAX_ITER 1000
#define AUTO_ARIMA_BIG 1e100

typedef struct AUTO_ARIMA_RESULT
{
  int p;
  int d;
  int q;
  var c;
  var sse;
  var aicc;
  var forecast;
  var* ar;
  var* ma;
} AUTO_ARIMA_RESULT;

// ---------- memory helpers ----------

var* aa_alloc_vars(int Count)
{
  if(Count <= 0)
    Count = 1;

  var* Ptr = malloc(Count*sizeof(var));
  if(Ptr == 0) {
    printf("\nAutoARIMA: memory allocation failed.");
    return 0;
  }

  memset(Ptr,0,Count*sizeof(var));
  return Ptr;
}

void aa_free_vars(var* Ptr)
{
  if(Ptr != 0)
    free(Ptr);
}

void init_auto_arima_result(AUTO_ARIMA_RESULT* R)
{
  if(R == 0)
    return;

  R->p = -1;
  R->d = 0;
  R->q = -1;
  R->c = 0.;
  R->sse = AUTO_ARIMA_BIG;
  R->aicc = AUTO_ARIMA_BIG;
  R->forecast = 0.;
  R->ar = 0;
  R->ma = 0;
}

void free_auto_arima_result(AUTO_ARIMA_RESULT* R)
{
  if(R == 0)
    return;

  aa_free_vars(R->ar);
  aa_free_vars(R->ma);
  R->ar = 0;
  R->ma = 0;
}

// ---------- small math helpers ----------

int aa_max_int(int A,int B)
{
  if(A > B)
    return A;
  return B;
}

var aa_mean(vars Series,int N)
{
  if(Series == 0)
    return 0.;
  if(N <= 0)
    return 0.;

  var Sum = 0.;
  int i;
  for(i = 0; i < N; i++)
    Sum += Series[i];

  return Sum / (var)N;
}

var aa_variance(vars Series,int N)
{
  if(Series == 0)
    return 0.;
  if(N < 2)
    return 0.;

  var Mu = aa_mean(Series,N);
  var Sum = 0.;
  int i;
  for(i = 0; i < N; i++) {
    var E = Series[i] - Mu;
    Sum += E * E;
  }

  return Sum / (var)N;
}

int aa_calculate_d(vars Series,int N)
{
  if(Series == 0)
    return 0;
  if(N < 3)
    return 0;

  var Var0 = aa_variance(Series,N);

  var* Diff1 = aa_alloc_vars(N-1);
  if(Diff1 == 0)
    return 0;

  int i;
  for(i = 0; i < N-1; i++)
    Diff1[i] = Series[i] - Series[i+1];

  var Var1 = aa_variance(Diff1,N-1);

  if(Var1 > 0.5 * Var0) {
    aa_free_vars(Diff1);
    return 0;
  }

  if(Var1 < 0.1 * Var0) {
    aa_free_vars(Diff1);
    return 1;
  }

  var* Diff2 = aa_alloc_vars(N-2);
  if(Diff2 == 0) {
    aa_free_vars(Diff1);
    return 1;
  }

  for(i = 0; i < N-2; i++)
    Diff2[i] = Diff1[i] - Diff1[i+1];

  var Var2 = aa_variance(Diff2,N-2);

  aa_free_vars(Diff1);
  aa_free_vars(Diff2);

  if(Var2 < 0.5 * Var1)
    return 2;

  return 1;
}

// Applies d-order differencing and writes chronological output:
//   Diff[0] = oldest transformed value
//   Diff[last] = newest transformed value
// Diff must be allocated by the caller with at least N elements.
// Returns the number of elements written, or 0 on error.
int aa_difference_series(int D,vars Close,int N,vars Diff)
{
  if(Close == 0)
    return 0;
  if(Diff == 0)
    return 0;
  if(N <= D)
    return 0;

  int i;

  if(D == 0) {
    for(i = 0; i < N; i++)
      Diff[i] = Close[N-1-i];
    return N;
  }

  if(D == 1) {
    int M = N-1;
    var* Tmp = aa_alloc_vars(M);
    if(Tmp == 0)
      return 0;

    for(i = 0; i < M; i++)
      Tmp[i] = Close[i] - Close[i+1];

    for(i = 0; i < M; i++)
      Diff[i] = Tmp[M-1-i];

    aa_free_vars(Tmp);
    return M;
  }

  if(D == 2) {
    int M1 = N-1;
    int M2 = N-2;

    var* First = aa_alloc_vars(M1);
    var* Tmp = aa_alloc_vars(M2);
    if(First == 0) {
      aa_free_vars(Tmp);
      return 0;
    }
    if(Tmp == 0) {
      aa_free_vars(First);
      return 0;
    }

    for(i = 0; i < M1; i++)
      First[i] = Close[i] - Close[i+1];

    for(i = 0; i < M2; i++)
      Tmp[i] = First[i] - First[i+1];

    for(i = 0; i < M2; i++)
      Diff[i] = Tmp[M2-1-i];

    aa_free_vars(First);
    aa_free_vars(Tmp);
    return M2;
  }

  return 0;
}

var aa_aicc_score(int N,var SSE,int P,int Q)
{
  int K = P + Q + 1;

  if(N - K - 1 <= 0)
    return AUTO_ARIMA_BIG;

  if(SSE <= 0.)
    return AUTO_ARIMA_BIG;

  var Sigma2 = SSE / (var)N;
  if(Sigma2 < 0.000000000000001)
    return AUTO_ARIMA_BIG;

  var AIC = (var)N * log(Sigma2) + 2. * (var)K;
  var AICc = AIC + (2. * (var)K * (var)(K + 1)) / (var)(N - K - 1);

  return AICc;
}

// Trains ARMA(P,Q) on chronological Series.
// ArOut must have at least max(1,P) elements.
// MaOut must have at least max(1,Q) elements.
int aa_arma_fit(int P,int Q,vars Series,int N,var* OutSSE,var* OutC,vars ArOut,vars MaOut)
{
  if(Series == 0)
    return 0;
  if(N < 2)
    return 0;
  if(OutSSE == 0)
    return 0;
  if(OutC == 0)
    return 0;
  if(ArOut == 0)
    return 0;
  if(MaOut == 0)
    return 0;

  *OutSSE = 0.;
  *OutC = aa_mean(Series,N);

  int PSize = 1;
  int QSize = 1;
  int ArGradSize = 1;
  int MaGradSize = 1;

  if(P > 0)
    PSize = P;
  if(Q > 0)
    QSize = Q;
  if(P > 0)
    ArGradSize = N * P;
  if(Q > 0)
    MaGradSize = N * Q;

  var* Ar = aa_alloc_vars(PSize);
  var* Ma = aa_alloc_vars(QSize);
  var* MAr = aa_alloc_vars(PSize);
  var* VAr = aa_alloc_vars(PSize);
  var* MMa = aa_alloc_vars(QSize);
  var* VMa = aa_alloc_vars(QSize);
  var* Errors = aa_alloc_vars(N);
  var* DEdC = aa_alloc_vars(N);
  var* DEdAR = aa_alloc_vars(ArGradSize);
  var* DEdMA = aa_alloc_vars(MaGradSize);
  var* GradAR = aa_alloc_vars(PSize);
  var* GradMA = aa_alloc_vars(QSize);

  if(Ar == 0 || Ma == 0 || MAr == 0 || VAr == 0 || MMa == 0 || VMa == 0 ||
     Errors == 0 || DEdC == 0 || DEdAR == 0 || DEdMA == 0 || GradAR == 0 || GradMA == 0) {
    aa_free_vars(Ar);
    aa_free_vars(Ma);
    aa_free_vars(MAr);
    aa_free_vars(VAr);
    aa_free_vars(MMa);
    aa_free_vars(VMa);
    aa_free_vars(Errors);
    aa_free_vars(DEdC);
    aa_free_vars(DEdAR);
    aa_free_vars(DEdMA);
    aa_free_vars(GradAR);
    aa_free_vars(GradMA);
    return 0;
  }

  var C = *OutC;

  var Beta1 = 0.9;
  var Beta2 = 0.999;
  var Eps = 0.00000001;
  var Eta = 0.001;

  var MC = 0.;
  var VC = 0.;
  var PrevSSE = AUTO_ARIMA_BIG;
  var Beta1Pow = 1.;
  var Beta2Pow = 1.;

  int Iter;
  int i;
  int j;
  int k;
  int t;

  for(Iter = 1; Iter < AUTO_ARIMA_MAX_ITER; Iter++) {
    var SSELocal = 0.;

    for(i = 0; i < N; i++) {
      Errors[i] = 0.;
      DEdC[i] = 0.;
    }

    for(i = 0; i < ArGradSize; i++)
      DEdAR[i] = 0.;

    for(i = 0; i < MaGradSize; i++)
      DEdMA[i] = 0.;

    for(i = 0; i < PSize; i++)
      GradAR[i] = 0.;

    for(i = 0; i < QSize; i++)
      GradMA[i] = 0.;

    Beta1Pow *= Beta1;
    Beta2Pow *= Beta2;

    // Forward pass: predictions and residuals.
    for(t = 1; t < N; t++) {
      var Pred = C;

      for(i = 0; i < P; i++) {
        if(t - 1 - i >= 0)
          Pred += Ar[i] * Series[t - 1 - i];
      }

      for(j = 0; j < Q; j++) {
        if(t - 1 - j >= 0)
          Pred += Ma[j] * Errors[t - 1 - j];
      }

      Errors[t] = Series[t] - Pred;
      SSELocal += Errors[t] * Errors[t];
    }

    if(Iter > 1) {
      if(abs(PrevSSE - SSELocal) < 0.0000000001) {
        *OutSSE = SSELocal;
        break;
      }
    }

    PrevSSE = SSELocal;

    // Backpropagation through the residual recursion.
    for(t = 1; t < N; t++) {
      DEdC[t] = -1.;

      for(j = 0; j < Q; j++) {
        if(t - 1 - j >= 0)
          DEdC[t] -= Ma[j] * DEdC[t - 1 - j];
      }

      for(i = 0; i < P; i++) {
        var DerAR = 0.;
        if(t - 1 - i >= 0)
          DerAR = -Series[t - 1 - i];

        for(j = 0; j < Q; j++) {
          if(t - 1 - j >= 0)
            DerAR -= Ma[j] * DEdAR[(t - 1 - j) * P + i];
        }

        DEdAR[t * P + i] = DerAR;
      }

      for(j = 0; j < Q; j++) {
        var DerMA = 0.;
        if(t - 1 - j >= 0)
          DerMA = -Errors[t - 1 - j];

        for(k = 0; k < Q; k++) {
          if(t - 1 - k >= 0)
            DerMA -= Ma[k] * DEdMA[(t - 1 - k) * Q + j];
        }

        DEdMA[t * Q + j] = DerMA;
      }
    }

    var GradC = 0.;

    for(t = 1; t < N; t++) {
      GradC += 2. * Errors[t] * DEdC[t];

      for(i = 0; i < P; i++)
        GradAR[i] += 2. * Errors[t] * DEdAR[t * P + i];

      for(j = 0; j < Q; j++)
        GradMA[j] += 2. * Errors[t] * DEdMA[t * Q + j];
    }

    GradC /= (var)N;

    for(i = 0; i < P; i++)
      GradAR[i] /= (var)N;

    for(j = 0; j < Q; j++)
      GradMA[j] /= (var)N;

    // Adam update for C.
    MC = Beta1 * MC + (1. - Beta1) * GradC;
    VC = Beta2 * VC + (1. - Beta2) * GradC * GradC;

    var MCHat = MC / (1. - Beta1Pow);
    var VCHat = VC / (1. - Beta2Pow);
    C -= Eta * MCHat / (sqrt(VCHat) + Eps);

    // Adam update for AR coefficients.
    for(i = 0; i < P; i++) {
      MAr[i] = Beta1 * MAr[i] + (1. - Beta1) * GradAR[i];
      VAr[i] = Beta2 * VAr[i] + (1. - Beta2) * GradAR[i] * GradAR[i];

      var MHatAR = MAr[i] / (1. - Beta1Pow);
      var VHatAR = VAr[i] / (1. - Beta2Pow);

      Ar[i] -= Eta * MHatAR / (sqrt(VHatAR) + Eps);
      Ar[i] = clamp(Ar[i],-0.99,0.99);
    }

    // Adam update for MA coefficients.
    for(j = 0; j < Q; j++) {
      MMa[j] = Beta1 * MMa[j] + (1. - Beta1) * GradMA[j];
      VMa[j] = Beta2 * VMa[j] + (1. - Beta2) * GradMA[j] * GradMA[j];

      var MHatMA = MMa[j] / (1. - Beta1Pow);
      var VHatMA = VMa[j] / (1. - Beta2Pow);

      Ma[j] -= Eta * MHatMA / (sqrt(VHatMA) + Eps);
      Ma[j] = clamp(Ma[j],-0.99,0.99);
    }
  }

  // Final SSE with optimized coefficients.
  for(i = 0; i < N; i++)
    Errors[i] = 0.;

  var FinalSSE = 0.;
  for(t = 1; t < N; t++) {
    var PredFinal = C;

    for(i = 0; i < P; i++) {
      if(t - 1 - i >= 0)
        PredFinal += Ar[i] * Series[t - 1 - i];
    }

    for(j = 0; j < Q; j++) {
      if(t - 1 - j >= 0)
        PredFinal += Ma[j] * Errors[t - 1 - j];
    }

    Errors[t] = Series[t] - PredFinal;
    FinalSSE += Errors[t] * Errors[t];
  }

  *OutSSE = FinalSSE;
  *OutC = C;

  for(i = 0; i < P; i++)
    ArOut[i] = Ar[i];

  for(j = 0; j < Q; j++)
    MaOut[j] = Ma[j];

  aa_free_vars(Ar);
  aa_free_vars(Ma);
  aa_free_vars(MAr);
  aa_free_vars(VAr);
  aa_free_vars(MMa);
  aa_free_vars(VMa);
  aa_free_vars(Errors);
  aa_free_vars(DEdC);
  aa_free_vars(DEdAR);
  aa_free_vars(DEdMA);
  aa_free_vars(GradAR);
  aa_free_vars(GradMA);

  return 1;
}

var aa_round_to_tick_size(var Price,var TickSize)
{
  if(TickSize <= 0.)
    return Price;

  return round(Price / TickSize) * TickSize;
}

// Main forecast function.
// max_p and max_q are exclusive upper bounds, same as the uploaded C version:
// max_p = 5 tests P = 0..4; max_q = 5 tests Q = 0..4.
// Result owns Result->ar and Result->ma after success; call free_auto_arima_result(&Result).
int auto_arima_forecast(vars Close,int N,var TickSize,int MaxP,int MaxQ,AUTO_ARIMA_RESULT* Result)
{
  if(Result == 0)
    return 0;

  init_auto_arima_result(Result);

  if(Close == 0)
    return 0;
  if(N < 5)
    return 0;
  if(MaxP <= 0)
    return 0;
  if(MaxQ <= 0)
    return 0;

  Result->forecast = Close[0];

  int D = aa_calculate_d(Close,N);
  Result->d = D;

  var* Diff = aa_alloc_vars(N);
  if(Diff == 0)
    return 0;

  int NDiff = aa_difference_series(D,Close,N,Diff);
  if(NDiff < 2) {
    aa_free_vars(Diff);
    return 0;
  }

  int MaxArSize = aa_max_int(MaxP,1);
  int MaxMaSize = aa_max_int(MaxQ,1);

  var* CandAR = aa_alloc_vars(MaxArSize);
  var* CandMA = aa_alloc_vars(MaxMaSize);
  var* BestAR = aa_alloc_vars(MaxArSize);
  var* BestMA = aa_alloc_vars(MaxMaSize);

  if(CandAR == 0 || CandMA == 0 || BestAR == 0 || BestMA == 0) {
    aa_free_vars(CandAR);
    aa_free_vars(CandMA);
    aa_free_vars(BestAR);
    aa_free_vars(BestMA);
    aa_free_vars(Diff);
    return 0;
  }

  int BestP = -1;
  int BestQ = -1;
  var BestC = 0.;
  var BestSSE = AUTO_ARIMA_BIG;
  var BestAICc = AUTO_ARIMA_BIG;

  int p;
  int q;
  int i;
  int j;

  for(p = 0; p < MaxP; p++) {
    for(q = 0; q < MaxQ; q++) {
      var SSE = 0.;
      var C = 0.;

      for(i = 0; i < MaxArSize; i++)
        CandAR[i] = 0.;
      for(j = 0; j < MaxMaSize; j++)
        CandMA[j] = 0.;

      if(!aa_arma_fit(p,q,Diff,NDiff,&SSE,&C,CandAR,CandMA))
        continue;

      var Score = aa_aicc_score(NDiff,SSE,p,q);

      if(Score < BestAICc) {
        BestAICc = Score;
        BestP = p;
        BestQ = q;
        BestC = C;
        BestSSE = SSE;

        for(i = 0; i < MaxArSize; i++)
          BestAR[i] = 0.;
        for(j = 0; j < MaxMaSize; j++)
          BestMA[j] = 0.;

        for(i = 0; i < p; i++)
          BestAR[i] = CandAR[i];
        for(j = 0; j < q; j++)
          BestMA[j] = CandMA[j];
      }
    }
  }

  aa_free_vars(CandAR);
  aa_free_vars(CandMA);

  if(BestP < 0 || BestQ < 0) {
    aa_free_vars(BestAR);
    aa_free_vars(BestMA);
    aa_free_vars(Diff);
    return 0;
  }

  var* Errors = aa_alloc_vars(NDiff);
  if(Errors == 0) {
    aa_free_vars(BestAR);
    aa_free_vars(BestMA);
    aa_free_vars(Diff);
    return 0;
  }

  // Reconstruct residual sequence with the best ARMA model.
  int t;
  for(t = 1; t < NDiff; t++) {
    var Pred = BestC;

    for(i = 0; i < BestP; i++) {
      if(t - 1 - i >= 0)
        Pred += BestAR[i] * Diff[t - 1 - i];
    }

    for(j = 0; j < BestQ; j++) {
      if(t - 1 - j >= 0)
        Pred += BestMA[j] * Errors[t - 1 - j];
    }

    Errors[t] = Diff[t] - Pred;
  }

  // One-step forecast on the differenced scale.
  var ForecastDiff = BestC;
  int LastT = NDiff - 1;

  for(i = 0; i < BestP; i++) {
    if(LastT - i >= 0)
      ForecastDiff += BestAR[i] * Diff[LastT - i];
  }

  for(j = 0; j < BestQ; j++) {
    if(LastT - j >= 1)
      ForecastDiff += BestMA[j] * Errors[LastT - j];
  }

  // Integrate forecast back to price space.
  var LastPrice = Close[0];
  var ForecastPrice = LastPrice;

  if(D == 0)
    ForecastPrice = ForecastDiff;
  else if(D == 1)
    ForecastPrice = LastPrice + ForecastDiff;
  else if(D == 2)
    ForecastPrice = LastPrice + ForecastDiff + (LastPrice - Close[1]);

  ForecastPrice = aa_round_to_tick_size(ForecastPrice,TickSize);

  Result->p = BestP;
  Result->d = D;
  Result->q = BestQ;
  Result->c = BestC;
  Result->sse = BestSSE;
  Result->aicc = BestAICc;
  Result->forecast = ForecastPrice;
  Result->ar = BestAR;
  Result->ma = BestMA;

  aa_free_vars(Errors);
  aa_free_vars(Diff);

  return 1;
}

// ---------- Standalone Zorro script demo ----------
// For strategy usage, replace this main() with your run() function.
// Zorro series are newest-first, so vars Close = series(priceClose()) matches the input convention.

void main()
{
  var Close[30] =
  {
    101.25, 101.10, 100.95, 100.70, 100.85,
    100.40, 100.20, 100.05, 99.90, 100.10,
    99.70, 99.55, 99.30, 99.50, 99.10,
    98.95, 99.20, 98.80, 98.60, 98.75,
    98.40, 98.20, 98.00, 97.85, 97.65,
    97.90, 97.55, 97.35, 97.20, 97.05
  };

  int N = 30;
  var TickSize = 0.01;
  int MaxP = 5;
  int MaxQ = 5;

  AUTO_ARIMA_RESULT R;
  init_auto_arima_result(&R);

  if(!auto_arima_forecast(Close,N,TickSize,MaxP,MaxQ,&R)) {
    printf("\nAutoARIMA forecast failed.");
    return;
  }

  printf("\nBest ARIMA(%i,%i,%i) AICc=%.10f SSE=%.10f",
    R.p,R.d,R.q,R.aicc,R.sse);

  printf("\nLast price: %.10f -> Forecast: %.10f",
    Close[0],R.forecast);

  printf("\nConstant c: %.10f",R.c);

  if(R.p > 0) {
    printf("\nAR coefficients:");
    int i;
    for(i = 0; i < R.p; i++)
      printf(" %.10f",R.ar[i]);
  }

  if(R.q > 0) {
    printf("\nMA coefficients:");
    int j;
    for(j = 0; j < R.q; j++)
      printf(" %.10f",R.ma[j]);
  }

  free_auto_arima_result(&R);
}

/*
// Example inside a Zorro strategy:

function run()
{
  BarPeriod = 1440;
  LookBack = 200;

  vars Close = series(priceClose());

  if(is(LOOKBACK))
    return;

  AUTO_ARIMA_RESULT R;
  init_auto_arima_result(&R);

  if(auto_arima_forecast(Close,100,PIP,5,5,&R)) {
    plot("ARIMA forecast",R.forecast,MAIN,BLUE);
    printf("\nARIMA(%i,%i,%i) forecast %.5f",R.p,R.d,R.q,R.forecast);
    free_auto_arima_result(&R);
  }
}
*/