#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 ranit_type, typename edge_eval, typename containers_reader>
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 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);

	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);
	std::vector<diff_type> minPenaltyNext(last);
	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)();

			for (bIdx = nextLowerBound, b = begin + nextLowerBound; furthest[aIdx] >= bIdx; ++b, ++bIdx) {
				assert(minPenaltyToEnd[bIdx] < (std::numeric_limits<penalty_type>::max)());
				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