#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

#define NDAYS 30
#define C0MIN 1000
#define MAXDAYS 512
#define MOVAVGDAYS 5
#define PLUSINFTY 1e300
#define WIDTH 15000
#define HEIGHT 10000

int readdata ( char reg[], char regname[], int c[], int C[], double ma[], double MA[], int madays, int maoffset, char date[][16] )
{
   FILE *fp;
   char fname[16], d[16];
   int n, nc, t;
   double manc;

   sprintf(fname, "db/%s.dat", reg);
   fp = fopen(fname, "r");
   if (fp == NULL) {
      fprintf(stderr, "*** Unable to read database \"%s\"\n", fname);
      exit(1);
   }
   fgets(regname, 128, fp); regname[strlen(regname)-1] = '\0';
   printf("+++ Reading data for %s\n", regname);
   n = -1;
   while (1) {
      fscanf(fp, "%s%d", d, &nc);
      if (feof(fp)) break;
      ++n;
      c[n] = nc; C[n] = C[n-1] + nc; strcpy(date[n], d);
   }
   fclose(fp);
   for (t=0; t<=maoffset; ++t)
      ma[t] = (double)C[t+maoffset] / (double)madays;
   manc = (double)C[madays-1];
   for (t=madays; t<=n; ++t) {
      manc -= c[t-madays];
      manc += c[t];
      ma[t-maoffset] = (double)manc / (double)madays;
   }
   MA[0] = ma[0];
   for (t=1; t<=n-maoffset; ++t) MA[t] = MA[t-1] + ma[t];
   return n;
}

void printdata ( char reg[], char regname[], int c[], int C[], double ma[], double MA[], int maoffset, char date[][16], int n )
{
   int i;

   printf("+++ Data for %s (%s)\n", reg, regname);
   for (i=0; i<=n-maoffset; ++i)
      printf("    %s %8d %10d   %8.1lf   %8.1lf\n", date[i], c[i], C[i], ma[i], MA[i]);
   for (; i<=n; ++i)
      printf("    %s %8d %10d         --         --\n", date[i], c[i], C[i]);
}

double f ( int t, double K, double r, double C0 )
{
   return K / (1 + ((K-C0)/C0)*exp(-r*t));
}

double partialK ( int t, double K, double r, double C0 )
{
   double ntr, dtr;
   ntr = 1 - exp(-r*t);
   dtr = 1 + ((K-C0)/C0)*exp(-r*t);
   dtr = dtr * dtr;
   return ntr / dtr;
}

double partialr ( int t, double K, double r, double C0 )
{
   double ntr, dtr;
   ntr = K * t * ((K-C0)/C0) * exp(-r*t);
   dtr = 1 + ((K-C0)/C0)*exp(-r*t);
   dtr = dtr * dtr;
   return ntr / dtr;
}

void estimateKr ( double C[], int n, int madays, int maoffset, double *Kp, double *rp )
{
   double K, r, delK, delr;
   double delfdelK[MAXDAYS], delfdelr[MAXDAYS], delC[MAXDAYS];
   double a1, a2, a3, a4, b1, b2;
   int i, t;

   K = 2 * C[n];
   r = log((2*C[n]-C[0])/C[0])/(double)n;
   // printf("    Iteration %3d: K = %lf, r = %lf\n", 0, K, r);

   i = 0;
   while (1) {
      for (t=0; t<=n; ++t) {
         delfdelK[t] = partialK(t,K,r,C[0]);
         delfdelr[t] = partialr(t,K,r,C[0]);
         delC[t] = C[t] - f(t,K,r,C[0]);
      }

      a1 = a2 = a3 = a4 = 0;
      b1 = b2 = 0;
      for (t=0; t<=n; ++t) {
         a1 += delfdelK[t] * delfdelK[t];
         a2 += delfdelK[t] * delfdelr[t];
         a3 += delfdelr[t] * delfdelK[t];
         a4 += delfdelr[t] * delfdelr[t];

         b1 += delfdelK[t] * delC[t];
         b2 += delfdelr[t] * delC[t];
      }

      delK = (a4 * b1 - a2 * b2) / (a1 * a4 - a2 * a3);
      delr = (a1 * b2 - a3 * b1) / (a1 * a4 - a2 * a3);

      K += delK;
      r += delr;

      ++i;
      // printf("    Iteration %3d: K = %lf, r = %lf\n", i, K, r);
      if ((abs(delK) < 0.5) && (abs(delr) < 0.0001)) break;
      if (isnan(K) || isnan(r) || (i > 100)) break;
   }

   *Kp = K; *rp = r;
}

