/* scaleswn.c */

/* $Name: fa34t26b4 $ - $Id: scaleswn.c,v 1.58 2006/06/24 02:34:14 wrp Exp $ */

/* as of 24 Sept, 2000 - scaleswn uses no global variables */

/*
	Provide statistical estimates using an extreme value distribution

	copyright (c) 1995, 1996, 2000 William R. Pearson

	This code provides multiple methods for scaling sequence
	similarity scores to correct for length effects.

	Currently, six methods are available:

	pst.zsflag = 0 - no scaling  (AVE_STATS)
	pst.zsflag = 1 - regression-scaled scores (REG_STATS)
	pst.zsflag = 2 - (revised) MLE Lmabda/K scaled scores (MLE_STATS)
	pst.zsflag = 3 - scaling using Altschul's parameters (AG_STATS)
	pst.zsflag = 4 - regression-scaled with iterative outlier removal (REGI_STATS)
	pst.zsflag = 5 = like 1, but length scaled variance (REG2_STATS)
	pst.zsflag = 6 = like 2, but uses lambda composition/scale (MLE2_STATS)
	pst.zsflag = 11 = 10 + 1 - use random shuffles, method 1
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>

#include <limits.h>

#include "defs.h"
#include "param.h"
#include "structs.h"
#ifndef PCOMPLIB
#include "mw.h"
#else
#include "p_mw.h"
#endif

#define MAXHIST 50
#define MAX_LLEN 200
#define LHISTC 5
#define VHISTC 5
#define MAX_SSCORE 300

#define LENGTH_CUTOFF 10 /* minimum database sequence length allowed, for fitting */

#define LN_FACT 10.0
#ifndef M_LN2
#define M_LN2 0.69314718055994530942
#endif
#define EULER_G 0.57721566490153286060
#define PI_SQRT6 1.28254983016186409554

#ifndef M_SQRT2
#define M_SQRT2 1.41421356237
#endif
#define LN200 5.2983173666
#define ZS_MAX 400.0	/* used to prevent underflow on some machines */
#define TOLERANCE 1.0e-12
#define TINY 1.0e-6

/* used by AVE_STATS, REG_STATS, REGI_STATS, REG2_STATS*/
struct rstat_str {
  double rho, rho_e, mu, mu_e, mean_var, var_e;  /* ?_e:std. error of ? */
/* used by REG2_STATS */
  double rho2, mu2, var_cutoff;
  int n_trimmed; /* excluded because of high z-score */
  int n1_trimmed, nb_trimmed, nb_tot; /* excluded because of bin */
};

/* used by AG_STATS, MLE_STATS */
struct ag_stat_str {
  double K, Lambda, H, a_n0f, a_n0;
};

/* used by MLE2_STATS */
struct mle2_stat_str {
  double a_n0;
  double mle2_a0, mle2_a1, mle2_a2, mle2_b1;
  double ave_comp, max_comp, ave_H;
};

struct pstat_str {
  double ngLambda, ngK, ngH;
  union {
    struct rstat_str rg;
    struct ag_stat_str ag;
    struct mle2_stat_str m2;
  } r_u;
};

#define AVE_STATS 0	/* no length effect, only mean/variance */
double find_zn(int score, double escore, int len, double comp, struct pstat_str *);

int proc_hist_n(struct stat_str *sptr, int n, 
		struct pstruct pst, struct hist_str *histp, int do_trim,
		struct pstat_str *);

#define REG_STATS 1	/* length-regression scaled */
#define REGI_STATS 4	/* length regression, iterative */
double find_zr(int score, double escore, int len, double comp, struct pstat_str *);
int proc_hist_r(struct stat_str *sptr, int n,
		struct pstruct pst, struct hist_str *histp,
		int do_trim, struct pstat_str *pu);

#define MLE_STATS 2	/* MLE for lambda, K */
double find_ze(int score, double escore, int len, double comp, struct pstat_str *);
int proc_hist_ml(struct stat_str *sptr, int n,
		 struct pstruct pst, struct hist_str *histp, int do_trim,
		 struct pstat_str *);

#define AG_STATS 3	/* Altschul-Gish parameters */
double find_za(int score, double escore, int len, double comp, struct pstat_str *);
int proc_hist_a(struct stat_str *sptr, int n,
		struct pstruct pst, struct hist_str *histp, int do_trim,
		struct pstat_str *);

#define REG2_STATS 5	/* length regression on mean + variance */
double find_zr2(int score, double escore, int len, double comp, struct pstat_str *);
int proc_hist_r2(struct stat_str *sptr, int n,
		 struct pstruct pst, struct hist_str *histp, int do_trim,
		 struct pstat_str *);

#define MLE2_STATS 6	/* MLE stats using comp(lambda) */
double find_ze2(int score, double escore, int length, double comp, struct pstat_str *);
int proc_hist_ml2(struct stat_str *sptr, int n,
		  struct pstruct pst, struct hist_str *histp, int do_trim,
		  struct pstat_str *);

#ifdef USE_LNSTATS
#define LN_STATS 2
double find_zl(int score, double escore, int len, double comp, struct pstat_str *);
int proc_hist_ln(struct stat_str *sptr, int n, 
		 struct pstruct pst, struct hist_str *histp, int do_trim,
		 struct pstat_str *);
#endif

/* scaleswn.c local variables that belong in their own structure */

double (*find_zp)(int score, double escore, int len, double comp, struct pstat_str *) = &find_zr;

/* void s_sort (double **ptr, int nbest); */
void ss_sort ( int *sptr, int n);

struct llen_str {
  int min, max;
  int max_score, min_score;
  int *hist;
  double *score_sums, *score2_sums;
  double *score_var;
  int max_length, min_length, zero_s;
  int fit_flag;
};

static void inithist(struct llen_str *, struct pstruct, int);
static void free_hist( struct llen_str *);
static void addhist(struct llen_str *, int, int, int);
static void prune_hist(struct llen_str *, int, int, int, long *);
void inithistz(int, struct hist_str *histp);
void addhistz(double zs, struct hist_str *histp);
void addhistzp(double zs, struct hist_str *histp);

static void fit_llen(struct llen_str *, struct rstat_str *);
static void fit_llen2(struct llen_str *, struct rstat_str *);
static void fit_llens(struct llen_str *, struct rstat_str *);

extern void sortbeste(struct beststr **bptr, int nbest);

/* void set_db_size(int, struct db_str *, struct hist_str *); */

#ifdef DEBUG
FILE *tmpf;
#endif

int
process_hist(struct stat_str *sptr, int nstats, 
	     struct mngmsg m_msg,
	     struct pstruct pst,
	     struct hist_str *histp,
	     struct pstat_str **ps_sp,
	     int do_hist)
{
  int zsflag, do_trim, i;
  struct pstat_str *ps_s;

  if (pst.zsflag < 0) {
    *ps_sp = NULL;
    return pst.zsflag;
  }

  if (*ps_sp == NULL) {
    if ((ps_s=(struct pstat_str *)calloc(1,sizeof(struct pstat_str)))==NULL) {
      fprintf(stderr," cannot allocate pstat_union: %ld\n",sizeof(struct pstat_str));
      exit(1);
    }
    else *ps_sp = ps_s;
  }
  else {
    ps_s = *ps_sp;
    memset(ps_s,0,sizeof(struct pstat_str));
  }

  ps_s->ngLambda = m_msg.Lambda;
  ps_s->ngK = m_msg.K;
  ps_s->ngH = m_msg.H;

  if (nstats < 10) pst.zsflag = AG_STATS;

  zsflag = pst.zsflag;

/*
#ifdef DEBUG
  if (pst.debug_lib) {
    tmpf=fopen("tmp_stats.res","w+");
    for (i=0; i<nstats; i++) fprintf(tmpf,"%d\t%d\n",sptr[i].score,sptr[i].n1);
    fclose(tmpf);
  }
#endif
*/

  if (zsflag >= 10) {
    zsflag -= 10;
    do_trim = 0;
  }
  else do_trim = 1;

#ifdef USE_LNSCALE
  if (zsflag==LN_STATS) {
    find_zp = &find_zl;
    pst.zsflag = proc_hist_ln(sptr, nstats, histp, do_trim, ps_s);
  }
#else
  if (zsflag==MLE_STATS) {
    find_zp = &find_ze;
    pst.zsflag = proc_hist_ml(sptr, nstats, pst, histp, do_trim, ps_s);
  }
#endif
  else if (zsflag==REG_STATS) {
    find_zp = &find_zr;
    pst.zsflag = proc_hist_r(sptr, nstats,pst, histp, do_trim,  ps_s);
  }
  else if (zsflag==AG_STATS) {
    find_zp = &find_za;
    pst.zsflag = proc_hist_a(sptr, nstats, pst, histp, do_trim,  ps_s);
  }
  else if (zsflag==REGI_STATS) {
    find_zp = &find_zr;
    pst.zsflag = proc_hist_r2(sptr,nstats, pst, histp, do_trim,  ps_s);
  }
  else if (zsflag==REG2_STATS) {
    find_zp = &find_zr2;
    pst.zsflag = proc_hist_r(sptr,nstats,pst, histp, do_trim,  ps_s);
  }
#if !defined(TFAST) && !defined(FASTX)
  else if (zsflag == MLE2_STATS) {
    find_zp = &find_ze2;
    pst.zsflag = proc_hist_ml2(sptr, nstats, pst, histp, do_trim,  ps_s);
  }
#endif
  else {	/* AVE_STATS */
    find_zp = &find_zn;
    pst.zsflag = proc_hist_n(sptr,nstats, pst, histp, do_trim,  ps_s);
  }

  if (!do_hist) {
    histp->entries = nstats; /* db->entries = 0; */
    inithistz(MAXHIST, histp);
    for (i = 0; i < nstats; i++) {
      if (sptr[i].n1 < 0) sptr[i].n1 = -sptr[i].n1;
      addhistz(find_zp(sptr[i].score,sptr[i].escore,sptr[i].n1,sptr[i].comp,ps_s),
	       histp);
    }
  }
  return pst.zsflag;
}

int
calc_thresh(struct pstruct pst, int nstats, 
	    double Lambda, double K, double H, double *zstrim)
{
  int max_hscore;
  double ave_n1, tmp_score, z, l_fact;

  if (pst.dnaseq == SEQT_DNA || pst.dnaseq == SEQT_RNA) {
    ave_n1 = 5000.0;
    l_fact = 1.0;
  }
  else {
    ave_n1 = 400.0;
    l_fact = 0.7;
  }

/*  max_hscore = MAX_SSCORE; */
/*  mean expected for pst.n0 * 400 for protein, 5000 for DNA */
/*  we want a number of offsets that is appropriate for the database size so
    far (nstats)
*/

/*
  the calculation below sets a high-score threshold using an
  ungapped lambda, but errs towards the high-score side by using
  E()=0.001 and calculating with 0.70*lambda, which is the correct for
  going from ungapped to -12/-2 gapped lambda with BLOSUM50
*/

#ifndef NORMAL_DIST
  tmp_score = 0.01/((double)nstats*K*(double)pst.n0*ave_n1);
  tmp_score = -log(tmp_score)/(Lambda*l_fact);
  max_hscore = (int)(tmp_score+0.5);

  z = 1.0/(double)nstats;
  z = (log(z)+EULER_G)/(- PI_SQRT6);
#else
  max_hscore = 100;
  z = 5.0;
#endif
  *zstrim = 10.0*z+50.0;
  return max_hscore;
}

