//
// This file is part of MorphoGraphX - http://www.MorphoGraphX.org
// Copyright (C) 2012-2015 Richard S. Smith and collaborators.
//
// If you use MorphoGraphX in your work, please cite:
//   http://dx.doi.org/10.7554/eLife.05864
//
// MorphoGraphX is free software, and is licensed under under the terms of the
// GNU General (GPL) Public License version 2.0, http://www.gnu.org/licenses.
//
#include <TypeRecognition.hpp>
#include <CellCluster.hpp>
#include <GraphUtils.hpp>
#include <MeshProcessHeatMap.hpp>
//#include <CellAtlasUtils.hpp>
#include <Triangulate.hpp>
#include <math.h>
#include <iostream>
#include <vector>
#include <fstream>
#include <limits>
#include <Progress.hpp>
#include <QFileDialog>
#include <QFileDialog>

namespace mgx
{

    QStringList getMeasureAttrList(Mesh *m)
    {
    Attributes *attributes;
    attributes = &m->attributes();
    QStringList attr = attributes->getAttrList();
    QStringList measureAttr;

    forall(const QString &text, attr) {
      QStringList list = text.split(" ");

      if(list.size() < 4) continue;
      if(list[0] != "Measure") continue;
      if(list[1] != "Label") continue;
      if(list[2] != "Double") continue;
      QString name;
      for(int i = 3; i<list.size(); i++){
        name = name + list[i] + " ";
      }
      measureAttr << name;
    }

    return measureAttr;
    }



    // parses the string used to save the selected attr map measures/features (String: processName1$processName2$...)
    // returns a set of QStrings with all attr map processNames as elements
    std::set<QString> parseMeasureString(QString parmsString)
    {
      //qDebug() << "parse\n";
      std::set<QString> selectedMeasuresNew;
      QStringList measures = parmsString.split("$");
      for(int i = 0; i < measures.size(); i++){
        if(measures[i] != "")
          selectedMeasuresNew.insert(measures[i]);
      }
      return selectedMeasuresNew;

    }


    // scale values of a heatmap to 0...1 [namrata: (value - mean)/stdev]
    bool scaleAttrMap(IntFloatAttr &measureMap, bool namrata, Point2d& minMax){    //Scale features
        double min = HUGE_VAL, max = -HUGE_VAL, mean = 0, standDev = 0;
        //Calculate mean
        forall(IntFloatP loop1, measureMap){
          if(loop1.first < 1) continue;
          mean += loop1.second;
          if(max < loop1.second)
            max = loop1.second;
          if(min > loop1.second)
            min = loop1.second;

          //std::cout << "value " << loop1.first << "/" << loop1.second << std::endl;
        }
        mean /= measureMap.size();

        if(minMax.x() < minMax.y()){
          min = minMax.x();
          max = minMax.y();
        }
        minMax.x() = min;
        minMax.y() = max;

        // calc stdev
        forall(IntFloatP loop1, measureMap){
          standDev += pow((loop1.second - mean), 2);
        }
        standDev /= measureMap.size();
        if(standDev > 0) standDev = pow(standDev, 0.5);

        std::cout << "scaled " << mean << "/ " << standDev << "/" << min << "/" << max << "/" << norm(max-min) << std::endl;

        for(auto it = measureMap.begin(); it!=measureMap.end(); it++){
        //forall(IntFloatP p, measureMap){ //Scale features with respect to Std dev]

          if(namrata){
            if(min == max or standDev < 1E-5)
              it->second = 1;
            else
              it->second = (it->second - mean)/standDev;
          } else {
            if(norm(max-min) < 1E-5) it->second = 1;
            else it->second = (it->second - min)/(max-min);
          }
          //std::cout << " new " << it->second;
        }

        return true;
    }



  // loop over all features, scale every one
  bool CellTypeRecognitionFeatureSelection::scaleData(){
      //forall(const QString &name, cc.SelectedFeaturesName){
      forall(const QString &name, selectedMeasures){
        Point2d minMax(0,0);
        scaleAttrMap(cc.measures[name],false, minMax);
          //qDebug()<<"\nscaled\t"<<name;
      }
      return true;
  }

  void CellTypeRecognitionFeatureSelection::populateMapTree(Mesh *mesh)
  {

    ui.mapTree->header()->resizeSection(1, 50);
    ui.mapTree->header()->resizeSection(0, 300);

    CellTypeRecognitionFeatureSelection *cellTypeMeasures;
    if(!getProcess("Mesh/Cell Types/Classification/A Select Measures", cellTypeMeasures))
      throw(QString("CellTypeRecognitionFeatureSelection:: Unable to make CellTypeRecognitionFeatureSelection process"));

    //QStringList processParms;
    //if(!getProcessParms<CellTypeRecognitionFeatureSelection>(this, processParms)) throw(QString("Unable to get Parameters"));


    QStringList attr = getMeasureAttrList(mesh);

    ui.mapTree->clear();
    //fillTreeWidgetWithMeasures(mesh,ui.mapTree);
    forall(const QString &text, attr) {
      QTreeWidgetItem *item = new QTreeWidgetItem(QStringList() << text);
      item->setFlags(Qt::ItemIsEnabled | Qt::ItemIsSelectable | Qt::ItemIsUserCheckable);
      item->setCheckState(1, Qt::Unchecked);
      item->setCheckState(2, Qt::Unchecked);
      ui.mapTree->addTopLevelItem(item);
    }

    for(int itemCount = 0; itemCount < ui.mapTree->topLevelItemCount(); ++itemCount){
      QString name = ui.mapTree->topLevelItem(itemCount)->text(0);
      bool found = false;

      forall(QString q, selectedMeasures){
        if(q==name.trimmed()) found = true;
      }
      if(found) ui.mapTree->topLevelItem(itemCount)->setCheckState(1, Qt::Checked);
    }

  }

  void CellTypeRecognitionFeatureSelection::changeMeshType(const QString& s)
  {
    if(s == "Area 2D") cc.folder = "Mesh/Heat Map/Measures/";
    else if(s == "Volume 3D") cc.folder = "Mesh/Heat Map/Measures 3D/";

    int sizeProcessTree = ui.processTree->topLevelItemCount();

    for(int i = 0; i < sizeProcessTree; i++){
      QTreeWidgetItem *item = ui.processTree->takeTopLevelItem(0);
      if(item)delete item;
    }

    findProcesses(currentMesh());
  }



    //Update and sort the feature box
    bool CellTypeRecognitionFeatureSelection::findProcesses(Mesh *mesh){
        QStringList procs = mgx::listProcesses();
        //int counter = 0;

        ui.mapTree->header()->resizeSection(1, 50);
        ui.mapTree->header()->resizeSection(0, 300);

        populateMapTree(mesh);

        forall(const QString& name, procs) {
          mgx::ProcessDefinition* def = mgx::getProcessDefinition(name);
          QStringList list = def->name.split("/");

          //if(ui.comboBoxMeshType)

          QStringList folder = cc.folder.split("/");

          QString subfolder = folder[2];

          if(list[0] == "Mesh" and list[1] == "Heat Map" and list[2] == subfolder){
            QTreeWidgetItem *item2 = new QTreeWidgetItem(QStringList() << list[3] + "/" + list[4]);
            item2->setFlags(Qt::ItemIsEnabled | Qt::ItemIsSelectable | Qt::ItemIsUserCheckable);
            ui.processTree->addTopLevelItem(item2);
          }
        }

        ui.mapTree->sortItems(0, Qt::AscendingOrder);
        ui.processTree->sortItems(0, Qt::AscendingOrder);
        for(int i = 0; i < 1; i++)
          ui.mapTree->resizeColumnToContents(i);

        return true;
    }

    //Populate the window with features
    bool CellTypeRecognitionFeatureSelection::initialize(QWidget *parent){
      //if(parms[0] == "2D")
      cc.folder = "Mesh/Heat Map/Measures/";
      //else if(parms[0] == "3D") cc.folder = "Mesh/Heat Map/Measures 3D";

      QString oldParms = parm("Selected Measures");

      Mesh *mesh = currentMesh();
      mesh->showHeat();
      QDialog dlg(parent);
      ui.setupUi(&dlg);

      // go through all available processes and fill the GUI with all relevant measures
      findProcesses(mesh);

      connect(ui.processTree, SIGNAL(itemClicked(QTreeWidgetItem *, int)),
                this, SLOT(on_svmTreeWidget_itemClicked(QTreeWidgetItem *, int)));
      connect(ui.mapTree, SIGNAL(itemClicked(QTreeWidgetItem *, int)),
                this, SLOT(on_svmMapTreeWidget_itemClicked(QTreeWidgetItem *, int)));
      connect(ui.selectAllButton, SIGNAL(clicked()), this, SLOT(on_selectAllButton_clicked()));
      connect(ui.unselectAllButton, SIGNAL(clicked()), this, SLOT(on_unselectAllButton_clicked()));
      connect(ui.comboBoxMeshType, SIGNAL(currentIndexChanged(const QString &)), this, SLOT(changeMeshType(const QString &)));
//qDebug() << "a " << parm("Selected Measures")<< "\n";
      if(dlg.exec() == QDialog::Accepted){
//qDebug() << "h " << parm("Selected Measures")<< "\n";
        // save the selected features
        setParm("Selected Measures", readSelectedMeasures());
        selectedMeasures = parseMeasureString(parm("Selected Measures"));
        //qDebug() << "again " << parm("Selected Measures")<< "\n";

        // save the heatmaps of the selected features
        //bool forceRecalc = ui.checkBoxRecalcMeasures->isChecked();
        bool forceRecalc = ui.checkBoxRecalcMeasures->isChecked();
        bool setScale = ui.checkBoxSetScale->isChecked();
        calcSelectedMeasures(mesh, forceRecalc,setScale);
      } else {
        setParm("Selected Measures", oldParms);
      }
      return true;
    }

