+++ /dev/null
-/* scaleswn.c */
-
-/* $Name: fa_34_26_5 $ - $Id: scaleswn.c,v 1.60 2007/04/26 18:32:48 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
-*/
-
-#define MAX_NIT 100
-
-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;
- int nit = 0;
- 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;
- nit++;
- } while (fabs((lambda - old_lambda)/lambda) > TINY && nit < MAX_NIT);
-
- /* fprintf(stderr," mle_cen nit: %d\n",nit); */
-
- if (nit >= MAX_NIT) return -1;
-
- 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 || pst.zsflag_f < 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
git://source.jalview.org / jabaws.git / blobdiff