int
proc_hist_r(struct stat_str *sptr, int nstats,
	    struct pstruct pst, struct hist_str *histp,
	    int do_trim, struct pstat_str *pu)
{
  int i, max_hscore;
  double zs, ztrim;
  char s_string[128];
  struct llen_str llen;
  char *f_string;
  llen.fit_flag=1;
  llen.hist=NULL;

  max_hscore = calc_thresh(pst, nstats, pu->ngLambda,
			   pu->ngK, pu->ngH, &ztrim);

  inithist(&llen,pst,max_hscore);
  f_string = &(histp->stat_info[0]);

  for (i = 0; i<nstats; i++)
    addhist(&llen,sptr[i].score,sptr[i].n1, max_hscore);

  if ((llen.max_score - llen.min_score) < 10) {
    free_hist(&llen);
    llen.fit_flag = 0;
    find_zp = &find_zn;
    return proc_hist_n(sptr, nstats, pst, histp, do_trim, pu);
  }

  fit_llen(&llen, &(pu->r_u.rg)); /* now we have rho, mu, rho2, mu2, mean_var
				 to set the parameters for the histogram */

  if (!llen.fit_flag) {	/* the fit failed, fall back to proc_hist_ml */
    free_hist(&llen);
    find_zp = &find_ze;
    return proc_hist_ml(sptr,nstats, pst, histp, do_trim, pu);
  }

  pu->r_u.rg.n_trimmed= pu->r_u.rg.n1_trimmed = pu->r_u.rg.nb_trimmed = 0;

  if (do_trim) {
    if (llen.fit_flag) {
      for (i = 0; i < nstats; i++) {
	zs = find_zr(sptr[i].score,sptr[i].escore,sptr[i].n1,sptr[i].comp, pu);
	if (zs < 20.0 || zs > ztrim) {
	  pu->r_u.rg.n_trimmed++;
	  prune_hist(&llen,sptr[i].score,sptr[i].n1, max_hscore,
		     &(histp->entries));
	}
      }
    }

  /*  fprintf(stderr,"Z-trimmed %d entries with z > 5.0\n", pu->r_u.rg.n_trimmed); */

    if (llen.fit_flag) fit_llens(&llen, &(pu->r_u.rg));

  /*   fprintf(stderr,"Bin-trimmed %d entries in %d bins\n", pu->r_u.rg.n1_trimmed,pu->r_u.rg.nb_trimmed); */
  }

  free_hist(&llen);

  /* put all the scores in the histogram */

  if (pst.zsflag < 10) s_string[0]='\0';
  else if (pst.zs_win > 0)
    sprintf(s_string,"(shuffled, win: %d)",pst.zs_win);
  else strncpy(s_string,"(shuffled)",sizeof(s_string));

  if (pst.zsflag == REG2_STATS || pst.zsflag == 10+REG2_STATS) 
    sprintf(f_string,"%s Expectation_v fit: rho(ln(x))= %6.4f+/-%6.3g; mu= %6.4f+/-%6.3f;\n rho2=%6.2f; mu2= %6.2f, 0's: %d Z-trim: %d  B-trim: %d in %d/%d",
	    s_string, pu->r_u.rg.rho*LN_FACT,sqrt(pu->r_u.rg.rho_e),pu->r_u.rg.mu,sqrt(pu->r_u.rg.mu_e),
	    pu->r_u.rg.rho2,pu->r_u.rg.mu2,llen.zero_s,
	    pu->r_u.rg.n_trimmed, pu->r_u.rg.n1_trimmed, pu->r_u.rg.nb_trimmed, pu->r_u.rg.nb_tot);
  else 
    sprintf(f_string,"%s Expectation_n fit: rho(ln(x))= %6.4f+/-%6.3g; mu= %6.4f+/-%6.3f\n mean_var=%6.4f+/-%6.3f, 0's: %d Z-trim: %d  B-trim: %d in %d/%d\n Lambda= %8.6f",
	    s_string,
	    pu->r_u.rg.rho*LN_FACT,sqrt(pu->r_u.rg.rho_e),pu->r_u.rg.mu,sqrt(pu->r_u.rg.mu_e), pu->r_u.rg.mean_var,sqrt(pu->r_u.rg.var_e),
	    llen.zero_s, pu->r_u.rg.n_trimmed, pu->r_u.rg.n1_trimmed, pu->r_u.rg.nb_trimmed, pu->r_u.rg.nb_tot,
	    PI_SQRT6/sqrt(pu->r_u.rg.mean_var));
  return REG_STATS;
}


int
proc_hist_r2(struct stat_str *sptr, int nstats,
	     struct pstruct pst, struct hist_str *histp,
	     int do_trim, struct pstat_str *pu)
{
  int i, nit, nprune, max_hscore;
  double zs, ztrim;
  char s_string[128];
  char *f_string;
  struct llen_str llen;

  llen.fit_flag=1;
  llen.hist=NULL;

  max_hscore = calc_thresh(pst, nstats, pu->ngLambda,
			   pu->ngK, pu->ngH, &ztrim);

  inithist(&llen, pst,max_hscore);
  f_string = &(histp->stat_info[0]);

  for (i = 0; i<nstats; i++)
    addhist(&llen,sptr[i].score,sptr[i].n1,max_hscore);

  pu->r_u.rg.n_trimmed= pu->r_u.rg.n1_trimmed = pu->r_u.rg.nb_trimmed = 0;
  if (do_trim) nit = 5;
  else nit = 0;

  while (nit-- > 0) {
    nprune = 0;
    fit_llen2(&llen, &(pu->r_u.rg));

    for (i = 0; i < nstats; i++) {
      if (sptr[i].n1 < 0) continue;
      zs = find_zr(sptr[i].score,sptr[i].escore,sptr[i].n1,sptr[i].comp,pu);
      if (zs < 20.0 || zs > ztrim ) {
	nprune++;
	pu->r_u.rg.n_trimmed++;
	prune_hist(&llen,sptr[i].score,sptr[i].n1,max_hscore,
		   &(histp->entries));
	sptr[i].n1 = -sptr[i].n1;
      }
    }
    /*    fprintf(stderr," %d Z-trimmed at %d\n",nprune,nit); */
    if (nprune < LHISTC) { break; }
  }

  fit_llen(&llen, &(pu->r_u.rg));

  free_hist(&llen);

  if (pst.zsflag < 10) s_string[0]='\0';
  else if (pst.zs_win > 0)
    sprintf(s_string,"(shuffled, win: %d)",pst.zs_win);
  else strncpy(s_string,"(shuffled)",sizeof(s_string));

  sprintf(f_string,"%s Expectation_i fit: rho(ln(x))= %6.4f+/-%6.3g; mu= %6.4f+/-%6.3f;\n mean_var=%6.4f+/-%6.3f 0's: %d Z-trim: %d N-it: %d\n Lambda= %8.6f",
	  s_string,
	  pu->r_u.rg.rho*LN_FACT,sqrt(pu->r_u.rg.rho_e),pu->r_u.rg.mu,sqrt(pu->r_u.rg.mu_e),
	  pu->r_u.rg.mean_var,sqrt(pu->r_u.rg.var_e),llen.zero_s,pu->r_u.rg.n_trimmed, nit,
	  PI_SQRT6/sqrt(pu->r_u.rg.mean_var));
  return REGI_STATS;
}

/* this procedure implements Altschul's pre-calculated values for lambda, K */

#include "alt_parms.h"

int
look_p(struct alt_p parm[], int gap, int ext, 
       double *K, double *Lambda, double *H);

int
proc_hist_a(struct stat_str *sptr, int nstats, 
	    struct pstruct pst, struct hist_str *histp,
	    int do_trim, struct pstat_str *pu)
{
  double Lambda, K, H;
  char *f_string;
  int r_v;
  int t_gdelval, t_ggapval;

#ifdef OLD_FASTA_GAP
  t_gdelval = pst.gdelval;
  t_ggapval = pst.ggapval;
#else
  t_gdelval = pst.gdelval+pst.ggapval;
  t_ggapval = pst.ggapval;
#endif

  f_string = &(histp->stat_info[0]);

  if (strcmp(pst.pamfile,"BL50")==0 || strcmp(pst.pamfile,"BLOSUM50")==0)
      r_v = look_p(bl50_p,t_gdelval,t_ggapval,&K,&Lambda,&H);
  else if (strcmp(pst.pamfile,"BL62")==0 || strcmp(pst.pamfile,"BLOSUM62")==0)
      r_v = look_p(bl62_p,t_gdelval,t_ggapval,&K,&Lambda,&H);
  else if (strcmp(pst.pamfile,"BL80")==0 || strcmp(pst.pamfile,"BLOSUM80")==0)
      r_v = look_p(bl80_p,t_gdelval,t_ggapval,&K,&Lambda,&H);
  else if (strcmp(pst.pamfile,"P250")==0)
      r_v = look_p(p250_p,t_gdelval,t_ggapval,&K,&Lambda,&H);
  else if (strcmp(pst.pamfile,"P120")==0)
      r_v = look_p(p120_p,t_gdelval,t_ggapval,&K,&Lambda,&H);
  else if (strcmp(pst.pamfile,"MD_10")==0)
      r_v = look_p(md10_p,t_gdelval,t_ggapval,&K,&Lambda,&H);
  else if (strcmp(pst.pamfile,"MD_20")==0)
      r_v = look_p(md20_p,t_gdelval,t_ggapval,&K,&Lambda,&H);
  else if (strcmp(pst.pamfile,"MD_40")==0)
      r_v = look_p(md40_p,t_gdelval,t_ggapval,&K,&Lambda,&H);
  else if (strcmp(pst.pamfile,"DNA")==0 || strcmp(pst.pamfile,"+5/-4")==0)
      r_v = look_p(nt54_p,t_gdelval,t_ggapval,&K,&Lambda,&H);
  else if (strcmp(pst.pamfile,"+3/-2")==0)
      r_v = look_p(nt32_p,t_gdelval,t_ggapval,&K,&Lambda,&H);
  else if (strcmp(pst.pamfile,"+1/-3")==0)
      r_v = look_p(nt13_p,t_gdelval,t_ggapval,&K,&Lambda,&H);
  else r_v = 0;

  pu->r_u.ag.Lambda = Lambda;
  pu->r_u.ag.K = K;
  pu->r_u.ag.H = H;

  if (r_v == 0) {
    fprintf(stderr,"Parameters not available for: %s: %d/%d\n",
	    pst.pamfile,t_gdelval,t_ggapval);

    find_zp = &find_zr;
    return proc_hist_r(sptr, nstats,pst, histp, do_trim, pu);
  }

  /*
    fprintf(stderr," the parameters are: Lambda: %5.3f K: %5.3f H: %5.3f\n",
	    Lambda, K, H);
	    */

    pu->r_u.ag.a_n0 = (double)pst.n0;
    pu->r_u.ag.a_n0f = log (K * pu->r_u.ag.a_n0)/H;

    sprintf(f_string,"Altschul/Gish params: n0: %d Lambda: %5.3f K: %5.3f H: %5.3f",
	    pst.n0,Lambda, K, H);
    return AG_STATS;
}

int 
ag_parm(char *pamfile, int gdelval, int ggapval, struct pstat_str *pu)
{
  double Lambda, K, H;
  int r_v;

  if (strcmp(pamfile,"BL50")==0)
    r_v = look_p(bl50_p,gdelval,ggapval,&K,&Lambda,&H);
  else if (strcmp(pamfile,"BL62")==0)
      r_v = look_p(bl62_p,gdelval,ggapval,&K,&Lambda,&H);
  else if (strcmp(pamfile,"P250")==0)
      r_v = look_p(p250_p,gdelval,ggapval,&K,&Lambda,&H);
  else if (strcmp(pamfile,"P120")==0)
      r_v = look_p(p120_p,gdelval,ggapval,&K,&Lambda,&H);
  else if (strcmp(pamfile,"MD_10")==0)
      r_v = look_p(md10_p,gdelval,ggapval,&K,&Lambda,&H);
  else if (strcmp(pamfile,"MD_20")==0)
      r_v = look_p(md20_p,gdelval,ggapval,&K,&Lambda,&H);
  else if (strcmp(pamfile,"MD_40")==0)
      r_v = look_p(md40_p,gdelval,ggapval,&K,&Lambda,&H);
  else if (strcmp(pamfile,"DNA")==0 || strcmp(pamfile,"+5/-4")==0)
      r_v = look_p(nt54_p,gdelval,ggapval, &K,&Lambda,&H);
  else if (strcmp(pamfile,"+3/-2")==0)
      r_v = look_p(nt32_p,gdelval,ggapval, &K,&Lambda,&H);
  else if (strcmp(pamfile,"+1/-3")==0)
      r_v = look_p(nt13_p,gdelval,ggapval, &K,&Lambda,&H);
  else r_v = 0;

  pu->r_u.ag.K = K;
  pu->r_u.ag.Lambda = Lambda;
  pu->r_u.ag.H = H;

  if (r_v == 0) {
    fprintf(stderr,"Parameters not available for: %s: %d/%d\n",
	    pamfile,gdelval,ggapval);
    }
  return r_v;
}

