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

/*
 * bondgraph.cpp
 *
 *  Created on: Oct 29, 2009
 *      Author: heber
 */

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

#include "CodePatterns/MemDebug.hpp"

#include <iostream>

#include "Atom/atom.hpp"
#include "Bond/bond.hpp"
#include "Graph/BondGraph.hpp"
#include "Box.hpp"
#include "Element/element.hpp"
#include "CodePatterns/Info.hpp"
#include "CodePatterns/Log.hpp"
#include "CodePatterns/Range.hpp"
#include "CodePatterns/Verbose.hpp"
#include "molecule.hpp"
#include "Element/periodentafel.hpp"
#include "Fragmentation/MatrixContainer.hpp"
#include "LinearAlgebra/Vector.hpp"
#include "World.hpp"
#include "WorldTime.hpp"

const double BondGraph::BondThreshold = 0.4;   //!< CSD threshold in bond check which is the width of the interval whose center is the sum of the covalent radii

BondGraph::BondGraph() :
    BondLengthMatrix(NULL),
    IsAngstroem(true)
{}

BondGraph::BondGraph(bool IsA) :
    BondLengthMatrix(NULL),
    IsAngstroem(IsA)
{}

BondGraph::~BondGraph()
{
  CleanupBondLengthTable();
}

void BondGraph::CleanupBondLengthTable()
{
  if (BondLengthMatrix != NULL) {
    delete(BondLengthMatrix);
  }
}

bool BondGraph::LoadBondLengthTable(
    std::istream &input)
{
  Info FunctionInfo(__func__);
  bool status = true;
  MatrixContainer *TempContainer = NULL;

  // allocate MatrixContainer
  if (BondLengthMatrix != NULL) {
    LOG(1, "MatrixContainer for Bond length already present, removing.");
    delete(BondLengthMatrix);
    BondLengthMatrix = NULL;
  }
  TempContainer = new MatrixContainer;

  // parse in matrix
  if ((input.good()) && (status = TempContainer->ParseMatrix(input, 0, 1, 0))) {
    LOG(1, "Parsing bond length matrix successful.");
  } else {
    ELOG(1, "Parsing bond length matrix failed.");
    status = false;
  }

  if (status) // set to not NULL only if matrix was parsed
    BondLengthMatrix = TempContainer;
  else {
    BondLengthMatrix = NULL;
    delete(TempContainer);
  }
  return status;
}

double BondGraph::GetBondLength(
    int firstZ,
    int secondZ) const
{
  double return_length;
  if ((firstZ < 0) || (firstZ >= (int)BondLengthMatrix->Matrix[0].size()))
    return -1.;
  if ((secondZ < 0) || (secondZ >= (int)BondLengthMatrix->Matrix[0][firstZ].size()))
    return -1.;
  if (BondLengthMatrix == NULL) {
    return_length = -1.;
  } else {
    return_length = BondLengthMatrix->Matrix[0][firstZ][secondZ];
  }
  LOG(4, "INFO: Request for length between " << firstZ << " and "
      << secondZ << ": " << return_length << ".");
  return return_length;
}

range<double> BondGraph::CovalentMinMaxDistance(
    const element * const Walker,
    const element * const OtherWalker) const
{
  range<double> MinMaxDistance(0.,0.);
  MinMaxDistance.first = OtherWalker->getCovalentRadius() + Walker->getCovalentRadius();
  MinMaxDistance.first *= (IsAngstroem) ? 1. : 1. / AtomicLengthToAngstroem;
  MinMaxDistance.last = MinMaxDistance.first + BondThreshold;
  MinMaxDistance.first -= BondThreshold;
  return MinMaxDistance;
}

range<double> BondGraph::BondLengthMatrixMinMaxDistance(
    const element * const Walker,
    const element * const OtherWalker) const
{
  ASSERT(BondLengthMatrix, "BondGraph::BondLengthMatrixMinMaxDistance() called without NULL BondLengthMatrix.");
  ASSERT(Walker, "BondGraph::BondLengthMatrixMinMaxDistance() - illegal element given.");
  ASSERT(OtherWalker, "BondGraph::BondLengthMatrixMinMaxDistance() - illegal other element given.");
  range<double> MinMaxDistance(0.,0.);
  MinMaxDistance.first = GetBondLength(Walker->getAtomicNumber()-1, OtherWalker->getAtomicNumber()-1);
  MinMaxDistance.first *= (IsAngstroem) ? 1. : 1. / AtomicLengthToAngstroem;
  MinMaxDistance.last = MinMaxDistance.first + BondThreshold;
  MinMaxDistance.first -= BondThreshold;
  return MinMaxDistance;
}

range<double> BondGraph::getMinMaxDistance(
    const element * const Walker,
    const element * const OtherWalker) const
{
  range<double> MinMaxDistance(0.,0.);
  if (BondLengthMatrix == NULL) {// safety measure if no matrix has been parsed yet
    LOG(5, "DEBUG: Using Covalent radii criterion for [min,max) distances.");
    MinMaxDistance = CovalentMinMaxDistance(Walker, OtherWalker);
  } else {
    LOG(5, "DEBUG: Using tabulated bond table criterion for [min,max) distances.");
    MinMaxDistance = BondLengthMatrixMinMaxDistance(Walker, OtherWalker);
  }
  return MinMaxDistance;
}