    bool CellTypeRecognitionFeatureSelection::run(){
      return true;
    }

    // generates or grabs an attribute map of a measure
    bool CellTypeRecognitionFeatureSelection::runMeasureProcess(Mesh* m, QString name, bool forceRecalc, bool justGrapAttrMap)
    {
        QString folder = cc.folder;
        QString proName = folder + name;
        QString proAttr = "Measure Label Double " + name;
        AttrMap<int, double>& measureMap = m->attributes().attrMap<int, double>(proAttr);
        qDebug() << "Found Attr Map: " << proAttr << " with size: " << measureMap.size() << "\n";

        // call the measure process if map empty
        if((measureMap.empty() or forceRecalc) and !justGrapAttrMap){
            QStringList parms;
            Process *pro = makeProcess(proName);

            // make sure Process exists
            if(pro){
              qDebug() << "Process call " << proName << "\n";
              //getLastParms(proName);
              pro->run();
            }
        // copy the map if it exitst
        } else {

          if(measureMap.empty()) return false;
          typedef std::pair<int,double> IntDouble;
          forall(IntDouble p, measureMap){
            m->labelHeat()[p.first] = p.second;
          }
        }

        m->heatMapBounds() = m->calcHeatMapBounds();
        m->showHeat();
        m->updateTriangles();

        return true;
    }

    //Update heat map for the selected map
    void CellTypeRecognitionFeatureSelection::on_svmTreeWidget_itemClicked(QTreeWidgetItem *item, int column){
      qDebug()<<" processTree ";
        Mesh *mesh = currentMesh();
        qDebug()<<ui.processTree->currentItem()->text(0);
        QString name = ui.processTree->currentItem()->text(0);
        name = name.trimmed();
        //QString prName = "Mesh/Heat Map/Measures/";
        //if(meshType3D) prName = "Mesh/Heat Map/Measures 3D/";
        bool forceRecalc = ui.checkBoxRecalcMeasures->isChecked();
        runMeasureProcess(mesh, name, forceRecalc, false);

        selectedMeasures = parseMeasureString(readSelectedMeasures());
        qDebug()<<" out ";
        populateMapTree(mesh);

        updateState();
        updateViewer();
    }

    //Update heat map for the selected map
    void CellTypeRecognitionFeatureSelection::on_svmMapTreeWidget_itemClicked(QTreeWidgetItem *item, int column){
        qDebug()<<" mapTree ";
        if(ui.mapTree->currentItem() == NULL) return; // in case there is no item selected but just the box ticked
        Mesh *mesh = currentMesh();
        qDebug()<<ui.mapTree->currentItem()->text(0);
        QString name = ui.mapTree->currentItem()->text(0);
        name = name.trimmed();
        //QString prName = "Mesh/Heat Map/Measures/";
        //if(meshType3D) prName = "Mesh/Heat Map/Measures 3D/";
        bool forceRecalc = ui.checkBoxRecalcMeasures->isChecked();
        runMeasureProcess(mesh, name, forceRecalc, false);

        updateState();
        updateViewer();
    }


    //Select all the features if Select all button is clicked
    void CellTypeRecognitionFeatureSelection::on_selectAllButton_clicked(){
        for(int itemCount = 0; itemCount < ui.mapTree->topLevelItemCount(); ++itemCount)
          ui.mapTree->topLevelItem(itemCount)->setCheckState(1, Qt::Checked);
    }
    //Unselect all the features if Select all button is clicked
    void CellTypeRecognitionFeatureSelection::on_unselectAllButton_clicked(){
        for(int itemCount = 0; itemCount < ui.mapTree->topLevelItemCount(); ++itemCount)
          ui.mapTree->topLevelItem(itemCount)->setCheckState(1, Qt::Unchecked);
    }


    bool CellTypeRecognitionFeatureSelection::calcSelectedMeasures(Mesh *mesh, bool forceRecalc, bool setScale){
        cc.measures.clear();
        // run all measure processes and populate cc data
        forall(const QString &name, selectedMeasures /*cc.SelectedFeaturesName*/){
          qDebug() << "runProcess " << name ;

          if(!runMeasureProcess(mesh, name, forceRecalc, true)){
            std::cout << "missing Attr Map: " << name.toUtf8().constData() << std::endl;
            return false;
          }

          cc.measures[name] = mesh->labelHeat();
          Point2d minMax(0,0);
          if(!setScale) minMax = cc.scales[name];
          std::cout << "scale " << minMax << std::endl;
          scaleAttrMap(cc.measures[name], false, minMax);
          qDebug() << "Scaled " << cc.measures[name].size();
          if(setScale){
            std::cout << "scale set " << minMax << std::endl;
            cc.scales[name] = minMax;
          }

        }
        //std::cout<<"\nDebug: "<<__LINE__<<"\n";
        return true;
    }

    // save the selected features and scale their values
    QString CellTypeRecognitionFeatureSelection::readSelectedMeasures(){
      qDebug() << "select features" << "\n";
      //Create a vector for selected features
      QString newFeatures;
      for(int itemCount = 0; itemCount < ui.mapTree->topLevelItemCount(); ++itemCount){
        if(Qt::Checked == ui.mapTree->topLevelItem(itemCount)->checkState(1)){
          QString name = ui.mapTree->topLevelItem(itemCount)->text(0);
          name = name.trimmed();
          newFeatures += name + "$";
          qDebug() << "selected Feature: " << name << "\n";
        }
        //qDebug() << newFeatures << "/" << parms.size() << "\n";

      }
      // overwrite selected Measures
      return newFeatures;
      //selectedMeasures = parseMeasureString(newFeatures);
    }

    void CellTypeRecognitionFeatureSelection::setSelectedMeasures(std::set<QString> newMeasures)
    {
      selectedMeasures = newMeasures;
      std::cout << "measures " << std::endl;
    }

    REGISTER_PROCESS(CellTypeRecognitionFeatureSelection);

    REGISTER_PROCESS(CellTypeRecognitionSpecification);


    //Clear cc data structure and/or model for new test data
    bool CellTypeRecognitionClear::run(bool clearInput, bool clearModel){
      if(clearInput){

	    CellTypeRecognitionFeatureSelection	*cellTypeMeasures;
        if(!getProcess("Mesh/Cell Types/Classification/A Select Measures", cellTypeMeasures))
          throw(QString("CellTypeRecognitionClear:: Unable to make CellTypeRecognitionFeatureSelection process"));
        cellTypeMeasures->setParm("Selected Measures", "");

         cc.measures.clear();
        //AttrMap<int, double>& measureMap = m->attributes().attrMap<int, double>(proAttr);
      }
      if(clearModel){
        if(model)
          svm_free_and_destroy_model(&model);     //Free the model space and param
      }
      return true;
    }
    REGISTER_PROCESS(CellTypeRecognitionClear);



    bool CellTypeRecognitionSimilarity::run(Mesh* m, Mesh* m2, QString distanceType, QString meshNr){

      // grab the selected measures from the parameter in the feature selection process
      std::set<QString> fileFeatures;

      CellTypeRecognitionFeatureSelection *cellTypeMeasures;
      if(!getProcess("Mesh/Cell Types/Classification/A Select Measures", cellTypeMeasures))
      throw(QString("CellTypeRecognitionSimilarity:: Unable to make CellTypeRecognitionFeatureSelection process"));

      fileFeatures = parseMeasureString(cellTypeMeasures->parm("Selected Measures"));

      // now create missing attribute maps using the feature selection process
      CellTypeRecognitionFeatureSelection fs(*this);
      fs.setSelectedMeasures(fileFeatures);
      std::cout << "simi " << std::endl;
      if(meshNr == "other") if(!fs.calcSelectedMeasures(m2, false, false)) return setErrorMessage("Not all required Attr Maps in other mesh present!");
      if(!fs.calcSelectedMeasures(m, false, false)) return setErrorMessage("Not all required Attr Maps present!");

      //auto &data = m->attributes().attrMap<int, QStringList>("Cell Type Recognition Selected Measures");
      //auto &data2 = m2->attributes().attrMap<int, QStringList>("Cell Type Recognition Selected Measures");

      AttrMap<int, QStringList> data, data2;

      data.clear();
      data2.clear();

      forall(QString s, fileFeatures){
        QString sExt = "Measure Label Double " + s;
        data[0] << sExt;
        if(meshNr == "other") data2[0] << sExt;

        qDebug() << "bla " << sExt << "\n";
      }

      std::set<int> selectedCells;
      std::set<int> allLabels;

      // find selected cells
      forall(const vertex& v, m->graph()){
        if(v->label == -1) continue;
        if(v->selected) selectedCells.insert(v->label);
        allLabels.insert(v->label);
      }

      if(meshNr == "other"){
        allLabels.clear();
        forall(const vertex& v, m2->graph()){
          if(v->label == -1) continue;
          allLabels.insert(v->label);
        }
      }

      if(selectedCells.empty()) return setErrorMessage("No cell selected!");

      std::map<QString, double> selectedValues;
      //std::map<int, double> cellDistance;

      // calc average values of selected cells
      forall(QString feat, data[0]){
        auto attrMapData = m->attributes().attrMap<int, double>(feat);
        forall(int l, selectedCells){
          selectedValues[feat] += attrMapData[l];//cc.measures[feat][l];
        }
        selectedValues[feat]/=selectedCells.size();
      }


      // calc distance to all other cells
      forall(int l, allLabels){
        double totalDis = 0;
        forall(QString feat, data[0]){

          auto attrMapData = m->attributes().attrMap<int, double>(feat);
          if(meshNr == "other"){
            attrMapData = m2->attributes().attrMap<int, double>(feat);
          }

          //double dis = norm(cc.measures[feat][l]-selectedValues[feat]);
          double dis = norm(attrMapData[l]-selectedValues[feat]);
          if(distanceType == "Euclidean") dis*=dis;
          totalDis+=dis;

        }
        if(distanceType == "Euclidean") totalDis=std::sqrt(totalDis);
        if(meshNr == "other"){
          m2->labelHeat()[l] = totalDis/fileFeatures.size();
        } else {
          m->labelHeat()[l] = totalDis/fileFeatures.size();
        }

      }
      if(meshNr == "other"){
        m2->setShowLabel("Label Heat");
        m2->heatMapBounds() = m->calcHeatMapBounds();
        m2->updateTriangles();
      } else {
        m->setShowLabel("Label Heat");
        m->heatMapBounds() = m->calcHeatMapBounds();
        m->updateTriangles();
      }

      return true;
    }

