#include <ncbi.h>
#include "blastapp.h"

static int		consistcmp PROTO((HSPPtr PNTR, HSPPtr PNTR));

/* consistn -- determine the consistent "n" parameter for each HSP */
void
consistn(hlp)
	HitListPtr	hlp;
{
	HSPPtr	hspstack[100];
	HSPPtr	PNTR list0;
	register HSPPtr	hp;
	register HSPPtr	PNTR list, PNTR listmax;
	size_t	cnts[3], cnt;
	register int	sgn, i;
	unsigned long	limit;

	if (hlp->hspcnt <= 1) {
		if (hlp->hspcnt == 1)
			hlp->hp->n = 1;
		return;
	}

	cnts[0] = cnts[1] = cnts[2] = 0;
	for (hp = hlp->hp; hp != NULL; hp = hp->next) {
		sgn = SIGN(hp->frame);
		/* hits on different strands are counted separately */
		++cnts[sgn+1];
		hp->n = 1;
	}
	for (i = 0; i < DIM(cnts); ++i) {
		if ((cnt = cnts[i]) < 2)
			continue;
		if (cnt <= DIM(hspstack))
			list0 = hspstack;
		else
			list0 = (HSPPtr PNTR) ckalloc(sizeof(*list0) * cnt);

		/* Convert to a linear list for sorting by qsort() */
		sgn = i - 1;
		for (hp = hlp->hp, list = list0; hp != NULL; hp = hp->next)
			if (SIGN(hp->frame) == sgn)
				*list++ = hp;

		hsort((CharPtr)list0, cnt, sizeof(*list0), (int (*)())consistcmp);

		for (list = list0, listmax = list0 + cnt; list < listmax; ++list) {
			limit = 2 * (*list)->s_pos + (*list)->len;
			(*list)->n =
				/* count the number that come before the current HSP */
				s_pos_less(list0, list-list0, limit, (*list)->score)
				/* count the number that come after the current HSP */
				+ s_pos_greater(list+1, listmax-list-1, limit, (*list)->score)
				/* count ourselves */
				+ 1;
		}

		if (list0 != hspstack)
			mem_free(list0);
	}
}

static int
consistcmp(hpp1, hpp2)
	register HSPPtr	PNTR hpp1, PNTR hpp2;
{
	register unsigned long	pos1, pos2;

	pos1 = 2 * (*hpp1)->q_pos + (*hpp1)->len;
	pos2 = 2 * (*hpp2)->q_pos + (*hpp2)->len;
	if (pos1 > pos2)
		return 1;
	if (pos1 < pos2)
		return -1;

	pos1 = 2 * (*hpp1)->s_pos + (*hpp1)->len;
	pos2 = 2 * (*hpp2)->s_pos + (*hpp2)->len;
	if (pos1 > pos2)
		return 1;
	if (pos1 < pos2)
		return -1;
	return 0;
}

/* s_pos_less -- return the longest increasing subsequence of s_pos values */
int
s_pos_less(list, len, poslimit, scorelimit)
	HSPPtr	PNTR list;
	unsigned	len;
	unsigned long	poslimit;
	Score_t	scorelimit;
{
	HSPPtr	PNTR y;
	HSPPtr	PNTR listmax;
	unsigned long	PNTR c, cstack[100];
	unsigned long	value;
	register int	i, j;

	listmax = list + len;

	if (len <= DIM(cstack))
		c = cstack;
	else
		c = (unsigned long PNTR) ckalloc(len * sizeof(*c));

	c[0] = 0;
	j = 0;
	for (y = list; y < listmax; ++y) {
		if ((*y)->score < scorelimit)
			continue;
		value = 2 * (*y)->s_pos + (*y)->len;
		if (value >= poslimit)
			continue;
		for (i = j; c[i] > value; --i)
			;
		if (i == j)
			++j;
		c[i+1] = value;
	}
	/* clean-up */
	if (c == cstack)
		return j;
	mem_free(c);
	return j;
}
/* s_pos_greater -- return the longest increasing subsequence of s_pos values */
int
s_pos_greater(list, len, poslimit, scorelimit)
	HSPPtr	PNTR list;
	unsigned	len;
	unsigned long	poslimit;
	Score_t	scorelimit;
{
	HSPPtr	PNTR y;
	HSPPtr	PNTR listmax;
	unsigned long	PNTR c, cstack[100];
	unsigned long	value;
	register int	i, j;

	listmax = list + len;

	if (len <= DIM(cstack))
		c = cstack;
	else
		c = (unsigned long PNTR) ckalloc(len * sizeof(*c));

	c[0] = 0;
	j = 0;
	for (y = list; y < listmax; ++y) {
		if ((*y)->score < scorelimit)
			continue;
		value = 2 * (*y)->s_pos + (*y)->len;
		if (value <= poslimit)
			continue;
		for (i = j; c[i] > value; --i)
			;
		if (i == j)
			++j;
		c[i+1] = value;
	}
	/* clean-up */
	if (c == cstack)
		return j;
	mem_free(c);
	return j;
}