//
// This file is part of 3DCellAtlas. 
// Copyright (C) 2015 George W. Bassel and collaborators.
//
// If you use 3DCellAtlas in your work, please cite:
//   http://dx.doi.org/10.1105/tpc.15.00175
//
// 3DCellAtlas is an AddOn for MorphoGraphX - http://www.MorphoGraphX.org
// Copyright (C) 2012-2015 Richard S. Smith and collaborators.
//
// 3DCellAtlas and MorphoGraphX are free software, and are licensed under under the terms of the 
// GNU General (GPL) Public License version 2.0, http://www.gnu.org/licenses.
// 
#include "DataAnalyzer.hpp"

#include <algorithm>

namespace mgx
{
//namespace process
//{
/*
General note on the code in DataAnalyzer:
This code was translated from the MatLab script "Data Analyzer" by Alexander Timothy Topham
into C++ for MGX
*/

// T_005 simply looks up the critical value form a t-distribution with alpha
// 0.05 for the given degrees of freedom, df, from a 2-tailed t.  Taken from
// values at
// http://easycalculation.com/statistics/t-distribution-critical-value-table.php
// It tops out at 200df.
double t005(double df)
  {

  double t [200] = {12.7065, 4.3026, 3.1824, 2.7764, 2.5706, 2.4469, 2.3646, 2.306, 2.2621, 2.2282, 
     2.201, 2.1788, 2.1604, 2.1448, 2.1314, 2.1199, 2.1098, 2.1009, 2.093, 2.086,
     2.0796, 2.0739, 2.0686, 2.0639, 2.0596, 2.0555, 2.0518, 2.0484, 2.0452, 2.0423, 
     2.0395, 2.0369, 2.0345, 2.0322, 2.0301, 2.0281, 2.0262, 2.0244, 2.0227, 2.0211,
     2.0196, 2.0181, 2.0167, 2.0154, 2.0141, 2.0129, 2.0117, 2.0106, 2.0096, 2.0086,
     2.0076, 2.0066, 2.0057, 2.0049, 2.0041, 2.0032, 2.0025, 2.0017, 2.001, 2.0003,
     1.9996, 1.999, 1.9983, 1.9977, 1.9971, 1.9966, 1.996, 1.9955, 1.995, 1.9944,
     1.9939, 1.9935, 1.993, 1.9925, 1.9921, 1.9917, 1.9913, 1.9909, 1.9904, 1.9901,
     1.9897, 1.9893, 1.9889, 1.9886, 1.9883, 1.9879, 1.9876, 1.9873, 1.987, 1.9867,
     1.9864, 1.9861, 1.9858, 1.9855, 1.9852, 1.985, 1.9847, 1.9845, 1.9842, 1.984,
     1.9837, 1.9835, 1.9833, 1.983, 1.9828, 1.9826, 1.9824, 1.9822, 1.982, 1.9818,
     1.9816, 1.9814, 1.9812, 1.981, 1.9808, 1.9806, 1.9805, 1.9803, 1.9801, 1.9799,
     1.9798, 1.9796, 1.9794, 1.9793, 1.9791, 1.979, 1.9788, 1.9787, 1.9785, 1.9784,
     1.9782, 1.9781, 1.9779, 1.9778, 1.9777, 1.9776, 1.9774, 1.9773, 1.9772, 1.9771,
     1.9769, 1.9768, 1.9767, 1.9766, 1.9765, 1.9764, 1.9762, 1.9761, 1.976, 1.9759,
     1.9758, 1.9757, 1.9756, 1.9755, 1.9754, 1.9753, 1.9752, 1.9751, 1.975, 1.9749,
     1.9748, 1.9747, 1.9746, 1.9745, 1.9744, 1.9744, 1.9743, 1.9742, 1.9741, 1.974,
     1.9739, 1.9739, 1.9738, 1.9737, 1.9736, 1.9735, 1.9735, 1.9734, 1.9733, 1.9732,
     1.9731, 1.9731, 1.973, 1.9729, 1.9729, 1.9728, 1.9727, 1.9727, 1.9726, 1.9725,
     1.9725, 1.9724, 1.9723, 1.9723, 1.9722, 1.9721, 1.9721, 1.972, 1.972, 1.9719};

  if(df < 201 && df > 0)
    return t[(int)df-1];  
  else if(df >= 201 && df <= 300)
    return 1.968;
  else if(df >= 301 && df <= 500)
    return 1.965;
  else if(df > 500)
    return 1.96;
  else
    //setErrorMessage("Error t_005, invalid degrees of freedom");
    return -1.0;
  //else if(df <= 0)
  //  error('stop feeding 0 degrees of freedom to t_005, Tiger');
  //else 
  //  error('what on earth are you feeding to t_005??');


  }
  //////////////////////////////////////////////////////////
  double tIntervalB(double t, std::vector<double> x, std::vector<double> y)
  {
/* T_INTERVAL returns the upper and lower confidence intervals based on the
% t-distribution. Defined as t*s/sqrt(n)
%
% varargin = cell array containing the value for t, dataA as a vector, and
% optionally dataB. With dataB specified, the confidence interval will be
% calculated using the additive variance of the two vectors.
*/
// dataB = y required here

  double meanX = 0.0, meanY = 0.0, stdX = 0.0, stdY = 0.0;

  for(double i = 0; i<x.size(); i++ ){
      meanX += x[i];
  }
  for(double i = 0; i<y.size(); i++ ){
      meanY += y[i];
  }
  if(x.size() > 0){
    meanX /= x.size();
    for(double i = 0; i<x.size(); i++ ){
      double fac = 1.0;
      if(x[i] - meanX < 0) fac = -1.0;
      stdX += (x[i] - meanX)*fac;
    }
    stdX /= x.size();
  }
  if(y.size() > 0){
    meanY /= y.size();
    for(double i = 0; i<y.size(); i++ ){
      double fac = 1.0;
      if(y[i] - meanY < 0) fac = -1.0;
      stdY += (y[i] - meanY)*fac;
    }
    stdY /= y.size();
  }

    //dataA = varargin{2};
    //dataB = varargin{3};
    
    double nX = x.size();
    double nY = y.size();
    
    return t * (stdX/sqrt(nX) + stdY/sqrt(nY)); 

  }


/*%LOGRATIO computes the log-ratio and 95% confidence intervals from x and y
%
% This function first takes the natural logarithm of the two column vetors
% x and y. These are then used to compute the log ratio of the mean of
% these two new vectors by subracting their means from one another.
% Confidence intervals are then computed using the t-distribution on these
% two log vectors. e is then raised to the power of these results to
% produce statistics on the same scale as the original data.
%
% x and y = column vectors of data. any '0's will be ignored, and if one or
% both of x and y are empty vectors, ratio will be returned as 'inf' and
% confidence intervals will both be 0.
%
% ratio = the log ratio of the means of x and y, end up as x/y
% upperCI = ratio plus the 95% confidence interval
% lowerCI = ratio minus the 95% confidence interval
%
% See also: T_INTERVAL, LOG, EXP
*/
bool DataAnalyzer::logRatio(double& ratio, double& upper, double& lower, std::vector<double> dataX, std::vector<double> dataY)
  {

  std::vector<double> processedX, processedY;
  double meanX = 0.0, meanY = 0.0;


  for(double i = 0; i<dataX.size(); i++ ){
    if(dataX[i] != 0){
      double logX = log(dataX[i]);
      processedX.push_back(logX);
      meanX += logX;
    }
  }
  for(double i = 0; i<dataY.size(); i++ ){
    if(dataY[i] != 0){
      double logY = log(dataY[i]);
      processedY.push_back(logY);
      meanY += logY;
    }
  }


  int df = processedX.size() + processedY.size() - 2;

  //if one or both of x and y are empty, there is no comparisson to make, so
  //set the ratio to inf and CIs to 0
  if(processedX.size() == 0 || processedY.size() == 0){
    ratio = -1;
    upper = 0;
    lower = 0;
    return true;
  } else {
    meanX /= processedX.size();
    meanY /= processedY.size();
    ratio = meanX - meanY;


    if(df > 1){
      double t = t005(df);
      double interval = tIntervalB(t, processedX, processedY);
      upper = exp(ratio + interval);
      lower = exp(ratio - interval);
    } else {
      upper = 0;
      lower = 0;
    }

    ratio = exp(ratio);

    return true;
  }
}

bool logRatioPerFile(double& ratio, double& upper, double& lower, std::vector<std::vector<double> > dataX, std::vector<std::vector<double> > dataY)
  {

  std::vector<double> processedX, processedY;
  std::vector<double> meanX(dataX.size()), meanY(dataY.size());

  int counter = 0, counter2 = 0;
  double meanNr = 0;

  double totalMeanX = 0, totalMeanY = 0;

  for(double i = 0; i<dataX.size(); i++ ){
    meanX[i] = 0;
    counter = 0;
    for(double j = 0; j<dataX[i].size(); j++ ){
      if(dataX[i][j] != 0){
        double logX = log(dataX[i][j]);
        processedX.push_back(logX);
        meanX[i] += logX;
        counter++;
      }
    }
    if(counter>0){ meanX[i]/=counter; counter2++; }
    meanNr+=counter;
    //cout << counter << "/" << counter2 << "//";
    totalMeanX+=meanX[i];
  }

  if(counter2>0){ totalMeanX/=counter2; meanNr/=counter2; }
  counter2 = 0;
  cout << meanNr << ",";
  meanNr = 0;



  for(double i = 0; i<dataY.size(); i++ ){
    meanY[i] = 0;
    counter = 0;
    for(double j = 0; j<dataY[i].size(); j++ ){
      if(dataY[i][j] != 0){
        double logY = log(dataY[i][j]);
        processedY.push_back(logY);
        meanY[i] += logY;
        counter++;
      }
    }
    if(counter>0){ meanY[i]/=counter; counter2++; }
    //cout << counter << "/" << counter2 << "//";
    meanNr+=counter;
    totalMeanY+=meanY[i];

  }
  //cout << endl;
  if(counter2>0){ totalMeanY/=counter2; meanNr/=counter2; }
  cout << meanNr << endl;

  int df = processedX.size() + processedY.size() - 2;

  //if one or both of x and y are empty, there is no comparisson to make, so
  //set the ratio to inf and CIs to 0
  if(processedX.size() == 0 || processedY.size() == 0){
    ratio = -1;
    upper = 0;
    lower = 0;
    return true;
  } else {
    
    ratio = totalMeanX - totalMeanY;

    if(df > 1){
      double t = t005(df);
      double interval = tIntervalB(t, processedX, processedY);
      upper = exp(ratio + interval);
      lower = exp(ratio - interval);
    } else {
      upper = 0;
      lower = 0;
    }

    ratio = exp(ratio);

    return true;
  }
}

