/******************************************************************************
* macro to test quantile normalization                                        *
*                                                                             *
* note: this code is under the license GPL>=2.0 since it is taken from:       *
*       http://www.bioconductor.org/packages/2.4/bioc/html/xps.html           *
*                                                                             *
* Author: Christian Stratowa, Vienna, Austria.                                *
******************************************************************************/

///////////////////////////
//
// in new root session(s):
//     > .L macroQuantileTest.C 
//     > QuantileTest("TestAB.txt", "TestAB_Quantile.txt", 126) 
//
///////////////////////////

#ifndef __CINT__
#include <Riostream.h>
#endif

#include "TFile.h"
#include "TMath.h"
#include "TString.h"
#include "TKey.h"
#include "TTree.h"

TFile    *fTmpFile = 0;
Double_t *fMean1   = 0;
Int_t     fNMean1  = 0;
Int_t     fNData   = 0;

const Int_t kBufSize = 1024;


//______________________________________________________________________________
//void TStat::Rank(Int_t n, Double_t *arr, Int_t *index, Int_t *rank, Bool_t down)
void Rank(Int_t n, Double_t *arr, Int_t *index, Int_t *rank, Bool_t down)
{
   if (n <= 0) return;
   if (n == 1) {
      index[0] = 0;
      rank[0] = 0;
      return;
   }//if

   TMath::Sort(n,arr,index,down);

   for (Int_t i=0; i<n; i++) {
      rank[index[i]] = i;
   }//for_i
}//Rank

//______________________________________________________________________________
//void TStat::Rank(Int_t n, Double_t *arr, Double_t *rank)
void Rank(Int_t n, Double_t *arr, Double_t *rank)
{
   Int_t i = 0;
   Int_t j = 0;
   
   while (i < n) {
      j = i;

      while ((j < n-1) && (arr[j] == arr[j+1])) j++;

      if (i != j) {
         for (Int_t k=i; k<=j; k++) rank[k] = (i + j + 2) / 2.0;
      } else {
         rank[i] = i + 1;
      }//if

      i = j + 1;
   }//while
}//Rank

//______________________________________________________________________________
//Double_t TStat::Mean(Int_t n, const Double_t *arr)
Double_t Mean(Int_t n, const Double_t *arr)
{
   if (n == 1) return arr[0];

   Double_t mean = 0.0;
   for (Int_t i=0; i<n; i++) mean += arr[i];

   return mean/n;
}//Mean

//______________________________________________________________________________
//Int_t XQuantileNormalizer::AddArray(Int_t n, Double_t *x, Short_t *msk,
Int_t AddArray(Int_t n, Double_t *x, Short_t *msk, const char *name)
{
   TDirectory *savedir = gDirectory;
   Int_t       err     = 0;

// Create temporary root file to store trees
   if (!fTmpFile) {
      fTmpFile = new TFile("tmp_rkq.root", "RECREATE");
   }//if

// Fill tree(s) with sorted data array and store in fTmpFile
   fTmpFile->cd();
   TString   tname = TString(name) + ".rkt";  //rkt - rank together
   Double_t  sort  = 0;
   Float_t   rank  = 0;
   Int_t     id    = 0;
   Int_t     i     = 0;

// Create tree with data and rank branches
   TTree *tree = new TTree(tname, "temporary tree");
   tree->Branch("sortBr", &sort, "sort/D");
   tree->Branch("rankBr", &rank, "rank/F");

// Get size of arrays for PM and MM
   for (i=0; i<n; i++) {if (msk[i] == 1) id++;}

// Create arrays
   Double_t *arr  = new Double_t[id];
   Int_t    *idx  = new Int_t[id];
   Double_t *ord  = new Double_t[id];
   Int_t    *rnki = new Int_t[id];
   Double_t *rnkd = new Double_t[id];

// Fill arr
   id = 0;
   for (i=0; i<n; i++) {
      if (msk[i] == 1) arr[id++] = x[i];
   }//for_i

// Sort array arr
   Rank(id, arr, idx, rnki, kFALSE);
   // get rank rnkd analogoulsy to rank used for RMA in Bioconductor
   for (Int_t i=0; i<id; i++) ord[i] = arr[idx[i]];
   Rank(id, ord, rnkd);

// Fill tree with sorted array 
   for (i=0; i<id; i++) {
      sort = arr[idx[i]];
      rank = rnkd[rnki[i]];
      tree->Fill();
   }//for_i

   tree->Write();
   fNData++;

   if (tree) {tree->Delete(""); tree = 0;}
   if (rnkd) {delete [] rnkd;   rnkd = 0;}
   if (rnki) {delete [] rnki;   rnki = 0;}
   if (ord)  {delete [] ord;    ord  = 0;}
   if (idx)  {delete [] idx;    idx  = 0;}
   if (arr)  {delete [] arr;    arr  = 0;}

   savedir->cd();

   return err;
}//AddArray