int
look_p(struct alt_p parm[], int gap, int ext,
       double *K, double *Lambda, double *H)
{
  int i;

  gap = -gap;
  ext = -ext;

  if (gap > parm[1].gap) {
    *K = parm[0].K;
    *Lambda = parm[0].Lambda;
    *H = parm[0].H;
    return 1;
  }

  for (i=1; parm[i].gap > 0; i++) {
    if (parm[i].gap > gap) continue;
    else if (parm[i].gap == gap && parm[i].ext > ext ) continue;
    else if (parm[i].gap == gap && parm[i].ext == ext) {
      *K = parm[i].K;
      *Lambda = parm[i].Lambda;
      *H = parm[i].H;
      return 1;
    }
    else break;
  }
  return 0;
}

/* uncensored and censored maximum likelihood estimates developed
   by Aaron Mackey based on a preprint from Sean Eddy */

int mle_cen  (struct stat_str *, int, int, double, double *, double *);

int
proc_hist_ml(struct stat_str *sptr, int nstats, 
	     struct pstruct pst, struct hist_str *histp,
	     int do_trim, struct pstat_str *pu)
{
  double f_cen;
  char s_string[128];
  char *f_string;

  f_string = &(histp->stat_info[0]);
  pu->r_u.ag.a_n0 = (double)pst.n0;

  if (pst.zsflag < 10) s_string[0]='\0';
  else if (pst.zs_win > 0)
    sprintf(s_string,"(shuffled, win: %d)",pst.zs_win);
  else strncpy(s_string,"(shuffled)",sizeof(s_string));

  if (!do_trim) {
    if (mle_cen(sptr, nstats, pst.n0, 0.0, &pu->r_u.ag.Lambda, &pu->r_u.ag.K) == -1)
      goto bad_mle;
    sprintf(f_string,"%s MLE statistics: Lambda= %6.4f;  K=%6.4g",
	    s_string,pu->r_u.ag.Lambda,pu->r_u.ag.K);
  }
  else {
    if (nstats/20 > 1000) f_cen = 1000.0/(double)nstats;
    else f_cen = 0.05;
    if (mle_cen(sptr, nstats, pst.n0, f_cen, &pu->r_u.ag.Lambda, &pu->r_u.ag.K) == -1)
      goto bad_mle;
    sprintf(f_string,"MLE_cen statistics: Lambda= %6.4f;  K=%6.4g (cen=%d)",
	    pu->r_u.ag.Lambda,pu->r_u.ag.K,(int)((double)nstats*f_cen));
  }    

  return MLE_STATS;
 bad_mle:
  find_zp = &find_zn;
  
  return proc_hist_n(sptr, nstats, pst, histp, do_trim, pu);
}

int
mle_cen2  (struct stat_str *, int, int, double, double *, double *, double *, double *);


int
proc_hist_ml2(struct stat_str *sptr, int nstats, 
	      struct pstruct pst, struct hist_str *histp,
	      int do_trim, struct pstat_str *pu)
{
  int i, ns=0, nneg=0;
  double f_cen, ave_lambda;
  char s_string[128], ex_string[64];
  char *f_string;

  f_string = &(histp->stat_info[0]);
  pu->r_u.m2.a_n0 = (double)pst.n0;

  if (pst.zsflag < 10) s_string[0]='\0';
  else if (pst.zs_win > 0)
    sprintf(s_string,"(shuffled, win: %d)",pst.zs_win);
  else strncpy(s_string,"(shuffled)",sizeof(s_string));

  pu->r_u.m2.ave_comp = 0.0;
  pu->r_u.m2.max_comp = -1.0;

  ns = nneg = 0;
  for (i=0; i<nstats; i++) {
    if (sptr[i].comp > pu->r_u.m2.max_comp) pu->r_u.m2.max_comp = sptr[i].comp;
    if (sptr[i].comp > 0.0) {
      pu->r_u.m2.ave_comp += log(sptr[i].comp);
      ns++;
    }
    else nneg++;
  }
  pu->r_u.m2.ave_comp /= (double)ns;
  pu->r_u.m2.ave_comp = exp(pu->r_u.m2.ave_comp);
  for (i=0; i<nstats; i++) if (sptr[i].comp < 0.0) {
    sptr[i].comp = pu->r_u.m2.ave_comp;
  }

  if (nneg > 0)
    sprintf(ex_string,"composition = -1 for %d sequences",nneg);
  else ex_string[0]='\0';

  if (!do_trim) {
    if (mle_cen2(sptr, nstats, pst.n0, 0.0,
	     &pu->r_u.m2.mle2_a0, &pu->r_u.m2.mle2_a1,
	     &pu->r_u.m2.mle2_a2, &pu->r_u.m2.mle2_b1) == -1) goto bad_mle2;
    ave_lambda = 1.0/(pu->r_u.m2.ave_comp*pu->r_u.m2.mle2_b1);

    sprintf(f_string,"%s MLE-2 statistics: a0= %6.4f;  a1=%6.4f; a2=%6.4f; b1=%6.4f\n  ave Lamdba: %6.4f",
	    s_string, pu->r_u.m2.mle2_a0, pu->r_u.m2.mle2_a1, pu->r_u.m2.mle2_a2, pu->r_u.m2.mle2_b1,ave_lambda);
  }
  else {
    if (nstats/20 > 500) f_cen = 500.0/(double)nstats;
    else f_cen = 0.05;
    if (mle_cen2(sptr, nstats, pst.n0, f_cen, &pu->r_u.m2.mle2_a0, &pu->r_u.m2.mle2_a1, &pu->r_u.m2.mle2_a2, &pu->r_u.m2.mle2_b1)== -1) goto bad_mle2;

    ave_lambda = 1.0/(pu->r_u.m2.ave_comp*pu->r_u.m2.mle2_b1);

    sprintf(f_string,"%s MLE-2-cen statistics: a0= %6.4f;  a1=%6.4f; a2=%6.4f; b1=%6.4f (cen=%d)\n  ave Lambda:%6.4f",
	    s_string, pu->r_u.m2.mle2_a0, pu->r_u.m2.mle2_a1, pu->r_u.m2.mle2_a2, pu->r_u.m2.mle2_b1, (int)((double)nstats*f_cen),ave_lambda);
  }    

  return MLE2_STATS;
 bad_mle2:
  find_zp = &find_zn;
  return proc_hist_n(sptr, nstats, pst, histp, do_trim, pu);
}

double first_deriv_cen(double lambda, struct stat_str *sptr, 
		       int start, int stop,
		       double sumlenL, double cenL,
		       double sumlenH, double cenH);

double second_deriv_cen(double lambda, struct stat_str *sptr,
			int start, int stop,
			double sumlenL, double cenL,
			double sumlenH, double cenH);

static void
st_sort (struct stat_str *v, int n) {
   int gap, i, j;
   int tmp;

   for (gap = 1; gap < n/3; gap = 3*gap +1) ;

   for (; gap > 0; gap = (gap-1)/3)
      for (i = gap; i < n; i++)
	 for (j = i - gap; j >= 0; j -= gap) {
	   if (v[j].score <= v[j + gap].score) break;

	   tmp = v[j].score;
	   v[j].score = v[j + gap].score;
	   v[j + gap].score = tmp;

	   tmp = v[j].n1;
	   v[j].n1 = v[j + gap].n1;
	   v[j + gap].n1 = tmp;
	 }
}

/* sptr[].score, sptr[].n1; sptr[] must be sorted
   int n = total number of samples
   int M = length of query
   double fn = fraction of scores to be censored fn/2.0 from top, bottom
   double *Lambda = Lambda estimate
   double *K = K estimate
*/
int
mle_cen(struct stat_str *sptr, int n, int M, double fc,
	double *Lambda, double *K) {

  double sumlenL, sumlenH, cenL, cenH;
  double sum_s, sum2_s, mean_s, var_s, dtmp;
  int start, stop;
  int i, nf;
  double deriv, deriv2, lambda, old_lambda, sum = 0.0;
  /*
   int sumlenL, int sumlenghtsR = sum of low (Left), right (High) seqs.
   int cenL, cenH = censoring score low, high 
  */

  nf = (fc/2.0) * n;
  start = nf;
  stop = n - nf;

  st_sort(sptr,n);

  sum_s = sum2_s = 0.0;
  for (i=start; i<stop; i++) {
    sum_s += sptr[i].score;
  }
  dtmp = (double)(stop-start);
  mean_s = sum_s/dtmp;

  for (i=start; i<stop; i++) {
    sum2_s += sptr[i].score * sptr[i].score;
  }
  var_s = sum2_s/(dtmp-1.0);

  sumlenL = sumlenH = 0.0;
  for (i=0; i<start; i++) sumlenL += (double)sptr[i].n1;
  for (i=stop; i<n; i++) sumlenH += (double)sptr[i].n1;

  if (nf > 0) {
    cenL = (double)sptr[start].score;
    cenH = (double)sptr[stop].score;
  }
  else {
    cenL = (double)sptr[start].score/2.0;
    cenH = (double)sptr[start].score*2.0;
  }

  if (cenL >= cenH) return -1;

  /* initial guess for lambda is 0.2 - this does not work for matrices
     with very different scales */
  /*  lambda = 0.2; */
  lambda = PI_SQRT6/sqrt(var_s);
  if (lambda > 1.0) {
    fprintf(stderr," Lambda initial estimate error: lambda: %6.4g; var_s: %6.4g\n",lambda,var_s);
    lambda = 0.2;
  }

  do {
    deriv = first_deriv_cen(lambda, sptr, start, stop,
			    sumlenL, cenL, sumlenH, cenH);
    /*   (uncensored version)
	 first_deriv(lambda, &sptr[start], stop - start))
    */

    /*  (uncensored version)
    deriv2 = second_deriv(lambda, &sptr[start], stop-start);
    */
    deriv2 = second_deriv_cen(lambda, sptr, start, stop,
			     sumlenL, cenL, sumlenH, cenH); 

    old_lambda = lambda;
    if (lambda - deriv/deriv2 > 0.0) lambda = lambda - deriv/deriv2;
    else lambda = lambda/2.0;
  } while (fabs((lambda - old_lambda)/lambda) > TINY);

  
  for(i = start; i < stop ; i++) {
    sum += (double) sptr[i].n1 * exp(- lambda * (double)sptr[i].score);
  }

  *Lambda = lambda;
  /* 
  *K = (double)(stop-start)/((double)M*sum);
  */
  *K = (double)n/((double)M*
		  (sum+sumlenL*exp(-lambda*cenL)-sumlenH*exp(-lambda*cenH)));
  return 0;
}

/*
double
first_deriv(double lambda, struct stat_str *sptr, int n) {

  int i;
  double sum = 0.0, sum1 = 0.0, sum2 = 0.0;
  double s, l, es;

  for(i = 0 ; i < n ; i++) {
    s = (double)sptr[i].score;
    l = (double)sptr[i].n1;
    es = exp(-lambda * s );
    sum += s;
    sum2 += l * es;
    sum1 += s * l * es;
  }

  return (1.0/lambda) - (sum/(double)n) + (sum1/sum2);
}
*/