  //////////////////////////////////////////////////////////
  // opens a single file and writes the information into tagData
  bool openFile(QString filename, tagDataMap& tagData, tagTypeVector& tagTypes, IntDoubleMap labelGUSMap, QString justFile)
  {

  QFile file(filename);
    if(!file.open(QIODevice::ReadOnly))
    {
      cout << "error  " << filename.toLocal8Bit().constData() << endl;
      //setErrorMessage(QString("File '%1' cannot be opened for reading").arg(filename));
      return false;
    }
    QTextStream ss(&file);
    QString line = ss.readLine();
    QStringList fields = line.split(",");

    // variables to be filled
    int counter = 0;

    bool ok;
    //cout << "start" << endl;
    while(ss.status() == QTextStream::Ok)
    {
      fields = ss.readLine().split(",");
      if(fields[0].isEmpty()) break;

      FileDataType line;

      int cellNr = fields[0].toInt(&ok);
      int cellType = fields[1].toInt(&ok);
      int asocCort = fields[2].toInt(&ok);
      double radial = fields[3].toDouble(&ok);
      double lenLong = fields[4].toDouble(&ok);
      double lenRad = fields[5].toDouble(&ok);
      double lenCirc = fields[6].toDouble(&ok);
      double wallArea = fields[7].toDouble(&ok);
      double cellVol = fields[8].toDouble(&ok);

      //if(counter != 0){
        line.cellNumber = cellNr;
        line.cellType = cellType;
        line.asocCortCell = asocCort;
        line.radial = radial;
        line.lenLong = lenLong;
        line.lenRad = lenRad;
        line.lenCirc = lenCirc;
        line.wallArea = wallArea;
        line.cellVol = cellVol;
        line.file = justFile;

        line.GUS = labelGUSMap[cellNr];

        //if(/*line.GUS !=0 && */cellType==7 && (asocCort==3))
        //  cout << "lll  " << cellNr << " t " << cellType << " acc " << asocCort << " gus " << line.GUS << "  " << justFile.toLocal8Bit().constData() << endl;

        //data.push_back(line);
        Vec2i currentTag (cellType, asocCort);
        bool exists = false;
        for(int j = 0; j < tagTypes.size(); j++){
          if(tagTypes[j] == currentTag)
            exists = true;
        }
        if(!exists)
          tagTypes.push_back(currentTag);
        //if(cellNr == 5425) cout << "Im in" << line.cellNumber << "  " << cellType*100 + asocCort << endl;
        tagData[cellType*100 + asocCort].push_back(line);
      //}
      //cout << "read file " << cellVol << endl;
      //cout << "lll  " << label << " parent" << parent << " long " << cellLengthLong << " rad " << cellLengthRad << " circ " << cellLengthCirc << endl;
      //uniqueLabels.push_back(label);
      //assocCC[label] = fields[2].toInt(&ok);
      //arclengths[label] = fields[3].toDouble(&ok);
      counter++;
    }

  return true;
  }