//______________________________________________________________________________
//Double_t *XQuantileNormalizer::GetArray(Int_t n, Double_t *x, Short_t *msk, 
Double_t *GetArray(Int_t n, Double_t *x, Short_t *msk, const char *name)
{
   TDirectory *savedir = gDirectory;
   fTmpFile->cd();

// Get tree with name
   TString tname = TString(name) + ".rkt";
   Float_t rank  = 0;
   Int_t   frank = 0;
   TTree  *tree  = (TTree*)fTmpFile->Get(tname.Data()); 
   tree->SetBranchAddress("rankBr", &rank);

// Sort mean in order of original tree data 
   Int_t id = 0;
   for (Int_t i=0; i<n; i++) {
      if (msk[i] == 1) {
         tree->GetEntry(id++);
         frank = (Int_t)floor(rank);
         x[i]  = fMean1[frank - 1];
      } else {
         x[i] = 0;  //not necessary?
      }//if
   }//for_i

   tree->Delete(""); tree = 0;
   savedir->cd();

   return x;
}//GetArray

//______________________________________________________________________________
//Int_t XQuantileNormalizer::Calculate(Int_t n, Double_t *x, Double_t *y, Short_t *msk)
Int_t Calculate(Int_t n, Double_t *x, Double_t *y, Short_t *msk)
{
   TDirectory *savedir = gDirectory;
   fTmpFile->cd();

// Calculate mean
   TTree   **treek = new TTree*[fNData];
   Double_t *sortk = new Double_t[fNData];
   for (Int_t k=0; k<fNData; k++) sortk[k] = 0;

// Init trees
   Int_t idx = 0;
   TKey *key = 0;
   TIter next(fTmpFile->GetListOfKeys());
   while ((key = (TKey*)next())) {
      TTree *tmptree = (TTree*)fTmpFile->Get(key->GetName());
      treek[idx] = tmptree;
      treek[idx]->SetBranchAddress("sortBr", &sortk[idx]);
      idx++;
   }//while

// Create data array
   Double_t *arr = new Double_t[fNData];

// Create mean array
   fNMean1 = treek[0]->GetEntries();
   fMean1 = new Double_t[fNMean1];

// Calculate mean values of all tree entries
   for (Int_t i=0; i<fNMean1; i++) {
      // read entry i from trees and fill array
      for (Int_t k=0; k<fNData; k++) {
         treek[k]->GetEntry(i);
         arr[k] = sortk[k];
      }//for_k

      // calculate mean of array
      fMean1[i] = Mean(fNData, arr);
      y[i]      = fMean1[i];
   }//for_i

// Cleanup
   if (arr) {delete [] arr; arr = 0;}
   delete [] sortk;
   for (Int_t k=0; k<fNData; k++) {delete treek[k]; treek[k] = 0;}
   delete [] treek;

   savedir->cd();

   return 0;
}//Calculate

//______________________________________________________________________________
void QuantileTest(const char *infile, const char *outfile, Int_t numrows = 126)
{
   Int_t err = 0;

// Reset global variables
   fNData  = 0;
   fNMean1 = 0;
   if(fTmpFile) {delete fTmpFile; fTmpFile = 0;}
   if(fMean1)   {delete fMean1;   fMean1   = 0;}

// Set  mask array to 1
   Short_t  *arrMask = new Short_t[numrows];
   for (Int_t i=0; i<numrows; i++) arrMask[i] = 1;

// Init arrays
   Double_t *arr1 = new Double_t[numrows];
   Double_t *arr2 = new Double_t[numrows];
   Double_t *arr3 = new Double_t[numrows];
   Double_t *arr4 = new Double_t[numrows];
   Double_t *out1 = new Double_t[numrows];
   Double_t *out2 = new Double_t[numrows];
   Double_t *out3 = new Double_t[numrows];
   Double_t *out4 = new Double_t[numrows];

// Import infile
   char     nextline[kBufSize];
   ifstream input;
   input.open(infile, ios::in);
   // read header line containing column names
   input.getline(nextline, kBufSize, '\n');
   // read data
   Int_t x, y;
   for (Int_t i=0; i<numrows; i++) {
      input.getline(nextline, kBufSize, '\n');
      sscanf(nextline,"%i %i %lf %lf %lf %lf", &x, &y, &arr1[i], &arr2[i], &arr3[i], &arr4[i]);
      if (input.eof()) break;
   }//for
   input.close();

// Sort arrays for quantile normalization
   err = AddArray(numrows, arr1, arrMask, "rnktree1");
   err = AddArray(numrows, arr2, arrMask, "rnktree2");
   err = AddArray(numrows, arr3, arrMask, "rnktree3");
   err = AddArray(numrows, arr4, arrMask, "rnktree4");

// Quantile normalization: array out1 not used!
   err = Calculate(numrows, 0, out1, arrMask);

// Get normalized arrays
   out1 = GetArray(numrows, out1, arrMask, "rnktree1");
   out2 = GetArray(numrows, out2, arrMask, "rnktree2");
   out3 = GetArray(numrows, out3, arrMask, "rnktree3");
   out4 = GetArray(numrows, out4, arrMask, "rnktree4");

// Export normalized arrays to outfile
   ofstream output(outfile, ios::out);
   output << "TestA1" << "\t" << "TestA2" << "\t" << "TestB1" << "\t" << "TestB2" << endl;
   for (Int_t i=0; i<numrows; i++) {
      output << out1[i] << "\t" << out2[i] << "\t" << out3[i] << "\t" << out4[i] << endl;
   }//for
   output.close();

// Cleanup
   delete [] arrMask;
   delete [] out4; delete [] out3; delete [] out2; delete [] out1;
   delete [] arr4; delete [] arr3; delete [] arr2; delete [] arr1;
}//QuantileTest