/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  2012 University of Bonn. All rights reserved.
 * Copyright (C)  2013 Frederik Heber. All rights reserved.
 * Please see the COPYING file or "Copyright notice" in builder.cpp for details.
 *
 *
 *   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/>.
 */

/*
 * TrainingData.cpp
 *
 *  Created on: 15.10.2012
 *      Author: heber
 */

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

//#include "CodePatterns/MemDebug.hpp"

#include "TrainingData.hpp"

#include <algorithm>
#include <boost/bind.hpp>
#include <boost/foreach.hpp>
#include <boost/lambda/lambda.hpp>
#include <iostream>
#include <sstream>

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

#include "Fragmentation/EdgesPerFragment.hpp"
#include "Fragmentation/Summation/SetValues/Fragment.hpp"
#include "FunctionApproximation/FunctionArgument.hpp"
#include "FunctionApproximation/FunctionModel.hpp"
#include "FunctionApproximation/Extractors.hpp"

void TrainingData::operator()(const range_t &range) {
  for (HomologyContainer::const_iterator iter = range.first; iter != range.second; ++iter) {
    const HomologyGraph &graph = iter->first;
    const Fragment &fragment = iter->second.fragment;
    const FragmentationEdges::edges_t &edges = iter->second.edges;
    FunctionModel::arguments_t all_args = Extractors::gatherAllSymmetricDistances(
        fragment.getPositions(),
        fragment.getAtomicNumbers(),
        edges,
        DistanceVector.size()
        );
    DistanceVector.push_back( all_args );
    const double &energy = iter->second.energy;
    EnergyVector.push_back( FunctionModel::results_t(1, energy) );
    // filter distances out of list of all arguments
    FunctionModel::list_of_arguments_t args = filter(graph, all_args);
    LOG(3, "DEBUG: Filtered arguments are " << args << ".");
    ArgumentVector.push_back( args );
  }
}

const double TrainingData::getL2Error(const FunctionModel &model) const
{
  double L2sum = 0.;

  FilteredInputVector_t::const_iterator initer = ArgumentVector.begin();
  OutputVector_t::const_iterator outiter = EnergyVector.begin();
  for (; initer != ArgumentVector.end(); ++initer, ++outiter) {
    const FunctionModel::results_t result = model((*initer));
    const double temp = fabs((*outiter)[0] - result[0]);
    L2sum += temp*temp;
  }
  return L2sum;
}

const double TrainingData::getLMaxError(const FunctionModel &model) const
{
  double Lmax = 0.;
//  size_t maxindex = -1;
  FilteredInputVector_t::const_iterator initer = ArgumentVector.begin();
  OutputVector_t::const_iterator outiter = EnergyVector.begin();
  for (; initer != ArgumentVector.end(); ++initer, ++outiter) {
    const FunctionModel::results_t result = model((*initer));
    const double temp = fabs((*outiter)[0] - result[0]);
    if (temp > Lmax) {
      Lmax = temp;
//      maxindex = std::distance(
//          const_cast<const FunctionApproximation::inputs_t &>(ArgumentVector).begin(),
//          initer
//          );
    }
  }
  return Lmax;
}

const TrainingData::L2ErrorConfigurationIndexMap_t
TrainingData::getWorstFragmentMap(
      const FunctionModel &model,
      const range_t &range) const
{
  L2ErrorConfigurationIndexMap_t WorseFragmentMap;
  // fragments make it into the container in reversed order, hence count from top down
  size_t index= std::distance(range.first, range.second)-1;
  InputVector_t::const_iterator distanceiter = DistanceVector.begin();
  FilteredInputVector_t::const_iterator initer = ArgumentVector.begin();
  OutputVector_t::const_iterator outiter = EnergyVector.begin();
  for (; initer != ArgumentVector.end(); ++initer, ++outiter, ++distanceiter) {
    // calculate value from potential
    const FunctionModel::list_of_arguments_t &args = *initer;
    const FunctionModel::results_t result = model(args);
    const double energy = (*outiter)[0];

    // insert difference into map
    const double error = fabs(energy - result[0]);
    WorseFragmentMap.insert( std::make_pair( error, index-- ) );

    {
      // give only the distances in the debugging text
      std::stringstream streamargs;
      BOOST_FOREACH (argument_t arg, *distanceiter) {
        streamargs << " " << arg.distance;
      }
      LOG(2, "DEBUG: frag.#" << index+1 << "'s error is |" << energy << " - " << result[0]
          << "| = " << error << " for args " << streamargs.str() << ".");
    }
  }

  return WorseFragmentMap;
}

const TrainingData::DistanceEnergyTable_t TrainingData::getDistanceEnergyTable() const
{
  TrainingData::DistanceEnergyTable_t table;

  /// extract distance member variable from argument_t and first value from results_t
  OutputVector_t::const_iterator ergiter = EnergyVector.begin();
  for (InputVector_t::const_iterator iter = DistanceVector.begin();
      iter != DistanceVector.end(); ++iter, ++ergiter) {
    ASSERT( ergiter != EnergyVector.end(),
        "TrainingData::getDistanceEnergyTable() - less output than input values.");
    std::vector< double > values(iter->size(), 0.);
    // transform all distances
    const FunctionModel::arguments_t &args = *iter;
    std::transform(
        args.begin(), args.end(),
        values.begin(),
        boost::bind(&argument_t::distance, _1));

    // get first energy value
    values.push_back((*ergiter)[0]);

    // push as table row
    table.push_back(values);
  }

  return table;
}

const FunctionModel::results_t TrainingData::getTrainingOutputAverage() const
{
  if (EnergyVector.size() != 0) {
    FunctionApproximation::outputs_t::const_iterator outiter = EnergyVector.begin();
    FunctionModel::results_t result(*outiter);
    for (++outiter; outiter != EnergyVector.end(); ++outiter)
      for (size_t index = 0; index < (*outiter).size(); ++index)
        result[index] += (*outiter)[index];
    LOG(2, "DEBUG: Sum of EnergyVector is " << result << ".");
    const double factor = 1./EnergyVector.size();
    std::transform(result.begin(), result.end(), result.begin(),
        boost::lambda::_1 * factor);
    LOG(2, "DEBUG: Average EnergyVector is " << result << ".");
    return result;
  }
  return FunctionModel::results_t();
}

std::ostream &operator<<(std::ostream &out, const TrainingData &data)
{
  const TrainingData::InputVector_t &DistanceVector = data.getAllArguments();
  const TrainingData::OutputVector_t &EnergyVector = data.getTrainingOutputs();
  out << "(" << DistanceVector.size()
      << "," << EnergyVector.size() << ") data pairs: " << std::endl;
  FunctionApproximation::inputs_t::const_iterator initer = DistanceVector.begin();
  FunctionApproximation::outputs_t::const_iterator outiter = EnergyVector.begin();
  for (; initer != DistanceVector.end(); ++initer, ++outiter) {
    for (size_t index = 0; index < (*initer).size(); ++index)
       out << "(" << (*initer)[index].indices.first << "," << (*initer)[index].indices.second
          << ") " << (*initer)[index].distance;
    out << " with energy ";
    out << (*outiter);
    out << std::endl;
  }
  return out;
}