#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