void nextdate ( char olddate[], char newdate[] )
{
   int dd, mm, yyyy;
   char date[32], *d, *D;

   strcpy(date,olddate);
   D = date;
   d = strchr(D, '-'); *d = '\0'; dd = atoi(D); D = d+1;
   d = strchr(D, '-'); *d = '\0'; mm = atoi(D); D = d+1;
   yyyy = atoi(D); D = d+1;
   if (mm == 2) {
     if (yyyy == 2020) {
        if (dd == 29) { mm = 3; dd = 1; } else ++dd;
     } else {
        if (dd == 28) { mm = 3; dd = 1; } else ++dd;
     }
   } else if ((mm == 4) || (mm == 6) || (mm == 9) || (mm == 11)) {
     if (dd == 30) { ++mm; dd = 1; } else ++dd;
   } else {
     if (dd == 31) {
        if (mm == 12) { ++yyyy; mm = 1; dd = 1; }
        else { ++mm; dd = 1; }
     } else ++dd;
   }
   sprintf(newdate, "%02d-%02d-%d", dd, mm, yyyy);
}

void printdate ( FILE *fp, int t, double xres, char date[] )
{
   int dd, mm;
   char M[16];

   date[2] = '\0';
   dd = atoi(date);
   date[5] = '\0';
   mm = atoi(date+3);
   date[2] = date[5] = '-';
   if (dd == 1) {
      switch (mm) {
         case 1: strcpy(M, "Jan"); break;
         case 2: strcpy(M, "Feb"); break;
         case 3: strcpy(M, "Mar"); break;
         case 4: strcpy(M, "Apr"); break;
         case 5: strcpy(M, "May"); break;
         case 6: strcpy(M, "Jun"); break;
         case 7: strcpy(M, "Jul"); break;
         case 8: strcpy(M, "Aug"); break;
         case 9: strcpy(M, "Sep"); break;
         case 10: strcpy(M, "Oct"); break;
         case 11: strcpy(M, "Nov"); break;
         case 12: strcpy(M, "Dec"); break;
      }
      fprintf(fp, "2 1 0 3 0 7 30 -1 -1 0.000 0 0 -1 0 0 2\n");
      fprintf(fp, "\t %d %d %d %d\n", (int)round(t * xres), HEIGHT-190, (int)round(t * xres), HEIGHT+210);
      fprintf(fp, "4 1 0 30 -1 2 20 0.0000 4 300 900 %d %d 1 %s\\001\n", (int)round(t * xres), HEIGHT + 600, M);
   }
}