  void switchTags(Vec2i& i, Vec2i& j)
  {

  Vec2i h = i;
  i = j;
  j = h;

  }

  //////////////////////////////////////////////////////////
  // goes through the folder and opens all files
  bool DataAnalyzer::readCellData(QString folder, tagDataMap& tagData, tagTypeVector& tagTypes, FilenameVector& filenames)
  {
  // go to folder, open all files in it and write data while taking care of the headers
    DIR *pDIR;
    struct dirent *entry;
    if( pDIR=opendir(folder.toLocal8Bit().constData()) ){

      while(entry = readdir(pDIR)){
        if( strcmp(entry->d_name, ".") != 0 && strcmp(entry->d_name, "..") != 0 ){
          // get data from current file and add data to folder data structure
          IntDoubleMap dummyMap;
          openFile(folder + "/" + entry->d_name, tagData, tagTypes, dummyMap, entry->d_name);
          filenames.push_back(entry->d_name);
          
        }
      }
      closedir(pDIR);
    }

  // sort the tagtypes (this will make things easier later)
  for(int i = 0; i < tagTypes.size(); i++){
    for(int j = 0; j < tagTypes.size(); j++){
      if(tagTypes[i][0] < tagTypes[j][0])
      {
        switchTags(tagTypes[i], tagTypes[j]);
      } else if(tagTypes[i][0] == tagTypes[j][0])
      {
        if(tagTypes[i][1] < tagTypes[j][1]) switchTags(tagTypes[i], tagTypes[j]);
      }
    }

  }

  return true;
  }

  // extract the right "column" from the data "table"
  std::vector<double> getDataVec(FileDataVector allData, int col)
  {
    std::vector<double> dataVector;

    for(int j = 0; j < allData.size(); j++){
      if(col==1)
        dataVector.push_back(allData[j].radial);
      else if(col==2)
        dataVector.push_back(allData[j].lenLong);
      else if(col==3)
        dataVector.push_back(allData[j].lenRad);
      else if(col==4)
        dataVector.push_back(allData[j].lenCirc);
      else if(col==5)
        dataVector.push_back(allData[j].wallArea);
      else if(col==6)
        dataVector.push_back(allData[j].cellVol);
    }

    return dataVector;
  }

 // extract the right "column" from the data "table"
  std::vector<std::vector<double> > getDataVecPerFile(FileDataVector allData, int col)
  {
    std::vector<std::vector<double> > dataVector;
    std::map<QString, int> fileIndexMap;

    for(int j = 0; j < allData.size(); j++){
      int currentIndex;
      if(fileIndexMap[allData[j].file] == 0){
        fileIndexMap[allData[j].file] = dataVector.size()+1;
        std::vector<double> emtpyVector;
        dataVector.push_back(emtpyVector);
      }
      currentIndex = fileIndexMap[allData[j].file] - 1;

      if(col==1)
        dataVector[currentIndex].push_back(allData[j].radial);
      else if(col==2)
        dataVector[currentIndex].push_back(allData[j].lenLong);
      else if(col==3)
        dataVector[currentIndex].push_back(allData[j].lenRad);
      else if(col==4)
        dataVector[currentIndex].push_back(allData[j].lenCirc);
      else if(col==5)
        dataVector[currentIndex].push_back(allData[j].wallArea);
      else if(col==6)
        dataVector[currentIndex].push_back(allData[j].cellVol);
    }

    return dataVector;
  }

/*
% SORTDATA takes two data sets (of the format output by SortData.m), and
% for each cell type at a given associated corticle cell (acc) applies the 
% matlab function found in analyseFunc to the data.
%
% the method found in analyseFunc is accessed using str2func(file) where
% file is the name of the .m file containing the function, which must be
% placed in the same directory as sortData. For each cell tag, if a
% matching tag cannot be found in the other data set, that point is egnored
% for the analysis.
%
% data1 and data2
% [cell#, type, acc, [data columns], tag, file]
%
% sortedData = a cell array with the first column containing the name of
% the metric (eg radius), and second column containing the results of this
% analysis in the form [cell type, acc, result, upperCI, lowerCI, tag]
*/

  //////////////////////////////////////////////////////////
  // sort the input data according to cell type and acc and calculate foreach tag the ratios
  bool DataAnalyzer::sortData(tagDataMap data1, tagDataMap data2, tagTypeVector tagTypes)
  {

  // loop through measures x tags and run the statistics
  // for all measures (rad / lenL / lenR / lenC / wallArea / vol)
  for(int i = 1; i <= 6; i++){
    // for all tags (= cell type x assoc cort cell)
    SortResultVector resVec;
    for(int j = 0; j < tagTypes.size(); j++){
      SortResultData res;

      int tag = tagTypes[j][0]*100 + tagTypes[j][1]; // current tag
      FileDataVector currentDataX = data1[tag];
      FileDataVector currentDataY = data2[tag];

      res.cellType = tagTypes[j][0];
      res.asocCortCell = tagTypes[j][1];

      if(currentDataX.size() < 1 || currentDataY.size() < 1){ // give error values if too small
        res.ratio = -1;
        res.upper = -1;
        res.lower = -1;
        res.measure = i;
      } else { // everything fine, do analysis
        double ratio, upper, lower;

        if(perCellAvg){
          // put the current measure in a vector
          std::vector<double> data1Vec = getDataVec(currentDataX,i);
          std::vector<double> data2Vec = getDataVec(currentDataY,i);

          // run the logratio
          logRatio(ratio, upper, lower, data1Vec, data2Vec);
        } else {
          // put the current measure in a vector
          std::vector<std::vector<double> > data1Vec = getDataVecPerFile(currentDataX,i);
          std::vector<std::vector<double> > data2Vec = getDataVecPerFile(currentDataY,i);

          cout << " tag " << tag << ",";
          // run the logratio
          logRatioPerFile(ratio, upper, lower, data1Vec, data2Vec);


        }
        //if(tag==320)
        //  cout << "ratio " << ratio << "/" << data1Vec.size() << "/" << data2Vec.size() << endl;

        // save the results
        res.ratio = ratio;
        res.upper = upper;
        res.lower = lower;
        res.measure = i;

      }

      // save results
      resVec.push_back(res);
    }
    resVecAll.push_back(resVec);
  }

  return true;
  }

