/*	initfa.c	*/

/* $Name: fa34t25d7 $ - $Id: initfa.c,v 1.133 2006/05/12 14:20:21 wrp Exp $ */

/* copyright (c) 1996, 1997, 1998  William R. Pearson and the U. of Virginia */

/* init??.c files provide function specific initializations */

/* h_init()	- called from comp_lib.c, comp_thr.c to initialize pstruct ppst
   		  which includes the alphabet, and pam matrix

   alloc_pam()	- allocate pam matrix space
   init_pam2()	- convert from 1D to 2D pam

   init_pamx()	- convert from 1D to 2D pam

   f_initenv()	- set up mngmsg and pstruct defaults
   f_getopt()	- read fasta specific command line options
   f_getarg()	- read ktup

   resetp()	- reset the parameters, scoring matrix for DNA-DNA/DNA-prot

   query_parm()	- ask for ktup
   last_init()	- some things must be done last

   f_initpam()	- set some parameters based on the pam matrix
*/

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

#ifdef UNIX
#include <sys/types.h>
#include <sys/stat.h>
#endif

#include "defs.h"
#include "structs.h"
#include "param.h"

#ifndef PCOMPLIB
#include "mw.h"
#else
#include "p_mw.h"
#endif

#define XTERNAL
#include "upam.h"
#include "uascii.h"
#undef XTERNAL

#define MAXWINDOW 32

int initpam(char *, struct pstruct *);
void extend_pssm(unsigned char *aa0, int n0, struct pstruct *ppst);
void build_xascii(int *qascii, char *save_str);
void ann_ascii(int *qascii, char *ann_arr);
extern int nrand(int);

/*  at some point, all the defaults should be driven from this table */
/*
#pgm	q_seq	l_seq	p_seq	matrix	g_open	g_ext	fr_shft	e_cut	ktup
#	-n/-p		-s	-e	-f	-h/-j	-E	argv[3]
fasta	prot(0)	prot(0)	prot(0)	bl50	-10	-2	-	10.0	2
fasta	dna(1)	dna(1)	dna(1)	+5/-4	-14	-4	-	2.0	6
ssearch	prot(0)	prot(0)	prot(0)	bl50	-10	-2	-	10.0	-
ssearch	dna(1)	dna(1)	dna(1)	+5/-4	-14	-4	-	2.0	-
fastx	dna(1)	prot(0)	prot(0)	BL50	-12	-2	-20	5.0	2
fasty	dna(1)	prot(0)	prot(0)	BL50	-12	-2	-20/-24	5.0	2
tfastx	dna(1)	prot(0)	prot(0)	BL50	-14	-2	-20	5.0	2
tfasty	dna(1)	prot(0)	prot(0)	BL50	-14	-2	-20/-24	5.0	2
fasts	prot(0)	prot(0)	prot(0)	MD20-MS	-	-	-	5.0	-
fasts	dna(1)	dna(1)	dna(1)	+2/-4	-	-	-	5.0	1
tfasts	prot(0)	dna(1)	prot(0)	MD10-MS	-	-	-	2.0	1
fastf	prot(0)	prot(0)	prot(0)	MD20	-	-	-	2.0	1
tfastf	prot(0)	dna(1)	prot(0)	MD10	-	-	-	1.0	1
fastm	prot(0)	prot(0)	prot(0)	MD20	-	-	-	5.0	1
fastm	dna(1)	dna(1)	dna(1)	+2/-4	-	-	-	2.0	1
tfastm	prot(0)	dna(1)	prot(0)	MD10	-	-	-	2.0	1
*/

struct pgm_def_str {
  int pgm_id;
  char *prog_func;
  char *pgm_abbr;
  char *iprompt0;
  char *ref_str;
  int PgmDID;
  char *smstr;
  int g_open_mod;
  int gshift;
  int hshift;
  int e_cut;
  int ktup;
};

char *ref_str_a[]={
  "\nPlease cite:\n W.R. Pearson & D.J. Lipman PNAS (1988) 85:2444-2448\n",
  "\nPlease cite:\n T. F. Smith and M. S. Waterman, (1981) J. Mol. Biol. 147:195-197; \n W.R. Pearson (1991) Genomics 11:635-650\n",
 "\nPlease cite:\n Pearson et al, Genomics (1997) 46:24-36\n",
 "\nPlease cite:\n Mackey et al. Mol. Cell. Proteomics  (2002) 1:139-147\n",
  "\nPlease cite:\n W.R. Pearson (1996) Meth. Enzymol. 266:227-258\n"
};

#define FA_PID	1
#define SS_PID	2
#define FX_PID	3
#define FY_PID	4
#define FS_PID	5
#define FF_PID	6
#define FM_PID	7
#define RSS_PID	8
#define RFX_PID	9
#define SSS_PID 10	/* old (slow) non-PG Smith-Waterman */
#define TFA_PID	FA_PID+10
#define TFX_PID	FX_PID+10
#define TFY_PID	FY_PID+10
#define TFS_PID	FS_PID+10
#define TFF_PID	FF_PID+10
#define TFM_PID FM_PID+10

struct pgm_def_str
pgm_def_arr[20] = {
  {0, "", "", "", NULL, 400, "", 0, 0, 0, 1.0, 0 },  /* 0 */
  {FA_PID, "FASTA", "fa",
   "FASTA searches a protein or DNA sequence data bank",
   NULL, 401, "BL50", 0, 0, 0, 10.0, 2}, /* 1 - FASTA */
  {SS_PID, "SSEARCH","gsw","SSEARCH searches a sequence data bank",
   NULL, 404, "BL50", 0, 0, 0, 10.0, 0}, /* 2 - SSEARCH */
  {FX_PID, "FASTX","fx",
   "FASTX compares a DNA sequence to a protein sequence data bank",
   NULL, 405, "BL50", -2, -20, 0, 5.0, 2}, /* 3 - FASTX */
  {FY_PID, "FASTY", "fy",
   "FASTY compares a DNA sequence to a protein sequence data bank",
   NULL, 405, "BL50", -2, -20, -24, 5.0, 2}, /* 4 - FASTY */
  {FS_PID, "FASTS", "fs",
   "FASTS compares linked peptides to a protein data bank",
   NULL, 400, "MD20-MS", 0, 0, 0, 5.0, 1}, /* 5 - FASTS */
  {FF_PID, "FASTF", "ff",
   "FASTF compares mixed peptides to a protein databank",
   NULL, 400, "MD20", 0, 0, 0, 2.0, 1 }, /* 6 - FASTF */
  {FM_PID, "FASTM", "fm",
   "FASTM compares ordered peptides to a protein data bank",
     NULL, 400, "MD20", 0, 0, 0, 5.0, 1 }, /* 7 - FASTM */
  {RSS_PID, "PRSS", "rss",
   "PRSS evaluates statistical signficance using Smith-Waterman",
   NULL, 401, "BL50", 0, 0, 0, 1000.0, 0 }, /* 8 - PRSS */
  {RFX_PID,"PRFX", "rfx",
   "PRFX evaluates statistical signficance using FASTX",
   NULL, 401, "BL50", -2, -20, -24, 1000.0, 2 }, /* 9 - PRFX */
  {SSS_PID, "OSEARCH","ssw","OSEARCH searches a sequence data bank",
   NULL, 404, "BL50", 0, 0, 0, 10.0, 0}, /* 2 - OSEARCH */
  {TFA_PID, "TFASTA", "tfa",
   "TFASTA compares a protein  to a translated DNA data bank",
   NULL, 402, "BL50", -2, 0, 0, 5.0, 2 },
  {0, "", "", "", NULL, 400, "", 0, 0, 0, 1.0, 0 },  /* 0 */
  {TFX_PID, "TFASTX", "tfx",
   "TFASTX compares a protein to a translated DNA data bank",
   NULL, 406, "BL50", -2, -20, 0, 2.0, 2},
  {TFY_PID, "TFASTY", "tfy",
   "TFASTY compares a protein to a translated DNA data bank",
   NULL, 406, "BL50", -2, -20, -24, 2.0, 2},
  {TFS_PID, "TFASTS", "tfs",
   "TFASTS compares linked peptides to a translated DNA data bank",
   NULL, 400, "MD10-MS", 0, 0, 0, 2.0, 2 },
  {TFF_PID, "TFASTF", "tff",
   "TFASTF compares mixed peptides to a protein databank",
   NULL, 400, "MD10", 0, 0, 0, 1.0, 1 },
  {TFM_PID, "TFASTM", "tfm",
   "TFASTM compares ordered peptides to a translated DNA databank",
   NULL, 400, "MD10", 0, 0, 0, 1.0, 1 }
};

struct msg_def_str {
  int pgm_id;
  int q_seqt;
  int l_seqt;
  int p_seqt;
  int sw_flag;
  int stages;
  int qframe;
  int nframe;
  int nrelv, srelv, arelv;
  char *f_id0, *f_id1, *label;
};

