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

static int cmp_hsp_by_midpt PROTO((HSPPtr PNTR, HSPPtr PNTR));
#define MIDPT_Q(A)	((A)->q_pos + (A)->len/2)
#define MIDPT_S(A)	((A)->s_pos + (A)->len/2)

void
pcnt_a(hlp)
	HitListPtr	hlp;
{
	HSPPtr	hlist[100], *hpp0, *hpp_by_pos;
	unsigned	iray[DIM(hlist)], PNTR ip;
	HSPPtr	hp, hp1, hp2;
	Score_t	score;
	unsigned	cnt, level, n, i, j;

	if (hlp->hp == NULL)
		return;

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

	if (cnt > DIM(hlist)) {
		hpp0 = (HSPPtr PNTR)ckalloc(sizeof(HSPPtr)*cnt);
		ip = (unsigned PNTR)ckalloc(sizeof(*ip)*cnt);
	}
	else {
		hpp0 = hlist;
		ip = iray;
	}

	hpp_by_pos = hpp0;

	for (hp = hlp->hp; hp != NULL; hp = hp->next) {
		*hpp_by_pos++ = hp;
		hp->n = 1;
	}
	hpp_by_pos = hpp0;

	(void) hsort((CharPtr)hpp_by_pos, (size_t)cnt, sizeof(*hpp_by_pos),
			(int (*)PROTO(()))&cmp_hsp_by_midpt);

	level = 0;
	ip[level] = 0;
	for (;;) {
		++level;
		for (ip[level] = ip[level-1]; ++ip[level] < cnt &&
				MIDPT_S(hpp_by_pos[ip[level]]) > MIDPT_S(hpp_by_pos[ip[level-1]);
				) {
		}
		if (ip[level] < cnt)
			continue;

		for (i=0; i < level; ++i) {
			score = hpp_by_pos[ip[i]]->score;
			for (j = n = 0; j < level; ++j) {
				if (hpp_by_pos[ip[j]]->score >= score)
					++n;
			}
			hpp_by_pos[ip[i]]->n = MAX(hpp_by_pos[ip[i]]->n, n);
		}

		for (; level > 0; ) {
			--level;
			for (; ++ip[level] < cnt; ) {
			}
		}
		if (level == 0)
			break;
	}


	if (hpp0 != hlist) {
		mem_free((VoidPtr)hpp0);
		mem_free((VoidPtr)ip);
	}
}


static int
cmp_hsp_by_midpt(h1, h2)
	HSPPtr	PNTR h1, PNTR h2;
{
	Coord_t	mid1, mid2;

	mid1 = MIDPT_Q(*h1);
	mid2 = MIDPT_Q(*h2);
	if (mid1 > mid2)
		return 1;
	if (mid1 < mid2)
		return -1;
	mid1 = MIDPT_S(*h1);
	mid2 = MIDPT_S(*h2);
	if (mid1 < mid2)
		return 1;
	if (mid1 > mid2)
		return -1;
	if ((*h1)->q_pos < (*h2)->q_pos)
		return 1;
	if ((*h1)->q_pos > (*h2)->q_pos)
		return -1;
	if ((*h1)->len < (*h2)->len)
		return 1;
	if ((*h2)->len > (*h2)->len)
		return -1;
	return 0;
}