/*
double
second_deriv(double lambda, struct stat_str *sptr, int n) {
  double sum1 = 0.0, sum2 = 0.0, sum3 = 0.0;
  double s, l, es;
  int i;

  for(i = 0 ; i < n ; i++) {
    l = (double)sptr[i].n1;
    s = (double)sptr[i].score;
    es = exp(-lambda * s);
    sum2 += l * es;
    sum1 += l * s * es;
    sum3 += l * s * s * es;
  }

  return ((sum1*sum1)/(sum2*sum2)) - (sum3/sum2) -  (1.0/(lambda*lambda));
}
*/

double
first_deriv_cen(double lambda, struct stat_str *sptr, int start, int stop,
		double sumlenL, double cenL, double sumlenH, double cenH) {
  int i;
  double sum = 0.0, sum1 = 0.0, sum2 = 0.0;
  double s, l, es;

  for(i = start ; i < stop ; i++) {
    s = (double)sptr[i].score;
    l = (double)sptr[i].n1;
    es = exp(-lambda * s );
    sum += s;
    sum2 += l * es;
    sum1 += s * l * es;
  }

  sum1 += sumlenL*cenL*exp(-lambda*cenL) - sumlenH*cenH*exp(-lambda*cenH);
  sum2 += sumlenL*exp(-lambda*cenL) - sumlenH*exp(-lambda*cenH);

  return (1.0 / lambda) - (sum /(double)(stop-start)) + (sum1 / sum2);
}

double
second_deriv_cen(double lambda, struct stat_str *sptr, int start, int stop,
		 double sumlenL, double cenL, double sumlenH, double cenH) {

  double sum1 = 0.0, sum2 = 0.0, sum3 = 0.0;
  double s, l, es;
  int i;

  for(i = start ; i < stop ; i++) {
    s = (double)sptr[i].score;
    l = (double)sptr[i].n1;
    es = exp(-lambda * s);
    sum2 += l * es;
    sum1 += l * s * es;
    sum3 += l * s * s * es;
  }

  sum1 += sumlenL*cenL*exp(-lambda*cenL) - sumlenH*cenH*exp(-lambda*cenH);
  sum2 += sumlenL*exp(-lambda * cenL) -  sumlenH*exp(-lambda * cenH);
  sum3 += sumlenL*cenL*cenL * exp(-lambda * cenL) -
    sumlenH*cenH*cenH * exp(-lambda * cenH);
  return ((sum1 * sum1) / (sum2 * sum2)) - (sum3 / sum2)
    - (1.0 / (lambda * lambda));
}

double mle2_func(double *params,
		 double *consts,
		 struct stat_str *values,
		 int n, int start, int stop);

void simplex(double *fitparams,
	     double *lambda,
	     int nparam,
	     double (*minfunc) (double *tryparams, double *consts, 
				struct stat_str *data, int ndata,
				int start, int stop),
	     double *consts,
	     void *data,
	     int ndata, int start, int stop
	     );

int
mle_cen2(struct stat_str *sptr, int n, int M, double fc,
	double *a0, double *a1, double *a2, double *b1) {

  double params[4], lambdas[4], consts[9];
  double avglenL, avglenH, avgcompL, avgcompH, cenL, cenH;
  int start, stop;
  int i, nf;

  nf = (fc/2.0) * n;
  start = nf;
  stop = n - nf;

  st_sort(sptr,n);

  /* choose arithmetic or geometic mean for compositions by appropriate commenting */

  if (nf > 0) {
    avglenL = avglenH = 0.0;
    avgcompL = avgcompH = 0.0;
    /* avgcompL = avgcompH = 1.0 */
    for (i=0; i<start; i++) {
      avglenL += (double)sptr[i].n1;
      avgcompL += (double)sptr[i].comp;
      /* avgcompL *= (double) sptr[i].comp; */
    }
    avglenL /= (double) start;
    avgcompL /= (double) start;
    /* avgcompL = pow(avgcompL, 1.0/(double) start); */
  
    for (i=stop; i<n; i++) {
      avglenH += (double)sptr[i].n1;
      avgcompH += (double)sptr[i].comp;
      /* avgcompH *= (double) sptr[i].comp; */
    }
    avglenH /= (double) (n - stop);
    avgcompH /= (double) (n - stop);
    /* avgcompL = pow(avgcompL, 1.0/(double) (n - stop)); */

    cenL = (double)sptr[start].score;
    cenH = (double)sptr[stop].score;
    if (cenL >= cenH) return -1;
  }
  else {
    avglenL = avglenH = cenL = cenH = 0.0;
    avgcompL = avgcompH = 1.0;
  }

  params[0] = 10.0;
  params[1] = -10.0;
  params[2] = 1.0;
  params[3] = 1.0;

  lambdas[0] = 1.0;
  lambdas[1] = 0.5;
  lambdas[2] = 0.1;
  lambdas[3] = 0.01;

  consts[0] = M;
  consts[1] = (double) start;
  consts[2] = (double) stop;
  consts[3] = cenL;
  consts[4] = cenH;
  consts[5] = avglenL;
  consts[6] = avglenH;
  consts[7] = avgcompL;
  consts[8] = avgcompH;

  simplex(params, lambdas, 4,
	  (double (*) (double *, double *, struct stat_str *, int, int, int) )mle2_func,
	  consts, sptr, n, start, stop);

  *a0 = params[0];
  *a1 = params[1];
  *a2 = params[2];
  *b1 = params[3];

  return 0;
}

double mle2_func(double *params,
		 double *consts,
		 struct stat_str *values,
		 int n, int start, int stop
		 ) {

  double a0, a1, a2, b1, M;
  double score, length, comp;
  double cenL, cenH, avglenL, avglenH, avgcompL, avgcompH;
  double L, y;

  int i;

  a0 = params[0];
  a1 = params[1];
  a2 = params[2];
  b1 = params[3];

  M = consts[0];
  /*
  start = (int) consts[1];
  stop = (int) consts[2];
  */
  cenL = consts[3];
  cenH = consts[4];
  avglenL = consts[5];
  avglenH = consts[6];
  avgcompL = consts[7];
  avgcompH = consts[8];

  L = 0;
  y = 0;

  if (start > 0) {
    y = -(cenL - (a0 + a1*avgcompL +a2*avgcompL*log(M*avglenL)))/(b1*avgcompL);
    L += (double) start * exp(y);
  }

  for(i = start ; i < stop ; i++) {
    score = (double) values[i].score;
    length = (double) values[i].n1;
    comp = (double) values[i].comp;

    y = - (score - (a0 + a1*comp + a2 * comp * log(M*length))) / (b1*comp);

    L += -y + exp(y) + log(b1 * comp);
  }

  if (stop < n) {
    y = -(cenH -(a0 + a1*avgcompH + a2*avgcompH*log(M*avglenH)))/(b1*avgcompH);
    L -= (double) (n - stop) * exp(y);
  }
  return L;
}

/* Begin Nelder-Mead simplex code: */

double evalfunc(double **param,
		double *vals,
		double *psums,
		double *ptry,
		int nparam,
		double (*minfunc) (double *params, double *consts,
				   struct stat_str *data, int ndata,
				   int start, int stop),
		double *consts,
		void *data,
		int ndata, int start, int stop,
		int ihi,
		double factor);

void simplex(double *fitparams,
	     double *lambda,
	     int nparam,
	     double (*minfunc) (double *tryparams, double *consts,
				struct stat_str *data, int ndata, 
				int start, int stop),
	     double *consts,
	     void *data,
	     int ndata,
	     int start,
	     int stop
	     )
{

  int i, j, ilo, ihi, inhi;
  double rtol, sum, tmp, ysave, ytry;
  double *psum, *vals, *ptry, **param;


  psum = (double *) calloc(nparam, sizeof(double));
  ptry = (double *) calloc(nparam, sizeof(double));

  vals = (double *) calloc(nparam + 1, sizeof(double));

  param = (double **) calloc(nparam + 1, sizeof(double *));
  param[0] = (double *) calloc((nparam + 1) * nparam, sizeof(double));
  for( i = 1 ; i < (nparam + 1) ; i++ ) {
    param[i] = param[0] + i * nparam;
  }

  /* Get our N+1 initial parameter values for the simplex */

  for( i = 0 ; i < nparam ; i++ ) {
    param[0][i] = fitparams[i];
  }

  for( i = 1 ; i < (nparam + 1) ; i++ ) {
    for( j = 0 ; j < nparam ; j++ ) {
      param[i][j] = fitparams[j] + lambda[j] * ( (i - 1) == j ? 1 : 0 );
    }
  }

  /* calculate initial values at the simplex nodes */

  for( i = 0 ; i < (nparam + 1) ; i++ ) {
    vals[i] = minfunc(param[i], consts, data, ndata, start, stop);
  }

  /* Begin Nelder-Mead simplex algorithm from Numerical Recipes in C */

  for( j = 0 ; j < nparam ; j++ ) {
    for( sum = 0.0, i = 0 ; i < nparam + 1 ; i++ ) {
      sum += param[i][j];
    }
    psum[j] = sum;
  }


  while( 1 ) {
/*
      determine which point is highest (ihi), next highest (inhi) and
      lowest (ilo) by looping over the points in the simplex
*/
    ilo = 0;

/*  ihi = vals[0] > vals[1] ? (inhi = 1, 0) : (inhi = 0, 1); */
    if(vals[0] > vals[1]) { ihi = 0; inhi = 1; }
    else { ihi = 1; inhi = 0; }

    for( i = 0 ; i < nparam + 1 ; i++) {
      if( vals[i] <= vals[ilo] ) ilo = i;
      if( vals[i] > vals[ihi] ) {
	inhi = ihi;
	ihi = i;
      } else if ( vals[i] > vals[inhi] && i != ihi ) inhi = i;
    }

    /* Are we finished? */

    rtol = 2.0 * fabs(vals[ihi] - vals[ilo]) / 
      (fabs(vals[ihi]) + fabs(vals[ilo]) + TINY);

    if( rtol < TOLERANCE ) {

/* put the best value and best parameters into the first index */

      tmp = vals[0];
      vals[0] = vals[ilo];
      vals[ilo] = tmp;

      for( i = 0 ; i < nparam ; i++ ) {
	tmp = param[0][i];
	param[0][i] = param[ilo][i];
	param[ilo][i] = tmp;
      }

      /* et voila, c'est finis */
      break;
    }

    /* Begin a new iteration */

    /* first, extrapolate by -1 through the face of the simplex across from ihi */

    ytry = evalfunc(param, vals, psum, ptry, nparam, minfunc, consts,
		    data, ndata, start, stop, ihi, -1.0);

    if( ytry <= vals[ilo] ) {

      /* Good result, try additional extrapolation by 2 */

      ytry = evalfunc(param, vals, psum, ptry, nparam, minfunc, consts, 
		      data, ndata, start, stop, ihi, 2.0);

    } else if ( ytry >= vals[inhi] ) {

      /* no good, look for an intermediate lower point by contracting */

      ysave = vals[ihi];
      ytry = evalfunc(param, vals, psum, ptry, nparam, minfunc, consts,
		      data, ndata, start, stop, ihi, 0.5);

      if( ytry >= ysave ) {

	/* Still no good.  Contract around lowest (best) point. */

	for( i = 0 ; i < nparam + 1 ; i++ ) {
	  if( i != ilo ) {
	    for ( j = 0 ; j < nparam ; j++ ) {
	      param[i][j] = psum[j] = 0.5 * (param[i][j] + param[ilo][j]);
	    }
	    vals[i] = minfunc(psum, consts, data, ndata, start, stop);
	  }
	}


	for( j = 0 ; j < nparam ; j++ ) {
	  for( sum = 0.0, i = 0 ; i < nparam + 1 ; i++ ) {
	    sum += param[i][j];
	  }
	  psum[j] = sum;
	}

      }
    }
  }
			   
  for( i = 0 ; i < nparam ; i++ ) {
    fitparams[i] = param[0][i];
  }

  if (ptry!=NULL) {
    free(ptry);
    ptry=NULL;
  }
  free(param[0]);
  free(param);
  free(vals);
  free(psum);
}