/* pgm_id    q_seqt     l_seqt   p_seqt sw_f st qf nf nrv srv arv s_ix */
struct msg_def_str msg_def_arr[20] = {
  {0, 0, 0,  0, 0, 0, 0, 0, 0, 0, 0, "", "", ""},	/* ID=0 */
  {FA_PID, SEQT_UNK, SEQT_PROT, SEQT_PROT, 1, 1, 1, -1, 3, 1, 3,
   "fa","sw", "opt"},
  {SS_PID, SEQT_UNK, SEQT_PROT, SEQT_PROT, 1, 1, 1, -1, 1, 1, 1,
   "sw","sw", "s-w"},
  {FX_PID, SEQT_DNA, SEQT_PROT, SEQT_PROT, 1, 1, 2, -1, 3, 1, 3,
   "fx","sx", "opt"},
  {FY_PID, SEQT_DNA, SEQT_PROT, SEQT_PROT, 1, 1, 2, -1, 3, 1, 3,
   "fy","sy", "opt"},
  {FS_PID, SEQT_UNK, SEQT_PROT, SEQT_PROT, 1, 1, 1, -1, 3, 2, 3,
   "fs","fs", "initn init1"},
  {FF_PID, SEQT_PROT,SEQT_PROT, SEQT_PROT, 1, 1, 1, -1, 3, 2, 3,
   "ff","ff", "initn init1"},
  {FM_PID, SEQT_PROT,SEQT_PROT, SEQT_PROT, 1, 1, 1, -1, 3, 2, 3,
   "fm","fm","initn init1"},
  {RSS_PID, SEQT_UNK,SEQT_PROT, SEQT_PROT, 0, 1, 1, -1, 1, 1, 1,
   "rss","sw","s-w"},
  {RFX_PID, SEQT_DNA,SEQT_PROT, SEQT_PROT, 0, 1, 2, -1, 3, 1, 3,
   "rfx","sx","opt"},
  {SSS_PID, SEQT_UNK,SEQT_PROT, SEQT_PROT, 1, 1, 1, -1, 1, 1, 1,
   "sw","sw", "s-w"},
  {TFA_PID, SEQT_PROT,SEQT_DNA, SEQT_PROT, 0, 1, 1, 6, 3, 1, 3,
   "tfa","fa","initn init1"},
  {0, 0, 0,  0, 0, 0, 0, 0, 0, 0, 0, "", "", ""},	/* ID=12 */
  {TFX_PID, SEQT_PROT,SEQT_DNA, SEQT_PROT, 1, 1, 1, 2, 3, 2, 3,
   "tfx","sx","initn opt"},
  {TFY_PID, SEQT_PROT,SEQT_DNA, SEQT_PROT, 1, 1, 1, 2, 3, 2, 3,
   "tfy","sy","initn opt"},
  {TFS_PID, SEQT_PROT,SEQT_DNA, SEQT_PROT, 1, 1, 1, 6, 3, 2, 3,
   "tfs","fs","initn init1"},
  {TFF_PID, SEQT_PROT,SEQT_DNA, SEQT_PROT, 1, 1, 1, 6, 3, 2, 3,
   "tff","ff","initn init1"},
  {TFM_PID, SEQT_PROT,SEQT_DNA, SEQT_PROT, 1, 1, 1, 6, 3, 2, 3,
   "tfm","fm","initn init1"}
};

int
get_pgm_id() {

  int rval;

#ifdef FASTA
#ifndef TFAST
  pgm_def_arr[FA_PID].ref_str = ref_str_a[0];
  rval=FA_PID;
#else
  pgm_def_arr[TFA_PID].ref_str = ref_str_a[0];
  rval=TFA_PID;
#endif
#endif

#ifdef FASTX
#ifndef TFAST
#ifndef PRSS
  pgm_def_arr[FX_PID].ref_str = ref_str_a[2];
  rval=FX_PID;
#else
  pgm_def_arr[RFX_PID].ref_str = ref_str_a[2];
  rval=RFX_PID;
#endif
#else
  pgm_def_arr[TFX_PID].ref_str = ref_str_a[2];
  rval=TFX_PID;
#endif
#endif

#ifdef FASTY
#ifndef TFAST
  pgm_def_arr[FY_PID].ref_str = ref_str_a[2];
  rval=FY_PID;
#else
  pgm_def_arr[TFY_PID].ref_str = ref_str_a[2];
  rval=TFY_PID;
#endif
#endif

#ifdef FASTS
#ifndef TFAST
  pgm_def_arr[FS_PID].ref_str = ref_str_a[3];
  rval=FS_PID;
#else
  pgm_def_arr[TFS_PID].ref_str = ref_str_a[3];
  rval=TFS_PID;
#endif
#endif

#ifdef FASTF
#ifndef TFAST
  pgm_def_arr[FF_PID].ref_str = ref_str_a[3];
  rval=FF_PID;
#else
  pgm_def_arr[TFF_PID].ref_str = ref_str_a[3];
  rval=TFF_PID;
#endif
#endif

#ifdef FASTM
#ifndef TFAST
  pgm_def_arr[FM_PID].ref_str = ref_str_a[3];
  rval=FM_PID;
#else
  pgm_def_arr[TFM_PID].ref_str = ref_str_a[3];
  rval=TFM_PID;
#endif
#endif

#ifdef SSEARCH
  pgm_def_arr[SS_PID].ref_str = ref_str_a[1];
  rval=SS_PID;
#endif

#ifdef OSEARCH
  pgm_def_arr[SSS_PID].ref_str = ref_str_a[1];
  rval=SSS_PID;
#endif

#ifdef PRSS
#ifndef FASTX
  pgm_def_arr[RSS_PID].ref_str = ref_str_a[4];
  rval=RSS_PID;
#endif
#endif

  return rval;
}

char *iprompt1=" test sequence file name: ";
char *iprompt2=" database file name: ";

char *verstr="version 3.4t25 Sept 2, 2005";

char   *s_optstr = "13Ac:f:g:h:j:k:oP:r:s:St:x:y:";

static int mktup=2;
static int ktup_set = 0;
static int gap_set=0;
static int del_set=0;
static int mshuff_set = 0;

extern int max_workers;

extern void s_abort(char *, char *);
extern void init_ascii(int ext_sq, int *sascii, int dnaseq);
extern int standard_pam(char *smstr, struct pstruct *ppst,
			int del_set, int gap_set);
extern void mk_n_pam(int *arr,int siz, int mat, int mis);
extern int karlin(int , int, double *, double *, double *);
extern void init_karlin_a(struct pstruct *, double *, double **);
extern int do_karlin_a(int **, struct pstruct *, double *,
		       double *, double *, double *, double *);

#if defined(TFAST) || defined(FASTX) || defined(FASTY)
extern void aainit(int tr_type, int debug);
#endif

char *iprompt0, *prog_func, *refstr;


/* Sets defaults assuming a protein sequence */
void h_init (struct pstruct *ppst, struct mngmsg *m_msp, char *pgm_abbr)
{
  struct pgm_def_str pgm_def;
  int i, pgm_id;

  ppst->pgm_id  = pgm_id =  get_pgm_id();
  pgm_def = pgm_def_arr[pgm_id];

  /* check that pgm_def_arr[] is valid */
  if (pgm_def.pgm_id != pgm_id) {
    fprintf(stderr,
	    "**pgm_def integrity failure: def.pgm_id %d != pgm_id %d**\n",
	    pgm_def.pgm_id, pgm_id);
    exit(1);
  }

  /* check that msg_def_arr[] is valid */
  if (msg_def_arr[pgm_id].pgm_id != pgm_id) {
    fprintf(stderr,
	    "**msg_def integrity failure: def.pgm_id %d != pgm_id %d**\n",
	    msg_def_arr[pgm_id].pgm_id, pgm_id);
    exit(1);
  }

  strncpy(pgm_abbr,pgm_def.pgm_abbr,MAX_SSTR);
  iprompt0 = pgm_def.iprompt0;
  refstr = pgm_def.ref_str;
  prog_func = pgm_def.prog_func;

  /* set up DNA query sequence if required*/
  if (msg_def_arr[pgm_id].q_seqt == SEQT_DNA) {
    memcpy(qascii,nascii,sizeof(qascii));
    m_msp->qdnaseq = SEQT_DNA;
  }
  else { 	/* when SEQT_UNK, start with protein */
    memcpy(qascii,aascii,sizeof(qascii));
    m_msp->qdnaseq = msg_def_arr[pgm_id].q_seqt;
  }

#if defined(FASTF) || defined(FASTS) || defined(FASTM)
  qascii[','] = ESS;
  /* also initialize aascii, nascii for databases */
  qascii['*'] = NA;
#endif

  /* initialize a pam matrix */
  strncpy(ppst->pamfile,pgm_def.smstr,MAX_FN);
  standard_pam(ppst->pamfile,ppst,del_set,gap_set);
  ppst->have_pam2 = 0;

  /* this is always protein by default */
  ppst->nsq = naa;
  ppst->nsqx = naax;
  for (i=0; i<=ppst->nsqx; i++) {
    ppst->sq[i] = aa[i];
    ppst->hsq[i] = haa[i];
    ppst->sqx[i]=aax[i];	/* sq = aa */
    ppst->hsqx[i]=haax[i];	/* hsq = haa */
  }
  ppst->sq[ppst->nsqx+1] = ppst->sqx[ppst->nsqx+1] = '\0';

  /* set up the c_nt[] mapping */

#if defined(FASTS) || defined(FASTF) || defined(FASTM)
  ppst->c_nt[ESS] = ESS;
#endif
  ppst->c_nt[0]=0;
  for (i=1; i<=nnt; i++) {
    ppst->c_nt[i]=gc_nt[i];
    ppst->c_nt[i+nnt]=gc_nt[i]+nnt;
  }
}

/*
 * alloc_pam(): allocates memory for the 2D pam matrix as well
 * as for the integer array used to transmit the pam matrix
 */
