/*
 *            Copyright 2009-2020 The VOTCA Development Team
 *                       (http://www.votca.org)
 *
 *      Licensed under the Apache License, Version 2.0 (the "License")
 *
 * You may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *              http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 *
 */

// Local VOTCA includes
#include "votca/xtp/orbreorder.h"
#include "votca/xtp/basisset.h"

namespace votca {
namespace xtp {

std::vector<Transposition> OrbReorder::computeTranspositions(
    std::vector<Index> vStart, std::vector<Index> vTarget) const {
  if (vStart.size() != vTarget.size()) {
    throw std::runtime_error(
        "Can't compute transpositions, reorder vectors have different "
        "size!\n");
  }
  std::vector<Transposition> transpositions;
  for (Index i = 0; i < static_cast<Index>(vStart.size()); i++) {
    std::vector<Index>::iterator it;
    it = std::find(vStart.begin(), vStart.end(), vTarget[i]);
    Index pos = std::distance(vStart.begin(), it);
    std::swap(vStart[i], vStart[pos]);
    if (i != pos) {
      transpositions.push_back(Transposition{i, pos});
    }
  }
  return transpositions;
}

std::vector<Index> OrbReorder::copySegment(const std::array<Index, 49>& input,
                                           Index start, Index size) const {
  return std::vector<Index>{input.begin() + start,
                            input.begin() + start + size};
}

OrbReorder::OrbReorder(std::array<Index, 49> reorder,
                       std::array<Index, 49> multipliers, bool reverse)
    : _multipliers(multipliers), _reorder(reorder), _reverse(reverse) {

  // Compute transpositions for every individual shell
  Index currentFunction = 0;
  for (int l = 0; l < 7; l++) {
    Index nrOfFunctions = NumFuncShell(static_cast<L>(l));
    if (!_reverse) {  // ordering from external to votca order
      _transpositions[l] = computeTranspositions(
          copySegment(_reorder, currentFunction, nrOfFunctions),
          copySegment(_votcaOrder, currentFunction, nrOfFunctions));
    } else {  // votca order to external order
      _transpositions[l] = computeTranspositions(
          copySegment(_votcaOrder, currentFunction, nrOfFunctions),
          copySegment(_reorder, currentFunction, nrOfFunctions));
    }
    currentFunction += nrOfFunctions;
  }
}

void OrbReorder::reorderOrbitals(Eigen::MatrixXd& moCoefficients,
                                 const AOBasis& basis) {

  for (const AOShell& shell : basis) {
    Index currentFunction = shell.getStartIndex();

    if (_reverse) {  // multiply first before reversing ordering
      // Get multiplier vector for shell
      std::vector<Index> shellmultiplier =
          copySegment(_multipliers, shell.getOffset(), shell.getNumFunc());

      // multiply shell
      for (Index i = 0; i < shell.getNumFunc(); i++) {
        moCoefficients.row(currentFunction + i) *= double(shellmultiplier[i]);
      }
    }

    // reorder shell
    Index l = static_cast<Index>(shell.getL());
    for (const Transposition& transposition : _transpositions[l]) {
      moCoefficients.row(currentFunction + transposition.first)
          .swap(moCoefficients.row(currentFunction + transposition.second));
    }

    if (!_reverse) {  // multiply after reordering
      // Get multiplier vector for shell
      std::vector<Index> shellmultiplier =
          copySegment(_multipliers, shell.getOffset(), shell.getNumFunc());

      // multiply shell
      for (Index i = 0; i < shell.getNumFunc(); i++) {
        moCoefficients.row(currentFunction + i) *= double(shellmultiplier[i]);
      }
    }
  }
}
void OrbReorder::reorderOperator(Eigen::MatrixXd& Matrixoperator,
                                 const AOBasis& basis) {
  // reorder rows first
  reorderOrbitals(Matrixoperator, basis);
  Matrixoperator.transposeInPlace();
  reorderOrbitals(Matrixoperator, basis);
  Matrixoperator.transposeInPlace();
}

}  // namespace xtp
}  // namespace votca