/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  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/>.
 */

/*
 * FragmentationChargeDensity.cpp
 *
 *  Created on: Aug 31, 2012
 *      Author: heber
 */


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

#include "CodePatterns/MemDebug.hpp"

#include "FragmentationChargeDensity.hpp"

#include <boost/mpl/for_each.hpp>
#include <boost/mpl/remove.hpp>
#include <sstream>

#include "CodePatterns/Assert.hpp"
#include "CodePatterns/Log.hpp"

#include "Fragmentation/KeySetsContainer.hpp"
#include "Fragmentation/Automation/createMatrixNrLookup.hpp"
#include "Fragmentation/Automation/extractJobIds.hpp"
#include "Fragmentation/Automation/parseKeySetFile.hpp"
#include "Fragmentation/Summation/IndexSetContainer.hpp"
#include "Fragmentation/Summation/OrthogonalSumUpPerLevel.hpp"
#include "Fragmentation/Summation/SubsetMap.hpp"

#include "Helpers/defs.hpp"

FragmentationChargeDensity::FragmentationChargeDensity(
    const std::map<JobId_t,MPQCData> &fragmentData,
    const std::string &KeySetFilename)
{
  // 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;
  parseKeySetFile(KeySet, KeySetFilename, FragmentCounter, NonHydrogenKeySets);

  /// 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
  OrthogonalSumUpPerLevel<MPQCDataGridMap_t, MPQCData, MPQCDataGridVector_t>(
      fragmentData, MatrixNrLookup, container, subsetmap,
      Result_Grid_fused, Result_perIndexSet_Grid);
  OrthogonalSumUpPerLevel<MPQCDataFragmentMap_t, MPQCData, MPQCDataFragmentVector_t>(
      fragmentData, MatrixNrLookup, container, subsetmap,
      Result_Fragment_fused, Result_perIndexSet_Fragment);
}

std::vector<SamplingGrid> FragmentationChargeDensity::getFullSampledGrid()
{
  std::vector<SamplingGrid> full_sample;
  full_sample.reserve(Result_Grid_fused.size()-1);
  // skip level 1, as the saturation hydrogen do not yet cancel each other
  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)));
  return full_sample;
}

const Fragment &FragmentationChargeDensity::getFragment()
{
  return boost::fusion::at_key<MPQCDataFused::fragment>(Result_Fragment_fused.back());
}