    REGISTER_PROCESS(CellTypeRecognitionSimilarity);



    bool CorrespondingVertices::run(Mesh* m, Mesh* m2){

      std::vector<vertex> selectedVtxs = m->selectedVertices();

      forall(vertex v, m2->graph()){
        v->selected = false;
      }


      forall(vertex v, selectedVtxs){

      }


      return true;
    }

    REGISTER_PROCESS(CorrespondingVertices);



    bool CellTypeRecognitionDataFile::writeDataFile(Mesh *mesh, QString filename, QString choice){

      vvGraph &S = mesh->graph();
      std::vector<int> allLabels = findAllLabels(S);


      std::set<QString> fileFeatures;

      if(filename.isEmpty())
          return false;
      if(!filename.endsWith(".txt", Qt::CaseInsensitive))
          filename += ".txt";
      QFile file(filename);
      // first get the selected measures
      if(choice == "New"){   //Depending on choice create or append file
          if(!file.open(QIODevice::WriteOnly))
              std::cout<<"Cannot Create a new file";

        // grab the selected measures from the parameter in the feature selection process
        CellTypeRecognitionFeatureSelection *cellTypeMeasures;
        if(!getProcess("Mesh/Cell Types/Classification/A Select Measures", cellTypeMeasures))
          throw(QString("MeshProcessHeatMap:: Unable to make CellTypeRecognitionFeatureSelection process"));
        fileFeatures = cellTypeMeasures->selectedMeasures ; //parseMeasureString(cellTypeMeasures->parm("Selected Measures"));
        qDebug() << "selected measures " << fileFeatures.size();
      } else { // Append
          if(!file.open(QIODevice::ReadOnly))
              std::cout<<"Cannot append to the file";
          //cc.SelectedFeaturesName.clear();
          // grab the selected measures from the existing file
          QTextStream in(&file);
          int size = in.readLine().toInt();
          //std::cout<<"\nsize of features\t" << size;
          if(size){
            while(size > 0){     //Define the number of rows in the training data file
                QString line = in.readLine();
                //cc.SelectedFeaturesName.insert(line);
                fileFeatures.insert(line);
                size--;
            }
            forall(const QString &name, fileFeatures/*cc.SelectedFeaturesName*/){
                qDebug() << "\n Append feature " << name << "\n";
            }

          }
          file.close();
          if(!file.open(QIODevice::WriteOnly | QIODevice::Append))
              std::cout<<"Cannot append to the file";

      }

      // if no features exit
      if(fileFeatures.empty()){
        file.close();
        return false;
      }

      // now create missing attribute maps using the feature selection process
      CellTypeRecognitionFeatureSelection fs(*this);
      fs.setSelectedMeasures(fileFeatures);
      if(!fs.calcSelectedMeasures(mesh, false, false)) return setErrorMessage("Not all required Attr Maps present!");

      // if new file write the selected measure names at the top
      QTextStream out(&file);
      if(choice == "New"){
          out << fileFeatures.size() << endl;      //Write the row size as the first line
          forall(const QString &name, fileFeatures){
              out << name << endl;
          }
      }


      // now go through all selected cells and write the values
      forall(int label, allLabels){
        int i = 0;
        if(label < 1) continue;
        //if(!QString::compare(trainTest, "Train", Qt::CaseSensitive)){       //Training data format
            if(mesh->parents()[label]){
              out << mesh->parents()[label] <<" ";
              forall(const QString &name, fileFeatures){
                //qDebug() << "write " << name << "/" << cc.measures[name].size() << "\n";
                forall(IntFloatP loop, cc.measures[name]){
                  //qDebug() << "values " << loop.first << "/" << loop.second << "\n";
                  if(loop.first == label and !std::isnan(loop.second)){
                    out << i << ":"<< loop.second <<" ";
                    i++;
                  }
                }
              }
            out << i << ":"<< "-1";
            out << endl;
            }
          //} else {      //Testing data format
          /*  out << v->label <<" ";
            forall(const QString &name, cc.SelectedFeaturesName){
              forall(IntFloatP loop, cc.measures[name]){
                if(loop.first == v->label and !isnan(loop.second)){
                  out << i << ":"<< loop.second <<" ";
                  i++;
                }
              }
            }
            out << i << ":"<< "-1" <<" ";
            out << endl;   */
          //}
      }
      file.close();

      return true;
    }


  bool CellTypeRecognitionDataFile::initialize(QWidget* parent)
  {
    QString filename = parm("Filename");
    if(parent)
      filename = QFileDialog::getSaveFileName(0, "Choose txt file to save or append", QDir::currentPath(), "txt files (*.txt)");
    if(filename.isEmpty())
      return false;
    if(!filename.endsWith(".txt", Qt::CaseInsensitive))
      filename += ".txt";
    setParm("Filename", filename);
    return true;
  }


    bool CellTypeRecognitionDataFile::run(Mesh *mesh, QString filename, QString choice){
        writeDataFile(mesh, filename, choice);   //Prepare data in libsvm format
        return true;
    }
    REGISTER_PROCESS(CellTypeRecognitionDataFile);

/*
    double CellTypeRecognitionTrain::do_cross_validation() {
        double cross_validation = 0.;
        int total_correct = 0;
        double total_error = 0;
        double sumv = 0, sumy = 0, sumvv = 0, sumyy = 0, sumvy = 0;
        double *target = Malloc(double, prob.l);
        int nr_fold = 10;
        svm_cross_validation(&prob, &param, nr_fold, target);   //Call to cross validation function in SVM.cpp. Refer SVM.cpp in libsvm
        if(param.svm_type == EPSILON_SVR || param.svm_type == NU_SVR) {
            for(int i=0;i<prob.l;i++) {
                double y = prob.y[i];
                double v = target[i];
                total_error += (v-y)*(v-y);
                sumv += v;
                sumy += y;
                sumvv += v*v;
                sumyy += y*y;
                sumvy += v*y;
            }
            std::cout<<"Cross Validation Mean squared error = " << total_error/prob.l << std::endl;
            std::cout<<"Cross Validation Squared correlation coefficient = " <<  ((prob.l*sumvy-sumv*sumy)*(prob.l*sumvy-sumv*sumy))/
            ((prob.l*sumvv-sumv*sumv)*(prob.l*sumyy-sumy*sumy)) << std::endl;
        } else {   //SVM type is C_SVC so this part is processed. Refer SVM.cpp in libsvm
            for(int i=0;i<prob.l;i++){
                if(target[i] == prob.y[i])
                    ++total_correct;
            }
            cross_validation = (double)total_correct/double(prob.l);
            std::cout<<"Cross Validation Accuracy = "<<(100.0*total_correct/prob.l) << " for gamma= " << param.gamma << " /C= " << param.C << std::endl;
        }
        free(target);
        return cross_validation;    //Return the cross validation accruacy
    }
*/
    /* TODO check this function
    bool CellTypeRecognitionTrain::backwardSelection(QString filename){
        double cv = 0;
        QFile file(filename);
        if(!file.open(QIODevice::ReadOnly))
        {
            return false;  //TODO: I need to exit with a proper msg
        }
        QTextStream in(&file);
        int rowSize = 0;
        QStringList fieldsx;

        for(uint i = 0; i < fileFeatures.size(); i++){
            QString line = in.readLine();
        }
        while(!in.atEnd()){     //Define the number of rows in the training data file
            QString line = in.readLine();
            fieldsx = line.split(" ");
            rowSize++;
        }
        file.close();
        std::cout<<__LINE__<<"\n";
        int trainFeatureSize = fieldsx.size() - 3;      //You want to remove -1 at the end as you do not want to check it
        std::cout<<"\ntrainFeatureSize\t"<<trainFeatureSize <<"\t"<< fieldsx.size();
        qDebug()<<fieldsx;
        while(trainFeatureSize >= 0){
            prob.l = rowSize;
            prob.y = Malloc(double, prob.l);
            prob.x = Malloc(struct svm_node *, prob.l);
            //svm_node** x = Malloc(svm_node*, prob.l);
            int i = 0;
            bool ok;
            if(!file.open(QIODevice::ReadOnly))
            {
                return false;
            }
            while(!in.atEnd()){     //Read the training data file
                QString line = in.readLine();
                QStringList fields = line.split(" ");
                if(!line.size())
                    continue;
                int j = 0, flag = 0;
                x_space = Malloc(struct svm_node, prob.l);
                prob.x[i] = &x_space[j];        //Prepare training data vectors
                prob.y[i] = fields[0].toInt(&ok);       //Mention labels of the training data

                for(int p = 1; p < fields.size() - 1; p++){
                    if(p - 1 == trainFeatureSize){
                        std::cout<<"\nNot taking\t"<<p - 1<<"\t"<<trainFeatureSize;
                        continue;
                    }
                    QStringList val = fields[p].split(":");
                    x_space[j].index = val[0].toInt(&ok);
                    x_space[j].value = val[1].toDouble(&ok);
                    ++j;
                    flag = 1;
                }
                if(flag == 1){
                    x_space[j++].index = -1;
                    i++;
                }
            }
            double cross_validation = do_cross_validation();       //Use cross validation to find the best parameters
            if(!cv){
                cv = cross_validation;
                std::cout<<"\nTotal features\t"<<trainFeatureSize <<"\tcv"<< cv;
            }
            else if(cv < cross_validation){
                cv = cross_validation;
                cc.eliFtrs.insert(trainFeatureSize);
                std::cout<<"\nEliminated feature number\t"<<trainFeatureSize <<"\tcv"<< cv;
            }
            file.close();       //Close the file
            std::cout<<"\ncross_validation cv cross_validation\t"<<cv<<"\t" << cross_validation;
            std::cout<<"\ntrainFeatureSize\t"<<trainFeatureSize;
            trainFeatureSize--;     //trainFeatureSize is fields.size - 1 which will not be reached by p in 1st iteration.


        }
        std::cout<<"\n\n\n\n\n";
        for (std::set<int>::iterator it = cc.eliFtrs.begin(); it != cc.eliFtrs.end(); ++it){
                std::cout<<"\nftr\t"<<*it;
        }
        std::cout<<"\n\n\n\n\n";
        return 0;

    }*/