double evalfunc(double **param,
		double *vals,
		double *psum,
		double *ptry,
		int nparam,
		double (*minfunc)(double *tryparam, double *consts,
				  struct stat_str *data, int ndata,
				  int start, int stop),
		double *consts,
		void *data,
		int ndata, int start, int stop,
		int ihi,
		double factor) {

  int j;
  double fac1, fac2, ytry;


  fac1 = (1.0 - factor) / nparam;
  fac2 = fac1 - factor;

  for( j = 0 ; j < nparam ; j++ ) {
    ptry[j] = psum[j] * fac1 - param[ihi][j] * fac2;
  }

  ytry = minfunc(ptry, consts, data, ndata, start, stop);

  if( ytry < vals[ihi] ) {
    vals[ihi] = ytry;
    for( j = 0 ; j < nparam ; j++ ) {
      psum[j] += ptry[j] - param[ihi][j];
      param[ihi][j] = ptry[j];
    }
  }

  return ytry;
}

/* end of Nelder-Mead simplex code */

int
proc_hist_n(struct stat_str *sptr, int nstats,
	    struct pstruct pst, struct hist_str *histp,
	    int do_trim, struct pstat_str *pu)
{
  int i, j;
  double s_score, s2_score, ssd, ztrim;
  int nit, max_hscore;
  char s_string[128];
  char *f_string;

  f_string = &(histp->stat_info[0]);

  max_hscore = calc_thresh(pst, nstats, pu->ngLambda,
			   pu->ngK, pu->ngH, &ztrim);

  s_score = s2_score = 0.0;

  for ( j = 0, i = 0; i < nstats; i++) {
    if (sptr[i].score > 0 && sptr[i].score <= max_hscore) {
      s_score += (ssd=(double)sptr[i].score);
      s2_score += ssd * ssd;
      j++;
    }
  }

  if (j > 1 ) {
    pu->r_u.rg.mu = s_score/(double)j;
    pu->r_u.rg.mean_var = s2_score - (double)j * pu->r_u.rg.mu * pu->r_u.rg.mu;
    pu->r_u.rg.mean_var /= (double)(j-1);
  }
  else {
    pu->r_u.rg.mu = 50.0;
    pu->r_u.rg.mean_var = 10.0;
  }
  
  if (pu->r_u.rg.mean_var < 0.01) {
    pu->r_u.rg.mean_var = (pu->r_u.rg.mu > 1.0) ? pu->r_u.rg.mu: 1.0;
  }

  /* now remove some scores */

  nit = 5;
  while (nit-- > 0) {
    pu->r_u.rg.n_trimmed = 0;

    for (i=0; i< nstats; i++) {
      if (sptr[i].n1 < 0) continue;
      ssd = find_zn(sptr[i].score,sptr[i].escore,sptr[i].n1,sptr[i].comp, pu);
      if (ssd > ztrim || ssd < 20.0) {
	/*      fprintf(stderr,"removing %3d %3d %4.1f\n",
		sptr[i].score, sptr[i].n1,ssd); */
	ssd = sptr[i].score;
	s_score -= ssd;
	s2_score -= ssd*ssd;
	j--;
	pu->r_u.rg.n_trimmed++;
	histp->entries--;
	sptr[i].n1 = -sptr[i].n1;
      }
    }

    if (j > 1 ) {
      pu->r_u.rg.mu = s_score/(double)j;
      pu->r_u.rg.mean_var = s2_score - (double)j * pu->r_u.rg.mu * pu->r_u.rg.mu;
      pu->r_u.rg.mean_var /= (double)(j-1);
    }
    else {
      pu->r_u.rg.mu = 50.0;
      pu->r_u.rg.mean_var = 10.0;
    }

    if (pu->r_u.rg.mean_var < 0.01) {
      pu->r_u.rg.mean_var = (pu->r_u.rg.mu > 1.0) ? pu->r_u.rg.mu: 1.0;
    }

    if (pu->r_u.rg.n_trimmed < LHISTC) {
      /*
	fprintf(stderr,"nprune %d at %d\n",nprune,nit);
	*/
      break;
    }
  }

  if (pst.zsflag < 10) s_string[0]='\0';
  else if (pst.zs_win > 0)
    sprintf(s_string,"(shuffled, win: %d)",pst.zs_win);
  else strncpy(s_string,"(shuffled)",sizeof(s_string));

  sprintf(f_string,"%s unscaled statistics: mu= %6.4f  var=%6.4f; Lambda= %6.4f",
	  s_string, pu->r_u.rg.mu,pu->r_u.rg.mean_var,PI_SQRT6/sqrt(pu->r_u.rg.mean_var));
  return AVE_STATS;
}

/*
This routine calculates the maximum likelihood estimates for the
extreme value distribution exp(-exp(-(-x-a)/b)) using the formula

	<lambda> = x_m - sum{ x[i] * exp (-x[i]<lambda>)}/sum{exp (-x[i]<lambda>)}
	<a> = -<1/lambda> log ( (1/nlib) sum { exp(-x[i]/<lambda> } )

	The <a> parameter can be transformed into and K
	of the formula: 1 - exp ( - K m n exp ( - lambda S ))
	using the transformation: 1 - exp ( -exp -(lambda S + log(K m n) ))
			1 - exp ( -exp( - lambda ( S + log(K m n) / lambda))

			a = log(K m n) / lambda
			a lambda = log (K m n)
			exp(a lambda)  = K m n 
	 but from above: a lambda = log (1/nlib sum{exp( -x[i]*lambda)})
	 so:            K m n = (1/n sum{ exp( -x[i] *lambda)})
			K = sum{}/(nlib m n )

*/

void
alloc_hist(struct llen_str *llen)
{
  int max_llen, i;
  max_llen = llen->max;

  if (llen->hist == NULL) {
    llen->hist = (int *)calloc((size_t)(max_llen+1),sizeof(int));
    llen->score_sums = (double *)calloc((size_t)(max_llen + 1),sizeof(double));
    llen->score2_sums =(double *)calloc((size_t)(max_llen + 1),sizeof(double));
    llen->score_var = (double *)calloc((size_t)(max_llen + 1),sizeof(double));
  }

  for (i=0; i< max_llen+1; i++) {
      llen->hist[i] = 0;
      llen->score_var[i] = llen->score_sums[i] = llen->score2_sums[i] = 0.0;
  }
}
  
void
free_hist(struct llen_str *llen)
{
  if (llen->hist!=NULL) {
    free(llen->score_var);
    free(llen->score2_sums);
    free(llen->score_sums);
    free(llen->hist);
    llen->hist=NULL;
  }
}

void
inithist(struct llen_str *llen, struct pstruct pst, int max_hscore)
{
  llen->max = MAX_LLEN;

  llen->max_score = -1;
  llen->min_score=10000;

  alloc_hist(llen);

  llen->zero_s = 0;
  llen->min_length = 10000;
  llen->max_length = 0;
}

void
addhist(struct llen_str *llen, int score, int length, int max_hscore)
{
  int llength; 
  double dscore;

  if ( score<=0 || length < LENGTH_CUTOFF) {
    llen->min_score = 0;
    llen->zero_s++;
    return;
  }

  if (score < llen->min_score) llen->min_score = score;
  if (score > llen->max_score) llen->max_score = score;

  if (length > llen->max_length) llen->max_length = length;
  if (length < llen->min_length) llen->min_length = length;
  if (score > max_hscore) score = max_hscore;

  llength = (int)(LN_FACT*log((double)length)+0.5);

  if (llength < 0 ) llength = 0;
  if (llength > llen->max) llength = llen->max;
  llen->hist[llength]++;
  dscore = (double)score;
  llen->score_sums[llength] += dscore;
  llen->score2_sums[llength] += dscore * dscore;
}

/* histogram will go from z-scores of 20 .. 100 with mean 50 and z=10 */

void
inithistz(int mh, struct hist_str *histp )
{
  int i;

  histp->z_calls = 0;

  histp->min_hist = 20;
  histp->max_hist = 120;

  histp->histint = (int)
    ((double)(histp->max_hist - histp->min_hist + 2)/(double)mh+0.5);
  histp->maxh = (int)
    ((double)(histp->max_hist - histp->min_hist + 2)/(double)histp->histint+0.5);

  if (histp->hist_a==NULL) {
    if ((histp->hist_a=(int *)calloc((size_t)histp->maxh,sizeof(int)))==
	NULL) {
      fprintf(stderr," cannot allocate %d for histogram\n",histp->maxh);
      histp->histflg = 0;
    }
    else histp->histflg = 1;
  }
  else {
    for (i=0; i<histp->maxh; i++) histp->hist_a[i]=0;
  }
  histp->entries = 0;
}

static double nrv[100]={
  0.3098900570,-0.0313400923, 0.1131975903,-0.2832547606, 0.0073672659,
  0.2914489107, 0.4209306311,-0.4630181404, 0.3326537896, 0.0050140359,
 -0.1117435426,-0.2835630301, 0.2302997065,-0.3102716394, 0.0819894916,
 -0.1676455701,-0.3782225018,-0.3204509938,-0.3594969187,-0.0308950398,
  0.2922813812, 0.1337170751, 0.4666577031,-0.2917784349,-0.2438179916,
  0.3002301394, 0.0231147123, 0.5687927366,-0.2318208709,-0.1476839273,
 -0.0385043851,-0.1213476523, 0.1486341995, 0.1027917167, 0.1409192644,
 -0.3280652579, 0.4232041455, 0.0775993309, 0.1159071787, 0.2769424442,
  0.3197284751, 0.1507346903, 0.0028580909, 0.4825103412,-0.0496843610,
 -0.2754357656, 0.6021881753,-0.0816123956,-0.0899148991, 0.4847183201,
  0.2151621865,-0.4542246220, 0.0690709102, 0.2461894193, 0.2126042295,
 -0.0767060668, 0.4819746149, 0.3323031326, 0.0177600676, 0.1143185210,
  0.2653977455, 0.0921872958,-0.1330986718, 0.0412287716,-0.1691604748,
 -0.0529679078,-0.0194157955,-0.6117493924, 0.1199067932, 0.0210243193,
 -0.5832259838,-0.1685528664, 0.0008591271,-0.1120347822, 0.0839125069,
 -0.2787486831,-0.1937017962,-0.1915733940,-0.7888453635,-0.3316745163,
  0.1180885226,-0.3347001067,-0.2477492636,-0.2445697600, 0.0001342482,
 -0.0015759812,-0.1516473992,-0.5202267615, 0.2136975210, 0.2500423188,
 -0.2402926401,-0.1094186280,-0.0618869933,-0.0815221188, 0.2623337275,
  0.0219427302 -0.1774469919, 0.0828245026,-0.3271952808,-0.0632898028};

void
addhistz(double zs, struct hist_str *histp)
{
  int ih, zi;
  double rv;

   rv = nrv[histp->z_calls++ % 100];
  zi = (int)(zs + 0.5+rv );

  if ((zi >= 0) && (zi <= 120)) histp->entries++;

  if (zi < histp->min_hist) zi = histp->min_hist;
  if (zi > histp->max_hist) zi = histp->max_hist;
  
  ih = (zi - histp->min_hist)/histp->histint;

  histp->hist_a[ih]++;
}

/* addhistzp() does not increase histp->entries since addhist did it already */
/*
void
addhistzp(double zs, struct hist_str *histp)
{
  int ih, zi;
  double rv;

  rv = nrv[histp->z_calls++ %100];
  zi = (int)(zs + 0.5 + rv);

  if (zi < histp->min_hist) zi = histp->min_hist;
  if (zi > histp->max_hist) zi = histp->max_hist;
  
  ih = (zi - histp->min_hist)/histp->histint;

  histp->hist_a[ih]++;
}
*/