  //////////////////////////////////////////////////////////
  // calculates the average value of a vector (from start to end)
  double avgVector(std::vector<double> data, int start, int end)
  {

  if(start<0 || start>end || end>data.size())
    return -1E20;

  double sum = 0;
  for(int i = start; i <= end; i++){
    sum+=data[i];
  }
  return (sum/(end-start+1));

  }

  //////////////////////////////////////////////////////////
  // calculates the sliding avg of a vector
  bool slidingAvg(std::vector<double>& data, int window)
  {

  int side = window/2; // default input 3 gives a side = 1
  std::vector<double> dataNew;

  if(side*2 >= data.size()) // data too small or window too large, do nothing
    return true;

  // OK now the beginning...
  for(int i = 0; i < side; i++){
    dataNew.push_back(avgVector(data,0,i+side));
  }

  // Lets do the easy bits
  for(int i = side; i < data.size()-side-1; i++){
    dataNew.push_back(avgVector(data,i-side,i+side));
  }

  // And the tail...
  for(int i = data.size()-side-1; i < data.size(); i++){
    dataNew.push_back(avgVector(data,i-side,data.size()-1));
  }

  data = dataNew;
  return true;
  }

  //////////////////////////////////////////////////////////
  // smooth the sortResVec with a sliding avg with window
  bool DataAnalyzer::smoothAnalyzed(SortResultVectorAll& sortResVec, double window, tagTypeVector tagTypes)
  {

  IntIntMap cellTypeCounter;

  // count how many asocCortCell types of each cellType
  for(int i = 0; i < tagTypes.size(); i++){
    //cout << "types " <<  tagTypes[i] << endl;
    cellTypeCounter[tagTypes[i][0]]++;
  }

  for(int t = 0; t < 6; t++){ // all different data sets
    for(int i = 0; i < cellTypeCounter.size(); i++){ // all cell types
      std::vector<double> dataRad;//, dataRadOrdered, indexData;
      for(int k = 0; k < sortResVec[0].size(); k++){
        if(sortResVec[t][k].cellType == i){
          dataRad.push_back(sortResVec[t][k].ratio); // fill dataRad in the right order
          //indexData.push_back(sortResVec[t][k].asocCortCell);
        }
      }
      // perform slidingavg
      slidingAvg(dataRad, window);

      // copy values back in the same order
      int counter = 0;
      for(int k = 0; k < sortResVec[t].size(); k++){
        //cout << "kkk   " << sortResVec[0][k].cellType << "  " << sortResVec[0][k].asocCortCell << endl;
        if(sortResVec[t][k].cellType == i){
          sortResVec[t][k].ratio = dataRad[counter];
          counter++;
        }
      }
    }

  }

  return true;
  }

  //////////////////////////////////////////////////////////
  // save the result in a file
  bool saveData(QString filename, QString folderOut, FileDataVector joinedData, SortResultVectorAll sortResVec, tagTypeVector tagTypes)
  {
    QFile file(folderOut + "/" + filename);
    if(!file.open(QIODevice::WriteOnly))
    {
      //setErrorMessage(QString("File '%1' cannot be opened for writing").arg(filename));
      return false;
    }
    // write the file
    QTextStream out(&file);

    // first the header
    out << "Cell Number,Radius,Longitudinal Length,Radial Length,Circumferential Length,Cell Wall Area,Cell Volume" << endl;

    // now the data
    for(int i = 0; i<joinedData.size(); i++){
      out << joinedData[i].cellNumber << "," << /*joinedData[i].cellType << "," << joinedData[i].asocCortCell << "," << */joinedData[i].radial << "," << joinedData[i].lenLong << "," << joinedData[i].lenRad << "," << joinedData[i].lenCirc << "," << joinedData[i].wallArea << "," << joinedData[i].cellVol << endl;
    }

  return true;
  }
  //////////////////////////////////////////////////////////
  // save cell type data with CI
  bool saveDataCI(QString filename, QString folderOut, SortResultVectorAll sortResVec, tagTypeVector tagTypes)
  {
  QFile file2(folderOut + "/" + filename);
    if(!file2.open(QIODevice::WriteOnly))
    {
      //setErrorMessage(QString("File '%1' cannot be opened for writing").arg(filename));
      return false;
    }
    // write the file
    QTextStream out2(&file2);

    // first the header
    out2 << "Cell Type,Associated Cortical Cell,Radial,Radial Upper CI,Radial Lower CI,Longitudinal Length,Longitudinal Upper CI,Longitudinal Lower CI,Radial Length,Radial Upper CI,Radial Lower CI,Circumferential Length,Circumferential Upper CI,Circumferential Lower CI,Cell Wall Area,Cell Wall Upper CI,Cell Wall Lower CI,Cell Volume,Cell Volume Upper CI,Cell Volume Lower CI" << endl;

    // now the data
    for(int i = 0; i<tagTypes.size(); i++){
      int tag = tagTypes[i][0]*100 + tagTypes[i][1];
      for(int k = 0; k < sortResVec[0].size(); k++){
        if(sortResVec[0][k].cellType*100 + sortResVec[0][k].asocCortCell == tag){
          out2 << sortResVec[0][k].cellType << "," << sortResVec[0][k].asocCortCell << ","
               << sortResVec[0][k].ratio << "," << sortResVec[0][k].upper << "," << sortResVec[0][k].lower << ","
               << sortResVec[1][k].ratio << "," << sortResVec[1][k].upper << "," << sortResVec[1][k].lower << ","
               << sortResVec[2][k].ratio << "," << sortResVec[2][k].upper << "," << sortResVec[2][k].lower << ","
               << sortResVec[3][k].ratio << "," << sortResVec[3][k].upper << "," << sortResVec[3][k].lower << ","
               << sortResVec[4][k].ratio << "," << sortResVec[4][k].upper << "," << sortResVec[4][k].lower << ","
               << sortResVec[5][k].ratio << "," << sortResVec[5][k].upper << "," << sortResVec[5][k].lower << ","
               << endl;
        }
      }
    }

  return true;
  }
  //////////////////////////////////////////////////////////
  // save the e factor files
  bool DataAnalyzer::saveEFactors(FileEFactorDataVector eFData, tagDataMap data, tagTypeVector tagTypes, FilenameVector filenames, QString folderOut)
  {

  // foreach input file there will be one new e factor file
  for(int k = 0; k<filenames.size(); k++){

    // assemble the filename
    QString filename = filenames[k];
    filename.remove(".csv", Qt::CaseInsensitive);
    filename += "_eFactor";
    filename += ".csv";

    QFile file(folderOut + "/" + filename );
    if(!file.open(QIODevice::WriteOnly))
    {
      //setErrorMessage(QString("File '%1' cannot be opened for writing").arg(filename));
      return false;
    }
    // write the file
    QTextStream out(&file);

    // first the header
    out << "Cell Number,Radial eFactor,Longitudinal eFactor,Circumferential eFactor" << endl;

    // now the data (a bit more tricky here as the cells have to be matched with the values of their type/acc)
    for(int i = 0; i<tagTypes.size(); i++){
      int tag = tagTypes[i][0]*100 + tagTypes[i][1];
      for(int j = 0; j<data[tag].size(); j++){
        for(int m = 0; m<eFData.size(); m++){
          if(eFData[m].cellType == tagTypes[i][0] && eFData[m].asocCortCell == tagTypes[i][1] && filenames[k] == data[tag][j].file)
            out << data[tag][j].cellNumber << "," << eFData[m].facRad << "," << eFData[m].facLong << "," << eFData[m].facCirc << endl;
        }
      }
    }

  }

  return true;
  }

