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

void
pvals(hlp, qlen, tlen, dblen)
	HitListPtr	hlp;
	unsigned long	qlen;
	unsigned long	tlen;
	unsigned long	dblen;
{
	extern double	Lambda, K, H;
	double	xtlen = tlen;
	double	xdblen = dblen;
	register HSPPtr	hp, minhp;
	register double	p;
	Score_t	highscore = SCORE_MIN, scoreray[3];
	int		sgn;
	double	e;

	scoreray[0] = scoreray[1] = scoreray[2] = 0;
	for (hp=minhp= hlp->hp; hp != NULL; hp = hp->next) {
		if (hp->score > highscore)
			highscore = hp->score;

		sgn = SIGN(hp->frame) + 1;
		scoreray[sgn] += hp->score;

		e = stoe(hp->score, Lambda, K, H, tlen, qlen);
		p = etopoissonp(hp->n, e);
		p = gapdecay(hp->n, p, gapdecayrate);

		/*
		Calculate the p-value adjusted for the length of the database
		sequence versus the total length of the database
		*/
		hp->P_val = etop(hp->expect = (p * xdblen)/xtlen);

		/* Remember which HSP had the smallest Poisson probability */
		if (hp->expect < minhp->expect)
			minhp = hp;
	}
	hlp->signifhsp = minhp;
	hlp->highscore = highscore;
	hlp->totalscore = MAX(scoreray[0], MAX(scoreray[1], scoreray[2]));
}


void
consist_pvals(hlp, qlen, tlen, dblen)
	HitListPtr	hlp;
	unsigned long	qlen;
	unsigned long	tlen;
	unsigned long	dblen;
{
	extern double	Lambda, K, H;
	double	xtlen = tlen;
	double	xdblen = dblen;
	register HSPPtr	hp, minhp;
	register double	p;
	Score_t	highscore = SCORE_MIN, scoreray[3];
	int		sgn;
	double	e;

	scoreray[0] = scoreray[1] = scoreray[2] = 0;
	for (hp=minhp= hlp->hp; hp != NULL; hp = hp->next) {
		if (hp->score > highscore)
			highscore = hp->score;

		sgn = SIGN(hp->frame) + 1;
		scoreray[sgn] += hp->score;

		e = stoe(hp->score, Lambda, K, H, tlen, qlen);
		p = etoconsistp(hp->n, e);
		p = gapdecay(hp->n, p, gapdecayrate);

		/*
		Calculate the p-value adjusted for the length of the database
		sequence versus the total length of the database
		*/
		hp->P_val = etop(hp->expect = (p * xdblen)/xtlen);

		/* Remember which HSP had the smallest Poisson probability */
		if (hp->expect < minhp->expect)
			minhp = hp;
	}
	hlp->signifhsp = minhp;
	hlp->highscore = highscore;
	hlp->totalscore = MAX(scoreray[0], MAX(scoreray[1], scoreray[2]));
}