source: src/Actions/FragmentationAction/FragmentationAutomationAction.cpp@ 0cd8cf

/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  2010-2012 University of Bonn. All rights reserved.
 * 
 *
 *   This file is part of MoleCuilder.
 *
 *    MoleCuilder is free software: you can redistribute it and/or modify
 *    it under the terms of the GNU General Public License as published by
 *    the Free Software Foundation, either version 2 of the License, or
 *    (at your option) any later version.
 *
 *    MoleCuilder is distributed in the hope that it will be useful,
 *    but WITHOUT ANY WARRANTY; without even the implied warranty of
 *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *    GNU General Public License for more details.
 *
 *    You should have received a copy of the GNU General Public License
 *    along with MoleCuilder.  If not, see <http://www.gnu.org/licenses/>.
 */

/*
 * FragmentationAutomationAction.cpp
 *
 *  Created on: May 18, 2012
 *      Author: heber
 */

// include config.h
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <boost/archive/text_iarchive.hpp>
// boost asio needs specific operator new
#include <boost/asio.hpp>

#include "CodePatterns/MemDebug.hpp"

// include headers that implement a archive in simple text format
#include <boost/archive/text_oarchive.hpp>
#include <boost/archive/text_iarchive.hpp>


#include <boost/mpl/remove.hpp>
#include <boost/lambda/lambda.hpp>

#include <iostream>

#include "CodePatterns/Assert.hpp"
#include "CodePatterns/Info.hpp"
#include "CodePatterns/Log.hpp"
#include "JobMarket/Jobs/FragmentJob.hpp"

#include "Fragmentation/Automation/createMatrixNrLookup.hpp"
#include "Fragmentation/Automation/extractJobIds.hpp"
#include "Fragmentation/Automation/FragmentationChargeDensity.hpp"
#include "Fragmentation/Automation/FragmentJobQueue.hpp"
#include "Fragmentation/Automation/FragmentationLongRangeResults.hpp"
#include "Fragmentation/Automation/FragmentationShortRangeResults.hpp"
#include "Fragmentation/Automation/MPQCFragmentController.hpp"
#include "Fragmentation/Automation/VMGDebugGridFragmentController.hpp"
#include "Fragmentation/Automation/VMGFragmentController.hpp"
#include "Fragmentation/EnergyMatrix.hpp"
#include "Fragmentation/ForceMatrix.hpp"
#include "Fragmentation/Fragmentation.hpp"
#include "Fragmentation/Homology/HomologyContainer.hpp"
#include "Fragmentation/Homology/HomologyGraph.hpp"
#include "Fragmentation/KeySet.hpp"
#include "Fragmentation/KeySetsContainer.hpp"
#include "Fragmentation/SetValues/Fragment.hpp"
#include "Fragmentation/SetValues/Histogram.hpp"
#include "Fragmentation/SetValues/IndexedVectors.hpp"
#include "Fragmentation/Summation/IndexSetContainer.hpp"
#include "Fragmentation/Summation/writeTable.hpp"
#include "Graph/DepthFirstSearchAnalysis.hpp"
#include "Jobs/MPQCJob.hpp"
#include "Jobs/MPQCData.hpp"
#include "Jobs/MPQCData_printKeyNames.hpp"
#ifdef HAVE_VMG
#include "Jobs/VMGDebugGridJob.hpp"
#include "Jobs/VMGJob.hpp"
#include "Jobs/VMGData.hpp"
#include "Jobs/VMGDataFused.hpp"
#include "Jobs/VMGDataMap.hpp"
#include "Jobs/VMGData_printKeyNames.hpp"
#endif
#include "World.hpp"

#include <fstream>
#include <iostream>
#include <string>
#include <vector>

#include <boost/mpl/for_each.hpp>

#include "Actions/FragmentationAction/FragmentationAutomationAction.hpp"

using namespace MoleCuilder;

// and construct the stuff
#include "FragmentationAutomationAction.def"
#include "Action_impl_pre.hpp"
/** =========== define the function ====================== */

class controller_AddOn;

// needs to be defined for using the FragmentController
controller_AddOn *getAddOn()
{
  return NULL;
}

void writeToFile(const std::string &filename, const std::string contents)
{
  std::ofstream tablefile(filename.c_str());
  tablefile << contents;
  tablefile.close();
}

/** Print (short range) energy, forces, and timings from received results.
 *
 * @param results summed up results container
 */