void exportfig ( double c[], double C[], int n, int ndays, int madays, double K, double r, char reg[], char regname[], char d[], double E, char winstart[], char winend[] )
{
   FILE *fp;
   char fname[32], date[16], newdate[16], peakdate[16], enddate[16];
   double A, p, q, pred, ppred;
   int peakday, maxdays, maxdays2, maxcases, ymax, t, fontparam;
   double xres, yres;

   printf("\n+++ Exporting to fig file\n");
   #ifdef _BESTFIT
   sprintf(fname, "ma%d/%sBF.fig", madays, reg);
   #else
   sprintf(fname, "ma%d/%s%d.fig", madays, reg, ndays);
   #endif
   fp = (FILE *)fopen(fname, "w");
   if (fp == NULL) {
      fprintf(stderr, "*** Unable to write to the result file \"%s\"\n", fname);
      exit(2);
   }

   fprintf(fp, "#FIG 3.2  Produced by predictor\nLandscape\nCenter\nMetric\nA4\n100.00\nSingle\n-2\n1200 2\n");

   A = (K - C[n-ndays]) / C[n-ndays];
   maxdays = (n - ndays) + round(log(1000 * A) / r);
   printf("maxdays = %d\n", maxdays);
   maxdays2 = (maxdays >= n) ? maxdays : n;
   printf("maxdays2 = %d\n", maxdays2);
   xres = (double)WIDTH / (double)maxdays2;
   ymax = maxcases = round(r * K / 4);
   for (t=0; t<=n; ++t) if (round(c[t]) > ymax) ymax = round(c[t]);
   printf("maxcases = %d\n", maxcases);
   printf("ymax = %d\n", ymax);
   yres = (double)HEIGHT / (double)ymax;

   fprintf(fp, "2 1 0 4 0 7 50 -1 -1 0.000 0 0 -1 0 0 3\n\t 0 0 0 %d %d %d\n", HEIGHT+10, WIDTH+200, HEIGHT+10);
   for (t=0; t<=n; ++t) {
      fprintf(fp, "2 1 0 2 4 7 40 -1 -1 0.000 0 0 7 0 0 2\n");
      fprintf(fp, "\t %d %d %d %d\n", (int)round(xres * t), HEIGHT, (int)round(xres * t), HEIGHT - (int)round(c[t] * yres));
   }

   fprintf(fp, "2 1 0 6 1 7 30 -1 -1 0.000 0 0 -1 0 0 %d\n", 1 + maxdays2);
   ppred = -1; peakday = -1;
   q = f(-(n-ndays+1),K,r,C[n-ndays]);
   for (t=0; t<=maxdays2; ++t) {
      p = f(t-(n-ndays),K,r,C[n-ndays]);
      pred = p - q;
      q = p;
      fprintf(fp, "\t %d %d\n", (int)round(xres * t), HEIGHT - (int)round(pred * yres));
      if (peakday < 0) {
         if (ppred > pred) peakday = t-1;
         ppred = pred;
      }
   }

   fprintf(fp, "2 1 0 4 0 7 30 -1 -1 0.000 0 0 -1 0 0 2\n");
   fprintf(fp, "\t -200 %d 200 %d\n", HEIGHT - (int)round(maxcases * yres), HEIGHT - (int)round(maxcases * yres));

   if (maxcases >= 100000) fontparam = 1530;
   else if (maxcases >= 10000) fontparam = 1275;
   else if (maxcases >= 1000) fontparam = 1020;
   else fontparam = 765;


   fprintf(fp, "4 0 0 20 -1 2 40 0.0000 4 585 4875 800 1000 %s\\001\n", regname);
   fprintf(fp, "4 2 0 20 -1 2 30 0.0000 4 330 %d -300 %d %d\\001\n", fontparam, HEIGHT + 150 - (int)round(maxcases * yres), maxcases);
   fprintf(fp, "4 0 0 20 -1 2 30 0.0000 4 435 7000 800 1800 Using %d-day moving averages\\001\n", madays);
   #ifdef _BESTFIT
   fprintf(fp, "4 0 0 20 -1 2 30 0.0000 4 435 5640 800 2400 Based on a %d-day window\\001\n", ndays);
   #else
   fprintf(fp, "4 0 0 20 -1 2 30 0.0000 4 435 4515 800 2400 Based on last %d days\\001\n", ndays);
   #endif
   fprintf(fp, "4 0 0 20 -1 2 30 0.0000 4 360 4725 800 3000 (%s to %s)\\001\n", winstart, winend);
   fprintf(fp, "4 0 0 20 -1 2 20 0.0000 4 225 4515 %d 1000 Expected total cases = %d\\001\n", WIDTH-4000, (int)round(K));
   fprintf(fp, "4 0 0 20 -1 2 20 0.0000 4 225 4515 %d 1500 Current total = %d\\001\n", WIDTH-4000, (int)round(C[n]));
   fprintf(fp, "4 0 0 20 -1 2 20 0.0000 4 225 4515 %d 2000 Estimated infection rate = %.3lf\\001\n", WIDTH-4000, r);

   strcpy(date, d);
   for (t=0; t<=maxdays2; ++t) {
      if (t == peakday) strcpy(peakdate,date);
      if (t == maxdays) strcpy(enddate, date);
      printdate(fp,t,xres,date);
      nextdate(date, newdate);
      strcpy(date, newdate);
   }

   printf("Peak date: %s\n", peakdate);
   printf("End date: %s\n", enddate);

   fprintf(fp, "4 0 0 20 -1 2 20 0.0000 4 225 4515 %d 2500 Expected peak on %s\\001\n", WIDTH-4000, peakdate);
   fprintf(fp, "4 0 0 20 -1 2 20 0.0000 4 225 4515 %d 3000 Expected end on %s\\001\n", WIDTH-4000, enddate);
   fprintf(fp, "4 0 0 20 -1 2 20 0.0000 4 225 4515 %d 3500 RMS error = %.3lf\\001\n", WIDTH-4000, E);

   fclose(fp);
}