range<double> BondGraph::getMinMaxDistance(
    const BondedParticle * const Walker,
    const BondedParticle * const OtherWalker) const
{
  return getMinMaxDistance(Walker->getType(), OtherWalker->getType());
}

range<double> BondGraph::getMinMaxDistanceSquared(
    const BondedParticle * const Walker,
    const BondedParticle * const OtherWalker) const
{
  // use non-squared version
  range<double> MinMaxDistance(getMinMaxDistance(Walker, OtherWalker));
  // and square
  MinMaxDistance.first *= MinMaxDistance.first;
  MinMaxDistance.last *= MinMaxDistance.last;
  return MinMaxDistance;
}

LinkedCell::LinkedCell_View BondGraph::getLinkedCell(const double max_distance) const
{
  return World::getInstance().getLinkedCell(max_distance);
}

std::set< const element *> BondGraph::getElementSetFromNumbers(const std::set<atomicNumber_t> &Set) const
{
  std::set< const element *> PresentElements;
  for(std::set<atomicNumber_t>::const_iterator iter = Set.begin(); iter != Set.end(); ++iter)
    PresentElements.insert( World::getInstance().getPeriode()->FindElement(*iter) );
  return PresentElements;
}

Box &BondGraph::getDomain() const
{
  return World::getInstance().getDomain();
}

unsigned int BondGraph::getTime() const
{
  return WorldTime::getTime();
}

bool BondGraph::operator==(const BondGraph &other) const
{
  if (IsAngstroem != other.IsAngstroem)
    return false;
  if (((BondLengthMatrix != NULL) && (other.BondLengthMatrix == NULL))
      || ((BondLengthMatrix == NULL) && (other.BondLengthMatrix != NULL)))
    return false;
  if ((BondLengthMatrix != NULL) && (other.BondLengthMatrix != NULL))
    if (*BondLengthMatrix != *other.BondLengthMatrix)
      return false;
  return true;
}

/** Corrects the bond degree by one at most if necessary.
 * \return number of corrections done
 */
int BondGraph::CorrectBondDegree(atom *_atom) const
{
  int NoBonds = 0;
  int OtherNoBonds = 0;
  int FalseBondDegree = 0;
  int CandidateDeltaNoBonds = -1;
  atom *OtherWalker = NULL;
  bond::ptr CandidateBond;

  NoBonds = _atom->CountBonds();
  LOG(3, "Walker " << *_atom << ": " << (int)_atom->getType()->getNoValenceOrbitals() << " > " << NoBonds << "?");
  const int DeltaNoBonds = (int)(_atom->getType()->getNoValenceOrbitals()) - NoBonds;
  if (DeltaNoBonds > 0) { // we have a mismatch, check all bonding partners for mismatch
    const BondList& ListOfBonds = _atom->getListOfBonds();
    for (BondList::const_iterator Runner = ListOfBonds.begin(); Runner != ListOfBonds.end(); (++Runner)) {
      OtherWalker = (*Runner)->GetOtherAtom(_atom);
      OtherNoBonds = OtherWalker->CountBonds();
      const int OtherDeltaNoBonds = (int)(OtherWalker->getType()->getNoValenceOrbitals()) - OtherNoBonds;
      LOG(3, "OtherWalker " << *OtherWalker << ": " << (int)OtherWalker->getType()->getNoValenceOrbitals() << " > " << OtherNoBonds << "?");
      if (OtherDeltaNoBonds > 0) { // check if possible candidate
        const BondList& OtherListOfBonds = OtherWalker->getListOfBonds();
        if ((CandidateBond == NULL)
            || (ListOfBonds.size() > OtherListOfBonds.size()) // pick the one with fewer number of bonds first
            || ((CandidateDeltaNoBonds < 0) || (CandidateDeltaNoBonds > OtherDeltaNoBonds)) ) // pick the one closer to fulfilling its valence first
        {
          CandidateDeltaNoBonds = OtherDeltaNoBonds;
          CandidateBond = (*Runner);
          LOG(3, "New candidate is " << *CandidateBond << ".");
        }
      }
    }
    if ((CandidateBond != NULL)) {
      CandidateBond->setDegree(CandidateBond->getDegree()+1);
      LOG(2, "Increased bond degree for bond " << *CandidateBond << ".");
    } else {
      ELOG(2, "Could not find correct degree for atom " << *_atom << ".");
      FalseBondDegree++;
    }
  }
  return FalseBondDegree;
};

/** Sets the weight of all connected bonds to one.
 */
void BondGraph::resetBondDegree(atom *_atom) const
{
  const BondList &ListOfBonds = _atom->getListOfBonds();
  for (BondList::const_iterator BondRunner = ListOfBonds.begin();
      BondRunner != ListOfBonds.end();
      ++BondRunner)
    (*BondRunner)->setDegree(1);
};