  bool CellTypeRecognitionTrain::initialize(QWidget* parent)
  {
    QString filename = parm("Filename");
    if(parent)
      filename = QFileDialog::getOpenFileName(0, "Choose txt file to load", QDir::currentPath(), "txt files (*.txt)");
    if(filename.isEmpty())
      return false;
    if(!filename.endsWith(".txt", Qt::CaseInsensitive))
      filename += ".txt";
    setParm("Filename", filename);
    return true;
  }
/*
  void CellTypeRecognitionTrain::optimizeSVMParameter()
  {

        // Tune the SVM model by changing C and gamma
        double cross_validation = 0;
        double Optimized_C = HUGE_VAL, max_cv = 0, Optimized_gamma = 0;
        //if(optimizeParms){
          for(double c = 1; c < 10000; c*= 4){
            for(double gamma = 0.1; gamma < 5; gamma*= 3){
              param.C = c;
              param.gamma = gamma;
              std::cout << "cross val " << std::endl;
              cross_validation = do_cross_validation();       //Use cross validation to find the best parameters
              if(cross_validation > max_cv){
                max_cv = cross_validation;
                Optimized_C = param.C;      //Select the optimized C and gamma
                Optimized_gamma = param.gamma;
              }
            }
          }
          if(Optimized_C > 0){        //If the model is optimized
            param.C = Optimized_C;
            param.gamma = Optimized_gamma;
          }
        //} else { // dont optimize, use user specified values
        //}
    std::cout << "Final C/Gamma//Cross validation "<< param.C << " / "<< param.gamma <<" // "<< max_cv << std::endl;
    if(max_cv < 0.6){
      std::cout << "Warning: the cross validation accuracy is very low. Consider re-training!" << std::endl;
    }
  }
*/
    bool CellTypeRecognitionTrain::run(Mesh *mesh, QString filename, bool optimizeParms, double pGamma, double pC,
      QString svmType, QString kernelType, bool featureReduction, double featureThreshold, int optGridSize, bool fullParameterSearch,
      bool featureInclusion, bool genetic, int kFold) {

        std::cout << "Train Debug: " << __LINE__ << std::endl;

        SVMClassifier svmcNew;

        svmc = svmcNew;


        svmc.initialize(svmType, kernelType, pGamma, pC, featureThreshold, kFold);

        const char *error_msg;

        CellTypeRecognitionFeatureSelection *cellTypeMeasures;
        if(!getProcess("Mesh/Cell Types/Classification/A Select Measures", cellTypeMeasures))
        throw(QString("MeshProcessHeatMap:: Unable to make CellTypeRecognitionFeatureSelection process"));
        fileFeatures = cellTypeMeasures->selectedMeasures; //parseMeasureString(cellTypeMeasures->parm("Selected Measures"));

        FileHandlerSVM fhSVM;
        FileData fd;

        fhSVM.readTrainingFile(fd,filename);
        svmc.allData = fd;
        svmc.optGridSize = optGridSize;

        svmc.generateDataSVM(fd, svmc.prob, fileFeatures, svmc.x_space);

        //error_msg = svm_check_parameter(&prob, &param);
        error_msg = svm_check_parameter(&svmc.prob, &svmc.param);
        if(error_msg) std::cout<<"ERROR: \n" << *error_msg;
        std::cout << "Parms checked" << __LINE__ << std::endl;

        param.C = pC;
        param.gamma = pGamma;

        if(featureReduction){
          // do feature reduction
          std::cout << "REDUCTION " << std::endl;

          std::set<QString> reducedFeatures = svmc.featureReduction(fileFeatures, fullParameterSearch);
          std::cout << "done nr features " << reducedFeatures.size();

          svm_problem probTest;
          svm_node *x_spaceTest;

          svmc.generateDataSVM(fd, probTest, reducedFeatures, x_spaceTest);

          svmc.prob = probTest;
          svmc.x_space = x_spaceTest;



          // if nr of features has changed to parms optimization again
          if(optimizeParms){
            svmc.optimizeSVMParameter(svmc.prob, svmc.param);
          }
          qDebug() << "final features:";
          forall(QString q, reducedFeatures){
            qDebug() << q;
          }
        } else if(featureInclusion){
          // do feature inclusion
          std::cout << "INCLUSION " << std::endl;

          std::set<QString> includedFeatures = svmc.featureInclusion(fileFeatures, fullParameterSearch);
          std::cout << "done nr features " << includedFeatures.size();

          svm_problem probTest;
          svm_node *x_spaceTest;

          svmc.generateDataSVM(fd, probTest, includedFeatures, x_spaceTest);

          svmc.prob = probTest;
          svmc.x_space = x_spaceTest;

          // if nr of features has changed to parms optimization again
          if(optimizeParms){
            svmc.optimizeSVMParameter(svmc.prob, svmc.param);
          }
          qDebug() << "final features:";
          forall(QString q, includedFeatures){
            qDebug() << q;
          }
        } else if(genetic){
          std::cout << "GENETIC " << std::endl;
          std::set<QString> includedFeatures = svmc.geneticFeatureSelection(fileFeatures);
          std::cout << "done nr features " << includedFeatures.size();

          svm_problem probTest;
          svm_node *x_spaceTest;

          svmc.generateDataSVM(fd, probTest, includedFeatures, x_spaceTest);

          svmc.prob = probTest;
          svmc.x_space = x_spaceTest;
          if(optimizeParms){
            svmc.optimizeSVMParameter(svmc.prob, svmc.param);
          }
          qDebug() << "final features:";
          forall(QString q, includedFeatures){
            qDebug() << q;
          }
        } else {
          std::cout << "SELECTED " << std::endl;
          if(optimizeParms){
            svmc.optimizeSVMParameter(svmc.prob, svmc.param);
          }
        }




std::cout << "Parms optimized" << __LINE__ << std::endl;
            //backwardSelection(cc, filename);
            //std::cout<<"\n\n\n\n\n";
        //for (std::set<int>::iterator it = cc.eliFtrs.begin(); it != cc.eliFtrs.end(); ++it){
        //  std::cout<<"\nftr\t"<<*it;
        //}
          //std::cout<<"\n\n\n\n\n";
          //model = svm_train(&prob, &param);       //Train the model with optimized parameters
          model = svm_train(&svmc.prob, &svmc.param);
          svm_save_model((filename + ".model").toUtf8().constData(), model);      //Save the model as a file to use later separately for test data, if needed
        //}
          std::cout << "Trained" << __LINE__ << std::endl;

        //for (std::set<int>::iterator it = cc.eliFtrs.begin(); it != cc.eliFtrs.end(); ++it){
        //        std::cout<<"\nftr\t"<<*it;
        //}

        // free the memory
        //svmc.freeMemory();
   /*     if(prob.y)
          free(prob.y);
        if(prob.x)
          free(prob.x);
        if(x_space)
          free(x_space);
        if(svmc.param.C){
          //svm_destroy_param(&param);
          svm_destroy_param(&svmc.param);
        }*/
std::cout << "Freed" << __LINE__ << std::endl;
        return true;
    }
    REGISTER_PROCESS(CellTypeRecognitionTrain);

