source: src/Actions/FragmentationAction/FragmentationAutomationAction.cpp@ 376a3b

/*
 * 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 <boost/mpl/remove.hpp>
#include <boost/lambda/lambda.hpp>

#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/FragmentationResults.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/SetValues/Fragment.hpp"
#include "Fragmentation/SetValues/Histogram.hpp"
#include "Fragmentation/SetValues/IndexedVectors.hpp"
#include "Fragmentation/HydrogenSaturation_enum.hpp"
#include "Fragmentation/KeySet.hpp"
#include "Fragmentation/KeySetsContainer.hpp"
#include "Fragmentation/Summation/OrthogonalSumUpPerLevel.hpp"
#include "Fragmentation/Summation/writeTable.hpp"
#include "Graph/DepthFirstSearchAnalysis.hpp"
#include "Helpers/defs.hpp"
#include "Jobs/MPQCJob.hpp"
#include "Jobs/MPQCData.hpp"
#include "Jobs/MPQCData_printKeyNames.hpp"
#include "Jobs/Grid/SamplingGrid.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 <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;
}

/** Helper function to get number of atoms somehow.
 *
 * Here, we just parse the number of lines in the adjacency file as
 * it should correspond to the number of atoms, except when some atoms
 * are not bonded, but then fragmentation makes no sense.
 *
 * @param path path to the adjacency file
 */
size_t getNoAtomsFromAdjacencyFile(const std::string &path)
{
  size_t NoAtoms = 0;

  // parse in special file to get atom count (from line count)
  std::string filename(path);
  filename += FRAGMENTPREFIX;
  filename += ADJACENCYFILE;
  std::ifstream adjacency(filename.c_str());
  if (adjacency.fail()) {
    LOG(0, endl << "getNoAtomsFromAdjacencyFile() - Unable to open " << filename << ", is the directory correct?");
    return false;
  }
  std::string buffer;
  while (getline(adjacency, buffer))
    NoAtoms++;
  LOG(1, "INFO: There are " << NoAtoms << " atoms.");

  return NoAtoms;
}



/** Place results from FragmentResult into EnergyMatrix and ForceMatrix.
 *
 * @param fragmentData MPQCData resulting from the jobs
 * @param MatrixNrLookup Lookup up-map from job id to fragment number
 * @param FragmentCounter total number of fragments
 * @param NoAtoms total number of atoms
 * @param Energy energy matrix to be filled on return
 * @param Force force matrix to be filled on return
 * @return true - everything ok, false - else
 */
bool putResultsintoMatrices(
    const std::map<JobId_t, MPQCData> &fragmentData,
    const std::map< JobId_t, size_t > &MatrixNrLookup,
    const size_t FragmentCounter,
    const size_t NoAtoms,
    EnergyMatrix &Energy,
    ForceMatrix &Force)
{
  for (std::map<JobId_t, MPQCData>::const_iterator dataiter = fragmentData.begin();
      dataiter != fragmentData.end(); ++dataiter) {
    const MPQCData &extractedData = dataiter->second;
    const JobId_t &jobid = dataiter->first;
    std::map< JobId_t, size_t >::const_iterator nriter = MatrixNrLookup.find(jobid);
    ASSERT( nriter != MatrixNrLookup.end(),
        "putResultsintoMatrices() - MatrixNrLookup does not contain id "
        +toString(jobid)+".");
    // place results into EnergyMatrix ...
    {
      MatrixContainer::MatrixArray matrix;
      matrix.resize(1);
      matrix[0].resize(1, extractedData.energies.total);
      if (!Energy.AddMatrix(
          std::string("MPQCJob ")+toString(jobid),
          matrix,
          nriter->second)) {
        ELOG(1, "Adding energy matrix failed.");
        return false;
      }
    }
    // ... and ForceMatrix (with two empty columns in front)
    {
      MatrixContainer::MatrixArray matrix;
      const size_t rows = extractedData.forces.size();
      matrix.resize(rows);
      for (size_t i=0;i<rows;++i) {
        const size_t columns = 2+extractedData.forces[i].size();
        matrix[i].resize(columns, 0.);
  //      for (size_t j=0;j<2;++j)
  //        matrix[i][j] = 0.;
        for (size_t j=2;j<columns;++j)
          matrix[i][j] = extractedData.forces[i][j-2];
      }
      if (!Force.AddMatrix(
          std::string("MPQCJob ")+toString(jobid),
          matrix,
          nriter->second)) {
        ELOG(1, "Adding force matrix failed.");
        return false;
      }
    }
  }
  // add one more matrix (not required for energy)
  MatrixContainer::MatrixArray matrix;
  matrix.resize(1);
  matrix[0].resize(1, 0.);
  if (!Energy.AddMatrix(std::string("MPQCJob total"), matrix, FragmentCounter))
    return false;
  // but for energy because we need to know total number of atoms
  matrix.resize(NoAtoms);
  for (size_t i = 0; i< NoAtoms; ++i)
    matrix[i].resize(2+NDIM, 0.);
  if (!Force.AddMatrix(std::string("MPQCJob total"), matrix, FragmentCounter))
    return false;

  return true;
}
/** Print MPQCData from received results.
 *
 * @param fragmentData MPQCData resulting from the jobs, each associated to a job
 * @param KeySetFilename filename with keysets to associate forces correctly
 * @param NoAtoms total number of atoms
 * @param full_sample summed up charge density of electrons from fragments on return
 * @param full_fragment summed up positions and charges of nuclei from fragments on return
 */