void printReceivedShortResults(
    const FragmentationShortRangeResults &results)
{
  // print tables (without eigenvalues, they go extra)
  {
    typedef boost::mpl::remove<MPQCDataEnergyVector_t, MPQCDataFused::energy_eigenvalues>::type
      MPQCDataEnergyVector_noeigenvalues_t;
    const std::string energyresult =
        writeTable<MPQCDataEnergyMap_t, MPQCDataEnergyVector_noeigenvalues_t >()(
            results.Result_Energy_fused, results.getMaxLevel());
    LOG(0, "Energy table is \n" << energyresult);
    std::string filename;
    filename += FRAGMENTPREFIX + std::string("_Energy.dat");
    writeToFile(filename, energyresult);
  }

  {
    const std::string eigenvalueresult;
    LOG(0, "Eigenvalue table is \n" << eigenvalueresult);
    std::string filename;
    filename += FRAGMENTPREFIX + std::string("_Eigenvalues.dat");
    writeToFile(filename, eigenvalueresult);
  }

  {
    const std::string forceresult =
        writeTable<MPQCDataForceMap_t, MPQCDataForceVector_t>()(
            results.Result_Force_fused, results.getMaxLevel());
    LOG(0, "Force table is \n" << forceresult);
    std::string filename;
    filename += FRAGMENTPREFIX + std::string("_Forces.dat");
    writeToFile(filename, forceresult);
  }
  // we don't want to print grid to a table
  {
    // print times (without flops for now)
    typedef boost::mpl::remove<
        boost::mpl::remove<MPQCDataTimeVector_t, MPQCDataFused::times_total_flops>::type,
        MPQCDataFused::times_gather_flops>::type
        MPQCDataTimeVector_noflops_t;
    const std::string timesresult =
        writeTable<MPQCDataTimeMap_t, MPQCDataTimeVector_noflops_t >()(
            results.Result_Time_fused, results.getMaxLevel());
    LOG(0, "Times table is \n" << timesresult);
    std::string filename;
    filename += FRAGMENTPREFIX + std::string("_Times.dat");
    writeToFile(filename, timesresult);
  }
}


/** Print long range energy from received results.
 *
 * @param results summed up results container
 */
void printReceivedFullResults(
    const FragmentationLongRangeResults &results)
{
  // print tables (without eigenvalues, they go extra)

  if (results.Result_LongRange_fused.size() >= (results.getMaxLevel()-2))
  {
    const std::string gridresult =
        writeTable<VMGDataMap_t, VMGDataVector_t >()(
            results.Result_LongRange_fused, results.getMaxLevel(), 2);
    LOG(0, "VMG table is \n" << gridresult);
    std::string filename;
    filename += FRAGMENTPREFIX + std::string("_VMGEnergy.dat");
    writeToFile(filename, gridresult);
  }

  if (results.Result_LongRangeIntegrated_fused.size() >= (results.getMaxLevel()-2))
  {
    const std::string gridresult =
        writeTable<VMGDataLongRangeMap_t, VMGDataLongRangeVector_t >()(
            results.Result_LongRangeIntegrated_fused, results.getMaxLevel(), 2);
    LOG(0, "LongRange table is \n" << gridresult);
    std::string filename;
    filename += FRAGMENTPREFIX + std::string("_LongRangeEnergy.dat");
    writeToFile(filename, gridresult);
  }
}