  bool LoadModelSVM::initialize(QWidget* parent)
  {
    QString filename = parm("Filename");
    if(filename.isEmpty() and parent)
      filename = QFileDialog::getOpenFileName(0, "Choose txt file to load", QDir::currentPath(), "libsvm model files (*.model)");
    if(filename.isEmpty())
      return false;
    if(!filename.endsWith(".model", Qt::CaseInsensitive))
      filename += ".model";
    setParm("Filename", filename);
    return true;
  }


    bool LoadModelSVM::run(QString filename) {

    model = svm_load_model(filename.toUtf8().constData());
    return true;
    }
    REGISTER_PROCESS(LoadModelSVM);

    bool CellTypeRecognitionClassification::run(Mesh *mesh, double classifyThreshold, bool forceRecalc) {

      if(!model){
        return setErrorMessage("No SVM model initialized! Load an existing model or train a new one!");
      }
//std::cout << "oha  " << param.C << "/" << param.cache_size << std::endl;
        vvGraph &S = mesh->graph();
        std::vector<int> allLabels = findAllLabels(S);
        //TODO: I need to change parent labels on GUI and save the entire data as training data
        // I need a pop up window to specify test data and training model

        // grab the selected measures from the parameter in the feature selection process
        CellTypeRecognitionFeatureSelection *cellTypeMeasures;
        if(!getProcess("Mesh/Cell Types/Classification/A Select Measures", cellTypeMeasures))
        throw(QString("MeshProcessHeatMap:: Unable to make CellTypeRecognitionFeatureSelection process"));
        std::set<QString>fileFeatures = cellTypeMeasures->selectedMeasures;//parseMeasureString(cellTypeMeasures->parm("Selected Measures"));

        // make sure the attribute maps are all there
        CellTypeRecognitionFeatureSelection fs(*this);
        fs.setSelectedMeasures(fileFeatures);
        if(!fs.calcSelectedMeasures(mesh, forceRecalc,false)) return setErrorMessage("Not all required Attr Maps present!");
        std::cout<<"Debug: "<<__LINE__<<"\n";
        int svm_type = -1, nr_class = -1, nr_features = -1;
        if(model){
            svm_type = svm_get_svm_type(model);     //Create test data      //Ask
            nr_class = svm_get_nr_class(model);
            nr_features = svm_get_nr_sv(model);
        }
        std::cout<<"Found SVM of type: "<<svm_type<< " with " << nr_class << " classes and " << nr_features << " support vectors" << std::endl;
        std::cout<<"Debug: "<<nr_class << "/" << __LINE__<<"\n";
        //double *prob_estimates = NULL;
        //prob_estimates = (double *) malloc(nr_class*sizeof(double));
        static const int nrOfClasses = nr_class;
        double prob_estimates[nrOfClasses];
        static const int nrOfFeatures = fileFeatures.size();
        static const int nrOfFeaturesSVM = svmc.modelFeatures.size();
        std::cout<<"Debug: "<< nrOfFeatures << "/" << nrOfFeaturesSVM << " \n";

        forall(int label, allLabels){
            if(label < 1) continue;
            int j = 0, cnt = 0;
            //std::cout<<"Debug: "<<__LINE__<<"\n";
            //svm_node *testnode = Malloc(svm_node, prob.l);
            //static const int testNodeSize = prob.l;// = nr_class;
            //std::cout<<"Debug: "<<testNodeSize << "/" << __LINE__<<"\n";
            svm_node *testnode = new svm_node[nrOfFeaturesSVM];
            //std::cout<<"Debug: "<<__LINE__<< fileFeatures.size() << "\n";
            forall(const QString &name, svmc.modelFeatures){
              forall(IntFloatP loop, cc.measures[name]){
                  if(loop.first == label){
                      //std::cout<<"\nDebug\t"<<__LINE__<<"\t" <<v->label <<"\t" << loop.second;
                      testnode[j].index = j;
                      testnode[j].value = loop.second;
                      ++j;
                  }
              }
              cnt++;
            }
            testnode[j].index = -1;   //Each row of properties should be terminated with a -1 according to the readme
            int predict_label = 0;
            std::cout << "classify " << label << "/" << j << std::endl;
            if(model) predict_label = svm_predict_probability(model,testnode,prob_estimates);    //Get the probabilty of predicted label
            std::cout << "done " << predict_label << std::endl;

            double maxEstimate  = -1;
            double sum = 0;
            std::cout << "class " << nr_class << std::endl;

            for(int p = 0; p < nr_class; p++){
              std::cout << "prob " << prob_estimates[p] << std::endl;
                sum+=prob_estimates[p];
                if(prob_estimates[p]>maxEstimate){     //If probabilty is above 0.8, then mark the label with parent label
                  maxEstimate = prob_estimates[p];
                  if(prob_estimates[p] > classifyThreshold) mesh->parents()[label] = predict_label;
                  mesh->labelHeat()[label] = maxEstimate;
                }
                //std::cout<<"probest "<<label << "/" << prob_estimates[p] << "/" << predict_label << std::endl;
            }
            //delete[] testnode;
            //std::cout<<"Debug: "<<sum<<"\n";
                   // if(testnode) free(testnode);
        }
        mesh->updateTriangles();    //Update the mesh with new labels
        //std::cout<<"\nDebug: "<<__LINE__<<"\n";
        //if(choice == "No"){

std::cout<<"Debug: "<<__LINE__<<"\n";
        //}
        return true;
    }
    REGISTER_PROCESS(CellTypeRecognitionClassification);

  bool KMeansTraining::run(Mesh *m, int steps, double threshold, bool useScaledMeasures)
  {

std::cout << "kmeans start" << std::endl;

    // get selected measures
        CellTypeRecognitionFeatureSelection *cellTypeMeasures;
        if(!getProcess("Mesh/Cell Types/Classification/A Select Measures", cellTypeMeasures))
        throw(QString("MeshProcessHeatMap:: Unable to make CellTypeRecognitionFeatureSelection process"));
        std::set<QString>fileFeatures = parseMeasureString(cellTypeMeasures->parm("Selected Measures"));

    CellTypeRecognitionFeatureSelection fs(*this);
    fs.setSelectedMeasures(fileFeatures);

    if(fileFeatures.size() == 0)
      return setErrorMessage("No Features Selected. Please run the Process Select Measures first.");

    // data structures for kMeans
    std::vector<AttrMap<int, double> > data;
    std::vector<std::vector<double> > inputClusters;

    std::map<int, int> parentLabel_ClusterIdxMap; // map for cell types -> which cluster idx
    std::vector<ClusterData> cluster;
    //cluster.clear();

    std::vector<int> labels = findAllLabels(m->graph());
std::cout << "kmeans init" << std::endl;
    forall(int label, labels){
      int parent = m->parents()[label];

      if(parent==0) continue; // only consider cells with assigned type
      int idx = parentLabel_ClusterIdxMap[parent];

      if(idx == 0){ // create new cluster
        std::cout << "new cluster " << std::endl;
        ClusterData newCluster;
        newCluster.associatedParent = parent;
        std::vector<double> posVec(fileFeatures.size());
        newCluster.position = posVec;
        newCluster.trainingCount = 0;
        cluster.push_back(newCluster);
        idx = cluster.size();
        parentLabel_ClusterIdxMap[parent] = idx;
      } //else { // cluster already exists
  std::cout << "idx " << idx << std::endl;
      int feature = 0;
      forall(const QString &name, fileFeatures){
        cluster[idx-1].position[feature] += cc.measures[name][label];
        std::cout << "f " << feature << "/" << cc.measures[name][label] << "/" << cluster[idx-1].position[feature] << std::endl;
        feature++;
      }
      cluster[idx-1].trainingCount++;
      //}

    }

    if(cluster.size() == 0)
      return setErrorMessage("No Parent labels found. Please specify some sample training data by assigning parent labels to some cells.");

    forall(ClusterData c, cluster){
      std::cout << "inp cluster " << std::endl;
      for(size_t i = 0; i<c.position.size(); i++){
      //forall(const double& d, c.position){
        c.position[i]/=c.trainingCount;
        std::cout << c.position[i] << "/" << c.trainingCount << std::endl;
      }
      inputClusters.push_back(c.position);
    }

    forall(const QString &name, fileFeatures){
      AttrMap<int, double> newMap;
      forall(IntFloatP p, cc.measures[name]){
        newMap[p.first] = p.second;
      }
      data.push_back(newMap);
    }

std::cout << "kmeans lloyd" << data.size() << "/" << inputClusters.size() << "/" << inputClusters[0].size() << "//" << inputClusters[1].size() << std::endl;
    if(!lloydKMeans(inputClusters, data, steps, threshold)) std::cout << "ERROR" << std::endl;

    kMeansClusters = inputClusters;
    clusterD = cluster;

    for(size_t c = 0; c<kMeansClusters.size(); c++){
      int feature = 0;
      forall(const QString &name, fileFeatures){
        std::cout << kMeansClusters[c][feature] << "/" << clusterD[c].associatedParent << std::endl;
        feature++;
      }
    }

    return true;
  }
  REGISTER_PROCESS(KMeansTraining);