void
prune_hist(struct llen_str *llen, int score, int length, int max_hscore,
	   long *entries)
{
  int llength;
  double dscore;

  if (score <= 0 || length < LENGTH_CUTOFF) return;

  if (score > max_hscore) score = max_hscore;

  llength = (int)(LN_FACT*log((double)length)+0.5);

  if (llength < 0 ) llength = 0;
  if (llength > llen->max) llength = llen->max;
  llen->hist[llength]--;
  dscore = (double)score;
  llen->score_sums[llength] -= dscore;
  llen->score2_sums[llength] -= dscore * dscore;

/*  (*entries)--; histp->entries is not yet initialized */
}  

/* fit_llen: no trimming
   (1) regress scores vs log(n) using weighted variance
   (2) calculate mean variance after length regression
*/

void
fit_llen(struct llen_str *llen, struct rstat_str *pr)
{
  int j;
  int n;
  int n_size;
  double x, y2, u, z;
  double mean_x, mean_y, var_x, var_y, covar_xy;
  double mean_y2, covar_xy2, var_y2, dllj;

  double sum_x, sum_y, sum_x2, sum_xy, sum_v, det, n_w;
  
/* now fit scores to best linear function of log(n), using
   simple linear regression */
  
  for (llen->min=0; llen->min < llen->max; llen->min++)
    if (llen->hist[llen->min]) break;
  llen->min--;

  for (n_size=0,j = llen->min; j < llen->max; j++) {
    if (llen->hist[j] > 1) {
      dllj = (double)llen->hist[j];
      llen->score_var[j] = llen->score2_sums[j]/dllj
	- (llen->score_sums[j]/dllj)*(llen->score_sums[j]/dllj);
      llen->score_var[j] /= (double)(llen->hist[j]-1);
      if (llen->score_var[j] <= 0.1 ) llen->score_var[j] = 0.1;
      n_size++;
    }
  }

  pr->nb_tot = n_size;

  n_w = 0.0;
  sum_x = sum_y = sum_x2 = sum_xy = sum_v = 0;
  for (j = llen->min; j < llen->max; j++)
    if (llen->hist[j] > 1) {
      x = j + 0.5;
      dllj = (double)llen->hist[j];
      n_w += dllj/llen->score_var[j];
      sum_x +=   dllj * x / llen->score_var[j] ;
      sum_y += llen->score_sums[j] / llen->score_var[j];
      sum_x2 +=  dllj * x * x /llen->score_var[j];
      sum_xy +=  x * llen->score_sums[j]/llen->score_var[j];
    }

  if (n_size < 5 ) {
    llen->fit_flag=0;
    pr->rho = 0;
    pr->mu = sum_y/n_w;
    return;
  }
  else {
    det = n_w * sum_x2 - sum_x * sum_x;
    if (det > 0.001) {
      pr->rho = (n_w * sum_xy  - sum_x * sum_y)/det;
      pr->rho_e = n_w/det;
      pr->mu = (sum_x2 * sum_y - sum_x * sum_xy)/det;
      pr->mu_e = sum_x2/det;
    }
    else {
      llen->fit_flag = 0;
      pr->rho = 0;
      pr->mu = sum_y/n_w;
      return;
    }
  }

  det = n_w * sum_x2 - sum_x * sum_x;
  pr->rho = (n_w * sum_xy  - sum_x * sum_y)/det;
  pr->mu = (sum_x2 * sum_y - sum_x * sum_xy)/det;

  n = 0;
  mean_x = mean_y = mean_y2 = 0.0;
  var_x = var_y = 0.0;
  covar_xy = covar_xy2 = 0.0;

  for (j = llen->min; j <= llen->max; j++) 
   if (llen->hist[j] > 1 ) {
      n += llen->hist[j];
      x = (double)j + 0.5;
      mean_x += (double)llen->hist[j] * x;
      mean_y += llen->score_sums[j];
      var_x += (double)llen->hist[j] * x * x;
      var_y += llen->score2_sums[j];
      covar_xy += x * llen->score_sums[j];
    }
  mean_x /= n; mean_y /= n;
  var_x = var_x / n - mean_x * mean_x;
  var_y = var_y / n - mean_y * mean_y;
  
  covar_xy = covar_xy / n - mean_x * mean_y;
/*
  pr->rho = covar_xy / var_x;
  pr->mu = mean_y - pr->rho * mean_x;
*/
  mean_y2 = covar_xy2 = var_y2 = 0.0;
  for (j = llen->min; j <= llen->max; j++) 
    if (llen->hist[j] > 1) {
      x = (double)j + 0.5;
      u = pr->rho * x + pr->mu;
      y2 = llen->score2_sums[j] - 2.0 * llen->score_sums[j] * u + llen->hist[j] * u * u;
/*
      dllj = (double)llen->hist[j];
      fprintf(stderr,"%.2f\t%d\t%g\t%g\n",x/LN_FACT,llen->hist[j],
	      llen->score_sums[j]/dllj,y2/dllj);
*/
      mean_y2 += y2;
      var_y2 += y2 * y2;
      covar_xy2 += x * y2;
      /*      fprintf(stderr,"%6.1f %4d %8d %8d %7.2f %8.2f\n",
	      x,llen->hist[j],llen->score_sums[j],llen->score2_sums[j],u,y2); */
    }
  
  pr->mean_var = mean_y2 /= (double)n;
  covar_xy2 = covar_xy2 / (double)n - mean_x * mean_y2;

  if (pr->mean_var <= 0.01) {
    llen->fit_flag = 0;
    pr->mean_var = (pr->mu > 1.0) ? pr->mu: 1.0;
  }

  /*
  fprintf(stderr," rho1/mu1: %.4f/%.4f mean_var %.4f\n",
	  pr->rho*LN_FACT,pr->mu,pr->mean_var);
  */
  if (n > 1) pr->var_e = (var_y2/n - mean_y2 * mean_y2)/(n-1);
  else pr->var_e = 0.0;

  if (llen->fit_flag) {
    pr->rho2 = covar_xy2 / var_x;
    pr->mu2 = pr->mean_var - pr->rho2 * mean_x;
  }
  else {
    pr->rho2 = 0;
    pr->mu2 = pr->mean_var;
  }

  if (pr->rho2 < 0.0 )
    z = (pr->rho2 * LN_FACT*log((double)llen->max_length) + pr->mu2 > 0.0) ? llen->max_length : exp((-1.0 - pr->mu2 / pr->rho2)/LN_FACT);
  else z =  pr->rho2 ? exp((1.0 - pr->mu2 / pr->rho2)/LN_FACT) : LENGTH_CUTOFF;
  if (z < 2*LENGTH_CUTOFF) z = 2*LENGTH_CUTOFF;

  pr->var_cutoff = pr->rho2 * LN_FACT*log(z) + pr->mu2;
}

/* fit_llens: trim high variance bins
   (1) regress scores vs log(n) using weighted variance
   (2) regress residuals vs log(n)
   (3) remove high variance bins
   (4) calculate mean variance after length regression
*/

void
fit_llens(struct llen_str *llen, struct rstat_str *pr)
{
  int j;
  int n, n_u2;
  double x, y, y2, u, u2, v, z;
  double mean_x, mean_y, var_x, var_y, covar_xy;
  double mean_y2, covar_xy2;
  double mean_u2, mean_3u2, dllj;
  double sum_x, sum_y, sum_x2, sum_xy, sum_v, det, n_w;

/* now fit scores to best linear function of log(n), using
   simple linear regression */
  
  for (llen->min=0; llen->min < llen->max; llen->min++)
    if (llen->hist[llen->min]) break;
  llen->min--;

  for (j = llen->min; j < llen->max; j++) {
    if (llen->hist[j] > 1) {
      dllj = (double)llen->hist[j];
      llen->score_var[j] = (double)llen->score2_sums[j]/dllj
	- (llen->score_sums[j]/dllj)*(llen->score_sums[j]/dllj);
      llen->score_var[j] /= (double)(llen->hist[j]-1);
      if (llen->score_var[j] <= 1.0 ) llen->score_var[j] = 1.0;
    }
  }
	  
  n_w = 0.0;
  sum_x = sum_y = sum_x2 = sum_xy = sum_v = 0;
  for (j = llen->min; j < llen->max; j++)
    if (llen->hist[j] > 1) {
      x = j + 0.5;
      dllj = (double)llen->hist[j];
      n_w += dllj/llen->score_var[j];
      sum_x +=   dllj * x / llen->score_var[j] ;
      sum_y += llen->score_sums[j] / llen->score_var[j];
      sum_x2 +=  dllj * x * x /llen->score_var[j];
      sum_xy +=  x * llen->score_sums[j]/llen->score_var[j];
    }

  det = n_w * sum_x2 - sum_x * sum_x;
  pr->rho = (n_w * sum_xy  - sum_x * sum_y)/det;
  pr->mu = (sum_x2 * sum_y - sum_x * sum_xy)/det;

/* printf(" rho1/mu1: %.2f/%.2f\n",pr->rho*LN_FACT,pr->mu); */

  n = 0;
  mean_x = mean_y = mean_y2 = 0.0;
  var_x = var_y = 0.0;
  covar_xy = covar_xy2 = 0.0;

  for (j = llen->min; j <= llen->max; j++) 
    if (llen->hist[j] > 1 ) {
      n += llen->hist[j];
      x = (double)j + 0.5;
      dllj = (double)llen->hist[j];
      mean_x += dllj * x;
      mean_y += llen->score_sums[j];
      var_x += dllj * x * x;
      var_y += llen->score2_sums[j];
      covar_xy += x * llen->score_sums[j];
    }
  mean_x /= n; mean_y /= n;
  var_x = var_x / n - mean_x * mean_x;
  var_y = var_y / n - mean_y * mean_y;
  
  covar_xy = covar_xy / n - mean_x * mean_y;
/*  pr->rho = covar_xy / var_x;
  pr->mu = mean_y - pr->rho * mean_x;
*/

  mean_y2 = covar_xy2 = 0.0;
  for (j = llen->min; j <= llen->max; j++) 
    if (llen->hist[j] > 1) {
      x = (double)j + 0.5;
      u = pr->rho * x + pr->mu;
      y2 = llen->score2_sums[j] - 2 * llen->score_sums[j] * u + llen->hist[j] * u * u;
      mean_y2 += y2;
      covar_xy2 += x * y2;
    }
  
  mean_y2 /= n;
  covar_xy2 = covar_xy2 / n - mean_x * mean_y2;
  pr->rho2 = covar_xy2 / var_x;
  pr->mu2 = mean_y2 - pr->rho2 * mean_x;

  if (pr->rho2 < 0.0 )
    z = (pr->rho2 * LN_FACT*log((double)llen->max_length) + pr->mu2 > 0.0) ? llen->max_length : exp((-1.0 - pr->mu2 / pr->rho2)/LN_FACT);
  else z =  pr->rho2 ? exp((1.0 - pr->mu2 / pr->rho2)/LN_FACT) : LENGTH_CUTOFF;
  if (z < 2* LENGTH_CUTOFF) z = 2*LENGTH_CUTOFF;

  pr->var_cutoff = pr->rho2*LN_FACT*log(z) + pr->mu2;

/*  fprintf(stderr,"\nminimum allowed predicted variance (%0.2f) at n = %.0f\n",
	 pr->var_cutoff,z);
*/
  mean_u2 = 0.0;
  n_u2 = 0;
  for ( j = llen->min; j < llen->max; j++) {
    y = j+0.5;
    dllj = (double)llen->hist[j];
    x = pr->rho * y + pr->mu;
    v = pr->rho2 * y + pr->mu2;
    if (v < pr->var_cutoff) v = pr->var_cutoff;
    if (llen->hist[j]> 1) {
      u2 =  (llen->score2_sums[j] - 2 * x * llen->score_sums[j] + dllj * x * x) - v*dllj;
      mean_u2 += llen->score_var[j] = u2*u2/(llen->hist[j]-1);
      n_u2++;
      /*      fprintf(stderr," %d (%d) u2: %.2f v*ll: %.2f %.2f\n",
	      j,llen->hist[j],u2,v*dllj,sqrt(llen->score_var[j])); */
    }
    else llen->score_var[j] = -1.0;
  }

  mean_u2 = sqrt(mean_u2/(double)n_u2);
  /* fprintf(stderr," mean s.d.: %.2f\n",mean_u2); */

  mean_3u2 = mean_u2*3.0;

  for (j = llen->min; j < llen->max; j++) {
    if (llen->hist[j] <= 1) continue;
    if (sqrt(llen->score_var[j]) > mean_3u2) {
      /*      fprintf(stderr," removing %d %d %.2f\n",
	     j, (int)(exp((double)j/LN_FACT)-0.5),
	     sqrt(llen->score_var[j]));
	     */
      pr->nb_trimmed++;
      pr->n1_trimmed += llen->hist[j];
      llen->hist[j] = 0;
    }
  }
  fit_llen(llen, pr);
}