  //////////////////////////////////////////////////////////
  // calculate the e factors
  bool DataAnalyzer::calcEFactors(SortResultVectorAll sortResVec, bool avg, int window, tagTypeVector tagTypes, FileEFactorDataVector& eFData)
  {

  IntIntMap cellTypeCounter;
  // count how many asocCortCell types of each cellType
  for(int i = 0; i < tagTypes.size(); i++){
    cellTypeCounter[tagTypes[i][0]]++;
  }

  int counter = 0;
  // for all cell types
  for(int i = 0; i < cellTypeCounter.size(); i++){
    std::vector<double> dataLong, dataRad, dataCirc;
    FileEFactorDataVector temp;
    for(int k = 0; k < sortResVec[0].size(); k++){
      if(sortResVec[0][k].cellType == i){
        FileEFactorDataType currentEFactors;
        currentEFactors.cellType = sortResVec[0][k].cellType;
        currentEFactors.asocCortCell = sortResVec[0][k].asocCortCell;
        double lon = sortResVec[1][k].ratio;
        double rad = sortResVec[2][k].ratio;
        double circ = sortResVec[3][k].ratio;
        double tot = lon + rad + circ;
        currentEFactors.facLong = lon/tot;
        currentEFactors.facRad = rad/tot;
        currentEFactors.facCirc = circ/tot;
        counter++;
        temp.push_back(currentEFactors); // save here as there might be smoothing necessary
        dataLong.push_back(currentEFactors.facLong);
        dataRad.push_back(currentEFactors.facRad);
        dataCirc.push_back(currentEFactors.facCirc);
      } // if
    } // for

    // sliding avg?
    if(avg){

      // do the avg
      slidingAvg(dataLong, window);
      slidingAvg(dataRad, window);
      slidingAvg(dataCirc, window);

      for(int k = 0; k < temp.size(); k++){
        // update the values
        temp[k].facLong = dataLong[k];
        temp[k].facRad = dataRad[k];
        temp[k].facCirc = dataCirc[k];
        eFData.push_back(temp[k]); // now save the data in the final structure
      }    

    } else {
      for(int k = 0; k < temp.size(); k++){
        eFData.push_back(temp[k]); // now save the data in the final structure
      } // for
    } // if
  } // for
  return true;
  }

  //////////////////////////////////////////////////////////
  // takes the result of sortData, joins it back with the original data to make it ready to get saved
  bool DataAnalyzer::joinData(SortResultVectorAll sortResVec, tagDataMap data, tagTypeVector tagTypes, FilenameVector& filenames, QString folderOut)
  {

  int fileCol = tagTypes.size();

  // OK now for each file name, we need to read down all the tags in column 10
  // and pull them out of each type from sortedData. First iterate over all
  // the file names and get the relevent data each time

  for(int i = 0; i < filenames.size(); i++){ // go through files & build data structure only with info from current file
    FileDataVector joinedData;
    FileDataVector currentData;

    for(int i2 = 0; i2 < tagTypes.size(); i2++){
      int tag = tagTypes[i2][0]*100 + tagTypes[i2][1]; // current tag
      for(int i3 = 0; i3 < data[tag].size(); i3++){
        if(filenames[i] == data[tag][i3].file) // save it if it also belongs to the right file
          currentData.push_back(data[tag][i3]);
      }
    }

    // Iterate over all the rows in currentData
    for(int j = 0; j < currentData.size(); j++){
      int currentTag = currentData[j].cellType*100 + currentData[j].asocCortCell;

      // Drop all data for air, root cap, vasculature and columella; these
      // are types 0, 5, 6 and 7. Converting the data to numerical matrix
      // when it's saved removes the empty rows automatically...
      if(currentData[j].cellType == 0 || currentData[j].cellType == 5 || currentData[j].cellType == 6 || currentData[j].cellType == 7){
        // ignore
      } else {
        // Now go thorugh each data type in sortedData and find the mean for
        // currentTag (yaaaay an n^3 loop! Could I sink any lower?). The tag
        // in each sortedData is in column 6.
        int foundC = 0;
        for(int k = 0; k < sortResVec[0].size(); k++){

          if(foundC > 1)
            cout << "WARNING: Problem in joinData" << endl;
          // find tag in sorted data
          if(sortResVec[0][k].cellType == currentData[j].cellType && sortResVec[0][k].asocCortCell == currentData[j].asocCortCell){
            // (rad / lenL / lenR / lenC / wallArea / vol)
            currentData[j].radial = sortResVec[0][k].ratio;
            currentData[j].lenLong = sortResVec[1][k].ratio;
            currentData[j].lenRad = sortResVec[2][k].ratio;
            currentData[j].lenCirc = sortResVec[3][k].ratio;
            currentData[j].wallArea = sortResVec[4][k].ratio;
            currentData[j].cellVol = sortResVec[5][k].ratio;
            foundC++;
            joinedData.push_back(currentData[j]);
          }
        }
        // not really needed
        /*if(!found){
          currentData[j].radial = -1;
          currentData[j].lenLong = -1;
          currentData[j].lenRad = -1;
          currentData[j].lenCirc = -1;
          currentData[j].wallArea = -1;
          currentData[j].cellVol = -1;
        }*/
        //joinedData.push_back(currentData[j]);
      }
    }

    saveData(filenames[i], folderOut, joinedData, sortResVec, tagTypes);
  }
  saveDataCI("Data_CI.csv", folderOut, sortResVec, tagTypes);
  return true;
  }