  bool KMeansClassification::run(Mesh *m)
  {

    // get selected measures
    CellTypeRecognitionFeatureSelection *cellTypeMeasures;

     if(!getProcess("Mesh/Cell Types/Classification/A Select Measures", cellTypeMeasures))
      throw(QString("MeshProcessHeatMap:: Unable to make CellTypeRecognitionFeatureSelection process"));


    std::set<QString> fileFeatures = parseMeasureString(cellTypeMeasures->parm("Selected Measures"));

    cellTypeMeasures->setSelectedMeasures(fileFeatures);

    if(kMeansClusters.size() == 0) return setErrorMessage("No Trained Clusters present!");

    std::vector<int> labels = findAllLabels(m->graph());

    std::cout << "kmeans lloyd done, now find cell types" << std::endl;
    std::map<int, int> parentCellCount;

    forall(int label, labels){
      double minDis = 1E20;
      int minClusterIdx = 0;

      for(size_t c = 0; c<kMeansClusters.size(); c++){

        double dis = 0;
        int feature = 0;
        forall(const QString &name, fileFeatures){
          double d = cc.measures[name][label] - kMeansClusters[c][feature];
          dis+=d*d;
          feature++;
        }
        if(dis<minDis){
          minDis = dis;
          minClusterIdx = c;
        }
      }
      m->parents()[label] = clusterD[minClusterIdx].associatedParent;
      parentCellCount[m->parents()[label]]++;

    }
      std::cout << "parent label/number of cells" << std::endl;
    forall(IntInt p, parentCellCount){
      std::cout << p.first << "/" << p.second << std::endl;
    }

std::cout << "kmeans done" << std::endl;
    m->updateTriangles();
    return true;

  }
  REGISTER_PROCESS(KMeansClassification);

  bool AugmentTrainingFile::initialize(QWidget* parent)
  {
    QString filename = parm("Filename");
    if(parent)
      filename = QFileDialog::getOpenFileName(0, "Choose txt file to load", QDir::currentPath(), "txt files (*.txt)");
    if(filename.isEmpty())
      return false;
    if(!filename.endsWith(".txt", Qt::CaseInsensitive))
      filename += ".txt";
    setParm("Filename", filename);
    return true;
  }

  bool AugmentTrainingFile::run(Mesh *m, QString filename, int repeats, double noise, bool balance, bool shuffle, int maxEntries)
  {

    FileHandlerSVM fhSVM;

    FileData fd;

    fhSVM.readTrainingFile(fd, filename);

    if(repeats > 0)
      fhSVM.augmentTrainingFile(fd, repeats, noise);

    if(balance)
      fhSVM.balanceTrainingFile(fd);

    if(shuffle)
      fhSVM.shuffleTrainingFile(fd);

    if(maxEntries > 0)
      fhSVM.limitTrainingFile(fd, maxEntries);

    fhSVM.writeTrainingFile(fd, filename, false);

    return true;

  }
  REGISTER_PROCESS(AugmentTrainingFile);


//duplicated function, also in CellAtlas.cpp
// set a number to a label
void setLabelProperties(QLabel* label, int x, int y, int number, bool reset, QRgb col)
  {
    int xOffset = 10+10;
    int yOffset = 100+7;
    int xSize = 10;
    int ySize = 10;

    int r = qRed(col);
    int g = qGreen(col);
    int b = qBlue(col);

    label->setGeometry(QRect(x+xOffset, y+yOffset, xSize, ySize));
    if(!reset)
    label->setText(QString::number(number+1));
    else
    label->setText("");
    QString s = "QLabel { color : rgb(" + QString::number(r) + ", " + QString::number(g) + ", " + QString::number(b) + "); }";
    label->setStyleSheet(s);

  }



// for returning the right UI label
 QLabel* getLabel(Ui_HeatMapClusterDialog& ui, int label)
  {
    if(label == 0){
      return ui.labelCluster1;
    } else if(label == 1){
      return ui.labelCluster2;
    } else if(label == 2){
      return ui.labelCluster3;
    } else if(label == 3){
      return ui.labelCluster4;
    } else if(label == 4){
      return ui.labelCluster5;
    } else if(label == 5){
      return ui.labelCluster6;
    } else if(label == 6){
      return ui.labelCluster7;
    } else if(label == 7){
      return ui.labelCluster8;
    } else if(label == 8){
      return ui.labelCluster9;
    } else if(label == 9){
      return ui.labelCluster10;
    } else if(label == 10){
      return ui.labelCluster11;
    } else if(label == 11){
      return ui.labelCluster12;
    } else if(label == 12){
      return ui.labelCluster13;
    } else if(label == 13){
      return ui.labelCluster14;
    } else if(label == 14){
      return ui.labelCluster15;
    } else if(label == 15){
      return ui.labelCluster16;
    } else if(label == 16){
      return ui.labelCluster17;
    } else if(label == 17){
      return ui.labelCluster18;
    } else if(label == 18){
      return ui.labelCluster19;
    } else if(label == 19){
      return ui.labelCluster20;
    }
    return ui.labelCluster20;
  }
// find the label that has to be changed
// requires ui
void drawNumber(Ui_HeatMapClusterDialog& ui, int x, int y, int label, int number)
  {
  QRgb col = qRgba(255,255,255,255);
  setLabelProperties(getLabel(ui, label), x, y, number,false,col);

  }

// reset the label number
void resetNumber(Ui_HeatMapClusterDialog& ui, int label)
  {
  QRgb col = qRgba(255,255,255,255);
  setLabelProperties(getLabel(ui, label), 0, 0, 0,true,col);

  }

// for calculating one kMeans Step for 2/n-dimensional cluster
// be aware of the scaling of the data as this can change results
void ClusterCells::kMeansN()
  {
    std::cout << "kMeans " << clMap.mode2D << std::endl;

    std::vector<AttrMap<int, double> > dataPoints;


    if(clMap.mode2D){
      // just take the selected X/Y
      dataPoints.push_back(clMap.cellFeatures[clMap.attrIdxMap[selectedX]-1]);
      dataPoints.push_back(clMap.cellFeatures[clMap.attrIdxMap[selectedY]-1]);
    } else {
      // take all
      dataPoints = clMap.cellFeatures;
    }

    std::vector<std::vector<double> > inputClusters;

    for(size_t i =0; i<clMap.cluster.size(); i++){
      ManualCluster& c = clMap.cluster[i];

      std::vector<double> cl;

      if(clMap.mode2D){
        // just take the selected X/Y
        cl.push_back(c.val[clMap.attrIdxMap[selectedX]-1]);
        cl.push_back(c.val[clMap.attrIdxMap[selectedY]-1]);
        std::cout << "cb " << c.val[clMap.attrIdxMap[selectedX]-1] << "/" << c.val[clMap.attrIdxMap[selectedY]-1] << std::endl;

      } else {
        // take all
        for(size_t j = 0; j<clMap.cellFeatures.size(); j++){
          cl.push_back(c.val[j]);
        }
      }

      inputClusters.push_back(cl);

    }

    //std::vector<Point2d> outputClusters;
    if(!lloydKMeans(inputClusters, dataPoints, 1)) cout << "ERROR" << endl;
std::cout << "run" << std::endl;
    int counter = 0;
    //for(std::map<int, Point2d>::iterator it = maxVec.begin(); it!= maxVec.end(); it++){
    for(size_t i =0; i<clMap.cluster.size(); i++){
      ManualCluster& c = clMap.cluster[i];

      if(clMap.mode2D){
        // just take the selected X/Y
        c.val[clMap.attrIdxMap[selectedX]-1] = inputClusters[i][0];
        c.val[clMap.attrIdxMap[selectedY]-1] = inputClusters[i][1];
        std::cout << "c " << counter << "/" << inputClusters[i][0] << "/" << inputClusters[i][1] << std::endl;
      } else {
        // take all
        for(size_t j = 0; j<clMap.cellFeatures.size(); j++){
          c.val[j] = inputClusters[i][j];
        }
      }
      counter++;
    }
    std::cout << "done" << std::endl;
    // update GUI
    setImage();

  }

// average the positions of cells with the same parent label and generate clusters from them
void ClusterCells::clustersFromParents()
{
  // go through all cells, get values
  std::map<int, int> parentCountMap;

  std::set<int> parentLabels;

  std::map<int, double> parentX, parentY;

  Mesh* m1 = currentMesh();

  forall(const vertex& v, m1->graph()){
    if(v->label < 1) continue;
  //std::cout << "cell" << v->label << std::endl;
    int pLabel = m1->parents()[v->label];
    if(parentLabels.find(pLabel) == parentLabels.end()){
      // new parent label
      parentLabels.insert(pLabel);
      parentCountMap[pLabel] = 1;
      //parentX[pLabel] = clMap.getValue(selectedX, v->label);//cellFeatures[clMap.attrIdxMap[selectedX]-1][v->label];
      //parentY[pLabel] = clMap.getValue(selectedY, v->label);//.cellFeatures[clMap.attrIdxMap[selectedY]-1][v->label];
    } else {
      // existing
      parentCountMap[pLabel] += 1;
    }
      parentX[pLabel] += clMap.getValue(selectedX, v->label);//.cellFeatures[clMap.attrIdxMap[selectedX]-1][v->label];
      parentY[pLabel] += clMap.getValue(selectedY, v->label);//.cellFeatures[clMap.attrIdxMap[selectedY]-1][v->label];
  }
    std::cout << "done" << std::endl;
    clMap.resetClusters();
    int counter = 0;
    forall(int p, parentLabels){

      std::cout << "pp " << parentX[p] << "/" << parentY[p] << "/" << parentCountMap[p] << std::endl;

      double xVal = parentX[p]/(double)parentCountMap[p];
      double yVal = parentY[p]/(double)parentCountMap[p];

      clMap.addCluster(selectedX, selectedY, realToImageCoord(clMap.xMinMax, xVal, clMap.gridSize), realToImageCoord(clMap.yMinMax, yVal, clMap.gridSize));
      clMap.cluster[counter].label = p % 16 > 9? 8 : p % 16 - 1;
      counter++;
    }

  changeHeatmap();
}

// generate the heatmap image
void ClusterCells::setImage()
  {
  if(ui.showParentsCheckBox->isChecked()) clMap.showParentLabels = true;
  else clMap.showParentLabels = false;

  if(ui.largerDotsCheckBox->isChecked()) clMap.largeDots = true;
  else clMap.largeDots = false;

  //if(ui.radioButton2D->isChecked()){
    clMap.mode2D = true;
  //} else {
  //  clMap.mode2D = false;
  //}

  QImage image;

  QRgb white = qRgba(255,255,255,255);

  selectedX = ui.xDimCmb->currentText();
  selectedY = ui.yDimCmb->currentText();
  std::cout << "selected " << selectedX.toLocal8Bit().constData() << " and " << selectedY.toLocal8Bit().constData() << std::endl;
  redraw = true;

  int background = 0;
  if(ui.background->isChecked()) background++;

  //std::cout << "setImage" << std::endl;
  clMap.setActiveMeasures(selectedX, selectedY, selectedHeat);
  if(ui.imageHeatButton->isChecked() and redraw){
    //std::cout << "heat" << std::endl;
    clMap.createHeatMap2D();
    //clMap.setHeatMax();
    image = clMap.createImage(true, ui.showParents0->isChecked(), false, background);
    std::cout << "img created" << std::endl;

  } else if (ui.imageCellsButton->isChecked() and redraw){

    clMap.createCellView();
    image = clMap.createImage(false, ui.showParents0->isChecked(), false, background);
  }
  redraw = false;

    // update min and max of the image window in the GUI
    ui.spinBoxSigmaXMin->setValue(clMap.xMinMax.x());
    ui.spinBoxSigmaXMax->setValue(clMap.xMinMax.y());
    ui.spinBoxSigmaYMin->setValue(clMap.yMinMax.x());
    ui.spinBoxSigmaYMax->setValue(clMap.yMinMax.y());
    ui.spinBoxScaleHeat->setValue(clMap.heatMax);

    int size = clMap.gridSize;
    int factor = clMap.gridFactor;
    int yOffset = size*factor-1;
    // reset the numbers
    for(int i = 0; i < 20; i++) {
      resetNumber(ui, i);
    }

    // draw the user defined cluster center with circle and big cross
    for(size_t i = 0; i < clMap.cluster.size(); i++){
      ManualCluster c = clMap.cluster[i];
      double valX = c.val[clMap.attrIdxMap[selectedX]-1];
      double valY = c.val[clMap.attrIdxMap[selectedY]-1];

      double xImage = interpolateArrayIndex(valX,clMap.xMinMax.x(),clMap.xMinMax.y(),clMap.gridSize-1);
      double yImage = interpolateArrayIndex(valY,clMap.yMinMax.x(),clMap.yMinMax.y(),clMap.gridSize-1);

      //std::cout << "drawCluster " << xImage << "/" << yImage << std::endl;

      Point2d currentPoint(xImage,yImage);//c.val;//maxVec2[i];

      drawCross(image, currentPoint[0]*factor+factor/2, yOffset-currentPoint[1]*factor-factor/2, size*factor, 4, white);
      drawCircle(image, currentPoint[0]*factor+factor/2, yOffset-currentPoint[1]*factor-factor/2, size*factor, white);
      if(i<20) // there are only 20 labels in the GUI
        drawNumber(ui, currentPoint[0]*factor+factor/2-2, yOffset-currentPoint[1]*factor-factor/2+2, i, c.label);
    }

    ui.imageHolder->setPixmap(QPixmap::fromImage(image));
    ui.imageHolder->setFocus();
    std::cout << "done setImage " << std::endl;
  }

// return the attr map for a specified measure
void ClusterCells::getAttrMap(AttrMap<int, double>& data, QString measure)
{
  Mesh* m1 = currentMesh();
  QString prefix = "Measure Label Double ";
  QString proName = "Mesh/Heat Map/Measures/";
  data = m1->attributes().attrMap<int, double>(prefix+measure);

  if(data.size() == 0){ // calculate it
    qDebug() << "Attr Map \t" << prefix+measure << " doesnt exist\n";
    QStringList parms;
    Process *pro = makeProcess(proName+measure);
    if(pro){
    qDebug() << "Process call\t" << proName+measure << "\n";
    //getLastParms(proName+measure);
    pro->run();
    data = m1->attributes().attrMap<int, double>(prefix+measure);
    }
  }
  clMap.addCellFeature(measure, data, selectedCells);
}