bool sumUpChargeDensity(
    const std::map<JobId_t,MPQCData> &fragmentData,
    const std::string &KeySetFilename,
    std::vector<SamplingGrid> &full_sample,
    Fragment &full_fragment)
{
  // create a vector of all job ids
  std::vector<JobId_t> jobids;
  std::transform(fragmentData.begin(),fragmentData.end(),
      std::back_inserter(jobids),
      boost::bind( &std::map<JobId_t,MPQCData>::value_type::first, boost::lambda::_1 )
  );

  // create lookup from job nr to fragment number
  size_t FragmentCounter = 0;
  const std::map< JobId_t, size_t > MatrixNrLookup =
      createMatrixNrLookup(jobids, FragmentCounter);

  // initialise keysets
  KeySetsContainer KeySet;
  {
    // else needs keysets without hydrogens
    std::stringstream filename;
    filename << FRAGMENTPREFIX << KEYSETFILE;
    if (!KeySet.ParseKeySets(KeySetFilename, filename.str(), FragmentCounter)) return false;
  }

  /// prepare for OrthogonalSummation

  // convert KeySetContainer to IndexSetContainer
  IndexSetContainer::ptr container(new IndexSetContainer(KeySet));
  // create the map of all keysets
  SubsetMap::ptr subsetmap(new SubsetMap(*container));

  /// convert all MPQCData to MPQCDataMap_t
  std::vector<MPQCDataGridMap_t> Result_Grid_fused(
    OrthogonalSumUpPerLevel<MPQCDataGridMap_t, MPQCData, MPQCDataGridVector_t>(
        fragmentData, MatrixNrLookup, container, subsetmap));
  std::vector<MPQCDataFragmentMap_t> Result_Fragment_fused(
    OrthogonalSumUpPerLevel<MPQCDataFragmentMap_t, MPQCData, MPQCDataFragmentVector_t>(
        fragmentData, MatrixNrLookup, container, subsetmap));
  // obtain full grid
  full_sample.clear();
  full_sample.reserve(Result_Grid_fused.size());
  for (std::vector<MPQCDataGridMap_t>::const_iterator iter = ++Result_Grid_fused.begin();
      iter !=Result_Grid_fused.end();
      ++iter)
    full_sample.push_back(boost::fusion::at_key<MPQCDataFused::sampled_grid>((*iter)));
  full_fragment = boost::fusion::at_key<MPQCDataFused::fragment>(Result_Fragment_fused.back());

  return true;
}