void exportfigc ( double C[], int n, int ndays, int madays, double K, double r, char reg[], char regname[], char d[], double E, char winstart[], char winend[] )
{
   FILE *fp;
   char fname[32], date[16], newdate[16], enddate[16];
   double A, P;
   int maxdays, maxdays2, t, fontparam;
   double xres, yres;

   printf("\n+++ Exporting cumulative data to fig file\n");
   sprintf(fname, "ma%d/%sCM.fig", madays, reg);
   fp = (FILE *)fopen(fname, "w");
   if (fp == NULL) {
      fprintf(stderr, "*** Unable to write to the result file \"%s\"\n", fname);
      exit(2);
   }

   fprintf(fp, "#FIG 3.2  Produced by predictor\nLandscape\nCenter\nMetric\nA4\n100.00\nSingle\n-2\n1200 2\n");

   A = (K - C[n-ndays]) / C[n-ndays];
   maxdays = (n - ndays) + round(log(1000 * A) / r);
   printf("maxdays = %d\n", maxdays);
   maxdays2 = (maxdays >= n) ? maxdays : n;
   printf("maxdays2 = %d\n", maxdays2);
   xres = (double)WIDTH / (double)maxdays2;
   yres = (double)HEIGHT / K;

   if (K >= 100000) fontparam = 1530;
   else if (K >= 10000) fontparam = 1275;
   else if (K >= 1000) fontparam = 1020;
   else fontparam = 765;

   fprintf(fp, "2 1 0 4 0 7 50 -1 -1 0.000 0 0 -1 0 0 3\n\t 0 0 0 %d %d %d\n", HEIGHT+10, WIDTH+200, HEIGHT+10);
   for (t=0; t<=n; ++t) {
      fprintf(fp, "2 1 0 3 4 7 30 -1 -1 0.000 0 0 7 0 0 2\n");
      fprintf(fp, "\t %d %d %d %d\n", (int)round(xres * t), HEIGHT, (int)round(xres * t), HEIGHT - (int)round(C[t] * yres));
   }

   fprintf(fp, "2 1 0 6 1 7 40 -1 -1 0.000 0 0 -1 0 0 %d\n", 1 + maxdays2);
   for (t=0; t<=maxdays2; ++t) {
      P = f(t-(n-ndays),K,r,C[n-ndays]);
      fprintf(fp, "\t %d %d\n", (int)round(xres * t), HEIGHT - (int)round(P * yres));
   }

   fprintf(fp, "2 1 0 3 0 7 30 -1 -1 0.000 0 0 -1 0 0 2\n");
   fprintf(fp, "\t -200 %d 200 %d\n", HEIGHT - (int)round(K * yres), HEIGHT - (int)round(K * yres));
   fprintf(fp, "4 1 0 20 -1 2 30 0.0000 4 330 %d 0 %d %d\\001\n", fontparam, HEIGHT - 200 - (int)round(K * yres), (int)round(K));

   fprintf(fp, "2 1 0 3 0 7 30 -1 -1 0.000 0 0 -1 0 0 2\n");
   fprintf(fp, "\t -200 %d 200 %d\n", HEIGHT - (int)round(C[n] * yres), HEIGHT - (int)round(C[n] * yres));
   // fprintf(fp, "\t 0 %d %d %d\n", HEIGHT - (int)round(C[n] * yres), (int)round(n * xres), HEIGHT - (int)round(C[n] * yres));
   fprintf(fp, "4 2 0 20 -1 2 30 0.0000 4 330 %d -300 %d %d\\001\n", fontparam, HEIGHT + 150 - (int)round(C[n] * yres), (int)round(C[n]));

   fprintf(fp, "4 0 0 20 -1 2 40 0.0000 4 585 4875 800 1000 %s\\001\n", regname);
   fprintf(fp, "4 0 0 20 -1 2 30 0.0000 4 435 7000 800 1800 Using %d-day moving averages\\001\n", madays);
   fprintf(fp, "4 0 0 20 -1 2 30 0.0000 4 375 5640 800 2350 Based on a %d-day window\\001\n", ndays);
   fprintf(fp, "4 0 0 20 -1 2 30 0.0000 4 360 4725 800 2900 (%s to %s)\\001\n", winstart, winend);

   strcpy(date, d);
   for (t=0; t<=maxdays2; ++t) {
      if (t == maxdays) strcpy(enddate, date);
      printdate(fp,t,xres,date);
      nextdate(date, newdate);
      strcpy(date, newdate);
   }

   printf("End date: %s\n", enddate);

   fprintf(fp, "4 0 0 20 -1 2 20 0.0000 4 225 4515 800 3600 Estimated infection rate = %.3lf\\001\n", r);
   fprintf(fp, "4 0 0 20 -1 2 20 0.0000 4 300 3960 800 4100 Expected end on %s\\001\n", enddate);
   fprintf(fp, "4 0 0 20 -1 2 20 0.0000 4 225 3055 800 4600 RMS error = %.3lf\\001\n", E);

   fclose(fp);
}