  // called when the measures in the GUI change
 void ClusterCells::changeHeatmap()
  {
  if(selectedX == ui.xDimCmb->currentText() and selectedY == ui.yDimCmb->currentText() and selectedHeat == ui.heatDimCmb->currentText()){
    redraw = false;
  } else {
    redraw = true;
    clMap.customMinMax = false;
  }
  selectedX = ui.xDimCmb->currentText();
  selectedY = ui.yDimCmb->currentText();
  selectedHeat = ui.heatDimCmb->currentText();

  AttrMap<int, double> dataX, dataY, dataHeat;
  getAttrMap(dataX,selectedX);
  getAttrMap(dataY,selectedY);
  if(selectedHeat != "none") getAttrMap(dataY,selectedHeat);

  setImage();

  }

void ClusterCells::changeSigma(double sigma)
  {
    clMap.sigma = sigma;
    redraw = true;
    changeHeatmap();
  }

// upon mouse click: if a maximum is nearby: activate the maximum (for drag and drop)
// if no maximum is nearby create a new one
void ClusterCells::setPosition(const QPoint& p)
  {
    Point2d currentPoint;
    currentPoint[0] = p.x()/clMap.gridFactor;
    currentPoint[1] = (clMap.gridSize*clMap.gridFactor-1-p.y())/clMap.gridFactor;

    //int maxVecSize = rcp.mainBody.maximaHeatMapSelected.size();

    int idx;
    // check whether there already is a maximum at the clicked position
    if(clMap.nearbyManualCluster(selectedX, selectedY, currentPoint, idx)){
      // yes : activate this maximum
      clMap.activeCluster = idx;

    } else {
      // no : create a new maximum
      clMap.addCluster(selectedX, selectedY, currentPoint.x(), currentPoint.y());
      clMap.activeCluster = -1;
    }

    changeHeatmap();
    std::cout << "done setPosition " << clMap.activeCluster << std::endl;
  }

// upon keypress: check whether there is a cluster nearby and if yes: assign the label to it
void ClusterCells::setClusterLabel(QString label){

    if(label == "1" or label == "2" or label == "3" or label == "4" or label == "5" or label == "6" or label == "7" or label == "8" or label == "9"){
      int labelNr = label.toInt() - 1;
      int idx;

      if(clMap.nearbyManualCluster(selectedX, selectedY, mousePos, idx)){
        // assign label to cluster
        clMap.cluster[idx].label = labelNr;
      }
    } else {
      cout << "Enter a number between 1 and 9" << endl;
    }
    changeHeatmap();

  }

// upon mouse release: check whether there is an activated maximum (=nearby maximum) and move it to the new position
void ClusterCells::setReleasePosition(const QPoint& p)
  {
    std::cout << "setReleasePosition " << clMap.gridFactor << "/" << clMap.activeCluster << std::endl;
    Point2d currentPoint;
    currentPoint[0] = p.x()/clMap.gridFactor;
    currentPoint[1] = (clMap.gridSize*clMap.gridFactor-1-p.y())/clMap.gridFactor;

      if(clMap.activeCluster != -1){    // check if there is an activated maximum
        //yes : move it to the new position
        clMap.updateCluster(clMap.activeCluster, selectedX, currentPoint.x());
        clMap.updateCluster(clMap.activeCluster, selectedY, currentPoint.y(), false);
      }

    changeHeatmap();

    clMap.activeCluster = -1;
    std::cout << "done setReleasePosition" << std::endl;
  }

// update the position of the mouse within the heatmap window
void ClusterCells::setMousePosition(const QPoint& p){

    Point2d currentPoint;
    currentPoint[0] = p.x()/clMap.gridFactor;
    currentPoint[1] = (clMap.gridSize*clMap.gridFactor-1-p.y())/clMap.gridFactor;
    mousePos = currentPoint;
  }