  double resampleAvg(std::vector<double> data)
  {
  std::vector<double> permutation;
  double sum = 0;


  for(int i = 0; i < data.size(); i++){
    int rnd = std::rand() % data.size();
    permutation.push_back(data[rnd]);
    sum += data[rnd];
  }

  return sum/data.size();
  }


  bool calcCI(std::vector<double> data, double& upper, double& lower)
  {
  int it = 100;
  std:vector<double> results;

  for(int i = 0; i < it; i++){
    results.push_back(resampleAvg(data));
  }
  sort(results.begin(), results.end());

  upper = results[94];
  lower = results[4];

  return true;
  }


  //////////////////////////////////////////////////////////
  // processGUSData
  bool processGUSData(IntPoint3fMap& averages, tagDataMap& tagData, tagTypeVector& tagTypes)
  {
  //cout << "GUS" << endl;
  // calc means
  // go through all existing tags
  for(int i = 0; i < tagTypes.size(); i++){
    std::vector<double> currentData;
    int tag = tagTypes[i][0]*100 + tagTypes[i][1]; // current tag

    double mean = 0;
    int zeroCounter = 0;
    // go through all cells within each tag
    for(int j = 0; j < tagData[tag].size(); j++){
      //mean+=tagData[tag][j].GUS;
      if(tagData[tag][j].GUS == -1){ zeroCounter++; }
      else {
        mean+=tagData[tag][j].GUS;
        currentData.push_back(tagData[tag][j].GUS);
      }
      //if(tag==703)
      //  cout << "703   " << tagData[tag][j].GUS << "  " << mean  << "  " << zeroCounter << endl;
    }

    if(tagData[tag].size() - zeroCounter > 0) mean/=(tagData[tag].size()- zeroCounter);
    averages[tag][0] = mean;

    if(currentData.size() > 7){
      double upper, lower;
      calcCI(currentData, upper, lower);
      averages[tag][1] = upper;
      averages[tag][2] = lower;
    } else {
      averages[tag][1] = -1;
      averages[tag][2] = -1;
    }

     // if(tag==703)
     //   cout << "703   " << tagData[tag].size() << "  " <<  zeroCounter << "  " << mean << endl;
  }
/* taken from upper/lower
  // calc std
  for(int i = 0; i < tagTypes.size(); i++){
    int tag = tagTypes[i][0]*100 + tagTypes[i][1]; // current tag

    double std = 0;
    double mean = averages[tag][0];
    for(int j = 0; j < tagData[tag].size(); j++){
      std+=std::max(mean-tagData[tag][j].GUS, tagData[tag][j].GUS-mean);
    }
    std/=tagData[tag].size();
    averages[tag][1] = mean + std;
    averages[tag][2] = mean - std;
  }
*/
  // replace GUS values by means
  for(int i = 0; i < tagTypes.size(); i++){
    int tag = tagTypes[i][0]*100 + tagTypes[i][1]; // current tag

    for(int j = 0; j < tagData[tag].size(); j++){
      tagData[tag][j].GUS = averages[tag][0];
      tagData[tag][j].GUSUpper = averages[tag][1];
      tagData[tag][j].GUSLower = averages[tag][2];
      //if(tag==104) 
      //cout << "tag  " << tag << "   " << averages[tag][0] << "  " << averages[tag][1] << endl;
    }
  }

  return true;
  }

  //////////////////////////////////////////////////////////
  // save GUS data 
  bool saveFileGUSExternal(QString folderOut, QString filename,  tagDataMap tagData, tagTypeVector tagTypes, IntDoubleMap& labelGUSMap)
  {

  // extract cell nr, type and assoc cort cell
  FileDataVector fileInfo;

    QFile file(filename);
    if(!file.open(QIODevice::ReadOnly))
    {
      cout << "error  " << filename.toLocal8Bit().constData() << endl;
      //setErrorMessage(QString("File '%1' cannot be opened for reading").arg(filename));
      return false;
    }
    QTextStream ss(&file);
    QString line = ss.readLine();
    QStringList fields = line.split(",");

    // variables to be filled
    int counter = 0;

    bool ok;
    //cout << "start" << endl;
    while(ss.status() == QTextStream::Ok)
    {
      FileDataType currentLine;
      fields = ss.readLine().split(",");
      if(fields[0].isEmpty()) break;

      currentLine.cellNumber = fields[0].toInt(&ok);
      currentLine.cellType = fields[1].toInt(&ok);
      currentLine.asocCortCell = fields[2].toInt(&ok);

      fileInfo.push_back(currentLine);

    }

  // go through all cells, get GUS value and write it into new file
  filename.remove(".csv", Qt::CaseInsensitive);
  filename += "_output";
  filename += ".csv";

    QFile file2(filename);
    if(!file2.open(QIODevice::WriteOnly))
    {
      //setErrorMessage(QString("File '%1' cannot be opened for writing").arg(filename));
      return false;
    }
    // write the file
    QTextStream out(&file2);

    // first the header
    out << "Cell Number,Mean Reporter Signal" << endl;
    cout << "2nd file " << endl;
    // now the data
    for(int i = 0; i<fileInfo.size(); i++){

      int tag = fileInfo[i].cellType*100 + fileInfo[i].asocCortCell; // current tag
      bool exists = false;
    //cout << "2nd file " << i << " tag " << tag << endl;
      // check if that tag exists
      for(int j = 0; j<tagTypes.size(); j++){
        int tag2 = tagTypes[j][0]*100 +  tagTypes[j][1];
        if(tag2 == tag) exists = true;
      }
      double gus = 0;
      if(exists) gus = tagData[tag][0].GUS;
      if(fileInfo[i].cellType >=1 && fileInfo[i].cellType<=4)
        labelGUSMap[fileInfo[i].cellNumber] = gus;
      out << fileInfo[i].cellNumber << "," << gus << endl;
    }





  return true;
  }

