/* Istruzioni:
////// .x nn.C
////// (per eseguire in una volta sola la conversione dei dati (da ascii a root-tree), la visualizzazione delle variabili, e l'addestramento e utilizzo della rete.)

////// oppure

////// .L nn.C                 (per caricare in memoria la macro; farlo tutte e sole le volte che la si e` modificata)
////// converti("pippo",K)        (per convertire pippo.txt in pippo.root, su K giorni)
////// visualizza("pippo")     (per visualizzare le variabili ora contenute in pippo.root)
////// reteneurale("pippo",K,N)  (per addestrare la rete sui dati di pippo.root con N cicli di addestramento, con l'algoritmo K, e salvarla nel file pesi.txt)
////// usarete("pesi.txt","pippo")     (per utilizzare la rete gia' addestrata di cui avete salvato i pesi in pesi.txt, sul sample pippo.root)

*/

#if !defined(__CINT__) || defined(__MAKECINT__)

#include <assert.h>
#include <fstream>
#include <iostream>
#include <iomanip>
#include <string>
#include <deque>
#include <numeric>

#include "TSystem.h"
#include "TCanvas.h"
#include "TMath.h"
#include "TFile.h"
#include "TTree.h"
#include "TEventList.h"
#include "TMultiLayerPerceptron.h"
#include "TMLPAnalyzer.h"

#endif
//#include "Riostream.h"

void nn()
{
  gROOT->Reset(); // fa il reset di tutte le variabili
  converti("dati",365); // converte il file testo in un file di ROOT; la variabile numerica e' il numero di giorni descritti dal file
  // se avete piu' di un file di input, ripetete il rigo sopra piu' volte mettendo ogni volta il nome del file (senza ".txt")
  visualizza("dati"); // grafico delle distribuzioni di tutte le variabili
  reteneurale("dati",1,4,2,100); // addestra la rete con il primo algoritmo, 4 neuroni nel primo strato, 2 nel secondo strato, e dopo l'addestramento salva i pesi in un file di testo
  usarete("pesi_Stochastic_4_2_100","dati"); // ricrea la rete con i pesi dati dall'addestramento, e la mette a confronto con la realta'
  /* Algoritmi a disposizione:
  1) Stochastic
  2) Batch
  3) SteepestDescent
  4) RibierePolak
  5) FletcherReeves
  6) BFGS
  */
}

void converti(Char_t name[40],Int_t ndays)
{
  // definizione dei nomi dei file di input e output
  Char_t input[40];
  Char_t output[40];
  sprintf(input,"%s.txt",name);
  sprintf(output,"%s.root",name);

  // apertura file di input
  ifstream in;
  in.open(input);
  // apertura file di output
  TFile *f = new TFile(output,"RECREATE");

  // lettura del contenuto del file di input
  Double_t x[366][366];
  Int_t nlines = 0;
  while (1) {
    for (Int_t i=0; i<ndays; i++) {
      in >> x[nlines][i];
    }
    if (!in.good()) break;
    nlines++;
  }

  cout << "number of lines: " << nlines << endl;
  for (Int_t j=0; j<nlines; j++) {
      cout << endl << "line n." << j << ": ";
      for (Int_t i=0; i<ndays; i++) {
	in >> x[j][i];
	cout << " " << x[j][i];
      }
      cout << endl << endl;
  }

  // ora dobbiamo salvare tutto in un "tree" di root

  // in quanto segue assumo che il file di input contenga 4 righe, corrispondenti rispettivamente a pressione, temperatura, umidita' e precipitazione; qualora il contenuto del file sia diverso, modificare di conseguenza
  Double_t pressione;
  Double_t temperatura;
  Double_t umidita;
  Double_t precipitazione; 

  TTree *tree = new TTree("T","il nostro tree");
  tree->Branch("pressione",&pressione,"pressione/D");
  tree->Branch("temperatura",&temperatura,"temperatura/D");
  tree->Branch("umidita",&umidita,"umidita/D");
  tree->Branch("precipitazione",&precipitazione,"precipitazione/D");

  for (Int_t i=0; i<ndays; i++) {
    // qui potete operare per regolarizzare le variabili (ad es. togliere fattori costanti per imporre media vicina a 0, o dividere per un fattore costante per imporre range tra -1 e 1):
    pressione=x[0][i];
    temperatura=x[1][i];
    umidita=x[2][i];
    precipitazione=x[3][i];
    if (pressione != -9999 && temperatura != -9999 && umidita != -9999 && precipitazione != -9999) {
      if (i<50) cout << endl << "day n." << i << ":  P=" << pressione << ", T=" << temperatura << ", U=" << umidita << ", Pr=" << precipitazione << endl; 
      tree->Fill();
    }
  }

  // operazioni conclusive: chiusura file di input, scrittura degli output
  in.close();
  tree->Write();
  f->Write();
  f->Close();
}

