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

BLAST_KarlinBlkPtr
lambdak(sp, sbp)
	BLAST_StrPtr	sp;
	BLAST_ScoreBlkPtr	sbp;
{
	BLAST_KarlinBlkPtr	kbp = NULL;
	BLAST_ResCompPtr	rcp;
	BLAST_ResFreqPtr	rfp, stdrfp;
	BLAST_ScoreFreqPtr	sfp;
	BLAST_AlphaMapPtr	amp;
	BLAST_Letter	ambiguous_chr;

	amp = sp->ap->inmap;
	rcp = BlastResCompNew(sp->ap);
	BlastResCompStr(rcp, sp, TRUE);

	/* Only count unambiguous residues in the effective sequence length */
	ambiguous_chr = NULLB;
	if (sp->ap->alphatype == BLAST_ALPHATYPE_AMINO_ACID)
		ambiguous_chr = 'X';
	if (sp->ap->alphatype == BLAST_ALPHATYPE_NUCLEIC_ACID)
		ambiguous_chr = 'N';
	if (ambiguous_chr != NULLB && BlastAlphaMapTst(amp, ambiguous_chr))
		BlastResCompSet(rcp, BlastAlphaMapChr(amp, ambiguous_chr), 0);
	if (BlastAlphaMapTst(amp, '-'))
		BlastResCompSet(rcp, BlastAlphaMapChr(amp, '-'), 0);
	sp->efflen = BlastResCompSum(rcp);

	rfp = BlastResFreqNew(sp->ap);
	BlastResFreqResComp(rfp, rcp);
	sfp = BlastScoreFreqNew(sbp->loscore, sbp->hiscore);
	stdrfp = BlastResFreqFind(sbp->a2);
	BlastScoreFreqCalc(sfp, sbp, rfp, stdrfp);

	if (sfp->score_avg >= 0.) {
		if (sp->ap->alphatype == BLAST_ALPHATYPE_AMINO_ACID) {
			if (sp->frame == 0)
				warning("Invalid (non-negative) expected score of %lg with the %s matrix.  The Karlin-Altschul K, Lambda and H parameters could not be computed.",
				(double)sfp->score_avg, sbp->name);
			else
				warning("Invalid (non-negative) expected score of %lg for the %+d reading frame with the %s matrix.  The Karlin-Altschul K, Lambda and H parameters could not be computed.",
					(double)sfp->score_avg, sp->frame, sbp->name);
		}
		else
			warning("Invalid (non-negative) expected score of %lg for the %+d strand with the %s matrix.  The Karlin-Altschul K, Lambda and H parameters could not be computed.",
			(double)sfp->score_avg, sp->frame, sbp->name);
		goto Error;
	}

	kbp = BlastKarlinBlkNew();
	kbp->q_frame = sp->frame;
	kbp->sbp = sbp;

	if (BlastKarlinBlkCalc(kbp, sfp) != BLAST_ERR_NONE) {
		BlastKarlinBlkDestruct(kbp);
		kbp = NULL;
		if (sp->ap->alphatype == BLAST_ALPHATYPE_AMINO_ACID) {
			if (sp->frame == 0)
				warning("Could not calculate Karlin-Altschul K, Lambda and H parameters for the %s matrix.",
					sbp->name);
			else
				warning("Could not calculate Karlin-Altschul K, Lambda and H parameters in the %+d reading with the %s matrix.",
					sp->frame, sbp->name);
		}
		else
			warning("Could not calculate Karlin-Altschul K, Lambda and H parameters for %+d strand with the %s matrix.",
				sp->frame, sbp->name);
		goto Error;
	}

Error:
	BlastResCompDestruct(rcp);
	BlastResFreqDestruct(rfp);
	BlastScoreFreqDestruct(sfp);

	return kbp;
}