void
alloc_pam (int d1, int d2, struct pstruct *ppst)
{
  int     i, *d2p;
  char err_str[128];

  if ((ppst->pam2[0] = (int **) malloc (d1 * sizeof (int *))) == NULL) {
     sprintf(err_str,"Cannot allocate 2D pam matrix: %d",d1);
     s_abort (err_str,"");
  }

  if ((ppst->pam2[1] = (int **) malloc (d1 * sizeof (int *))) == NULL) {
     sprintf(err_str,"Cannot allocate 2D pam matrix: %d",d1);
     s_abort (err_str,"");
  }

  if ((d2p = pam12 = (int *) calloc (d1 * d2, sizeof (int))) == NULL) {
     sprintf(err_str,"Cannot allocate 2D pam matrix: %d",d1);
     s_abort (err_str,"");
   }

   for (i = 0; i < d1; i++, d2p += d2)
      ppst->pam2[0][i] = d2p;

   if ((d2p=pam12x= (int *) malloc (d1 * d2 * sizeof (int))) == NULL) {
     sprintf(err_str,"Cannot allocate 2d pam matrix: %d",d2);
     s_abort (err_str,"");
   }

   for (i = 0;  i < d1; i++, d2p += d2)
      ppst->pam2[1][i] = d2p;
}

/*
 *  init_pam2(struct pstruct pst): Converts 1-D pam matrix to 2-D
 */
void
init_pam2 (struct pstruct *ppst) {
  int     i, j, k, nsq;

  nsq = ppst->nsq;

  ppst->pam2[0][0][0] = -BIGNUM;
  ppst->pam_h = -1; ppst->pam_l = 1;

  k = 0;
  for (i = 1; i <= nsq; i++) {
    ppst->pam2[0][0][i] = ppst->pam2[0][i][0] = -BIGNUM;
    for (j = 1; j <= i; j++) {
      ppst->pam2[0][j][i] = ppst->pam2[0][i][j] = pam[k++] - ppst->pamoff;
      if (ppst->pam_l > ppst->pam2[0][i][j]) ppst->pam_l =ppst->pam2[0][i][j];
      if (ppst->pam_h < ppst->pam2[0][i][j]) ppst->pam_h =ppst->pam2[0][i][j];
    }
  }
}

void
init_pamx (struct pstruct *ppst) {
  int     i, j, k, nsq, pam_xx, pam_xm;
  int sa_x, sa_t, tmp;

  nsq = ppst->nsq;

  ppst->nt_align = (ppst->dnaseq== SEQT_DNA || ppst->dnaseq == SEQT_RNA);

  if (ppst->nt_align) {
    sa_x = pascii['N'];
    sa_t = sa_x;
  }
  else {
    sa_x = pascii['X'];
    sa_t = pascii['*'];
  }

  if (ppst->dnaseq == SEQT_RNA) {
    tmp = ppst->pam2[0][nascii['G']][nascii['G']] - 1;
    ppst->pam2[0][nascii['A']][nascii['G']] = 
      ppst->pam2[0][nascii['C']][nascii['T']] = 
      ppst->pam2[0][nascii['C']][nascii['U']] = tmp;
  }

  if (ppst->pam_x_set) {
    for (i=1; i<=nsq; i++) {
      ppst->pam2[0][sa_x][i] = ppst->pam2[0][i][sa_x]=ppst->pam_xm;
      ppst->pam2[0][sa_t][i] = ppst->pam2[0][i][sa_t]=ppst->pam_xm;
    }
    ppst->pam2[0][sa_x][sa_x]=ppst->pam_xx;
    ppst->pam2[0][sa_t][sa_t]=ppst->pam_xx;
  }
  else {
    ppst->pam_xx = ppst->pam2[0][sa_x][sa_x];
    ppst->pam_xm = ppst->pam2[0][1][sa_x];
  }

  pam_xx = ppst->pam_xx;
  pam_xm = ppst->pam_xm;

  if (ppst->ext_sq_set) {	/* using extended alphabet */
    /* fill in pam2[1] matrix */
    ppst->pam2[1][0][0] = -BIGNUM;
    /* fill in additional parts of the matrix */
    for (i = 1; i <= nsq; i++) {

      /* -BIGNUM to all matches vs 0 */
      ppst->pam2[0][0][i+nsq] = ppst->pam2[0][i+nsq][0] = 
	ppst->pam2[1][0][i+nsq] = ppst->pam2[1][i+nsq][0] = 
	ppst->pam2[1][0][i] = ppst->pam2[1][i][0] = -BIGNUM;

      for (j = 1; j <= nsq; j++) {

	/* replicate pam2[0] to i+nsq, j+nsq */
	ppst->pam2[0][i+nsq][j] = ppst->pam2[0][i][j+nsq] =
	  ppst->pam2[0][i+nsq][j+nsq] = ppst->pam2[1][i][j] =
	  ppst->pam2[0][i][j];

	/* set the high portion of pam2[1] to the corresponding value
	   of pam2[1][sa_x][j] */

	ppst->pam2[1][i+nsq][j] = ppst->pam2[1][i][j+nsq]=
	  ppst->pam2[1][i+nsq][j+nsq]=ppst->pam2[0][sa_x][j];
      }
    }
  }
}

/*  function specific initializations */
void
f_initenv (struct mngmsg *m_msp, struct pstruct *ppst, unsigned char **aa0) {
  struct msg_def_str m_msg_def;
  int pgm_id;

  pgm_id=ppst->pgm_id;
  m_msg_def = msg_def_arr[pgm_id];

  m_msp->last_calc_flg=0;

  strncpy(m_msp->f_id0,m_msg_def.f_id0,sizeof(m_msp->f_id0));
  strncpy(m_msp->f_id1,m_msg_def.f_id1,sizeof(m_msp->f_id1));
  strncpy (m_msp->label, m_msg_def.label, sizeof(m_msp->label));

#ifndef SSEARCH
  strncpy (m_msp->alab[0],"initn",20);
  strncpy (m_msp->alab[1],"init1",20);
  strncpy (m_msp->alab[2],"opt",20);
#else
  strncpy (m_msp->alab[0],"s-w opt",20);
#endif

  ppst->gdelval += pgm_def_arr[pgm_id].g_open_mod;
  ppst->sw_flag = m_msg_def.sw_flag;
  m_msp->e_cut=pgm_def_arr[pgm_id].e_cut;

  ppst->score_ix = 0;
  ppst->histint = 2;
  m_msp->qframe = m_msg_def.qframe;
  ppst->sw_flag = m_msg_def.sw_flag;
  m_msp->nframe = m_msg_def.nframe;
  m_msp->nrelv = m_msg_def.nrelv;
  m_msp->srelv = m_msg_def.srelv;
  m_msp->arelv = m_msg_def.arelv;
  m_msp->stages = m_msg_def.stages;
#if defined(PRSS)
  m_msp->shuff_wid = 0;
  m_msp->shuff_max = 200;
#endif

  /* see param.h for the definition of all these */

  m_msp->qshuffle = 0;
  m_msp->nm0 = 1;
  m_msp->escore_flg = 0;

  /* pam information */
  ppst->pam_pssm = 0;
#if defined(FASTS) || defined(FASTF) || defined(FASTM)
   ppst->pam_xx = ppst->pam_xm = 0;
#else
  ppst->pam_xx = 1;  /* set >0 to use pam['X']['X'] value */
  ppst->pam_xm = -1;  /* set >0 to use pam['X']['A-Z'] value */
#endif
  ppst->pam_x_set = 0;
  ppst->pam_set = 0;
  ppst->pam_pssm = 0;
  ppst->p_d_set = 0;
  ppst->pamoff = 0;
  ppst->ext_sq_set = 0;

  if (pgm_def_arr[ppst->pgm_id].ktup > 0) {
    mktup = 2;
    ppst->param_u.fa.bestscale = 300;
    ppst->param_u.fa.bestoff = 36;
    ppst->param_u.fa.bkfact = 6;
    ppst->param_u.fa.scfact = 3;
    ppst->param_u.fa.bktup = 2;
    ppst->param_u.fa.ktup = 0;
    ppst->param_u.fa.bestmax = 50;
    ppst->param_u.fa.pamfact = 1;
    ppst->param_u.fa.altflag = 0;
    ppst->param_u.fa.optflag = 1;
    ppst->param_u.fa.iniflag = 0;
    ppst->param_u.fa.optcut = 0;
    ppst->param_u.fa.optcut_set = 0;
    ppst->param_u.fa.cgap = 0;
    ppst->param_u.fa.optwid = MAXWINDOW;
  }

}

/*  switches for fasta only */

static int shift_set=0;
static int subs_set=0;
static int sw_flag_set=0;
static int nframe_set=0;
static int wid_set=0;

void
f_getopt (char copt, char *optarg,
	  struct mngmsg *m_msg, struct pstruct *ppst)
{
  int pgm_id;
  char *bp;

  pgm_id = ppst->pgm_id;