void visualizza(Char_t name[40])
{

  // definizione dei nomi dei file
  Char_t rootfile[40];
  sprintf(rootfile,"%s.root",name);

  // apertura del file e del suo tree
  TFile f(rootfile);  
  TTree *tree = (TTree*)gDirectory->Get("T"); // metto "T" perche' e' il nome del tree contenuto nel file

  c1 = new TCanvas("c1","Variabili",200,10,1000,900);
  c1->SetFillColor(18);
  pad1 = new TPad("pad1","Pressione",0.05,0.50,0.45,0.95,21);
  pad2 = new TPad("pad2","Temperatura",0.05,0.05,0.45,0.45,21);
  pad3 = new TPad("pad3","Umidita`",0.50,0.50,0.95,0.95,21);
  pad4 = new TPad("pad4","Precipitazione",0.50,0.05,0.95,0.45,21);
  pad1->Draw();
  pad2->Draw();
  pad3->Draw();
  pad4->Draw();
  pad1->cd();
  tree->Draw("pressione");
  pad2->cd();
  tree->Draw("temperatura");
  pad3->cd();
  tree->Draw("umidita");
  pad4->cd();
  tree->Draw("precipitazione");

}

void reteneurale(Char_t name[40], Int_t nAlgo, Int_t neuro1, Int_t neuro2, Int_t ntrain)
{
  if (!gROOT->GetClass("TMultiLayerPerceptron")) {
    gSystem->Load("libMLP");
    cout << endl << "E' stata caricata la libreria ROOT per le reti neurali." << endl << endl;
  }

  // qui mettere le stesse variabili messe sopra
  Double_t pressione;
  Double_t temperatura;
  Double_t umidita;
  Double_t precipitazione;
  
  // definizione dei nomi dei file
  Char_t rootfile[40];
  sprintf(rootfile,"%s.root",name);

  // apertura del file e del suo tree
  TFile f(rootfile); // questo sara' il file usato per il training della rete
  TTree *tree = (TTree*)gDirectory->Get("T"); // metto "T" perche' e' il nome del tree contenuto nel file
  tree->SetBranchAddress("pressione",&pressione);
  tree->SetBranchAddress("temperatura",&temperatura);
  tree->SetBranchAddress("umidita",&umidita);
  tree->SetBranchAddress("precipitazione",&precipitazione);

  // creiamo un oggetto della classe "rete neurale"
  // istruzioni qui: http://root.cern.ch/root/html310/TMultiLayerPerceptron.html
  TString struttura("pressione,temperatura,umidita"); // inputs
  struttura+=":";
  struttura+=neuro1;
  if (neuro2>0) { // due strati nascosti
    struttura+=":";
    struttura+=neuro2;
  }
  struttura+=":";
  struttura+="@precipitazione"; //output; la @ significa che la rete lo considera normalizzato

  TMultiLayerPerceptron *mlp = new TMultiLayerPerceptron(struttura,tree);

  TEventList train;
  TEventList test;
  for (int i = 0; i < tree->GetEntries(); i++) {
    if (i < 150)
      train.Enter(i);
    else
      test.Enter(i);
  }
  train.Print();
  test.Print();

  mlp->SetTrainingDataSet(&train);
  mlp->SetTestDataSet(&test);

  ///// algoritmo di addestramento:
  if (nAlgo == 1) mlp->SetLearningMethod(TMultiLayerPerceptron::kStochastic);
  if (nAlgo == 2) mlp->SetLearningMethod(TMultiLayerPerceptron::kBatch);
  if (nAlgo == 3) mlp->SetLearningMethod(TMultiLayerPerceptron::kSteepestDescent);
  if (nAlgo == 4) mlp->SetLearningMethod(TMultiLayerPerceptron::kRibierePolak);
  if (nAlgo == 5) mlp->SetLearningMethod(TMultiLayerPerceptron::kFletcherReeves);
  if (nAlgo == 6) mlp->SetLearningMethod(TMultiLayerPerceptron::kBFGS);
  if (nAlgo<1 || nAlgo>6) cout << "Algoritmo non valido." << endl;
  TString algo;
  if (nAlgo == 1) algo="Stochastic";
  if (nAlgo == 2) algo="Batch";
  if (nAlgo == 3) algo="SteepestDescent";
  if (nAlgo == 4) algo="RibierePolak";
  if (nAlgo == 5) algo="FletcherReeves";
  if (nAlgo == 6) algo="BFGS";
  cout << "Algoritmo scelto: " << algo << endl;

  // addestramento:
  mlp->Randomize(); // randomizza i pesi tra neuroni prima dell'addestramento
  mlp->Train(ntrain, "text,graph,update=10");

  // Use TMLPAnalyzer to see what it looks for
  TCanvas* mlpa_canvas = new TCanvas("mlpa_canvas",
				     "Network analysis");
  mlpa_canvas->SetWindowPosition(160,80);
  TMLPAnalyzer ana(mlp);
  // Initialisation
  ana.GatherInformations();
  // output to the console
  ana.CheckNetwork();
  mlp->Draw();

  //  mlp->DrawResult(0, "test"); // perche' non funziona???

  // salva i pesi di questa rete per poterli usare dopo:
  TString nomepesi("pesi_");
  nomepesi+=algo;
  nomepesi+="_";
  nomepesi+=neuro1;
  nomepesi+="_";
  nomepesi+=neuro2;
  nomepesi+="_";
  nomepesi+=ntrain;
  nomepesi+=".txt";
  mlp->DumpWeights(nomepesi);
  cout << "File dei pesi salvato col nome " << nomepesi << endl;

  // per scoprire come fissare l'algoritmo di addestramento, e in generale come usare la rete, leggere attentamente http://root.cern.ch/root/html310/TMultiLayerPerceptron.html
}
 
