#pragma once

#ifndef TCG_SEQUENCE_OPS_HPP

#define TCG_SEQUENCE_OPS_HPP

#include "../sequence_ops.h"

#ifdef min

#undef min

#endif

namespace tcg {

namespace sequence_ops {

//***********************************************************************************

//    Minimal Path Functions

//***********************************************************************************

template <typename containers_reader="" edge_eval,="" ranit_type,="" typename=""></typename>

bool minimalPath(ranit_type begin, ranit_type end, edge_eval &eval,

                 containers_reader &output) {

  typedef typename ranit_type::difference_type diff_type;

  typedef typename edge_eval::penalty_type penalty_type;

  ranit_type a, b;

  diff_type i, j, m, n = end - begin, last = n - 1;

  // Precache the longest edge possible from each vertex, imposing that furthest

  // nodes have increasing indices.

  std::vector<diff_type> furthest(n);</diff_type>

  diff_type currFurthest = furthest[last] = last;

  for (i = last - 1; i >= 0; --i) {

    currFurthest = furthest[i] =

        std::min(eval.furthestFrom(begin + i) - begin, currFurthest);

    if (currFurthest == i)

      return false;  // There exists no path from start to end - quit

  }

  // Iterate from begin to end, using the maximum step allowed. The number of

  // iterations is the number of output edges.

  for (i = 0, m = 0; i < last; i = furthest[i], ++m)

    ;

  // Also, build the iteration sequence. It will define the upper bounds for the

  // k-th vertex of the output.

  std::vector<diff_type> upperBound(m + 1);</diff_type>

  for (i = 0, j = 0; i <= m; j = furthest[j], ++i) upperBound[i] = j;

  // Now, the body of the algorithm

  std::vector<penalty_type> minPenaltyToEnd(n);</penalty_type>

  std::vector<diff_type> minPenaltyNext(last);</diff_type>

  diff_type aIdx, bIdx;

  penalty_type newPenalty;

  minPenaltyToEnd[last] = 0;

  diff_type nextLowerBound;

  for (j = m - 1, nextLowerBound = last; j >= 0; --j) {

    // Build the minimal penalty to end (also storing the next iterator

    // leading to it) from each vertex of the polygon, assuming the minimal

    // number of edges from the vertex to end.

    // The j-th polygon vertex must lie in the [lowerBound, upperBound[j]]

    // interval, whereas the (j+1)-th will be in [nextLowerBound,

    // upperBound[j+1]].

    // Please, observe that we always have upperBound[j] < nextLowerBound due

    // to the minimal edge count constraint.

    for (aIdx = upperBound[j]; aIdx >= 0 && furthest[aIdx] >= nextLowerBound;

         --aIdx) {

      a = begin + aIdx;

      // Search the vertex next to a which minimizes the penalty to end - and

      // store it.

      minPenaltyToEnd[aIdx] = (std::numeric_limits<penalty_type>::max)();</penalty_type>

      for (bIdx = nextLowerBound, b = begin + nextLowerBound;

           furthest[aIdx] >= bIdx; ++b, ++bIdx) {

        assert(minPenaltyToEnd[bIdx] <

               (std::numeric_limits<penalty_type>::max)());</penalty_type>

        newPenalty = eval.penalty(a, b) + minPenaltyToEnd[bIdx];

        if (newPenalty < minPenaltyToEnd[aIdx]) {

          minPenaltyToEnd[aIdx] = newPenalty;

          minPenaltyNext[aIdx]  = bIdx;

        }

      }

    }

    // Update

    nextLowerBound = aIdx;

    ++nextLowerBound;

  }

  // Finally, build the output sequence

  output.openContainer(begin);

  for (i = minPenaltyNext[0], j = 1; j < m; i = minPenaltyNext[i], ++j)

    output.addElement(begin + i);

  output.addElement(begin + last);

  output.closeContainer();

  return true;

}

}

}  // namespace tcg::sequence_ops

#endif  // TCG_SEQUENCE_OPS_HPP