  switch (copt) {
  case '1': 
    if (pgm_def_arr[pgm_id].ktup > 0) {
      ppst->param_u.fa.iniflag=1;
    }
     break;
  case '3':
    nframe_set = 1;
    if (pgm_id == TFA_PID) {
      m_msg->nframe = 3; break;
    }
    else {
      m_msg->nframe = 1;	/* for TFASTXY */
      m_msg->qframe = 1;  /* for FASTA, FASTX */
    }
    break;
  case 'A':
    ppst->sw_flag= 1;
    sw_flag_set = 1;
    break;
  case 'c':
    if (pgm_def_arr[pgm_id].ktup > 0) {
      sscanf (optarg, "%d", &ppst->param_u.fa.optcut);
      ppst->param_u.fa.optcut_set = 1;
    }
    break;
  case 'f':
    sscanf (optarg, "%d", &ppst->gdelval);
    if (ppst->gdelval > 0) ppst->gdelval = -ppst->gdelval;
    del_set = 1;
    break;
  case 'g':
    sscanf (optarg, "%d", &ppst->ggapval);
    if (ppst->ggapval > 0) ppst->ggapval = -ppst->ggapval;
    gap_set = 1;
    break;
  case 'h':
    sscanf (optarg, "%d", &ppst->gshift);
    if (ppst->gshift > 0) ppst->gshift = -ppst->gshift;
    shift_set = 1;
    break;
  case 'j':
    sscanf (optarg, "%d", &ppst->gsubs);
    subs_set = 1;
    break;
  case 'k':
    sscanf (optarg, "%d", &m_msg->shuff_max);
    mshuff_set = 1;
    break;
  case 'o':
    if (pgm_def_arr[pgm_id].ktup > 0) {
      ppst->param_u.fa.optflag = 0;
      msg_def_arr[pgm_id].nrelv = m_msg->nrelv = 2;
    }
    break;
  case 'P':
    strncpy(ppst->pgpfile,optarg,MAX_FN);
    if ((bp=strchr(ppst->pgpfile,' '))!=NULL) {
      *bp='\0';
      ppst->pgpfile_type = atoi(bp+1);
    }
    else ppst->pgpfile_type = 0;
    ppst->pgpfile[MAX_FN-1]='\0';
    ppst->pam_pssm = 1;
    break;
  case 'r':
    sscanf(optarg,"%d/%d",&ppst->p_d_mat,&ppst->p_d_mis);
    if (ppst->p_d_mat > 0 && ppst->p_d_mis < 0) {
      ppst->p_d_set = 1;
      strncpy(ppst->pamfile,optarg,40);
    }
    break;
  case 's':
    strncpy (ppst->pamfile, optarg, 120);
    ppst->pamfile[120-1]='\0';
    if (!standard_pam(ppst->pamfile,ppst,del_set, gap_set)) {
      initpam (ppst->pamfile, ppst);
    }
    ppst->pam_set=1;
    break;
  case 'S':	/* turn on extended alphabet for seg */
    ppst->ext_sq_set = 1;
    break;
  case 't':
    if (tolower(optarg[0])=='t') {
      m_msg->term_code = aascii['*']; optarg++;
    }
    if (*optarg) {sscanf (optarg, "%d", &ppst->tr_type);}
    break;
  case 'x':
    if (strchr(optarg,',')!=NULL) {
      sscanf (optarg,"%d,%d",&ppst->pam_xx, &ppst->pam_xm);
    }
    else {
      sscanf (optarg,"%d",&ppst->pam_xx);
      ppst->pam_xm = ppst->pam_xx;
    }
    ppst->pam_x_set=1;
    break;
  case 'y':
    if (pgm_def_arr[pgm_id].ktup > 0) {
      sscanf (optarg, "%d", &ppst->param_u.fa.optwid);
      wid_set = 1;
    }
    break;
  }
}

void
f_lastenv (struct mngmsg *m_msg, struct pstruct *ppst)
{
  char save_str[MAX_SSTR];

#if !defined(FASTM) && !defined(FASTS) && !defined(FASTF)
  strncpy(save_str,"*",sizeof(save_str));
#else
  strncpy(save_str,",",sizeof(save_str));
#endif

  if (m_msg->qdnaseq == SEQT_UNK) {
    build_xascii(qascii,save_str);
    if (m_msg->ann_flg) ann_ascii(qascii,m_msg->ann_arr);
  }  

/* this check allows lc DNA sequence queries with FASTX */
#if defined(FASTA) && !defined(FASTS) && !defined(FASTM) && !defined(FASTF)
  else
   init_ascii(ppst->ext_sq_set,qascii,m_msg->qdnaseq);
#endif
}

void
f_getarg (int argc, char **argv, int optind,
	  struct mngmsg *m_msg, struct pstruct *ppst)
{

  if (pgm_def_arr[ppst->pgm_id].ktup > 0) {
    if (argc - optind >= 4) {
      sscanf (argv[optind + 3], "%d", &ppst->param_u.fa.ktup);
      ktup_set = 1;
    }
    else
      ppst->param_u.fa.ktup = -ppst->param_u.fa.bktup;
  }
  
  if (ppst->pgm_id == RSS_PID && argc - optind > 3) {
    sscanf (argv[optind + 3], "%d", &m_msg->shuff_max);
  }

  if (ppst->pgm_id == RFX_PID && argc - optind > 4) {
    sscanf (argv[optind + 4], "%d", &m_msg->shuff_max);
  }
}

/* fills in the query ascii mapping from the parameter
   ascii mapping.
*/

void
re_ascii(int *qascii, int *pascii) {
  int i;

  for (i=0; i < 128; i++) {
    if (qascii[i] > '@' || qascii[i] < ESS) {
      qascii[i] = pascii[i];
    }
  }
}


/* recode has become function specific to accommodate FASTS/M */
/* modified 28-Dec-2004 to ensure that all mapped characters
   are valid */
int
recode(unsigned char *seq, int n, int *qascii, int nsqx) {
  int i,j;
  char save_c;

#if defined(FASTS) || defined(FASTM)
  qascii[',']=ESS;
#endif

  for (i=0; i < n; i++) {
    save_c = seq[i];
    if (seq[i] > '@') seq[i] = qascii[seq[i]];
    if (seq[i] > nsqx && seq[i]!=ESS) {
      fprintf(stderr, "*** Warning - unrecognized residue at %d:%c - %2d\n",
	      i,save_c,save_c);
      seq[i] = qascii['X'];
    }
  }
  seq[i]=EOSEQ;
  return i;
}

/* here we have the query sequence, all the command line options,
   but we need to set various parameter options based on the type
   of the query sequence (m_msg->qdnaseq = 0:protein/1:DNA) and
   the function (FASTA/FASTX/TFASTA)
*/

/* this resetp is for conventional a FASTA/TFASTXYZ search */
void
resetp (struct mngmsg *m_msg, struct pstruct *ppst) {
  int i, pgm_id;

  pgm_id = ppst->pgm_id;

#if defined(TFAST)
  if (m_msg->qdnaseq == SEQT_DNA || m_msg->qdnaseq == SEQT_RNA) {
    fprintf(stderr," %s compares a protein to a translated\n\
DNA sequence library.  Do not use a DNA query/scoring matrix.\n",prog_func);
    exit(1);
  }
#else
#if (defined(FASTX) || defined(FASTY))
  if (!(m_msg->qdnaseq == SEQT_DNA || m_msg->qdnaseq == SEQT_RNA)) {
    fprintf(stderr," FASTX/Y compares a DNA sequence to a protein database\n");
    fprintf(stderr," Use a DNA query\n");
    exit(1);
  }
#endif
#endif

  if (msg_def_arr[pgm_id].q_seqt == SEQT_UNK){
    if (m_msg->qdnaseq == SEQT_DNA || m_msg->qdnaseq == SEQT_RNA) {
      msg_def_arr[pgm_id].q_seqt = m_msg->qdnaseq;
      msg_def_arr[pgm_id].p_seqt = SEQT_DNA;
      msg_def_arr[pgm_id].l_seqt = SEQT_DNA;
      if (m_msg->qdnaseq == SEQT_DNA) msg_def_arr[pgm_id].qframe = 2;
    }
    else {
      msg_def_arr[pgm_id].q_seqt = SEQT_PROT;
    }
  }

  ppst->dnaseq = msg_def_arr[pgm_id].p_seqt;
  if (!sw_flag_set) ppst->sw_flag = msg_def_arr[pgm_id].sw_flag;
  if (!m_msg->e_cut_set) m_msg->e_cut=pgm_def_arr[pgm_id].e_cut;

  if (ppst->dnaseq == SEQT_DNA && m_msg->qdnaseq==SEQT_RNA) {
    ppst->dnaseq = SEQT_RNA;
    ppst->nt_align = 1;
  }
  if (ppst->dnaseq==SEQT_DNA) pascii = &nascii[0];
  else if (ppst->dnaseq==SEQT_RNA) {
    pascii = &nascii[0];
    ppst->sq[nascii['T']] = 'U';
  }
  else pascii = &aascii[0];
  m_msg->ldnaseq = msg_def_arr[pgm_id].l_seqt;
  if (m_msg->ldnaseq & SEQT_DNA) {
    memcpy(lascii,nascii,sizeof(lascii));
#ifndef TFAST
#ifdef DNALIB_LC
   init_ascii(ppst->ext_sq_set,lascii,m_msg->ldnaseq);
#endif
#else
  /* no init_ascii() because we translate lower case library sequences */
#endif
  }
  else {
    memcpy(lascii,aascii,sizeof(lascii));	/* initialize lib mapping */

#if defined(FASTF) || defined(FASTS) || defined(FASTM)
    lascii['*'] = NA;
#endif
    init_ascii(ppst->ext_sq_set,lascii,m_msg->ldnaseq);
  }

  if (!nframe_set) {
    m_msg->qframe = msg_def_arr[pgm_id].qframe;
    m_msg->nframe = msg_def_arr[pgm_id].nframe;
  }

  /* the possibilities:
  	     -i  -3	qframe	revcomp
   FA_D/FX   -    -        2       0	
   FA_D/FX   +    -        2       1	
   FA_D/FX   -    +        1       0	
   FA_D/FX   +    +        2       1	
  */

  if (m_msg->qdnaseq == SEQT_DNA) {
    m_msg->nframe = 1;
    if (m_msg->qframe == 1 && m_msg->revcomp==1) {
      m_msg->qframe = m_msg->revcomp+1;
    }
  }
  else if (m_msg->qdnaseq == SEQT_RNA) {
    m_msg->qframe = m_msg->revcomp+1;
    m_msg->nframe = 1;
  }

  /* change settings for DNA search */
  if (ppst->dnaseq == SEQT_DNA || ppst->dnaseq == SEQT_RNA) {
    ppst->histint = 4;

    if (!del_set) {
#ifdef OLD_FASTA_GAP
      ppst->gdelval = -16;	/* def. del penalty */
#else
      ppst->gdelval = -12;	/* def. open penalty */
#endif
    }
    if (!gap_set) ppst->ggapval = -4;	/* def. gap penalty */

    if (pgm_def_arr[pgm_id].ktup > 0) {
      /* these parameters are used to scale optcut, they should be replaced
	 by statistically based parameters */
      if (!wid_set) ppst->param_u.fa.optwid = 16;
      ppst->param_u.fa.bestscale = 80;
      ppst->param_u.fa.bkfact = 5;
      ppst->param_u.fa.scfact = 1;
      ppst->param_u.fa.bktup = 6;
      ppst->param_u.fa.bestmax = 80;
      ppst->param_u.fa.bestoff = 45;

      if (!sw_flag_set) {
	ppst->sw_flag = 0;
	strncpy(m_msg->f_id1,"bs",sizeof(m_msg->f_id1));
      }

      /* largest ktup */
      mktup = 6;
      
      if (ppst->param_u.fa.pamfact >= 0) ppst->param_u.fa.pamfact = 0;
      if (ppst->param_u.fa.ktup < 0)
	ppst->param_u.fa.ktup = -ppst->param_u.fa.bktup;
    }

    ppst->nsq = nnt;
    ppst->nsqx = nntx;
    for (i=0; i<=ppst->nsqx; i++) {
      ppst->hsq[i] = hnt[i];
      ppst->sq[i] = nt[i];
      ppst->hsqx[i] = hntx[i];
      ppst->sqx[i] = ntx[i];
    }
    ppst->sq[ppst->nsqx+1] = ppst->sqx[ppst->nsqx+1] = '\0';

    if (!ppst->pam_set) {
      if (ppst->p_d_set)
	mk_n_pam(npam,nnt,ppst->p_d_mat,ppst->p_d_mis);
#if !defined(FASTS) && !defined(FASTM)
      else if (ppst->pamfile[0]=='\0' || strncmp(ppst->pamfile,"BL50",4)==0) {
	strncpy (ppst->pamfile, "+5/-4", sizeof(ppst->pamfile));
      }
#else
      else if (strncmp(ppst->pamfile,"MD20",4)==0) {
	strncpy (ppst->pamfile, "+2/-2", sizeof(ppst->pamfile));
	ppst->p_d_mat = +2;
	ppst->p_d_mis = -2;
      	mk_n_pam(npam,nnt,ppst->p_d_mat,ppst->p_d_mis);
      }
#endif
      pam = npam;
    }

    strncpy (m_msg->sqnam, "nt",sizeof(m_msg->sqnam));
    strncpy (m_msg->sqtype, "DNA",sizeof(m_msg->sqtype));
  }	/* end DNA reset */

  else {  /* other parameters for protein comparison */
    if (pgm_def_arr[pgm_id].ktup > 0) {
      if (!wid_set) {
	if (ppst->param_u.fa.ktup==1) ppst->param_u.fa.optwid = 32;
	else ppst->param_u.fa.optwid = 16;
      }
    }
    if (!del_set) {ppst->gdelval += pgm_def_arr[pgm_id].g_open_mod;}
    if (!shift_set) {ppst->gshift = pgm_def_arr[pgm_id].gshift;}
    if (!subs_set) {ppst->gsubs = pgm_def_arr[pgm_id].hshift;}
  }

}

