#include "macrovsa.hpp"
#include "regex.hpp"
#include "file.hpp"
#include "time.hpp"
#include "stats.hpp"
#include <functional>

using namespace macrovsa;

/* Experiment corresponding to the (Viéville & Mercier 2024) draft
 */
int main()
{
  // Verification of symbolic derivations and reductions, dimension = 1000
  {
    Symbol::setDimension(1024);
    // Returns parsing and expression reduction
    auto get = [] (String s0) {
      // Escapes latex meta chars
      auto s2l = [] (String s)
      {
        return aidesys::regexReplace(aidesys::regexReplace(s, "([\\{}_])", "\\$1"), "~", "${}^\\sim$");
      };
      const Symbol& s1 = Symbol::fromJSON(s0), & s2 = algo::reduce(s1);
      return "{\\bf [I]: } {\\tt " + s2l(s0) + " } \\\\ \\hline {\\tt {\\bf [P]: } " + s2l(s1.asString()) + " } \\\\ \\hline {\\bf [R]: } {\\tt " + s2l(s2.asString()) + " }\\\\\\hline\\hline\n";
    };
    aidesys::save("../public/macrovsa_experiments/symbolic_reduction.tex", "\\begin{tabular}{|l|}\\hline\n" +
                  // "{\\bf [I]}nput in weak json syntax \\\\ \\hline {\\bf [P]}arsed form \\\\ \\hline {\\bf [R]}educed form \\\\\\hline\\hline\n" +
                  get("[ {b y: c x: [a b]} {b y: c x: [b a]}]") +
                  get("[ [{name: a tau: 0.5 sigma: 0.1} [ ]] {name: a tau: 0.5 sigma: 0.1}]") +
                  get("{b y: a x: {u y: a x: [ {name: c tau: 2 sigma: 0.1}]}") +
                  get("{b y: c x: {b y: c x: {b y: c x: {u y: c x: {u y: c x: {u y: c x: a}}}}}}") +
                  "\\end{tabular}\n");
  }
  // Verification of the noise and magnitude calculations w.r.t. the mesoscipoc level
  {
    Symbol y("y"), x("x");
    Binding Byx(y, x), By_Byx(y, Byx, false);
    Binding Byy(y, y), By_Byy(y, Byy, false);
    Bundling Bu;
    Bu.add(x), Bu.add(y);
    const Symbol *symbols[] = {
      &y, &x, &Byx, &By_Byx, &Byy, &By_Byy
    };
    const std::string strings[] = {
      "y", "x", "{b y: y x: x}", "{u y: y x: {b y: y x: x}}", "{b y: y x: y}", "{u y: y x: {b y: y x: y}}"
    };
    const unsigned int expressions_count = 8;
    const std::string labels[] = {
      "B_y x \\cdot x", "B_{y^\\sim} B_y x \\cdot x", "\\|B_y x\\|^2", "\\|B_{y^\\sim} B_y x\\|^2", "B_y y \\cdot y", "B_{y^\\sim} B_y y \\cdot y", "\\|B_y y\\|^2", "\\|B_{y^\\sim} B_y y\\|^2"
    };
    const unsigned int op1[] = {
      2, 3, 2, 3, 4, 5, 4, 5
    };
    const unsigned int op2[] = {
      0, 0, 2, 3, 1, 1, 4, 5
    };
    const unsigned int dimensions[] = { 10 * 10, 20 * 20, 32 * 32, 50 * 50, 64 * 64, 100 * 100, 1000 * 1000, 5000 * 5000};
    const unsigned int dimensions_count = 8;
    // Mesoscopic computations <label, dimension> => <data[N], stat>;
    const unsigned int N = 100;
    std::map < std::pair < unsigned int, unsigned int >, std::vector < double >> data;
    std::map < std::pair < unsigned int, unsigned int >, std::string > data_stats;
    std::map < std::pair < unsigned int, unsigned int >, Belief > beliefs;
    // Computation average durations
    double binding_computation_durations[dimensions_count];
    double bundling_computation_durations[dimensions_count];
    // Computation loops
    for(unsigned c = 0; c < dimensions_count; c++) {
      printf("macrovsa_experiments d=%d \n", dimensions[c]);
      binding_computation_durations[c] = 0;
      bundling_computation_durations[c] = 0;
      for(unsigned int n = 0; n < N; n++) {
	if(c > 5)
	  printf("{d : %d n: %d now: %.0f}\n", dimensions[c], n, 0.001*aidesys::now(false, false));
        Symbol::setDimension(dimensions[c]);
	x.getVector(), y.getVector();
       // Binding computation duration estimation
        {
          aidesys::now(false, false);
          Byx.getVector(), By_Byx.getVector(), Byy.getVector(), By_Byy.getVector();
          binding_computation_durations[c] += 0.25 * aidesys::now(false, true);
        }
         // Bundling computation duration estimation
        {
          aidesys::now(false, false);
	  Bu.getVector();
          bundling_computation_durations[c] += 0.5 * aidesys::now(false, true);
        }
         // Mesoscopic similarities computation
        for(unsigned int e = 0; e < expressions_count; e++) {
          double r = algo::msim(*symbols[op1[e]], *symbols[op2[e]]);
          data[std::pair < unsigned int, unsigned int > (e, dimensions[c])].push_back(r);
          // Macroscopic computation
          if(n == 0) {
            Belief b = algo::sim(strings[op1[e]], strings[op2[e]]);
            beliefs[std::pair < unsigned int, unsigned int > (e, dimensions[c])] = b;
          }
        }
      }
    }
    printf("macrovsa_experiments statistics\n");
    // Computes the related statistics
    std::map < unsigned int, unsigned int > models_counts;
    std::map < unsigned int, std::vector < double >> normal_divergences;
    std::map < std::string, std::string > normal_divergences_stats;
    for(unsigned c = 0; c < dimensions_count; c++) {
      binding_computation_durations[c] /= N;
      bundling_computation_durations[c] /= N;
      for(unsigned int e = 0; e < expressions_count; e++) {
        auto i = std::pair < unsigned int, unsigned int > (e, dimensions[c]);
        data_stats[i] = aidesys::getStat(data[i]);
        models_counts[(int) aidesys::getStatValue("best-model", data_stats[i])]++;
        normal_divergences[c].push_back(aidesys::getStatValue("normal-divergence", data_stats[i]));
      }
      normal_divergences_stats[aidesys::echo("%6d", dimensions[c])] = aidesys::getStat(normal_divergences.at(c));
    }
    printf("macrovsa_experiments output\n");
    // Output the data results for a given print and label
    {
      // Returns the concatenation of c times the string s
      auto repeatString = [](String s, unsigned int c) {
	std::string r;
	for (unsigned int i = 0; i < c ; i++)
	  r += s;
	return r;
      };
      std::string tabular_header = "\\begin{tabular}{|"+repeatString("l|", dimensions_count + 2) + "}\n\\hline\n";
      auto getLine = [labels,dimensions_count](std::function < std::string(unsigned int e, unsigned int c) > print, unsigned int e)
      {
	std::string r = "${\\bf " + labels[e] + "}$";
	for (unsigned int i = 0; i < dimensions_count ; i++)
	  r += "&" + print(e, i) ;
	return r + "\\\\\\hline\n";
      };
      // Returns mesoscopic means and stdev for a given label
      auto getM = [labels, data_stats, dimensions, getLine](unsigned int e)
      {
        // Returns stats means and stdev parameters
        auto print = [data_stats, dimensions](unsigned int e, unsigned int c)
        {
          String stat = data_stats.at(std::pair < unsigned int, unsigned int > (e, dimensions[c]));
          return aidesys::echo("%.2g±%.2g", aidesys::getStatValue("mean", stat), aidesys::getStatValue("stdev", stat));
        };
        return getLine(print, e);
      };
      // Noise distribution
      {
        printf("Noise-best-model: { ");
        for(auto it = models_counts.cbegin(); it != models_counts.cend(); it++) {
          unsigned int model_index = it->first;
          std::string model_name = model_index == (unsigned int) -1 ? "uniform" : model_index == 0 ? "normal" : aidesys::echo("gamma_%d", model_index);
          printf("%s: %.0f%% ", model_name.c_str(), (100.0 * it->second) / (dimensions_count * expressions_count));
        }
        printf("}\n");
        aidesys::plotStatBoxes("mesoscopic_noise_normal_divergence", normal_divergences_stats);
      }
      // Noise comparison output
      {
        // Returns belief for a given label
        auto getB = [labels, beliefs, dimensions, getLine](unsigned int e)
        {
          // Returns beliefs
          auto print = [beliefs, dimensions](unsigned int e, unsigned int c)
          {
	    return aidesys::echo("0±%.2g", beliefs.at(std::pair < unsigned int, unsigned int > (e, dimensions[c])).sigma);
          };
          return getLine(print, e);
        };
        auto getR = [labels, data_stats, beliefs, dimensions, getLine](unsigned int e)
        {
          // Returns beliefs
          auto print = [data_stats, beliefs, dimensions](unsigned int e, unsigned int c)
          {
            String stat = data_stats.at(std::pair < unsigned int, unsigned int > (e, dimensions[c]));
            double sigma_meso = aidesys::getStatValue("stdev", stat);
            double sigma_macro = beliefs.at(std::pair < unsigned int, unsigned int > (e, dimensions[c])).sigma;
            return aidesys::echo("%.2g", sigma_meso == 0 ? 0 : sigma_macro / sigma_meso);
          };
          return getLine(print, e);
        };
	// Returns the dimensions
	auto getDims = [dimensions, dimensions_count]()
	{
	  std::string r = "Dimension: ";
	  for (unsigned int i = 0; i < dimensions_count ; i++)
	    r += "&" + aidesys::echo("%d", dimensions[i]);
	  return r + "\\\\\\hline\n";
	};
	std::string dimp1 = aidesys::echo("%d", dimensions_count +1);
        aidesys::save("../public/macrovsa_experiments/meso_versus_macro_noise.tex",
                      tabular_header + getDims() +
                      "\\multicolumn{"+dimp1+"}{|l|}{Mesoscopic estimation over "+aidesys::echo("%d", N)+" samples}\\\\\\hline\n" +
                      getM(0) + getM(1) + getM(4) + getM(5) +
                      "\\multicolumn{"+dimp1+"}{|l|}{Macroscopic bias and standard-deviation}\\\\\\hline\n" +
                      getB(0) + getB(1) + getB(4) + getB(5) +
                      "\\multicolumn{"+dimp1+"}{|l|}{Standard-deviation macroscopic/microscopic ratio}\\\\\\hline\n" +
                      getR(0) + getR(1) + getR(4) + getR(5) +
                      "\\end{tabular}\n");
      }
      // Binding magnitude output
      {
        aidesys::save("../public/macrovsa_experiments/magnitudes.tex", tabular_header + getM(2) + getM(3) + getM(6) + getM(7) + "\\end{tabular}\n");
      }
      // Output computation time result
      {
        // Computation interpolation using maple
        {
          std::string X, Y, Xu, Yu;
          for(unsigned c = 0; c < dimensions_count; c++) {
            X += (c == 0 ? "" : ", ") + aidesys::echo("%d", dimensions[c]);
            Y += (c == 0 ? "" : ", ") + aidesys::echo("%g", binding_computation_durations[c]);
	    if (c < dimensions_count-1) {
	      Xu += (c == 0 ? "" : ", ") + aidesys::echo("%g", dimensions[c]);
	      Yu += (c == 0 ? "" : ", ") + aidesys::echo("%g", bundling_computation_durations[c]);
	    }
          }
          aidesys::save("durations_interpolation.mpl",
                        "# Generated by macrovsa_experiments.C do NOT edit\n"
                        "with(Statistics):with(plots):\n"
                        "X := Vector([" + X + "], datatype=float):\n"
                         "Y := Vector([" + Y + "], datatype=float):\n"
                         "Xu := Vector([" + Xu + "], datatype=float):\n"
                         "Yu := Vector([" + Yu + "], datatype=float):\n"
			"fit1 := Fit(a+b*d^c, X, Y, d, initialvalues = [a = 0.05, b = 0.00003, c = 1.35], output = [leastsquaresfunction, residualstandarddeviation]);\n"
			"f1 := unapply(fit1[1], d):\n"
			"fit2 := Fit(a+b*d^1.5, X, Y, d, output=[leastsquaresfunction, residualstandarddeviation], summarize=true);\n"
			"f2 := unapply(fit2[1], d):\n"
			"plotsetup(jpeg, plotoutput=\"../public/macrovsa_experiments/durations_interpolation.jpg\", plotoptions=\"width=600,height=600\"):\n"
			"display(plot(f1(d), d=100..10000), pointplot(X, Y), title=\"binding computation time in msec\");\n"
			"fit3 := Fit(a+b*d, Xu, Yu, d, initialvalues = [a = 0.00, b = 0.00003], output = [leastsquaresfunction, residualstandarddeviation]);\n"
			"f3 := unapply(fit3[1], d):\n"
			"plotsetup(jpeg, plotoutput=\"../public/macrovsa_experiments/durations_interpolation_2.jpg\", plotoptions=\"width=600,height=600\"):\n"
			"display(plot(f3(d), d=100..10000), pointplot(Xu, Yu), title=\"bundling computation time in msec\");\n"
			"quit:\n");
          // T_binding = 0.048 + 0.28 d^(1.35) / 10000 +- 0.05
        }
      }
    }
  }
  // Generates a pdf for control
  {
    aidesys::save("macrovsa_experiments.tex",
                  "% Generated by macrovsa_experiments.C do NOT edit\n"
                  "\\section*{Macrovsa experiments results}\n\n"
                  "\\subsection*{Sympolic reduction}\n\n"
                  "\\input{../public/macrovsa_experiments/symbolic_reduction.tex}\n\n"
                  "\\subsection*{Mesoscopic versus macroscopic noise comparison}\n\n"
                  "\\input{../public/macrovsa_experiments/meso_versus_macro_noise.tex}\n\n"
                  "\\subsection*{Binding magnitudes}\n\n"
                  "\\input{../public/macrovsa_experiments/magnitudes.tex}\n\n"
                  "\\subsection*{Binding computation durations}\n\n"
                  "\\includegraphics[width=0.5\\textwidth]{../public/macrovsa_experiments/durations_interpolation.jpg}\n"
                 "\\subsection*{Bundling computation durations}\n\n"
                  "\\includegraphics[width=0.5\\textwidth]{../public/macrovsa_experiments/durations_interpolation_2.jpg}\n");
  }
  printf("macrovsa_experiments done!\n");
  return 0;
}