struct s2str {double s; int n;};
void s2_sort ( struct s2str *sptr, int n);

void
fit_llen2(struct llen_str *llen, struct rstat_str *pr)
{
  int j;
  int n, n_y2, llen_delta, llen_del05;
  int n_size;
  double x, y2, u;
  double mean_x, mean_y, var_x, var_y, covar_xy;
  double mean_y2, covar_xy2;
  struct s2str *ss2;

  double sum_x, sum_y, sum_x2, sum_xy, sum_v, det, n_w;
  
/* now fit scores to best linear function of log(n), using
   simple linear regression */
  
  for (llen->min=0; llen->min < llen->max; llen->min++)
    if (llen->hist[llen->min]) break;

  for ( ; llen->max > llen->min; llen->max--)
    if (llen->hist[llen->max]) break;

  for (n_size=0,j = llen->min; j < llen->max; j++) {
    if (llen->hist[j] > 1) {
      llen->score_var[j] = llen->score2_sums[j]/(double)llen->hist[j]
	- (llen->score_sums[j]/(double)llen->hist[j])
	* (llen->score_sums[j]/(double)llen->hist[j]);
      llen->score_var[j] /= (double)(llen->hist[j]-1);
      if (llen->score_var[j] <= 1.0 ) llen->score_var[j] = 1.0;
      n_size++;
    }
  }
	  
  n_w = 0.0;
  sum_x = sum_y = sum_x2 = sum_xy = sum_v = 0;
  for (j = llen->min; j < llen->max; j++)
    if (llen->hist[j] > 1) {
      x = j + 0.5;
      n_w += (double)llen->hist[j]/llen->score_var[j];
      sum_x +=   (double)llen->hist[j] * x / llen->score_var[j] ;
      sum_y += llen->score_sums[j] / llen->score_var[j];
      sum_x2 +=  (double)llen->hist[j] * x * x /llen->score_var[j];
      sum_xy +=  x * llen->score_sums[j]/llen->score_var[j];
    }

  if (n_size < 5 ) {
    llen->fit_flag=0;
    pr->rho = 0;
    pr->mu = sum_y/n_w;
  }
  else {
    det = n_w * sum_x2 - sum_x * sum_x;
    if (det > 0.001) {
      pr->rho = (n_w * sum_xy  - sum_x * sum_y)/det;
      pr->rho_e = n_w/det;
      pr->mu = (sum_x2 * sum_y - sum_x * sum_xy)/det;
      pr->mu_e = sum_x2/det;
    }
    else {
      llen->fit_flag = 0;
      pr->rho = 0;
      pr->mu = sum_y/n_w;
    }
  }

  det = n_w * sum_x2 - sum_x * sum_x;
  pr->rho = (n_w * sum_xy  - sum_x * sum_y)/det;
  pr->mu = (sum_x2 * sum_y - sum_x * sum_xy)/det;

/*   fprintf(stderr," rho1/mu1: %.2f/%.2f\n",pr->rho*LN_FACT,pr->mu); */

  n = 0;
  mean_x = mean_y = mean_y2 = 0.0;
  var_x = var_y = 0.0;
  covar_xy = covar_xy2 = 0.0;

  for (j = llen->min; j <= llen->max; j++) 
    if (llen->hist[j] > 1 ) {
      n += llen->hist[j];
      x = (double)j + 0.5;
      mean_x += (double)llen->hist[j] * x;
      mean_y += llen->score_sums[j];
      var_x += (double)llen->hist[j] * x * x;
      var_y += llen->score2_sums[j];
      covar_xy += x * llen->score_sums[j];
    }
  mean_x /= n; mean_y /= n;
  var_x = var_x / n - mean_x * mean_x;
  var_y = var_y / n - mean_y * mean_y;
  
  covar_xy = covar_xy / n - mean_x * mean_y;
/*
  pr->rho = covar_xy / var_x;
  pr->mu = mean_y - pr->rho * mean_x;
*/

  if ((ss2=(struct s2str *)calloc(llen->max+1,sizeof(struct s2str)))==NULL) {
    fprintf(stderr," cannot allocate ss2\n");
    return;
  }

  mean_y2 = 0.0;
  n_y2 = n = 0;
  for (j = llen->min; j <= llen->max; j++) 
    if (llen->hist[j] > VHISTC) {
      n++;
      n_y2 += ss2[j].n = llen->hist[j];
      x = (double)j + 0.5;
      u = pr->rho * x + pr->mu;
      ss2[j].s = y2 = llen->score2_sums[j] - 2*llen->score_sums[j]*u + llen->hist[j]*u*u;
      mean_y2 += y2;
    }
  pr->mean_var = mean_y2/(double)n_y2;

  s2_sort(ss2+llen->min,llen->max-llen->min+1);
  
  /*  fprintf(stderr,"llen->min: %d, max: %d\n",llen->min,llen->max); */
  llen_delta = 0;
  for (j=llen->min; j<=llen->max; j++) {
    if (ss2[j].n > 1) {
      llen_delta++;
/*      fprintf(stderr,"%d\t%d\t%.2f\t%.4f\n",
	      j,ss2[j].n,ss2[j].s,ss2[j].s/ss2[j].n);
*/
    }
  }

  llen_del05 = llen_delta/20;
  mean_y2 = 0.0;
  n_y2 = 0;
  for (j = llen->min; j<llen->min+llen_del05; j++) {
    pr->n1_trimmed += ss2[j].n;
    pr->nb_trimmed++;
  }
  for (j = llen->min+llen_del05; j <= llen->min+llen_delta-llen_del05; j++) 
    if (ss2[j].n > 1) {
      mean_y2 += ss2[j].s;
      n_y2 += ss2[j].n;
    }
  for (j = llen->min+llen_delta-llen_del05+1; j< llen->max; j++) {
    pr->n1_trimmed += ss2[j].n;
    pr->nb_trimmed++;
  }
  
  free(ss2);
  if (n_y2 > 1) pr->mean_var = mean_y2/(double)n_y2;

  /*    fprintf(stderr," rho1/mu1: %.4f/%.4f mean_var: %.4f/%d\n",
	  pr->rho*LN_FACT,pr->mu,pr->mean_var,n); */

    pr->var_e = 0.0;
}

/* REG_STATS - Z() from rho/mu/mean_var */
double find_zr(int score, double escore, int length, double comp, struct pstat_str *pu)
{
  double log_len, z;
  
  if (score <= 0) return 0;
  if ( length < LENGTH_CUTOFF) return 0;

  log_len = LN_FACT*log((double)(length));
/*  var = pu->r_u.rg.rho2 * log_len + pu->r_u.rg.mu2;
  if (var < pu->r_u.rg.var_cutoff) var = pu->r_u.rg.var_cutoff;
*/

  z = ((double)score - pu->r_u.rg.rho * log_len - pu->r_u.rg.mu) / sqrt(pu->r_u.rg.mean_var);

  return (50.0 + z*10.0);
}

/* REG2_STATS Z() from rho/mu, rho2/mu2 */
double find_zr2(int score, double escore, int length, double comp, struct pstat_str *pu)
{
  double log_len, var;
  double z;
  
  if ( length < LENGTH_CUTOFF) return 0;

  log_len = LN_FACT*log((double)(length));

  var = pu->r_u.rg.rho2 * log_len + pu->r_u.rg.mu2;
  if (var < pu->r_u.rg.var_cutoff) var = pu->r_u.rg.mean_var;

  z = ((double)score - pu->r_u.rg.rho * log_len - pu->r_u.rg.mu) / sqrt(var);

  return (50.0 + z*10.0);
}

#ifdef USE_LNSTATS
/* LN_STATS - ln()-scaled mu, mean_var */
double find_zl(int score, int length, double comp, struct pstat_str *pu)
{
  double ls, z;
  
  ls = (double)score*LN200/log((double)length);

  z = (ls - pu->r_u.rg.mu) / sqrt(pu->r_u.rg.mean_var);

  return (50.0 + z*10.0);
}
#endif

/* MLE_STATS - Z() from MLE for lambda, K */
double
find_ze(int score, double escore, int length, double comp, struct pstat_str *pu)
{
  double z, mp, np, a_n1;
  
  a_n1 = (double)length; 

  mp = pu->r_u.ag.a_n0;
  np = a_n1;

  if (np < 1.0) np = 1.0;
  if (mp < 1.0) mp = 1.0;

  z = pu->r_u.ag.Lambda * score - log(pu->r_u.ag.K * np * mp);

  z = -z + EULER_G;
  z /= - PI_SQRT6;

  return (50.0 + z*10.0);
}

/* MLE2_STATS - Z() from MLE for mle_a0..2, mle_b1, length, comp */
double
find_ze2(int score, double escore, int length, double comp, struct pstat_str *pu)
{
  double z, mp, np, a_n1;
  
  a_n1 = (double)length; 

  if (comp <= 0.0) comp = pu->r_u.m2.ave_comp;

  /* avoid very biased comp estimates */
  /* comp = exp((4.0*log(comp)+log(pu->r_u.m2.ave_comp))/5.0); */

  mp = pu->r_u.m2.a_n0;
  np = a_n1;

  if (np < 1.0) np = 1.0;
  if (mp < 1.0) mp = 1.0;

  z = (-(pu->r_u.m2.mle2_a0 + pu->r_u.m2.mle2_a1 * comp + pu->r_u.m2.mle2_a2 * comp * log(np * mp)) + score) / (pu->r_u.m2.mle2_b1 * comp);

  z = -z + EULER_G;
  z /= - PI_SQRT6;

  return (50.0 + z*10.0);
}

/* AG_STATS - Altschul-Gish Lamdba, K */
double
find_za(int score, double escore, int length, double comp, struct pstat_str *pu)
{
  double z, mp, np, a_n1, a_n1f;
  
  a_n1 = (double)length; 
  a_n1f = log(a_n1)/pu->r_u.ag.H;

  mp = pu->r_u.ag.a_n0 - pu->r_u.ag.a_n0f - a_n1f;
  np = a_n1 - pu->r_u.ag.a_n0f - a_n1f;

  if (np < 1.0) np = 1.0;
  if (mp < 1.0) mp = 1.0;

  z = pu->r_u.ag.Lambda * score - log(pu->r_u.ag.K * np * mp);

  z = -z + EULER_G;
  z /= - PI_SQRT6;

  return (50.0 + z*10.0);
}

double find_zn(int score, double escore, int length, double comp, struct pstat_str *pu)
{
  double z;
  
  z = ((double)score - pu->r_u.rg.mu) / sqrt(pu->r_u.rg.mean_var);

  return (50.0 + z*10.0);
}

/* computes E value for a given z value, assuming extreme value distribution */
double
z_to_E(double zs, long entries, struct db_str db)
{
  double e, n;

  /*  if (db->entries < 5) return (double)db.entries; */
  if (entries < 1) { n = db.entries;}
  else {n = entries;}

  if (zs > ZS_MAX) return 0.0;

#ifndef NORMAL_DIST
  e = exp(- PI_SQRT6 * zs - .577216);
  return n * (e > .01 ? 1.0 - exp(-e) : e);
#else
  return n * erfc(zs/M_SQRT2)/2.0; 
#endif
}