/* query_parm()	this function asks for any additional parameters
	that have not been provided.  Could be null. */
void
query_parm (struct mngmsg *m_msp, struct pstruct *ppst)
{
   char    qline[40];

   if (pgm_def_arr[ppst->pgm_id].ktup > 0) {
     if (ppst->param_u.fa.ktup < 0)
       ppst->param_u.fa.ktup = -ppst->param_u.fa.ktup;

     if (ppst->param_u.fa.ktup == 0) {
       printf (" ktup? (1 to %d) [%d] ", mktup, ppst->param_u.fa.bktup);
       if (fgets (qline, sizeof(qline), stdin) == NULL) exit (0);
       else sscanf(qline,"%d",&ppst->param_u.fa.ktup);
     }
     if (ppst->param_u.fa.ktup == 0)
       ppst->param_u.fa.ktup = ppst->param_u.fa.bktup;
     else ktup_set = 1;
   }

#if defined(PRSS)
   if (m_msp->shuff_max < 10) m_msp->shuff_max = 200;

   if (!mshuff_set) {
     printf(" number of shuffles [%d]? ",m_msp->shuff_max);
     fflush(stdout);
     if (fgets (qline, sizeof(qline), stdin) == NULL) exit (0);
     else sscanf(qline,"%d",&m_msp->shuff_max);
   }

   if (ppst->zs_win == 0) {
     printf (" local (window) (w) or uniform (u) shuffle [u]? ");
     if (fgets (qline, sizeof(qline), stdin) == NULL) exit (0);
     else if (qline[0]=='w' || qline[0]=='W') {
       m_msp->shuff_wid = 20;
       printf(" local shuffle window size [%d]? ",m_msp->shuff_wid);
       if (fgets (qline, sizeof(qline), stdin) == NULL) exit (0);
       else sscanf(qline,"%d",&m_msp->shuff_wid);
     }
   }
#endif
}

/* last_init() cannot look at aa0, n0, because it is only run once,
   it is not run before each new aa0 search */
void
last_init (struct mngmsg *m_msg, struct pstruct *ppst
#ifdef PCOMPLIB
	   ,int nnodes
#endif
	   )
{
  int ix_l, ix_i, i, pgm_id;
  double *kar_p;
  double aa0_f[MAXSQ];

  pgm_id = ppst->pgm_id;

#if defined(FASTF) || defined(FASTS) || defined(FASTM)
  m_msg->nohist = 1;
  m_msg->shuff_max = 2000;
#ifndef PCOMPLIB
  ppst->shuff_node = m_msg->shuff_max/max_workers;
#else
  ppst->shuff_node = m_msg->shuff_max/nnodes;
#endif
#endif

  if (m_msg->aln.llen < 1) {
    m_msg->aln.llen = 60;
  }

#ifndef PCOMPLIB
#if defined(FASTX) || defined(FASTY) || defined(TFAST)
  /* set up translation tables: faatran.c */
  aainit(ppst->tr_type,ppst->debug_lib);
#endif
#endif

/* a sanity check */
#if !defined(TFAST)
   if (m_msg->revcomp && m_msg->qdnaseq!=SEQT_DNA && m_msg->qdnaseq!=SEQT_RNA) {
     fprintf(stderr," cannot reverse complement protein\n");
     m_msg->revcomp = 0;
   }
#endif

   if (pgm_def_arr[pgm_id].ktup > 0) {

     if (ppst->param_u.fa.ktup < 0)
       ppst->param_u.fa.ktup = -ppst->param_u.fa.ktup;

     if (ppst->param_u.fa.ktup < 1 || ppst->param_u.fa.ktup > mktup) {
       fprintf(stderr," warning ktup = %d out of range [1..%d], reset to %d\n",
	       ppst->param_u.fa.ktup, mktup, ppst->param_u.fa.bktup);
       ppst->param_u.fa.ktup = ppst->param_u.fa.bktup;
     }
   }

   if (pgm_id == TFA_PID) {
     m_msg->revcomp *= 3;
     if (m_msg->nframe == 3) m_msg->nframe += m_msg->revcomp;
   }
   else if (pgm_id == TFX_PID || pgm_id == TFY_PID) {
     if (m_msg->nframe == 1) m_msg->nframe += m_msg->revcomp;
   }

#if !defined(TFAST)
  /* for fasta/fastx searches, itt iterates the the query strand */
  m_msg->nitt1 = m_msg->qframe-1;
#else
  /* for tfasta/tfastxy searches, itt iterates the library frames */
  m_msg->nitt1 = m_msg->nframe-1;
#endif

  if (pgm_def_arr[pgm_id].ktup > 0) {
    if (ppst->param_u.fa.ktup>=2 && !wid_set) {
      ppst->param_u.fa.optwid=16;
      switch (pgm_id) {
      case FA_PID:
	m_msg->thr_fact = 16;
	break;
      case FX_PID:
      case FY_PID:
	m_msg->thr_fact = 8;
	break;
      case TFA_PID:
      case TFX_PID:
      case TFY_PID:
	m_msg->thr_fact = 4;
	break;
      default:
	m_msg->thr_fact = 2;
      }
    }
    else { m_msg->thr_fact = 4;}
  }
  else m_msg->thr_fact = 1;

#if defined(PRSS)
   if (m_msg->shuff_max < 10) m_msg->shuff_max = 200;
   if (ppst->zsflag < 10) ppst->zsflag += 10;
   if (ppst->zs_win > 0) {
     m_msg->shuff_wid = ppst->zs_win;
   }
#endif

   if (pgm_def_arr[ppst->pgm_id].ktup > 0) {
     if (ppst->param_u.fa.iniflag) {
       ppst->score_ix = 1;
       strncpy (m_msg->label, "initn init1", sizeof(m_msg->label));
     }
     else if (ppst->param_u.fa.optflag) {
       ppst->score_ix = 2;
       m_msg->stages = 1;
     }
   }

   if (!ppst->have_pam2) {
     alloc_pam (MAXSQ, MAXSQ, ppst);
     init_pam2(ppst);
   }
   init_pamx(ppst);

   if (ppst->pam_ms) {
     if (m_msg->qdnaseq == SEQT_PROT) {
       /* code to make 'L'/'I' identical scores */
       ix_l = pascii['L'];
       ix_i = pascii['I'];
       ppst->pam2[0][ix_l][ix_i] = ppst->pam2[0][ix_i][ix_l] =
	 ppst->pam2[0][ix_l][ix_l] = ppst->pam2[0][ix_i][ix_i] =
	 (ppst->pam2[0][ix_l][ix_l]+ppst->pam2[0][ix_i][ix_i]+1)/2;
       for (i=1; i<=ppst->nsq; i++) {
	 ppst->pam2[0][i][ix_i] = ppst->pam2[0][i][ix_l] =
	   (ppst->pam2[0][i][ix_l]+ppst->pam2[0][i][ix_i]+1)/2;
	 ppst->pam2[0][ix_i][i] = ppst->pam2[0][ix_l][i] =
	   (ppst->pam2[0][ix_i][i]+ppst->pam2[0][ix_l][i]+1)/2;
       }

       /* code to make 'Q'/'K' identical scores */
       if (!shift_set) {
	 ix_l = pascii['Q'];
	 ix_i = pascii['K'];
	 ppst->pam2[0][ix_l][ix_i] = ppst->pam2[0][ix_i][ix_l] =
	   ppst->pam2[0][ix_l][ix_l] = ppst->pam2[0][ix_i][ix_i] =
	   (ppst->pam2[0][ix_l][ix_l]+ppst->pam2[0][ix_i][ix_i]+1)/2;
	 for (i=1; i<=ppst->nsq; i++) {
	   ppst->pam2[0][i][ix_i] = ppst->pam2[0][i][ix_l] =
	     (ppst->pam2[0][i][ix_l]+ppst->pam2[0][i][ix_i]+1)/2;
	   ppst->pam2[0][ix_i][i] = ppst->pam2[0][ix_l][i] =
	     (ppst->pam2[0][ix_i][i]+ppst->pam2[0][ix_l][i]+1)/2;
	 }
       }
     }
   }

   /* once we have a complete pam matrix, we can calculate Lambda and K 
      for "average" sequences */
   kar_p = NULL;
   init_karlin_a(ppst, aa0_f, &kar_p);
   do_karlin_a(ppst->pam2[0], ppst, aa0_f,
	       kar_p, &m_msg->Lambda, &m_msg->K, &m_msg->H);
   free(kar_p);

#if defined(FASTF) || defined(FASTS) || defined(FASTM)
   if (ppst->ext_sq_set) {
     fprintf(stderr," -S not available on [t]fast[fs]\n");
     ppst->ext_sq_set = 0;

     /* reset sascii to ignore -S, map lc */
     init_ascii(0,lascii,0);
   }
#endif
}