bool appendToHomologyFile(
    const boost::filesystem::path &homology_file,
    const FragmentationShortRangeResults &shortrangeresults,
    const FragmentationLongRangeResults &longrangeresults
    )
{
  /// read homology container (if present)
  HomologyContainer homology_container;
  if (boost::filesystem::exists(homology_file)) {
    std::ifstream returnstream(homology_file.string().c_str());
    if (returnstream.good()) {
      boost::archive::text_iarchive ia(returnstream);
      ia >> homology_container;
    } else {
      ELOG(2, "Failed to parse from " << homology_file.string() << ".");
    }
    returnstream.close();
  } else {
    LOG(2, "Could not open " << homology_file.string()
        << ", creating empty container.");
  }

  /// append all fragments to a HomologyContainer
  HomologyContainer::container_t values;
  const size_t FragmentCounter = shortrangeresults.Result_perIndexSet_Energy.size();

  // convert KeySetContainer to IndexSetContainer
  IndexSetContainer::ptr ForceContainer(new IndexSetContainer(shortrangeresults.getForceKeySet()));
  const IndexSetContainer::Container_t &Indices = shortrangeresults.getContainer();
  const IndexSetContainer::Container_t &ForceIndices = ForceContainer->getContainer();
  IndexSetContainer::Container_t::const_iterator iter = Indices.begin();
  IndexSetContainer::Container_t::const_iterator forceiter = ForceIndices.begin();
  /// go through all fragments
  for (;iter != Indices.end(); ++iter, ++forceiter) // go through each IndexSet
  {
    /// create new graph entry in HomologyContainer which is (key,value) type
    LOG(1, "INFO: Creating new graph with " << **forceiter << ".");
    HomologyGraph graph(**forceiter);
    LOG(2, "DEBUG: Created graph " << graph << ".");
    const IndexSet::ptr &index = *iter;
    HomologyContainer::value_t value;

    // obtain fragment as key
    std::map<IndexSet::ptr, std::pair< MPQCDataFragmentMap_t, MPQCDataFragmentMap_t> >::const_iterator fragmentiter
        = longrangeresults.Result_perIndexSet_Fragment.find(index);
    ASSERT( fragmentiter != longrangeresults.Result_perIndexSet_Fragment.end(),
        "appendToHomologyFile() - cannot find index "+toString(*index)
        +" in FragmentResults.");
    value.first = boost::fusion::at_key<MPQCDataFused::fragment>(fragmentiter->second.first);

    // obtain energy as value
    std::map<IndexSet::ptr, std::pair< MPQCDataEnergyMap_t, MPQCDataEnergyMap_t> >::const_iterator energyiter
        = shortrangeresults.Result_perIndexSet_Energy.find(index);
    ASSERT( energyiter != shortrangeresults.Result_perIndexSet_Energy.end(),
        "appendToHomologyFile() - cannot find index "+toString(*index)
        +" in FragmentResults.");
    // value.second = boost::fusion::at_key<MPQCDataFused::energy_total>(energyiter->second.first); // values
    value.second = boost::fusion::at_key<MPQCDataFused::energy_total>(energyiter->second.second); // contributions
    values.insert( std::make_pair( graph, value) );
  }
  homology_container.insert(values);

  LOG(1, "INFO: Listing all present atomic ids ...");
  std::stringstream output;
  for (World::AtomIterator iter = World::getInstance().getAtomIter();
      iter != World::getInstance().atomEnd(); ++iter)
    output << (*iter)->getId() << " ";
  LOG(1, "INFO: { " << output.str() << "} .");

  // for debugging: print container
  LOG(1, "INFO: Listing all present homologies ...");
  for (HomologyContainer::container_t::const_iterator iter =
      homology_container.begin(); iter != homology_container.end(); ++iter) {
    LOG(1, "INFO: graph " << iter->first << " has Fragment "
        << iter->second.first << " and associated energy " << iter->second.second << ".");
  }

  /// store homology container again
  std::ofstream outputstream(homology_file.string().c_str());
  if (outputstream.good()) { // check if opened
    boost::archive::text_oarchive oa(outputstream);
    oa << homology_container;
    if (outputstream.fail()) { // check if correctly written
      LOG(1, "Failed to write to file " << homology_file.string() << ".");
      return false;
    } else
      outputstream.close();
  } else {
    LOG(1, "Failed to open file " << homology_file.string()
        << " for writing.");
    return false;
  }
  return true;
}