void usarete(Char_t filepesi[40],Char_t name[40])
{
  if (!gROOT->GetClass("TMultiLayerPerceptron")) {
    gSystem->Load("libMLP");
    cout << endl << "E' stata caricata la libreria ROOT per le reti neurali." << endl << endl;
  }

  // qui mettere le stesse variabili messe sopra
  Double_t pressione;
  Double_t temperatura;
  Double_t umidita;
  Double_t precipitazione;
  
  // definizione dei nomi dei file
  Char_t rootfile[40];
  sprintf(rootfile,"%s.root",name);

  // apertura del file e del suo tree
  TFile f(rootfile); // questo sara' il file usato per il training della rete
  TTree *tree = (TTree*)gDirectory->Get("T"); // metto "T" perche' e' il nome del tree contenuto nel file
  tree->SetBranchAddress("pressione",&pressione);
  tree->SetBranchAddress("temperatura",&temperatura);
  tree->SetBranchAddress("umidita",&umidita);
  tree->SetBranchAddress("precipitazione",&precipitazione);

  // definiamo la struttura della rete neurale (dev'essere identica a quella usata in reteneurale())
  TMultiLayerPerceptron *mlp =
    new TMultiLayerPerceptron("pressione,temperatura,umidita"
			      ":5:2:"
			      "@precipitazione",
			      tree);
  mlp->LoadWeights(filepesi);
  //  mlp->Draw(); // questo in realta' non ci serve

  // loop sul numero di entries nel sample
  Double_t inputs[3];
  Double_t output;
  Double_t target;
  TH1F *h1;
  TH1F *h2;
  h1 = new TH1F("h1","output",20,-1,10);
  h2 = new TH1F("h2","target",20,-1,10);
  Int_t pioveNN=0;
  Int_t pioveR=0;
  Int_t falsiPos=0;
  Int_t falsiNeg=0;
  Double_t soglia=1; // soglia sotto la quale non consideriamo piovosa la giornata
  cout << "Numero giorni totali: " << tree->GetEntries() << endl;
  for (int i = 0; i < tree->GetEntries(); i++) {
    tree->GetEntry(i);
    inputs[0]=pressione;
    inputs[1]=temperatura;
    inputs[2]=umidita;
    output=mlp->Evaluate(0,inputs);
    target=precipitazione;
    h1->Fill(target);
    h2->Fill(output);
    if (output > soglia) pioveNN++;
    if (target > soglia) pioveR++;
    if (output > soglia && target < soglia) falsiPos++;
    if (output < soglia && target > soglia) falsiNeg++;
  }

  // visualizza confronto tra output e target
  gStyle->SetOptStat(0000);
  TCanvas* c = new TCanvas("c","Output / Target comparison"); 
  c->SetWindowPosition(260,100);
  c->cd();
  h1->SetLineWidth(2);
  h2->SetLineWidth(2);
  h1->SetLineColor(kBlack);
  h2->SetLineColor(kRed);
  h2->DrawCopy();
  h1->DrawCopy("same");
  //legenda
  Double_t ymax=max(h1->GetMaximum(),h2->GetMaximum());
  TLegend *legend = new TLegend(7.5,0.75*ymax,9.5,0.95*ymax,"Precipitazioni","br");
  legend->AddEntry(h1,"previsione","l");
  legend->AddEntry(h2,"realta`","l");
  legend->Draw();

  cout << "Ha piovuto in " << pioveR << " giorni contro " << pioveNN << " previsti dalla rete." << endl << "Falsi positivi: " << falsiPos << ", falsi negativi: " << falsiNeg << endl;

}

Double_t max(Double_t a, Double_t b)
{
  Double_t res=a;
  if (b>a) res=b;

  return res;
}