/* this function is left over from the older FASTA format scoring
   matrices that allowed additional parameters (bktup, bkfact) to be
   set in the scoring matrix.  It is no longer used.  A modern version
   would set parameters based on lambda and K.
*/
/*
void
f_initpam (line, ppst)
char   *line;
struct pstruct *ppst;
{
   if (sscanf (line, " %d %d %d %d %d %d %d", &ppst->param_u.fa.scfact,
	       &ppst->param_u.fa.bestoff, &ppst->param_u.fa.bestscale,
	       &ppst->param_u.fa.bkfact, &ppst->param_u.fa.bktup,
	       &ppst->param_u.fa.bestmax, &ppst->histint) != 7)
   {
      printf ("  bestcut parameters - bad format\n");
      exit (1);
   }
}
*/

/* alloc_pam2 creates a profile structure */
int **
alloc_pam2p(int len, int nsq) {
  int i;
  int **pam2p;

  if ((pam2p = (int **)calloc(len,sizeof(int *)))==NULL) {
    fprintf(stderr," Cannot allocate pam2p: %d\n",len);
    return NULL;
  }

  if((pam2p[0] = (int *)calloc((nsq+1)*len,sizeof(int)))==NULL) {
    fprintf(stderr, "Cannot allocate pam2p[0]: %d\n", (nsq+1)*len);
    free(pam2p);
    return NULL;
  }

  for (i=1; i<len; i++) {
    pam2p[i] = pam2p[0] + (i*(nsq+1));
  }

  return pam2p;
}

void free_pam2p(int **pam2p) {
  if (pam2p) {
    free(pam2p[0]);
    free(pam2p);
  }
}

/* sortbest has now become comparison function specific so that we can use
   a different comparison for fasts/f 
*/
#if !defined(FASTS) && !defined (FASTF) && !defined(FASTM)
#ifndef PCOMPLIB
void
qshuffle() {}
#endif

int
last_calc(unsigned char *aa0, unsigned char *aa1, int maxn,
	  struct beststr **bestp_arr, int nbest,
	  struct mngmsg *m_msg, struct pstruct *pst,
	  void **f_str, void *rs_str)
{
  return nbest;
}