Action::state_ptr FragmentationFragmentationAutomationAction::performCall() {
  boost::asio::io_service io_service;

  // TODO: Have io_service run in second thread and merge with current again eventually

  size_t Exitflag = 0;
  std::map<JobId_t, MPQCData> fragmentData;
  {
    const size_t NumberJobs = FragmentJobQueue::getInstance().size();
    MPQCFragmentController mpqccontroller(io_service);
    mpqccontroller.setHost(params.host.get());
    mpqccontroller.setPort(params.port.get());
    // Phase One: obtain ids
    mpqccontroller.requestIds(NumberJobs);

    // Phase Two: add MPQCJobs and send
    const size_t NoJobs = mpqccontroller.addJobsFromQueue(
        params.DoLongrange.get() ? MPQCData::DoSampleDensity : MPQCData::DontSampleDensity);
    LOG(1, "INFO: Added " << NoJobs << " from FragmentJobsQueue.");
    mpqccontroller.run();

    // Phase Three: calculate result
    mpqccontroller.waitforResults(NumberJobs);
    mpqccontroller.getResults(fragmentData);

    Exitflag += mpqccontroller.getExitflag();
  }

  FragmentationShortRangeResults shortrangeresults(
      fragmentData,
      FragmentJobQueue::getInstance().getKeySets(),
      FragmentJobQueue::getInstance().getFullKeySets());
  shortrangeresults(fragmentData);
  printReceivedShortResults(shortrangeresults);

#ifdef HAVE_VMG
  if (params.DoLongrange.get()) {
  if ( World::getInstance().getAllAtoms().size() == 0) {
    ELOG(1, "Please load the full molecule into the world before starting this action.");
    return Action::failure;
  }

  // obtain combined charge density
  FragmentationChargeDensity summedChargeDensity(
      fragmentData,
      FragmentJobQueue::getInstance().getKeySets());
  const std::vector<SamplingGrid> full_sample = summedChargeDensity.getFullSampledGrid();
  LOG(1, "INFO: There are " << fragmentData.size() << " short-range and "
      << full_sample.size() << " level-wise long-range jobs.");

  // Phase Four: obtain more ids
  std::map<JobId_t, VMGData> longrangeData;
  {
    VMGFragmentController vmgcontroller(io_service);
    vmgcontroller.setHost(params.host.get());
    vmgcontroller.setPort(params.port.get());
    const size_t NoJobs = fragmentData.size()+full_sample.size();
    vmgcontroller.requestIds(NoJobs);

    // Phase Five: create VMGJobs
    const size_t near_field_cells = params.near_field_cells.get();
    const size_t interpolation_degree = params.interpolation_degree.get();
    if (!vmgcontroller.createLongRangeJobs(
        fragmentData,
        full_sample,
        summedChargeDensity.getFragment(),
        near_field_cells,
        interpolation_degree))
      return Action::failure;

    // Phase Six: calculate result
    vmgcontroller.waitforResults(NoJobs);
    vmgcontroller.getResults(longrangeData);
    ASSERT( NoJobs == longrangeData.size(),
        "FragmentationFragmentationAutomationAction::performCall() - number of MPQCresults+"
        +toString(full_sample.size())+"="+toString(NoJobs)
        +" and VMGresults "+toString(longrangeData.size())+" don't match.");
    Exitflag += vmgcontroller.getExitflag();
  }

  // remove full solution corresponding to full_sample from map (must be highest ids), has to be treated extra
  std::map<JobId_t, VMGData>::iterator iter = longrangeData.end();
  for (size_t i=0;i<full_sample.size();++i)
    --iter;
  std::map<JobId_t, VMGData>::iterator remove_iter = iter;
  std::vector<VMGData> fullsolutionData;
  for (; iter != longrangeData.end(); ++iter)
    fullsolutionData.push_back(iter->second);
  longrangeData.erase(remove_iter, longrangeData.end());

  // Final phase: sum up and print result
  FragmentationLongRangeResults longrangeresults(
      fragmentData,
      longrangeData,
      FragmentJobQueue::getInstance().getKeySets(),
      FragmentJobQueue::getInstance().getFullKeySets());
  longrangeresults(
      fragmentData,
      longrangeData,
      fullsolutionData,
      full_sample);
  printReceivedFullResults(longrangeresults);

  // append all keysets to homology file
  if ((Exitflag == 0) && (!params.homology_file.get().empty())) {
    const boost::filesystem::path &homology_file = params.homology_file.get();
    if (homology_file.string() != "") {
      LOG(1, "INFO: Appending HomologyGraphs to file " << homology_file.string() << ".");
      if (!appendToHomologyFile(homology_file, shortrangeresults, longrangeresults))
        Exitflag = 1;
    }
  }

  std::map<JobId_t, std::string> debugData;
  {
    if (!full_sample.empty()) {
      // create debug jobs for each level to print the summed-up potential to vtk files
      VMGDebugGridFragmentController debugcontroller(io_service);
      debugcontroller.setHost(params.host.get());
      debugcontroller.setPort(params.port.get());
      debugcontroller.requestIds(full_sample.size());
      if (!debugcontroller.createDebugJobs(full_sample))
        return Action::failure;
      debugcontroller.waitforResults(full_sample.size());
      debugcontroller.getResults(debugData);
      Exitflag += debugcontroller.getExitflag();
    }
  }
  }
#endif

  // now clear all present jobs as we are done
  FragmentJobQueue::getInstance().clear();

  return (Exitflag == 0) ? Action::success : Action::failure;
}

Action::state_ptr FragmentationFragmentationAutomationAction::performUndo(Action::state_ptr _state) {
  return Action::success;
}

Action::state_ptr FragmentationFragmentationAutomationAction::performRedo(Action::state_ptr _state){
  return Action::success;
}

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

bool FragmentationFragmentationAutomationAction::shouldUndo() {
  return false;
}
/** =========== end of function ====================== */