  //////////////////////////////////////////////////////////
  // save GUS
  bool saveFileGUS(QString filename, QString filename2, QString folderOut, tagDataMap tagData, tagTypeVector tagTypes)
  {
/* // not needed anymore, information is in the CI file
    QFile file(folderOut + "/" + filename);
    if(!file.open(QIODevice::WriteOnly))
    {
      //setErrorMessage(QString("File '%1' cannot be opened for writing").arg(filename));
      return false;
    }
    // write the file
    QTextStream out(&file);

    // first the header
    out << "Cell Type,Associated Cortical Cell,Metric" << endl;

    // now the data
    for(int i = 0; i<tagTypes.size(); i++){
      int tag = tagTypes[i][0]*100 + tagTypes[i][1]; // current tag
      out << tagData[tag][0].cellType << "," << tagData[tag][0].asocCortCell << "," << tagData[tag][0].GUS << endl;
    }
*/
  // save the file for the confidence intervalls
  QFile file2(folderOut + "/" + filename2);

  if(!file2.open(QIODevice::WriteOnly))
    {
      //setErrorMessage(QString("File '%1' cannot be opened for writing").arg(filename));
      return false;
    }
    // write the file
    QTextStream out2(&file2);

    // first the header
    out2 << "Cell Type,Associated Cortical Cell,Mean Reporter Signal, Upper CI, Lower CI" << endl;

    // now the data
    for(int i = 0; i<tagTypes.size(); i++){
      int tag = tagTypes[i][0]*100 + tagTypes[i][1]; // current tag
      QString upper, lower;
      if(tagData[tag][0].GUSUpper == -1){
        upper = "";
      } else {
        upper = QString::number(tagData[tag][0].GUSUpper);
      }

      if(tagData[tag][0].GUSLower == -1){
        lower = "";
      //cout << "aaaa " << lower.toLocal8Bit().constData() << "  " << tagData[tag][0].GUSLower << endl;
      } else {
        lower = QString::number(tagData[tag][0].GUSLower);
      //cout << "ffff " << lower.toLocal8Bit().constData() << "  " << tagData[tag][0].GUSLower << endl;
      }

      out2 << tagData[tag][0].cellType << "," << tagData[tag][0].asocCortCell << "," << tagData[tag][0].GUS << "," << upper << "," << lower << endl;
    }


  return true;
  }

  //////////////////////////////////////////////////////////
  // opens a single file and writes the information into tagData
  bool openFileGUS(QString filename, IntDoubleMap& labelGUSMap)
  {

  QFile file(filename);
    if(!file.open(QIODevice::ReadOnly))
    {
      cout << "error  " << filename.toLocal8Bit().constData() << endl;
      //setErrorMessage(QString("File '%1' cannot be opened for reading").arg(filename));
      return false;
    }
    QTextStream ss(&file);
    QString line = ss.readLine();
    QStringList fields = line.split(",");

    // variables to be filled
    int counter = 0;

    bool ok;
    //cout << "start" << endl;
    while(ss.status() == QTextStream::Ok)
    {
      fields = ss.readLine().split(",");
      if(fields[0].isEmpty()) break;

      int cellNr = fields[0].toInt(&ok);
      double gus = fields[1].toDouble(&ok);

      labelGUSMap[cellNr] = gus;

      
    }

  return true;
  }


//////////////////////////////////////////////////////////
  // goes through the folder, checks all files in the folder and check whether theres is filename+gusFileExt
  bool readCellDataGUS(QString folder, QString gusFileExt, tagDataMap& tagData, tagTypeVector& tagTypes, FilenameVector& filenames)
  {

  IntDoubleMap labelGUSMap;



  // go to folder, open all files in it and write data while taking care of the headers
    DIR *pDIR;
    struct dirent *entry;
    if( pDIR=opendir(folder.toLocal8Bit().constData()) ){

      while(entry = readdir(pDIR)){
        if( strcmp(entry->d_name, ".") != 0 && strcmp(entry->d_name, "..") != 0 ){
          // get data from current file and add data to folder data structure
          //openFile(folder + "/" + entry->d_name, tagData, tagTypes, entry->d_name);
          filenames.push_back(entry->d_name);    
        }
      }
      closedir(pDIR);
    }

  for(int i = 0; i < filenames.size(); i++){
  QString currentFile = filenames[i];
  QString searchFile = currentFile;
  searchFile.remove(".csv", Qt::CaseInsensitive);
  searchFile += gusFileExt;
  searchFile += ".csv";
  //cout << currentFile.toLocal8Bit().constData() << "  " << searchFile.toLocal8Bit().constData() << endl;
    for(int j = 0; j < filenames.size(); j++){
      if(filenames[j] == searchFile){
        cout << "found pair " << currentFile.toLocal8Bit().constData() << "  " << searchFile.toLocal8Bit().constData() << endl;
        // now open these two files
        // init labelGUSmap with invalid entries
        for(int i = 0; i<10000; i++){
          labelGUSMap[i] = -1;
        }
        // GUS file to an intintmap
        openFileGUS(folder + "/" + searchFile, labelGUSMap);

        int oldLength = tagData.size();
        openFile(folder + "/" + currentFile, tagData, tagTypes, labelGUSMap, currentFile);
        int newLength = tagData.size();

        // now merge the GUS and file data
        //for(int k = oldLength-1; k < newLength; k++){
          //tagData[k].GUS = labelGUSMap[tagData[k].cellNumber];
        //}

      }
    }
  }


  // sort the tagtypes (this will make things easier later)
  for(int i = 0; i < tagTypes.size(); i++){
    for(int j = 0; j < tagTypes.size(); j++){
      if(tagTypes[i][0] < tagTypes[j][0])
      {
        switchTags(tagTypes[i], tagTypes[j]);
      } else if(tagTypes[i][0] == tagTypes[j][0])
      {
        if(tagTypes[i][1] < tagTypes[j][1]) switchTags(tagTypes[i], tagTypes[j]);
      }
    }

  }

  return true;
  }