double
zs_to_p(double zs)
{
  double e, z;

  /*  if (db.entries < 5) return 0.0; */

  z = (zs - 50.0)/10.0;

  if (z > ZS_MAX) return 0.0;

#ifndef NORMAL_DIST
  e = exp(- PI_SQRT6 * z - EULER_G);
  return (e > .01 ? 1.0 - exp(-e) : e);
#else
  return erfc(zs/M_SQRT2)/2.0; 
#endif
}

double
zs_to_bit(double zs, int n0, int n1)
{
  double z, a_n0, a_n1;

  z = (zs - 50.0)/10.0;
  a_n0 = (double)n0;
  a_n1 = (double)n1;

  return (PI_SQRT6 * z + EULER_G + log(a_n0*a_n1))/M_LN2 ;
}

/* computes E-value for a given z value, assuming extreme value distribution */
double
zs_to_E(double zs,int n1, int dnaseq, long entries, struct db_str db)
{
  double e, z, k;

  /*  if (db->entries < 5) return 0.0; */

  z = (zs - 50.0)/10.0;

  if (z > ZS_MAX ) return 0.0;

  if (entries < 1) entries = db.entries;

  if (dnaseq == SEQT_DNA || dnaseq == SEQT_RNA) {
    k = (double)db.length /(double)n1;
    if (db.carry > 0) {
      k += ((double)db.carry * (double)LONG_MAX)/(double)n1;
    }
  }
  else k = (double)entries;

  if (k < 1.0) k = 1.0;

#ifndef NORMAL_DIST
  z *= PI_SQRT6;
  z += EULER_G;
  e = exp(-z);
  return k * (e > .01 ? 1.0 - exp(-e) : e);
#else
  return k * erfc(z/M_SQRT2)/2.0; 
#endif
}

#ifdef NORMAL_DIST
double np_to_z(double, int *);
#endif

/* computes E-value for a given z value, assuming extreme value distribution */
double
E_to_zs(double E, long entries)
{
  double e, z;
  int error;

  e = E/(double)entries;

#ifndef NORMAL_DIST
  z = (log(e)+EULER_G)/(- PI_SQRT6);
  return z*10.0+50.0;
#else
  z = np_to_z(1.0-e,&error);

  if (!error) return z*10.0+50.0;
  else return 0.0;
#endif
}

/* computes 1.0 - E value for a given z value, assuming extreme value
   distribution */
double
zs_to_Ec(double zs, long entries)
{
  double e, z;

  if (entries < 5) return 0.0;

  z = (zs - 50.0)/10.0;

  if (z > ZS_MAX) return 1.0;

#ifndef NORMAL_DIST
  e =  exp(- PI_SQRT6 * z - EULER_G);
  return (double)entries * (e > .01 ?  exp(-e) : 1.0 - e);
#else
  return (double)entries*erf(z/M_SQRT2)/2.0; 
#endif
}

/* calculate a threshold score, given an E() value and Lambda,K,H */

int
E1_to_s(double e_val, int n0, int n1, struct pstat_str *pu) {
  double mp, np, a_n0, a_n0f, a_n1;
  int score;

  a_n0 = (double)n0;
  a_n1 = (double)n1;
  a_n0f = log(pu->r_u.ag.K * a_n0 * a_n1)/pu->r_u.ag.H;

  mp = a_n0 - a_n0f;
  np = a_n1 - a_n0f;

  if (np < 1.0) np = 1.0;
  if (mp < 1.0) mp = 1.0;

  score = (int)((log( pu->r_u.ag.K * mp * np) - log(e_val))/pu->r_u.ag.Lambda +0.5);
  if (score < 0) score = 0;
  return score;
}

/* no longer used; stat_str returned by process_hist
void
summ_stats(char *s_str, struct pstat_str *pu)
{
  strcpy(s_str,f_string);
}
*/

void
vsort(v,s,n)
	double *v; int *s, n;
{
  int gap, i, j;
  double tmp;
  int itmp;
	
  for (gap=n/2; gap>0; gap/=2)
    for (i=gap; i<n; i++)
      for (j=i-gap; j>=0; j -= gap) {
	if (v[j] >= v[j+gap]) break;
	tmp = v[j]; v[j]=v[j+gap]; v[j+gap]=tmp;
	itmp = s[j]; s[j]=s[j+gap]; s[j+gap]=itmp;
      }
}

/*
void s_sort (double **ptr, int nbest)
{
  int gap, i, j;
  double *tmp;

  for (gap = nbest/2; gap > 0; gap /= 2)
    for (i = gap; i < nbest; i++)
      for (j = i - gap; j >= 0; j-= gap) {
	if (*ptr[j] >= *ptr[j + gap]) break;
	tmp = ptr[j];
	ptr[j] = ptr[j + gap];
	ptr[j + gap] = tmp;
      }
}
*/

void ss_sort (int *ptr, int n)
{
  int gap, i, j;
  int tmp;

  for (gap = n/2; gap > 0; gap /= 2)
    for (i = gap; i < n; i++)
      for (j = i - gap; j >= 0; j-= gap) {
	if (ptr[j] >= ptr[j + gap]) break;
	tmp = ptr[j];
	ptr[j] = ptr[j + gap];
	ptr[j + gap] = tmp;
      }
}


void s2_sort (struct s2str *ptr, int n)
{
  int gap, i, j;
  struct s2str tmp;

  for (gap = n/2; gap > 0; gap /= 2)
    for (i = gap; i < n; i++)
      for (j = i - gap; j >= 0; j-= gap) {
	if (ptr[j].s >= ptr[j + gap].s) break;
	tmp.s = ptr[j].s;
	tmp.n = ptr[j].n;
	ptr[j].s = ptr[j + gap].s;
	ptr[j].n = ptr[j + gap].n;
	ptr[j + gap].s = tmp.s;
	ptr[j + gap].n = tmp.n;
      }
}

void last_stats() {}

void
scale_scores(struct beststr **bptr, int nbest, struct db_str db,
	     struct pstruct pst, struct pstat_str *rs)
{
  int i;
  double zscore;

  if (pst.zsflag < 0 ) return;

  for (i=0; i<nbest; i++) {
    zscore = find_zp(bptr[i]->score[pst.score_ix], bptr[i]->escore,
		     bptr[i]->n1,bptr[i]->comp,rs);
    bptr[i]->zscore = zscore;
    bptr[i]->escore
      =zs_to_E(zscore,bptr[i]->n1,pst.dnaseq, pst.zdb_size,db);
  }
  sortbeste(bptr,nbest);
}

#ifdef NORMAL_DIST
/*     ALGORITHM AS241  APPL. STATIST. (1988) VOL. 37, NO. 3

       Produces the normal deviate Z corresponding to a given lower
       tail area of P; Z is accurate to about 1 part in 10**16.

       The hash sums below are the sums of the mantissas of the
       coefficients.   They are included for use in checking
       transcription.
*/

double np_to_z(double p, int *fault) {

  double q, r, ppnd16;

  double zero = 0.0, one = 1.0, half = 0.5;
  double split1 = 0.425, split2 = 5.0;
  double const1 = 0.180625, const2 = 1.6;

/*   Coefficients for P close to 0.5 */

  double a0 = 3.3871328727963666080e0;
  double a1 = 1.3314166789178437745e+2;
  double a2 = 1.9715909503065514427e+3;
  double a3 = 1.3731693765509461125e+4;
  double a4 = 4.5921953931549871457e+4;
  double a5 = 6.7265770927008700853e+4;
  double a6 = 3.3430575583588128105e+4;
  double a7 = 2.5090809287301226727e+3;
  double b1 = 4.2313330701600911252e+1;
  double b2 = 6.8718700749205790830e+2;
  double b3 = 5.3941960214247511077e+3;
  double b4 = 2.1213794301586595867e+4;
  double b5 = 3.9307895800092710610e+4;
  double b6 = 2.8729085735721942674e+4;
  double b7 = 5.2264952788528545610e+3;

  double sum_ab= 55.8831928806149014439;
/* 
  Coefficients for P not close to 0, 0.5 or 1.
*/

  double c0 = 1.42343711074968357734;
  double c1 = 4.63033784615654529590;
  double c2 = 5.76949722146069140550;
  double c3 = 3.64784832476320460504;
  double c4 = 1.27045825245236838258;
  double c5 = 2.41780725177450611770e-1;
  double c6 = 2.27238449892691845833e-2;
  double c7 = 7.74545014278341407640e-4;
  double d1 = 2.05319162663775882187;
  double d2 = 1.67638483018380384940;
  double d3 = 6.89767334985100004550e-1;
  double d4 = 1.48103976427480074590e-1;
  double d5 = 1.51986665636164571966e-2;
  double d6 = 5.47593808499534494600e-4;
  double d7 = 1.05075007164441684324e-9;

  double sum_cd=49.33206503301610289036;
/*
       Coefficients for P near 0 or 1.
*/
  double e0 = 6.65790464350110377720e0;
  double e1 = 5.46378491116411436990e0;
  double e2 = 1.78482653991729133580e0;
  double e3 = 2.96560571828504891230e-1;
  double e4 = 2.65321895265761230930e-2;
  double e5 = 1.24266094738807843860e-3;
  double e6 = 2.71155556874348757815e-5;
  double e7 = 2.01033439929228813265e-7;
  double f1 = 5.99832206555887937690e-1;
  double f2 = 1.36929880922735805310e-1;
  double f3 = 1.48753612908506148525e-2;
  double f4 = 7.86869131145613259100e-4;
  double f5 = 1.84631831751005468180e-5;
  double f6 = 1.42151175831644588870e-7;
  double f7 = 2.04426310338993978564e-15;

  double sum_ef=47.52583317549289671629;

  double sum_tmp = 0.0;

  /*
  sum_tmp = a0+a1+a2+a3+a4+a5+a6+a7+b1+b2+b3+b4+b5+b6+b7;
  if (fabs(sum_tmp - sum_ab) > 1e-12) {
    fprintf (stderr," sum_ab error: %lg %lg\n",sum_tmp,sum_ab);
    *fault = 1;
    return zero;
  }

  sum_tmp = c0+c1+c2+c3+c4+c5+c6+c7+d1+d2+d3+d4+d5+d6+d7;
  if (fabs(sum_tmp - sum_cd) > 1e-12) {
    fprintf (stderr," sum_cd error: %lg %lg\n",sum_tmp,sum_cd);
    *fault = 1;
    return zero;
  }
  sum_tmp = e0+e1+e2+e3+e4+e5+e6+e7+f1+f2+f3+f4+f5+f6+f7;
  if (fabs(sum_tmp - sum_ef) > 1e-12) {
    fprintf (stderr," sum_ef error: %lg %lg\n",sum_tmp,sum_ef);
    *fault = 1;
    return zero;
  }
  */

  *fault = 0;
  q = p - half;
  if (fabs(q) <= split1) {
    r = const1 - q * q;
    return q * (((((((a7 * r + a6) * r + a5) * r + a4) * r + a3)
		    * r + a2) * r + a1) * r + a0) /
      (((((((b7 * r + b6) * r + b5) * r + b4) * r + b3)
	 * r + b2) * r + b1) * r + one);
  }
  else {
    r = (q < zero) ?  p : one - p;
    if (r <= zero) {
      *fault = 1;
      return zero;
    }
    r = sqrt(-log(r));
    if (r <= split2) {
      r -= const2;
      ppnd16 = (((((((c7 * r + c6) * r + c5) * r + c4) * r + c3)
		  * r + c2) * r + c1) * r + c0) /
	(((((((d7 * r + d6) * r + d5) * r + d4) * r + d3)
	   * r + d2) * r + d1) * r + one);
    }
    else {
      r -= split2;
      ppnd16 = (((((((e7 * r + e6) * r + e5) * r + e4) * r + e3)
		  * r + e2) * r + e1) * r + e0) /
	(((((((f7 * r + f6) * r + f5) * r + f4) * r + f3)
	   * r + f2) * r + f1) * r + one);
    }
    if (q < zero) return -ppnd16;
    else return ppnd16;
  }
}
#endif