/** Print MPQCData from received results.
 *
 * @param fragmentData MPQCData resulting from the jobs, associated to job id
 * @param KeySetFilename filename with keysets to associate forces correctly
 * @param NoAtoms total number of atoms
 * @param full_sample summed up charge from fragments on return
 */
bool printReceivedMPQCResults(
    const std::map<JobId_t, MPQCData> &fragmentData,
    const std::string &KeySetFilename,
    size_t NoAtoms)
{
  // create a vector of all job ids
  std::vector<JobId_t> jobids;
  std::transform(fragmentData.begin(),fragmentData.end(),
      std::back_inserter(jobids),
      boost::bind( &std::map<JobId_t,MPQCData>::value_type::first, boost::lambda::_1 )
  );

  // create lookup from job nr to fragment number
  size_t FragmentCounter = 0;
  const std::map< JobId_t, size_t > MatrixNrLookup=
      createMatrixNrLookup(jobids, FragmentCounter);


  // place results into maps
  EnergyMatrix Energy;
  ForceMatrix Force;
  if (!putResultsintoMatrices(fragmentData, MatrixNrLookup, FragmentCounter, NoAtoms, Energy, Force))
    return false;

  // initialise keysets
  KeySetsContainer KeySet;
  KeySetsContainer ForceKeySet;
  if (!Energy.InitialiseIndices()) return false;

  if (!Force.ParseIndices(KeySetFilename.c_str())) return false;

  {
    // else needs keysets without hydrogens
    std::stringstream filename;
    filename << FRAGMENTPREFIX << KEYSETFILE;
    if (!KeySet.ParseKeySets(KeySetFilename, filename.str(), FragmentCounter)) return false;
  }

  {
    // forces need keysets including hydrogens
    std::stringstream filename;
    filename << FRAGMENTPREFIX << FORCESFILE;
    if (!ForceKeySet.ParseKeySets(KeySetFilename, filename.str(), FragmentCounter)) return false;
  }

  // combine all found data
  if (!KeySet.ParseManyBodyTerms()) return false;

  EnergyMatrix EnergyFragments;
  ForceMatrix ForceFragments;
  if (!EnergyFragments.AllocateMatrix(Energy.Header, Energy.MatrixCounter, Energy.RowCounter, Energy.ColumnCounter)) return false;
  if (!ForceFragments.AllocateMatrix(Force.Header, Force.MatrixCounter, Force.RowCounter, Force.ColumnCounter)) return false;

  if(!Energy.SetLastMatrix(0., 0)) return false;
  if(!Force.SetLastMatrix(0., 2)) return false;

  for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
    // --------- sum up energy --------------------
    LOG(1, "INFO: Summing energy of order " << BondOrder+1 << " ...");
    if (!EnergyFragments.SumSubManyBodyTerms(Energy, KeySet, BondOrder)) return false;
    if (!Energy.SumSubEnergy(EnergyFragments, NULL, KeySet, BondOrder, 1.)) return false;

    // --------- sum up Forces --------------------
    LOG(1, "INFO: Summing forces of order " << BondOrder+1 << " ...");
    if (!ForceFragments.SumSubManyBodyTerms(Force, KeySet, BondOrder)) return false;
    if (!Force.SumSubForces(ForceFragments, KeySet, BondOrder, 1.)) return false;
  }

  // for debugging print resulting energy and forces
  LOG(1, "INFO: Resulting energy is " << Energy.Matrix[ FragmentCounter ][0][0]);
  std::stringstream output;
  for (int i=0; i< Force.RowCounter[FragmentCounter]; ++i) {
    for (int j=0; j< Force.ColumnCounter[FragmentCounter]; ++j)
      output << Force.Matrix[ FragmentCounter ][i][j] << " ";
    output << "\n";
  }
  LOG(1, "INFO: Resulting forces are " << std::endl << output.str());

  return true;
}

/** Print MPQCData from received results.
 *
 * @param results summed up results container
 */
