source: src/Actions/AnalysisAction/PairCorrelationAction.cpp@ 192f6e

/*
 * PairCorrelationAction.cpp
 *
 *  Created on: May 9, 2010
 *      Author: heber
 */

#include "Helpers/MemDebug.hpp"

#include "Actions/AnalysisAction/PairCorrelationAction.hpp"
#include "analysis_correlation.hpp"
#include "boundary.hpp"
#include "linkedcell.hpp"
#include "log.hpp"
#include "element.hpp"
#include "molecule.hpp"
#include "periodentafel.hpp"
#include "vector.hpp"
#include "World.hpp"

#include <iostream>
#include <string>

using namespace std;

#include "UIElements/UIFactory.hpp"
#include "UIElements/Dialog.hpp"
#include "Actions/MapOfActions.hpp"

const char AnalysisPairCorrelationAction::NAME[] = "pair-correlation";

AnalysisPairCorrelationAction::AnalysisPairCorrelationAction() :
  Action(NAME)
{}

AnalysisPairCorrelationAction::~AnalysisPairCorrelationAction()
{}

Action::state_ptr AnalysisPairCorrelationAction::performCall() {
  Dialog *dialog = UIFactory::getInstance().makeDialog();
  int ranges[3] = {1, 1, 1};
  double BinEnd = 0.;
  double BinStart = 0.;
  double BinWidth = 0.;
  molecule *Boundary = NULL;
  string outputname;
  string binoutputname;
  bool periodic;
  ofstream output;
  ofstream binoutput;
  std::vector< element *> elements;
  string type;
  Vector Point;
  BinPairMap *binmap = NULL;
  MoleculeListClass *molecules = World::getInstance().getMolecules();

  // first dialog: Obtain which type of correlation
  dialog->queryString(NAME, &type, MapOfActions::getInstance().getDescription(NAME));
  if(dialog->display()) {
    delete dialog;
  } else {
    delete dialog;
    return Action::failure;
  }

  // second dialog: Obtain parameters specific to this type
  dialog = UIFactory::getInstance().makeDialog();
  if (type == "P")
    dialog->queryVector("position", &Point, World::getInstance().getDomain(), false, MapOfActions::getInstance().getDescription("position"));
  if (type == "S")
    dialog->queryMolecule("molecule-by-id", &Boundary, MapOfActions::getInstance().getDescription("molecule-by-id"));
  dialog->queryElement("elements", &elements, MapOfActions::getInstance().getDescription("elements"));
  dialog->queryDouble("bin-start", &BinStart, MapOfActions::getInstance().getDescription("bin-start"));
  dialog->queryDouble("bin-width", &BinWidth, MapOfActions::getInstance().getDescription("bin-width"));
  dialog->queryDouble("bin-end", &BinEnd, MapOfActions::getInstance().getDescription("bin-end"));
  dialog->queryString("output-file", &outputname, MapOfActions::getInstance().getDescription("output-file"));
  dialog->queryString("bin-output-file", &binoutputname, MapOfActions::getInstance().getDescription("bin-output-file"));
  dialog->queryBoolean("periodic", &periodic, MapOfActions::getInstance().getDescription("periodic"));

  if(dialog->display()) {
    output.open(outputname.c_str());
    binoutput.open(binoutputname.c_str());
    if (type == "E") {
      PairCorrelationMap *correlationmap = NULL;
      if (periodic)
        correlationmap = PeriodicPairCorrelation(World::getInstance().getMolecules(), elements, ranges);
      else
        correlationmap = PairCorrelation(World::getInstance().getMolecules(), elements);
      OutputPairCorrelation(&output, correlationmap);
      binmap = BinData( correlationmap, BinWidth, BinStart, BinEnd );
      OutputCorrelation ( &binoutput, binmap );
      delete(binmap);
      delete(correlationmap);
    } else if (type == "P")  {
      cout << "Point to correlate to is  " << Point << endl;
      CorrelationToPointMap *correlationmap = NULL;
      if (periodic)
        correlationmap  = PeriodicCorrelationToPoint(molecules, elements, &Point, ranges);
      else
        correlationmap = CorrelationToPoint(molecules, elements, &Point);
      OutputCorrelationToPoint(&output, correlationmap);
      binmap = BinData( correlationmap, BinWidth, BinStart, BinEnd );
      OutputCorrelation ( &binoutput, binmap );
      delete(binmap);
      delete(correlationmap);
    } else if (type == "S") {
      ASSERT(Boundary != NULL, "No molecule specified for SurfaceCorrelation.");
      const double radius = 4.;
      double LCWidth = 20.;
      if (BinEnd > 0) {
        if (BinEnd > 2.*radius)
          LCWidth = BinEnd;
        else
          LCWidth = 2.*radius;
      }

      // get the boundary
      class Tesselation *TesselStruct = NULL;
      const LinkedCell *LCList = NULL;
      // find biggest molecule
      int counter  = molecules->ListOfMolecules.size();
      bool *Actives = new bool[counter];
      counter = 0;
      for (MoleculeList::iterator BigFinder = molecules->ListOfMolecules.begin(); BigFinder != molecules->ListOfMolecules.end(); BigFinder++) {
        Actives[counter++] = (*BigFinder)->ActiveFlag;
        (*BigFinder)->ActiveFlag = (*BigFinder == Boundary) ? false : true;
      }
      LCList = new LinkedCell(Boundary, LCWidth);
      FindNonConvexBorder(Boundary, TesselStruct, LCList, radius, NULL);
      CorrelationToSurfaceMap *surfacemap = NULL;
      if (periodic)
        surfacemap = PeriodicCorrelationToSurface( molecules, elements, TesselStruct, LCList, ranges);
      else
        surfacemap = CorrelationToSurface( molecules, elements, TesselStruct, LCList);
      delete LCList;
      OutputCorrelationToSurface(&output, surfacemap);
      // re-set ActiveFlag
      counter = 0;
      for (MoleculeList::iterator BigFinder = molecules->ListOfMolecules.begin(); BigFinder != molecules->ListOfMolecules.end(); BigFinder++) {
        (*BigFinder)->ActiveFlag = Actives[counter++];
      }
      delete[] Actives;
      // check whether radius was appropriate
      {
        double start; double end;
        GetMinMax( surfacemap, start, end);
        if (LCWidth < end)
          DoeLog(1) && (eLog()<< Verbose(1) << "Linked Cell width is smaller than the found range of values! Bins can only be correct up to: " << radius << "." << endl);
      }
      binmap = BinData( surfacemap, BinWidth, BinStart, BinEnd );
      OutputCorrelation ( &binoutput, binmap );
      delete TesselStruct;  // surfacemap contains refs to triangles! delete here, not earlier!
      delete(binmap);
      delete(surfacemap);
    } else {
      return Action::failure;
    }
    output.close();
    binoutput.close();
    delete dialog;
    return Action::success;
  } else {
    delete dialog;
    return Action::failure;
  }
}

Action::state_ptr AnalysisPairCorrelationAction::performUndo(Action::state_ptr _state) {
//  ParserLoadXyzState *state = assert_cast<ParserLoadXyzState*>(_state.get());

  return Action::failure;
//  string newName = state->mol->getName();
//  state->mol->setName(state->lastName);
//
//  return Action::state_ptr(new ParserLoadXyzState(state->mol,newName));
}

Action::state_ptr AnalysisPairCorrelationAction::performRedo(Action::state_ptr _state){
  return Action::failure;
}

bool AnalysisPairCorrelationAction::canUndo() {
  return false;
}

bool AnalysisPairCorrelationAction::shouldUndo() {
  return false;
}

const string AnalysisPairCorrelationAction::getName() {
  return NAME;
}