double RMSerror ( double C[], int n, int ndays, double K, double r )
{
   int t;
   double p, E;

   E = 0;
   for (t=0; t<=n; ++t) {
      p = f(t-(n-ndays),K,r,C[n-ndays]);
      E += (C[t] - p) * (C[t] - p);
   }
   E /= n+1;
   E = sqrt(E);
   return E;
}

int bestfit ( double C[], int n, int madays, int maoffset, double *Ep, char date[][16], char winstart[], char winend[], double *Kp, double *rp )
{
   double K, r, E, Emin;
   int m, ndays, best;

   Emin = 1e300; best = -1;
   for (m=30; m<=n; ++m) {
      for (ndays=30; ndays<=m; ++ndays) {
         if (m < ndays) break;
         estimateKr(C+m-ndays,ndays,madays,maoffset,&K,&r);
         if (isnan(r) || isnan(K) || (r < 0) || (r > 1) || (K <= C[n]) || (K > 5 * C[n])) continue;
         E = RMSerror(C,n,ndays,K,r);
         // printf("N = %d, RMS error = %lf\n", ndays, E);
         if (E < Emin) {
            Emin = E;
            best = ndays;
            strcpy(winstart, date[m-ndays]);
            strcpy(winend, date[m]);
            *Kp = K;
            *rp = r;
         }
      }
   }
   if (best != -1) {
      printf("Best value of N = %d\n", best);
      printf("Start date: %s\n", winstart);
      printf("End date: %s\n", winend);
      if (Ep != NULL) *Ep = Emin;
   }
   return best;
}

int main ( int argc, char *argv[] )
{
   int n, ndays, madays, maoffset, c[MAXDAYS], C[MAXDAYS];
   char reg[16], regname[128], date[MAXDAYS][16], winstart[16], winend[16];
   double K, r, E, ma[MAXDAYS], MA[MAXDAYS];

   if (argc == 1) strcpy(reg, "IN"); else strcpy(reg, argv[1]);
   #ifdef _BESTFIT
   madays = (argc >= 3) ? atoi(argv[2]) : MOVAVGDAYS;
   #else
   ndays = (argc >= 3) ? atoi(argv[2]) : NDAYS;
   madays = (argc >= 4) ? atoi(argv[3]) : MOVAVGDAYS;
   #endif
   maoffset = madays / 2;

   n = readdata(reg, regname, c, C, ma, MA, madays, maoffset, date);
   printdata(reg, regname, c, C, ma, MA, maoffset, date, n);

   #ifdef _BESTFIT
   if (n < 30+maoffset) {
      fprintf(stderr, "*** Insufficient data\n");
      exit(3);
   }
   #else
   if (n < ndays+maoffset) {
      fprintf(stderr, "*** Insufficient data\n");
      exit(3);
   }
   #endif

   #ifdef _BESTFIT
   ndays = bestfit(MA,n-maoffset,madays,maoffset,&E,date,winstart,winend,&K,&r);
   #else
   estimateKr(MA+n-ndays-maoffset,ndays,madays,maoffset,&K,&r);
   strcpy(winstart,date[n-ndays-maoffset]);
   strcpy(winend,date[n-maoffset]);
   #endif
   if (isnan(r) || isnan(K)) {
      fprintf(stderr, "*** Estimation failed\n");
      exit(4);
   }
   printf("K = %d, C[n] = %d\n", (int)K, (int)round(MA[n-maoffset]));
   if ((r < 0) || (r > 1) || (K <= MA[n-maoffset]) || (K > 5 * MA[n-maoffset])) {
      fprintf(stderr, "*** Estimation failed\n");
      exit(5);
   }

   #ifndef _BESTFIT
   E = RMSerror(MA,n-maoffset,ndays,K,r);
   #endif

   printf("\nK = %lf, r = %.3lf, E = %.3lf\n", K, r, E);

   exportfig(ma,MA,n-maoffset,ndays,madays,K,r,reg,regname,date[0],E,winstart,winend);

   #ifdef _BESTFIT
   exportfigc(MA,n-maoffset,ndays,madays,K,r,reg,regname,date[0],E,winstart,winend);
   #endif

   exit(0);
}