void printReceivedFullResults(
    const FragmentationResults &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);

  const std::string gridresult =
      writeTable<VMGDataLongRangeMap_t, VMGDataLongRangeVector_t >()(
          results.Result_LongRangeIntegrated_fused, results.getMaxLevel());
  LOG(0, "LongRange table is \n" << gridresult);

  const std::string eigenvalueresult;
  LOG(0, "Eigenvalue table is \n" << eigenvalueresult);

  const std::string forceresult =
      writeTable<MPQCDataForceMap_t, MPQCDataForceVector_t>()(
          results.Result_Force_fused, results.getMaxLevel());
  LOG(0, "Force table is \n" << 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);
}


Action::state_ptr FragmentationFragmentationAutomationAction::performCall() {
  boost::asio::io_service io_service;
  MPQCFragmentController mpqccontroller(io_service);
  mpqccontroller.setHost(params.host.get());
  mpqccontroller.setPort(params.port.get());
  mpqccontroller.setLevel(params.level.get());
  VMGFragmentController vmgcontroller(io_service);
  vmgcontroller.setHost(params.host.get());
  vmgcontroller.setPort(params.port.get());
  VMGDebugGridFragmentController debugcontroller(io_service);
  debugcontroller.setHost(params.host.get());
  debugcontroller.setPort(params.port.get());

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

  // Phase One: obtain ids
  std::vector< boost::filesystem::path > jobfiles = params.jobfiles.get();
  mpqccontroller.requestIds(jobfiles.size());

  // Phase Two: create and add MPQCJobs
  if (!mpqccontroller.addJobsFromFiles(params.executable.get().string(), jobfiles))
    return Action::failure;

  // Phase Three: calculate result
  mpqccontroller.waitforResults(jobfiles.size());
  std::map<JobId_t, MPQCData> fragmentData;
  mpqccontroller.getResults(fragmentData);

#ifdef HAVE_VMG
  if (params.DoLongrange.get()) {
  ASSERT( World::getInstance().getAllAtoms().size() != 0,
      "FragmentationFragmentationAutomationAction::performCall() - please load the full molecule into the World before.");

  // obtain combined charge density
  LOG(1, "INFO: Parsing fragment files from " << params.path.get() << ".");
  std::vector<SamplingGrid> full_sample; // have charges from level 2 onward summed up
  Fragment full_fragment;
  sumUpChargeDensity(
      fragmentData,
      params.path.get(),
      full_sample,
      full_fragment);

  // Phase Four: obtain more ids
  vmgcontroller.requestIds(fragmentData.size()+full_sample.size());

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

  // Phase Six: calculate result
  vmgcontroller.waitforResults(fragmentData.size()+full_sample.size());
  std::map<JobId_t, VMGData> longrangeData;
  vmgcontroller.getResults(longrangeData);
  ASSERT( fragmentData.size()+full_sample.size() == longrangeData.size(),
      "FragmentationFragmentationAutomationAction::performCall() - number of MPQCresults+"
      +toString(full_sample.size())+" "+toString(fragmentData.size()+full_sample.size())
      +" and VMGresults "+toString(longrangeData.size())+" don't match.");

  // 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
  {
    FragmentationResults results(
        fragmentData,
        longrangeData,
        fullsolutionData,
        params.path.get(),
        getNoAtomsFromAdjacencyFile(params.path.get()),
        full_sample);

    LOG(1, "INFO: Parsing fragment files from " << params.path.get() << ".");
    printReceivedFullResults(results);

    if (!full_sample.empty()) {
      // create debug jobs for each level to print the summed-up potential to vtk files
      debugcontroller.requestIds(full_sample.size());
      if (!debugcontroller.createDebugJobs(full_sample))
        return Action::failure;
      debugcontroller.waitforResults(full_sample.size());
      std::map<JobId_t, std::string> debugData;
      debugcontroller.getResults(debugData);
    }
  }
  }
#else
  // Final phase: print result
  {
    LOG(1, "INFO: Parsing fragment files from " << params.path.get() << ".");
    printReceivedMPQCResults(
        fragmentData,
        params.path.get(),
        getNoAtomsFromAdjacencyFile(params.path.get()));
  }
#endif
  size_t Exitflag = vmgcontroller.getExitflag() + mpqccontroller.getExitflag();

  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 ====================== */