void sortbest (bptr, nbest, irelv)
struct beststr **bptr;
int nbest, irelv;
{
    int gap, i, j;
    struct beststr *tmp;

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

void show_aux(FILE *fp, struct beststr *bptr) {}
void header_aux(FILE *fp) {}

#else
void sortbest (bptr, nbest, irelv)
struct beststr **bptr;
int nbest, irelv;
{
    int gap, i, j;
    struct beststr *tmp;

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

#if defined(FASTS) || defined(FASTM)

#ifndef PCOMPLIB
/* this shuffle is for FASTS  */
/* convert ',' -> '\0', shuffle each of the substrings */
void
qshuffle(unsigned char *aa0, int n0, int nm0) {

  unsigned char **aa0start, *aap, tmp;
  int i,j,k, ns;

  if ((aa0start=(unsigned char **)calloc(nm0+1,
					 sizeof(unsigned char *)))==NULL) {
    fprintf(stderr,"cannot calloc for qshuffle %d\n",nm0);
    exit(1);
  }

  aa0start[0]=aa0;
  for (k=1,i=0; i<n0; i++) {
    if (aa0[i]==EOSEQ || aa0[i]==ESS) {
      aa0[i]='\0';
      aa0start[k++] = &aa0[i+1];
    }
  }  

  /* aa0start has the beginning of each substring */
  for (k=0; k<nm0; k++) {
    aap=aa0start[k];
    ns = strlen((char *)aap);
    for (i=ns; i>1; i--) {
      j = nrand(i);
      tmp = aap[j];
      aap[j] = aap[i-1];
      aap[i-1] = tmp;
    }
    aap[ns] = 0;
  }

  for (k=1; k<nm0; k++) {
/*  aap = aa0start[k];
    while (*aap) fputc(pst.sq[*aap++],stderr);
    fputc('\n',stderr);
*/
    aa0start[k][-1]=ESS;
  }

  free(aa0start);
}
#endif
#endif

#ifdef FASTF
#ifndef PCOMPLIB
void qshuffle(unsigned char *aa0, int n0, int nm0) {

  int i, j, k, nmpos;
  unsigned char tmp;
  int nmoff;
  
  nmoff = (n0 - nm0 - 1)/nm0 + 1;

  for (i = nmoff-1 ; i > 0 ; i--) {

    /* j = nrand(i); if (i == j) continue;*/       /* shuffle columns */ 
    j = (nmoff -1 ) - i; 
    if (i <= j) break; /* reverse columns */

    /* swap all i'th column residues for all j'th column residues */
    for(nmpos = 0, k = 0 ; k < nm0 ; k++, nmpos += nmoff+1 ) {
      tmp = aa0[nmpos + i];
      aa0[nmpos + i] = aa0[nmpos + j];
      aa0[nmpos + j] = tmp;
    }
  }
}
#endif
#endif


/* show additional best_str values */
void show_aux(FILE *fp, struct beststr *bptr) {
  fprintf(fp," %2d %3d",bptr->segnum,bptr->seglen);
}

void header_aux(FILE *fp) {
  fprintf(fp, " sn  sl");
}
#endif

void
fill_pam(int **pam2p, int n0, int nsq, double **freq2d, double scale) {
  int i, j;
  double freq;

  /* fprintf(stderr, "scale: %g\n", scale); */
  
  /* now fill in the pam matrix: */
  for (i = 0 ; i < n0 ; i++) {
    for (j = 1 ; j <=20 ; j++) {
      freq = scale * freq2d[i][j-1];
      if ( freq < 0.0) freq -= 0.5;
      else freq += 0.5;
      pam2p[i][j] = (int)(freq);
    }
  }
}

double
get_lambda(int **pam2p, int n0, int nsq, unsigned char *aa0) {
  double lambda, H;
  double *pr, tot, sum;
  int i, ioff, j, min, max;

  /* get min and max scores */
  min = BIGNUM;
  max = -BIGNUM;
  if(pam2p[0][1] == -BIGNUM) {
    ioff = 1;
    n0++;
  } else {
    ioff = 0;
  }

  for (i = ioff ; i < n0 ; i++) {
    for (j = 1; j <= nsq ; j++) {
      if (min > pam2p[i][j])
	min = pam2p[i][j];
      if (max < pam2p[i][j])
	max = pam2p[i][j];
    }
  }

  /*  fprintf(stderr, "min: %d\tmax:%d\n", min, max); */
  
  if ((pr = (double *) calloc(max - min + 1, sizeof(double))) == NULL) {
    fprintf(stderr, "Couldn't allocate memory for score probabilities: %d\n", max - min + 1);
    exit(1);
  }

  tot = (double) rrtotal * (double) rrtotal * (double) n0;
  for (i = ioff ; i < n0 ; i++) {
    for (j = 1; j <= nsq ; j++) {
      pr[pam2p[i][j] - min] +=
	(double) ((double) rrcounts[aascii[aa0[i]]] * (double) rrcounts[j]) / tot;
    }
  }

  sum = 0.0;
  for(i = 0 ; i <= max-min ; i++) { 
    sum += pr[i];
    /*    fprintf(stderr, "%3d: %g %g\n", i+min, pr[i], sum); */
  }
  /*  fprintf(stderr, "sum: %g\n", sum); */

  for(i = 0 ; i <= max-min ; i++) { pr[i] /= sum; }

  if (!karlin(min, max, pr, &lambda, &H)) {
    fprintf(stderr, "Karlin lambda estimation failed\n");
  }

  /*   fprintf(stderr, "lambda: %g\n", lambda); */
  free(pr);

  return lambda;
}

/*
   *aa0 - query sequence
   n0   - length
   pamscale - scaling for pam matrix - provided by apam.c, either
              0.346574 = ln(2)/2 (P120, BL62) or
	      0.231049 = ln(2)/3 (P250, BL50) 
*/

static unsigned char ustandard_aa[] ="\0ARNDCQEGHILKMFPSTWYV";

void
read_pssm(unsigned char *aa0, int n0, int nsq,
	  double pamscale, 
	  FILE *fp, int pgpf_type, struct pstruct *ppst) {
  int i, j, len, k;
  int qi, rj;	/* qi - index query; rj - index residues (1-20) */
  int **pam2p;
  int first, too_high;
  unsigned char *query, ctmp;
  char dline[512];
  double freq, **freq2d, lambda, new_lambda;
  double scale, scale_high, scale_low;

  pam2p = ppst->pam2p[0];

  if (pgpf_type == 0) {

    if(1 != fread(&len, sizeof(int), 1, fp)) {
      fprintf(stderr, "error reading from checkpoint file: %d\n", len);
      exit(1);
    }

    if(len != n0) {
      fprintf(stderr, "profile length (%d) and query length (%d) don't match!\n",
	      len,n0);
      exit(1);
    }

    /* read over query sequence stored in BLAST profile */
    if(NULL == (query = (unsigned char *) calloc(len+2, sizeof(char)))) {
      fprintf(stderr, "Couldn't allocate memory for query!\n");
      exit(1);
    }

    if(len != fread(query, sizeof(char), len, fp)) {
      fprintf(stderr, "Couldn't read query sequence from profile: %s\n", query);
      exit(1);
    }
  }
  else if (pgpf_type == 1) {

    if ((fgets(dline,sizeof(dline),fp) == NULL)  ||
	(1 != sscanf(dline, "%d",&len))) {
      fprintf(stderr, "error reading from checkpoint file: %d\n", len);
      exit(1);
    }

    if(len != n0) {
      fprintf(stderr, "profile length (%d) and query length (%d) don't match!\n",
	      len,n0);
      exit(1);
    }

    /* read over query sequence stored in BLAST profile */
    if(NULL == (query = (unsigned char *) calloc(len+2, sizeof(char)))) {
      fprintf(stderr, "Couldn't allocate memory for query!\n");
      exit(1);
    }

    if (fgets((char *)query,len+2,fp)==NULL) {
      fprintf(stderr, "Couldn't read query sequence from profile: %s\n", query);
      exit(1);
    }
  }  
  else {
    fprintf(stderr," Unrecognized PSSM file type: %d\n",pgpf_type);
    exit(1);
  }

  /* currently we don't do anything with query; ideally, we should
     check to see that it actually matches aa0 ... */

  /* quick 2d array alloc: */
  if((freq2d = (double **) calloc(n0, sizeof(double *))) == NULL) {
    fprintf(stderr, "Couldn't allocate memory for frequencies!\n");
    exit(1);
  }

  if((freq2d[0] = (double *) calloc(n0 * 20, sizeof(double))) == NULL) {
    fprintf(stderr, "Couldn't allocate memory for frequencies!\n");
    exit(1);
  }

  /* a little pointer arithmetic to fill out 2d array: */
  for (i = 1 ; i < n0 ; i++) {
    freq2d[i] = freq2d[i-1] + 20;
  }

  if (pgpf_type == 0) {
    for (qi = 0 ; qi < n0 ; qi++) {
      for (rj = 0 ; rj < 20 ; rj++) {
	if(1 != fread(&freq, sizeof(double), 1, fp)) {
	  fprintf(stderr, "Error while reading frequencies!\n");
	  exit(1);
	}
	freq2d[qi][rj] = freq;
      }
    }
  }
  else {
    for (qi = 0 ; qi < n0 ; qi++) {
      if ((fgets(dline,sizeof(dline),fp) ==NULL) ||
      (k = sscanf(dline,"%c %lg %lg %lg %lg %lg %lg %lg %lg %lg %lg %lg %lg %lg %lg %lg %lg %lg %lg %lg %lg\n",
		 &ctmp, &freq2d[qi][0], &freq2d[qi][1], &freq2d[qi][2], &freq2d[qi][3], &freq2d[qi][4], 
		 &freq2d[qi][5], &freq2d[qi][6], &freq2d[qi][7], &freq2d[qi][8], &freq2d[qi][9],
		 &freq2d[qi][10], &freq2d[qi][11], &freq2d[qi][12], &freq2d[qi][13], &freq2d[qi][14],
		      &freq2d[qi][15], &freq2d[qi][16], &freq2d[qi][17], &freq2d[qi][18], &freq2d[qi][19]))<1) {
	fprintf(stderr, "Error while reading frequencies: %d read!\n",k);
	exit(1);
      }
      for (rj=0; rj<20; rj++) { freq2d[qi][rj] /= 10.0; }	/* reverse scaling */
    }
  }

  for (qi = 0 ; qi < n0 ; qi++) {
    for (rj = 0 ; rj < 20 ; rj++) {
      if (freq2d[qi][rj] > 1e-20) {
	freq = log(freq2d[qi][rj] /((double) (rrcounts[rj+1])/(double) rrtotal));
	freq /= pamscale; /* this gets us close to originial pam scores */
	freq2d[qi][rj] = freq;
      }
      else {		
	/* when blastpgp decides to leave something out, it puts 0's in all the frequencies
	   in the binary checkpoint file.  In the ascii version, however, it uses BLOSUM62
	   values.  I will put in scoring matrix values as well */

	freq2d[qi][rj] = ppst->pam2[0][aa0[qi]][rj+1];
      }
    }
  }

  /* now figure out the right scale */
  scale = 1.0;
  lambda = get_lambda(ppst->pam2[0], 20, 20, ustandard_aa);

  /* should be near 1.0 because of our initial scaling by ppst->pamscale */
  /* fprintf(stderr, "real_lambda: %g\n", lambda); */

  /* get initial high/low scale values: */
  first = 1;
  while (1) {
    fill_pam(pam2p, n0, 20, freq2d, scale);
    new_lambda = get_lambda(pam2p, n0, 20, query); 

    if (new_lambda > lambda) {
      if (first) {
	first = 0;
	scale = scale_high = 1.0 + 0.05;
	scale_low = 1.0;
	too_high = 1;
      } else {
	if (!too_high) break;
	scale = (scale_high += scale_high - 1.0);
      }
    } else if (new_lambda > 0) {
      if (first) {
	first = 0;
	scale_high = 1.0;
	scale = scale_low = 1.0 - 0.05;
	too_high = 0;
      } else {
	if (too_high) break;
	scale = (scale_low += scale_low - 1.0);
      }
    } else {
      fprintf(stderr, "new_lambda (%g) <= 0; matrix has positive average score", new_lambda);
      exit(1);
    }
  }

  /* now do binary search between low and high */
  for (i = 0 ; i < 10 ; i++) {
    scale = 0.5 * (scale_high + scale_low);
    fill_pam(pam2p, n0, 20, freq2d, scale);
    new_lambda = get_lambda(pam2p, n0, 20, query);
    
    if (new_lambda > lambda) scale_low = scale;
    else scale_high = scale;
  }

  scale = 0.5 * (scale_high + scale_low);
  fill_pam(pam2p, n0, 20, freq2d, scale);

  /* 
  fprintf(stderr, "final scale: %g\n", scale);

  for (i = 0 ; i < n0 ; i++) {
    fprintf(stderr, "%4d %c:  ", i+1, query[i]);
    for (j = 1 ; j <= 20 ; j++) {
      fprintf(stderr, "%4d", pam2p[i][j]);
    }
    fprintf(stderr, "\n");
  }
  */

  free(freq2d[0]);
  free(freq2d);

  free(query);
}

#if defined(SSEARCH) || (defined(PRSS) && !defined(FASTX))
void
parse_pssm_asn_fa(FILE *afd, int *n_rows, int *n_cols,
		  unsigned char **query, double ***freqs);

void
read_asn_pssm(unsigned char *aa0, int n0, int nsq,
	      double pamscale, 
	      FILE *fp, struct pstruct *ppst) {

  int i, j, len, k;
  int qi, rj;	/* qi - index query; rj - index residues (1-20) */
  int **pam2p;
  int first, too_high;
  unsigned char *query, ctmp;
  char dline[512];
  double freq, **freq2d, lambda, new_lambda;
  double scale, scale_high, scale_low;
  int n_rows, n_cols;

  pam2p = ppst->pam2p[0];

  parse_pssm_asn_fa(fp, &n_rows, &n_cols, &query, &freq2d);

  if (n_cols != n0) { 
    fprintf(stderr, " query length: %d != n_cols: %d\n",n0, n_cols);
    exit(1);
  }

  for (qi = 0 ; qi < n0 ; qi++) {
    for (rj = 0 ; rj < 20 ; rj++) {
      if (freq2d[qi][rj] > 1e-20) {
	freq = log(freq2d[qi][rj] /((double) (rrcounts[rj+1])/(double) rrtotal));
	freq /= pamscale; /* this gets us close to originial pam scores */
	freq2d[qi][rj] = freq;
      }
      else {		
	/* when blastpgp decides to leave something out, it puts 0's in all the frequencies
	   in the binary checkpoint file.  In the ascii version, however, it uses BLOSUM62
	   values.  I will put in scoring matrix values as well */

	freq2d[qi][rj] = ppst->pam2[0][aa0[qi]][rj+1];
      }
    }
  }

  /* now figure out the right scale */
  scale = 1.0;
  lambda = get_lambda(ppst->pam2[0], 20, 20, ustandard_aa);

  /* should be near 1.0 because of our initial scaling by ppst->pamscale */
  /* fprintf(stderr, "real_lambda: %g\n", lambda); */

  /* get initial high/low scale values: */
  first = 1;
  while (1) {
    fill_pam(pam2p, n0, 20, freq2d, scale);
    new_lambda = get_lambda(pam2p, n0, 20, query); 

    if (new_lambda > lambda) {
      if (first) {
	first = 0;
	scale = scale_high = 1.0 + 0.05;
	scale_low = 1.0;
	too_high = 1;
      } else {
	if (!too_high) break;
	scale = (scale_high += scale_high - 1.0);
      }
    } else if (new_lambda > 0) {
      if (first) {
	first = 0;
	scale_high = 1.0;
	scale = scale_low = 1.0 - 0.05;
	too_high = 0;
      } else {
	if (too_high) break;
	scale = (scale_low += scale_low - 1.0);
      }
    } else {
      fprintf(stderr, "new_lambda (%g) <= 0; matrix has positive average score", new_lambda);
      exit(1);
    }
  }

  /* now do binary search between low and high */
  for (i = 0 ; i < 10 ; i++) {
    scale = 0.5 * (scale_high + scale_low);
    fill_pam(pam2p, n0, 20, freq2d, scale);
    new_lambda = get_lambda(pam2p, n0, 20, query);
    
    if (new_lambda > lambda) scale_low = scale;
    else scale_high = scale;
  }

  scale = 0.5 * (scale_high + scale_low);
  fill_pam(pam2p, n0, 20, freq2d, scale);

  /* 
  fprintf(stderr, "final scale: %g\n", scale);

  for (i = 0 ; i < n0 ; i++) {
    fprintf(stderr, "%4d %c:  ", i+1, query[i]);
    for (j = 1 ; j <= 20 ; j++) {
      fprintf(stderr, "%4d", pam2p[i][j]);
    }
    fprintf(stderr, "\n");
  }
  */

  free(freq2d[0]);
  free(freq2d);

  free(query);
}
#endif

void
last_params(unsigned char *aa0, int n0, 
	    struct mngmsg *m_msg,
	    struct pstruct *ppst
#ifdef PCOMPLIB
	    , struct qmng_str *qm_msg
#endif
	    ) {
  int i, nsq;
  FILE *fp;

  if (n0 < 0) { return;}

  ppst->n0 = m_msg->n0;

  if (ppst->ext_sq_set) { nsq = ppst->nsqx; }
  else {nsq = ppst->nsq;}

/* currently, profiles are only available for SSEARCH, PRSS */
#if defined(SSEARCH) || defined(PRSS)

  ppst->pam2p[0] = alloc_pam2p(n0,nsq);
  ppst->pam2p[1] = alloc_pam2p(n0,nsq);

  if (ppst->pam_pssm) {
    if ((ppst->pgpfile_type == 0) && (fp=fopen(ppst->pgpfile,"rb"))) {
      read_pssm(aa0, n0, ppst->nsq, ppst->pamscale, fp, 0, ppst);
      extend_pssm(aa0, n0, ppst);
    }
    else if ((ppst->pgpfile_type == 1) && (fp=fopen(ppst->pgpfile,"r"))) {
      read_pssm(aa0, n0, ppst->nsq, ppst->pamscale, fp, 1, ppst);
      extend_pssm(aa0, n0, ppst);
    }
#if defined(SSEARCH) || (defined(PRSS) && !defined(FASTX))
    else if ((ppst->pgpfile_type == 2) && (fp=fopen(ppst->pgpfile,"rb"))) {
      read_asn_pssm(aa0, n0, ppst->nsq, ppst->pamscale, fp, ppst);
      extend_pssm(aa0, n0, ppst);
    }
#endif
    else {
      fprintf(stderr," Could not open PSSM file: %s\n",ppst->pgpfile);
      ppst->pam_pssm = 0;
      return;
    }
  }
#endif

#if defined(FASTF) || defined(FASTS) || defined(FASTM)
  m_msg->nm0 = 1;
  for (i=0; i<n0; i++)
    if (aa0[i]==EOSEQ || aa0[i]==ESS) m_msg->nm0++;

/*
  for FASTS, we can do statistics in one of two different ways
  if there are <= 10 query fragments, then we calculate probabilistic
  scores for every library sequence.  If there are > 10 fragments, this
  takes much too long and too much memory, so we use the old fashioned
  raw score only z-score normalized method initially, and then calculate
  the probabilistic scores for the best hits.  To scale those scores, we
  also need a set of random probabilistic scores.  So we do the qshuffle
  to get them.

  For FASTF, precalculating probabilities is prohibitively expensive,
  so we never do it; FASTF always acts like FASTS with nfrags>10.

*/

#if defined(FASTS) || defined(FASTM)
  if (m_msg->nm0 > 10) m_msg->escore_flg = 0;
  else m_msg->escore_flg = 1;
#endif

  if (m_msg->escore_flg && (ppst->zsflag&1)) {
    m_msg->last_calc_flg = 0;
    m_msg->qshuffle = 0;
  }
  else {	/* need random query, second set of 2000 scores */
    m_msg->last_calc_flg = 1;
    m_msg->qshuffle = 1;
  }
#else
  m_msg->last_calc_flg = 0;
  m_msg->qshuffle = 0;
  m_msg->escore_flg = 0;
  m_msg->nm0 = 1;
#endif

/* adjust the ktup if appropriate */  

  if (!ktup_set && pgm_def_arr[ppst->pgm_id].ktup > 0) {
    if (m_msg->qdnaseq == SEQT_PROT) {
      ppst->param_u.fa.ktup = pgm_def_arr[ppst->pgm_id].ktup;
#if defined(FASTS) || defined(FASTM)
      if (n0 > 100) ppst->param_u.fa.ktup = 2;
#endif
      if (n0 < 40) ppst->param_u.fa.ktup = 1;
    }
    else if (m_msg->qdnaseq == SEQT_DNA || m_msg->qdnaseq == SEQT_RNA) {
      if (n0 < 20) ppst->param_u.fa.ktup = 1;
#if defined(FASTS) || defined(FASTM)
      /* with the current (April 12 2005) dropfs2.c - ktup cannot be > 2 */
      else ppst->param_u.fa.ktup = 2;
#else
      else if (n0 < 50) ppst->param_u.fa.ktup = 2;
      else if (n0 < 100)  ppst->param_u.fa.ktup = 3;
#endif
    }
  }

#ifdef PCOMPLIB
  qm_msg->nm0 = m_msg->nm0;
  qm_msg->escore_flg = m_msg->escore_flg;
  qm_msg->qshuffle = m_msg->qshuffle;
  qm_msg->pam_pssm = 0;
#endif
}

/* given a good profile in ppst->pam2p[0], make an extended profile
   in ppst->pam2p[1]
*/
void
extend_pssm(unsigned char *aa0, int n0, struct pstruct *ppst) {

  int i, j, nsq;
  int sa_x, sa_t, sa_b, sa_z;
  int **pam2p0, **pam2p1;

  nsq = ppst->nsq;

  pam2p0 = ppst->pam2p[0];
  pam2p1 = ppst->pam2p[1];

  sa_x = pascii['X'];
  sa_t = pascii['*'];
  sa_b = pascii['B'];
  sa_z = pascii['Z'];

  /* fill in boundaries, B, Z, *, X */
  for (i=0; i<n0; i++) {
    pam2p0[i][0] = -BIGNUM;
    pam2p0[i][sa_b] = (int)
      (((float)pam2p0[i][pascii['N']]+(float)pam2p0[i][pascii['D']]+0.5)/2.0);
    pam2p0[i][sa_z] = (int)
      (((float)pam2p0[i][pascii['Q']]+(float)pam2p0[i][pascii['E']]+0.5)/2.0);
    pam2p0[i][sa_x] = ppst->pam_xm;
    pam2p0[i][sa_t] = ppst->pam_xm;
  }

  /* copy pam2p0 into pam2p1 */
  for (i=0; i<n0; i++) {
    pam2p1[i][0] = -BIGNUM;
    for (j=1; j<=ppst->nsq; j++) {
      pam2p1[i][j] = pam2p0[i][j];
    }
  }

  /* then fill in extended characters, if necessary */
  if (ppst->ext_sq_set) {
    for (i=0; i<n0; i++) {
      for (j=1; j<=ppst->nsq; j++) {
	pam2p0[i][nsq+j] = pam2p0[i][j];
	pam2p1[i][nsq+j] = ppst->pam_xm;
      }
    }
  }
}