  // fill the combo boxes with all available measures and/or attr maps
  void ClusterCells::fillMeasures(){

    if(ui.comboMeshType->currentText() == "2D"){
      subfolder = "Measures";
    } else {
      subfolder = "Measures 3D";
    }

    ui.xDimCmb->clear();
    ui.yDimCmb->clear();
    ui.heatDimCmb->clear();
    QStringList defaultEntries;
    defaultEntries << "none";
    ui.heatDimCmb->addItems(defaultEntries);

    int idx = 0;
    int selIdxX = 0;
    int selIdxY = 0;

    // find measure processes
    QStringList procs = mgx::listProcesses();

    QStringList addProcs;
    std::set<QString> procsSet;
    QStringList addAttr;
    forall(const QString& name, procs) {
      mgx::ProcessDefinition* def = mgx::getProcessDefinition(name);
      QStringList list = def->name.split("/");

      if(list[0] == "Mesh" and list[1] == "Heat Map" and list[2] == subfolder){
        QString newProc = list[3] + "/" + list[4];
        addProcs << newProc;
        procsSet.insert(newProc);
        if(newProc == selectedX) selIdxX = idx;
        if(newProc == selectedY) selIdxY = idx;
        idx++;
      }

    }

    // now the attr maps
    Mesh* m = currentMesh();
    Attributes *attributes;
    attributes = &m->attributes();
    QStringList attr = attributes->getAttrList();
    QStringList measureAttr;

    forall(const QString &text, attr) {
      QStringList list = text.split(" ");

      if(list.size() < 4) continue;
      if(list[0] != "Measure") continue;
      QString name;
      for(int i = 3; i<list.size(); i++){
        name = name + list[i] + " ";
      }
      name = name.trimmed();
      if(procsSet.find(name) == procsSet.end()){
        addProcs << name;
          //if(name == selectedX) selIdxX = idx;
          //if(name == selectedY) selIdxY = idx;
          //idx++;
      }
    }

    addProcs.sort(Qt::CaseSensitive);

    forall(const QString &name, attr){
      if(name == selectedX) selIdxX = idx;
      if(name == selectedY) selIdxY = idx;
      idx++;
    }

    ui.xDimCmb->addItems(addProcs);
    ui.yDimCmb->addItems(addProcs);
    ui.heatDimCmb->addItems(addProcs);

    ui.xDimCmb->setCurrentIndex(selIdxX);
    ui.yDimCmb->setCurrentIndex(selIdxY);
    ui.heatDimCmb->setCurrentIndex(0);

  }

// parse the parms string
bool ClusterCells::parseParmsString(QString parmsString)
{

  //std::map<int,Point2d> selectedClustersNew;

  std::set<QString> selectedMeasuresNew;
  QStringList measures = parmsString.split("$");
  if(measures.size()<3) return false;

  selectedX = measures[0];
  selectedY = measures[1];
  sigma = measures[2].toDouble();
  gaussianMode = stringToBool(measures[3]);

  for(int i = 4; i < measures.size(); i=i+3){
    if(measures.size()<i+2) i = measures.size();
    else {
      //int label = measures[i].toInt();
      double posX = measures[i+1].toDouble();
      double posY = measures[i+2].toDouble();

      clMap.addCluster(selectedX, selectedY, posX, posY);

      //selectedClustersNew[label] = Point2d(posX, posY);
    }

  }
  //selectedClusters = selectedClustersNew;
  return true;
}

// create the parms string
QString ClusterCells::createParmsString()
{
  QString parmsNew;

  parmsNew += selectedX;
  parmsNew += "$" + selectedY;
  parmsNew += "$" + QString::number(sigma);
  if(gaussianMode)
    parmsNew += "$Yes";
  else
    parmsNew += "$No";

  //typedef std::pair<int, Point2d> IntP2d;
  forall(ManualCluster c, clMap.cluster){
    parmsNew += "$" + QString::number(c.label);
    parmsNew += "$" + QString::number(c.val[clMap.attrIdxMap[selectedX]]);
    parmsNew += "$" + QString::number(c.val[clMap.attrIdxMap[selectedY]]);
  }

  return parmsNew;
}

  // reset and delete all cluster
  void ClusterCells::resetClusters()
  {
  clMap.resetClusters();
  changeHeatmap();

  }

void ClusterCells::updateScreen()
{
std::cout << "update screen" << std::endl;
  Mesh* m1 = currentMesh();
  clMap.setActiveMeasures(selectedX, selectedY, selectedHeat);
  if(ui.imageHeatButton->isChecked()){
    clMap.setParentsHeatMapView();
  } else if (ui.imageCellsButton->isChecked()){
    clMap.setParentsCellView();
  }

  //forall(const vertex& v, m1->graph()){
  forall(const int& l, selectedCells){
    //if(v->label == -1) continue;
    m1->parents()[l] = clMap.parentLabels[l];
  }
      changeHeatmap();
  m1->setShowLabel("Label");
  m1->setUseParents(true);
  m1->updateTriangles();
  updateState();
  updateViewer();
}


void ClusterCells::changeMinMax(){

  clMap.customMinMax = true;

  clMap.xMinMax.x() = ui.spinBoxSigmaXMin->value();
  clMap.xMinMax.y() = ui.spinBoxSigmaXMax->value();
  clMap.yMinMax.x() = ui.spinBoxSigmaYMin->value();
  clMap.yMinMax.y() = ui.spinBoxSigmaYMax->value();

  clMap.heatMax = ui.spinBoxScaleHeat->value();

  changeHeatmap();
}

// restore the old parent labels (note that they can be overwritten if the GUI is closed using "OK")
void ClusterCells::restoreParents()
{
  Mesh* m1 = currentMesh();
  forall(const vertex& v, m1->graph()){
    if(v->label == -1) continue;
    clMap.parentLabels[v->label] = backupParents[v->label];
    m1->parents()[v->label] = backupParents[v->label];
  }
  changeHeatmap();
  m1->setShowLabel("Label");
  m1->setUseParents(true);
  m1->updateTriangles();
  updateState();
  updateViewer();
}


  bool ClusterCells::initialize(QWidget* parent)
  {
    ClusterMap clMapNew;
    clMap = clMapNew;
    clMap.sigma = sigma = 4;
    gaussianMode = false;
    clMap.customMinMax = false;
    clMap.gridSize = 480;
    clMap.gridFactor = 1;
    //bool parmsOk = parseParmsString(parms[0]);

    if(!parent){
      return true;
    }
    //cout << "sel11 " << endl;
    //if(parmsOk){
      //rcp.mainBody.maximaHeatMapSelected = selectedClusters;

    //}

    // find selected cell labels & backup parents
    std::set<int> selCells, allCells;
    Mesh* m1 = currentMesh();
    forall(const vertex& v, m1->graph()){
      if(v->label == -1) continue;
      clMap.parentLabels[v->label] = m1->parents()[v->label];
      backupParents[v->label] = m1->parents()[v->label];
      allCells.insert(v->label);
      if(v->selected) selCells.insert(v->label);
    }

    if(selCells.empty()) selCells = allCells;
    selectedCells = selCells;

    dlg = new QDialog(parent);
    ui.setupUi(dlg);
    this->dlg = dlg;
    //if(ui.radioButton2D->isChecked()) clMap.mode2D = true;
    //else clMap.mode2D = false;

    fillMeasures();

    //fillMeasures();

    if(gaussianMode) ui.imageCellsButton->setChecked(true);
    else ui.imageHeatButton->setChecked(true);

    //ui.preselectButton->setEnabled(false);
    ui.spinBoxSigma->setValue(sigma);

    changeHeatmap();

    // check which is selected
    connect(ui.yDimCmb, SIGNAL(currentIndexChanged(QString)), this, SLOT(changeHeatmap()));
    connect(ui.xDimCmb, SIGNAL(currentIndexChanged(QString)), this, SLOT(changeHeatmap()));
    connect(ui.heatDimCmb, SIGNAL(currentIndexChanged(QString)), this, SLOT(changeHeatmap()));

    connect(ui.kMeansButton, SIGNAL(released()), this, SLOT(kMeansN()));
    connect(ui.updateButton, SIGNAL(released()), this, SLOT(updateScreen()));
    connect(ui.buttonReset, SIGNAL(released()), this, SLOT(resetClusters()));
    connect(ui.clustersFromParentsButton, SIGNAL(released()), this, SLOT(clustersFromParents()));

    // calculate measure & update heatmap

    connect(ui.spinBoxSigma, SIGNAL(valueChanged(double)), this, SLOT(changeSigma(double)));

    connect(ui.showParentsCheckBox, SIGNAL(stateChanged(int)), this, SLOT(setImage()));
    connect(ui.largerDotsCheckBox, SIGNAL(stateChanged(int)), this, SLOT(setImage()));
    connect(ui.showParents0, SIGNAL(stateChanged(int)), this, SLOT(setImage()));
    connect(ui.showParents1, SIGNAL(stateChanged(int)), this, SLOT(setImage()));
    connect(ui.background, SIGNAL(stateChanged(int)), this, SLOT(setImage()));

    connect(ui.imageCellsButton,SIGNAL(toggled(bool)),this,SLOT(setImage()));
    connect(ui.imageHeatButton,SIGNAL(toggled(bool)),this,SLOT(setImage()));

    connect(ui.restoreParentsButton,SIGNAL(released()),this,SLOT(restoreParents()));
    connect(ui.rescaleButton,SIGNAL(released()),this,SLOT(changeMinMax()));

    connect(ui.comboMeshType, SIGNAL(currentIndexChanged(QString)), this, SLOT(fillMeasures()));

    // clustering
    connect(ui.imageHolder, SIGNAL(mousePressed( const QPoint& )), this, SLOT(setPosition(const QPoint&)));
    connect(ui.imageHolder, SIGNAL(mouseReleased( const QPoint& )), this, SLOT(setReleasePosition(const QPoint&)));
    connect(ui.imageHolder, SIGNAL(keyPressed( const QString& )), this, SLOT(setClusterLabel(QString)));
    connect(ui.imageHolder, SIGNAL(mouseMove( const QPoint& )), this, SLOT(setMousePosition(const QPoint&)));

    if(dlg->exec() == QDialog::Accepted){
      ///selectedClusters = rcp.mainBody.maximaHeatMapSelected;
      sigma = clMap.sigma;

      //parms[0] = createParmsString();

      return true;
    } else {
      return false;
    }

  }

  // CellAtlas GUI process for HeatMap clustering
  bool ClusterCells::run(Mesh *m)
  {
    updateScreen();

    return true;

  }
  REGISTER_PROCESS(ClusterCells);

  }