  //////////////////////////////////////////////////////////
  // smoothing the GUS data, this is not too different from the other smoothing
  bool smoothGUS(tagDataMap& tagData, tagTypeVector& tagTypes, int window)
  {

  IntIntMap cellTypeCounter;

  // count how many asocCortCell types of each cellType
  for(int i = 0; i < tagTypes.size(); i++){
    //cout << "types " <<  tagTypes[i] << endl;
    cellTypeCounter[tagTypes[i][0]]++;
  }

  //for(int t = 0; t < 6; t++){ // all different data sets
    for(int i = 0; i < cellTypeCounter.size(); i++){
      std::vector<double> dataVec;
      for(int k = 0; k < tagTypes.size(); k++){
        int tag = tagTypes[k][0]*100 + tagTypes[k][1]; // current tag
        if(tagData[tag][0].cellType == i){
          dataVec.push_back(tagData[tag][0].GUS);
          //cout << "b " << tag << "  " << tagData[tag][0].GUS << endl;
        }
      }
      //cout << "now avg " <<  window << endl;
      // perform slidingavg
      slidingAvg(dataVec, window);
      //cout << "avg done " << endl;

      //for(int k = 0; k < dataVec.size(); k++){
      //  cout << dataVec[k] << endl;
     // }


      // copy values back
      int counter = 0;
      for(int k = 0; k < tagTypes.size(); k++){
        int tag = tagTypes[k][0]*100 + tagTypes[k][1]; // current tag
        //cout << "kkk   " << sortResVec[0][k].cellType << "  " << sortResVec[0][k].asocCortCell << endl;
        if(tagData[tag][0].cellType == i){
          for(int j = 0; j < tagData[tag].size(); j++){
            tagData[tag][j].GUS = dataVec[counter];
          }
          //cout << "a " << tag << "  " << tagData[tag][0].GUS << endl;
          counter++;
        }

      }
    }

  //}



  return true;
  }

  // filters high and low values out
  bool filterGUS(bool upperBool, bool lowerBool, bool upperAbs, bool lowerAbs, double upperFilterValue, double lowerFilterValue, tagDataMap& tagData, tagTypeVector& tagTypes, QString filename)
  {

  //if(!filterAbs){
  for(int i = 0; i < tagTypes.size(); i++){ // go through all tag types
    int tag = tagTypes[i][0]*100 + tagTypes[i][1]; // current tag
    double maxGUS = -1E20;
    double minGUS = 1E20;
    for(int j = 0; j < tagData[tag].size(); j++){
      if(filename == tagData[tag][j].file){ // check if right file
        double gus = tagData[tag][j].GUS;
        if(maxGUS < gus && gus != -1) maxGUS = gus;
        if(minGUS > gus && gus != -1) minGUS = gus;
      }
    }

    double upper = upperFilterValue;
    double lower = lowerFilterValue;
    if(!upperAbs) upper = upperFilterValue/100.0 * (maxGUS-minGUS);
    if(!lowerAbs) lower = lowerFilterValue/100.0 * (maxGUS-minGUS);

    //cout << "filter " << tag << "  " << upper << "    "  << lower << "  " << maxGUS << "   " << minGUS << endl;

    for(int j = 0; j < tagData[tag].size(); j++){
      if(filename == tagData[tag][j].file){ // check if right file
        if(upperBool){
          if(tagData[tag][j].GUS > upper) { //cout << "filter " << tag << "  " << upper << "    "  << lower << "  " << maxGUS << "   " << minGUS << endl; cout << "filtered up " << tag << "  " << tagData[tag][j].GUS << endl;
            tagData[tag][j].GUS = -1; }
        }
        if(lowerBool){
          if(tagData[tag][j].GUS < lower) { //cout << "filtered low " << tag << "  " << tagData[tag][j].GUS << endl; 
            tagData[tag][j].GUS = -1; }
        }
      }
    }


  }

  return true;
  }


  //////////////////////////////////////////////////////////
  // main method GUS
  bool DataAnalyzer::analyzeGUS(QString folderInput, /*QString folderOut, */QString gusFileExt, bool mergeBool, QString fileToMerge, bool avg, int window, QString upperFilter, double upperFilterValue, QString lowerFilter, double lowerFilterValue, IntDoubleMap& GUSInfo)
  {
  IntPoint3fMap averages;

  QString folderOut = folderInput;

  // go through all files in the folder and check whether theres is filename+gusFileExt
  // if yes, open heatmap file & open GUS file
  readCellDataGUS(folderInput, gusFileExt, dataX, tagTypesX, filenamesX);

  for(int i = 0; i < filenamesX.size(); i++){
    bool upper = false, lower = false, upperAbs = false, lowerAbs = false;
    if(upperFilter == "No Filter"){
      // do nothing
    } else if(upperFilter == "Percentage"){
      upper = true;
    } else if(upperFilter == "Value"){
      upper = true;
      upperAbs = true;
    }

    if(lowerFilter == "No Filter"){
      // do nothing
    } else if(lowerFilter == "Percentage"){
      lower = true;
    } else if(lowerFilter == "Value"){
      lower = true;
      lowerAbs = true;
    }

    if(upper || lower)
      filterGUS(upper, lower, upperAbs, lowerAbs, upperFilterValue, lowerFilterValue, dataX, tagTypesX,filenamesX[i]);


  }
  // processing
  processGUSData(averages, dataX, tagTypesX);

  // smooth if required
  if(avg) smoothGUS(dataX, tagTypesX, window);

  QString filename = "output.csv";
  QString filename2 = "output_CI.csv";
  
  // save the result to one file
  saveFileGUS(filename, filename2, folderOut, dataX, tagTypesX);

  // if required save to another file
  if(mergeBool) saveFileGUSExternal(folderOut, fileToMerge, dataX, tagTypesX, GUSInfo);

  return true;
  }

  //////////////////////////////////////////////////////////
  // main method
  bool DataAnalyzer::analyze(QString folderX, QString folderY, QString folderOut, bool treatment, bool avg, int window, bool perCell)
  {

  perCellAvg = perCell;

  // read in treatment
  readCellData(folderX, dataX, tagTypesX, filenamesX);
  // read in control
  readCellData(folderY, dataY, tagTypesY, filenamesY);

  // calculate the log ratios
  sortData(dataX, dataY, tagTypesX);

  // join and smooth the data and save it into MGX heatmap csv's
  if(!treatment){ // join it with control data
    if(avg) smoothAnalyzed(resVecAll, window, tagTypesX);
    joinData(resVecAll, dataX, tagTypesX, filenamesX, folderOut);
  } else { // join it with treatment data
    if(avg) smoothAnalyzed(resVecAll, window, tagTypesY);
    joinData(resVecAll, dataY, tagTypesY, filenamesY, folderOut);
  }

  // now to the E-Factors
  // calculate them & save the result according to the user selection
  if(!treatment){
    calcEFactors(resVecAll, avg, window, tagTypesX, efData);
    saveEFactors(efData, dataX, tagTypesX, filenamesX, folderOut);
  } else {
    calcEFactors(resVecAll, avg, window, tagTypesY, efData);
    saveEFactors(efData, dataY, tagTypesY, filenamesY, folderOut);
  }



  return true;
}
//}
}
