/*** Program alcomp1.c for comparison of two alignments
Input:  alignments in CLUSTAL or STOCKHOLM format

Output:
line
#TAB<aln1 name>TAB<aln2 name>TAB<length_aln1>
TAB<length of "gap-filtered" aln1>TAB<length_aln2>TAB<length of "gap-filtered" aln2>
TAB<Nseqs in aln1>TAB<neff_aln1>TAB<Nseqs in aln2>TAB<neff_aln2>TAB<score>TAB<Evalue>\n

and optimal Smith-Waterman alignment.

Command line:
-i <infile1> -j <infile2>  -o <outfile>
[-g <threshold for "gap content filtering" of columns, e.g. 0.5>
-t <threshold for marking "moderately gapped" regions, e.g. 0.25 >
  -f <score formula> -s <path to substitution matrix> -q <path to frequency of pairs matrix>
-e <switch for adjustment of gap extension penalty according to gap content in opposite column>
-r <switch for reduction of column-column scores according to gap content>  ]


***/


#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <ctype.h>
#include <string.h>
#include <malloc.h>
#include <stddef.h>
#include "clustalw.h"
#include "blosum62.h"

#define NR_END 1
#define FREE_ARG char*

#define SQUARE(a) ((a)*(a))
#define NUM_METHOD 9
#define MAX_WINDOW 20
#define MAX_DELTASITE 20
#define MAXSTR   100001
#define INDI -100

#define JMAX 40
#define IA 16807
#define IM 2147483647
#define AM (1.0/IM)
#define IQ 127773
#define IR 2836
#define NTAB 32
#define NDIV (1+(IM-1)/NTAB)
#define EPS 1.2e-7
#define RNMX (1.0-EPS)

#define NRANSI
#define SWAP(a,b) temp=(a);(a)=(b);(b)=temp;
#define M 7
#define NSTACK 50

#define LN2 0.69314718055994528623
#define LAMB_UNG 0.3176

#define NCOLMAX 100000
#define NFILES_MAX 1500
#define NSEQ_OUT 1


char *digit="0123456789";
void nrerror(char error_text[]);
char *cvector(long nl, long nh);
int *ivector(long nl, long nh);
double *dvector(long nl, long nh);
char **cmatrix(long nrl, long nrh, long ncl, long nch);
int **imatrix(long nrl, long nrh, long ncl, long nch);
double **dmatrix(long nrl, long nrh, long ncl, long nch);
char **cmatrix(long nrl, long nrh, long ncl, long nch);
double ***d3tensor(long nrl,long nrh,long ncl,long nch,long ndl,long ndh);

void free_ivector(int *v, long nl, long nh);
void free_dvector(double *v, long nl, long nh);
void free_cvector(char *v, long nl, long nh);
void free_dmatrix(double **m, long nrl, long nrh, long ncl, long nch);
void free_imatrix(int **m, long nrl, long nrh, long ncl, long nch);
void free_cmatrix(char **m, long nrl, long nrh, long ncl, long nch);

int a3let2num(char *let);
int am2num_c(int c);
int am2num(int c);
int am2numBZX(int c);

char am2lower(char inchr);
char am2upper(char inchr);

static void *mymalloc(int size);
char *strsave(char *str);
char *strnsave(char *str, int l);
static char **incbuf(int n, char **was);
static int *incibuf(int n, int *was);

void err_readali(int err_num);
void readali(char *filename);
static void printali_ali(char *argt, int chunk, int n1, int n2, int len, char **aname1, char **aname2,
char **aseqGap1, char **aseqGap2, int *start1, int *start2, int *positive, int **col_score, int score);
int **ali_char2int(char **aseq,int start_num, int start_seq);
int **read_alignment2int(char *filename,int start_num,int start_seq);

char **traceback(char **aseq_mat1, char **aseq_mat2, int n1, int n2, int start_ali1, int start_ali2, int
end_ali1, int end_ali2, int **tracebackDir, int **flagNewGapQuery, int **flagNewGapDb, char **aseqGapTr1, char **aseqGapTr2);

void **traceback_outputPos(int start_ali1, int start_ali2, int end_ali1, int end_ali2, int **tracebackDir, int **flagNewGapQuery, int **flagNewGapDb, int *apos1, int *apos2);

void counter(int b);
double effective_number(int **ali, int *marks, int n, int start, int end);
double effective_number_nogaps(int **ali, int *marks, int n, int start, int end);
double effective_number_nogaps_expos(int **ali, int *marks, int n, int start, int end, int pos);

void **freqInt(int **ali,int nal, int alilen, int **f,int *num_gaps,int
*effindiarr,double gap_threshold, double *p_comp);
void **freqIntMaskGaps(int **ali,int nal, int alilen, int **f, double gap_threshold, double gapRegionMin, double *p_comp);

double *overall_freq(int **ali, int startp, int endp, int *mark);
double *overall_freq_wgt(int **ali,int startp,int endp,int *mark,double *wgt);
double *h_weight(int **ali, int ip);
double **h_freq(int **ali, double **f, double **hfr);
double *entro_conv(double **f, int **ali, double *econv);

double *variance_conv(double **f, int **ali, double **oaf, double *vconv);
double *pairs_conv(double **f,int **ali,int **matrix1,int indx,double *pconv);

 	   typedef struct _score_Vector{
		int *noGap, *gapExists, *noGapOld, *gapExistsOld, *prevScoreGapQueryOld, *noGapStore, *gapExistsStore, *prevScoreGapQueryOldStore;
	   } score_Vector;

	int *dbSequence, queryEnd, dbEnd, queryStart,dbStart;
	int gapOpen, gapExtend, dbLength, queryLength;
/*****	double **query, **matrix; *****/
 static int SmithWatermanScore( double **score_matrix, int queryLength, int dbLength, int gapOpen, int gapextend, int queryEnd, int dbEnd, int **tracebackDir,
int **flagNewGapQuery, int **flagNewGapDb);

static int SmithWatermanFindStart( double **score_matrix, int queryLength,  int dbLength,  int gapOpen, int gapExtend, int queryEnd, int dbEnd, int
score, int queryStart, int dbStart);

int score, End1, End2, Start1, Start2;

int gap__open, gap__extend;

int *Sequence2;
double lambda_al, score_scale;

int ScoreForTwoRows(double *subjectRow, double *queryRow);
double ScoreForTwoRows_Model(int *cntRow1, int *cntRow2, double *pseudoCntRow1, double *pseudoCntRow2);

double ScoreForTwoRows_Model2(int *cntRow1, int *cntRow2, double *pseudoCntRow1, double *pseudoCntRow2);
double ScoreForTwoRows_Model3(int *cntRow1, int *cntRow2, double *pseudoCntRow1, double *pseudoCntRow2, double n_effColumn1, double n_effColumn2);
double ScoreForTwoRows_Model4(int *cntRow1, int *cntRow2, double *pseudoCntRow1, double *pseudoCntRow2, double n_effColumn1, double n_effColumn2);
double ScoreForTwoRows_Model5(int *cntRow1, int *cntRow2, double *pseudoCntRow1, double *pseudoCntRow2, double n_effColumn1, double n_effColumn2, double score_scale, double b);
double ScoreForTwoRows_Model6(int pos1, int pos2, double score_scale, double b);

double ScoreForTwoRows_smat3_21(int pos1, int pos2);
double ScoreForTwoRows_smat3_22(int pos1, int pos2);
double ScoreForTwoRows_smat3_23(int pos1, int pos2);
double ScoreForTwoRows_smat3_27(int pos1, int pos2);
double ScoreForTwoRows_smat3_28(int pos1, int pos2);

/* JP */
double ScoreForTwoRows_smat1q_2d(int pos1, int pos2);

double Sgap6_smat(int pos1, int pos2, double b);
double Sgap6_smat_off(int pos1, int pos2, double b);
double (*sgapfcn)(int pos1, int pos2, double b);

double GapExtend1(int pos2, double b);
double GapExtend1_off(int pos2, double b);
double (*g_e1)(int pos2, double b);
double GapExtend2(int pos1, double b);
double GapExtend2_off(int pos1, double b);
double (*g_e2)(int pos1, double b);

int *ScoreOverColumn (int colScore, int flag1, int flag2, int flag3, int flag4, int flag5, int flag6, int *column_Score);

void *ReadRef (char *inputfile);
int CompareAlnVsReferenceAln (int *apos1, int *apos2, int *aposref1, int *aposref2, int start_ref1, int start_ref2, int end_ref1, int end_ref2 /* , int coverage1, int coverage2, int accuracy1, int accuracy2*/ );

double p_dayhoff[]={0, 0.0143859, 0.0384319, 0.0352742, 0.0195027, 0.086209, 0.048466, 0.0708956, 0.0866279, 0.0193078,
0.0832518, 0.0457631, 0.0610127, 0.0695179, 0.0390894, 0.0367281, 0.0570451, 0.0580589, 0.0244313, 0.043972, 0.0620286};

double p_rbnsn[]={0.013298, 0.038556, 0.032165, 0.022425, 0.090191, 0.05142, 0.064409, 0.078047, 0.019246, 0.073772, 0.052028, 0.058413,
0.071198, 0.044873, 0.042644, 0.05364, 0.062949, 0.021992, 0.051295, 0.057438};

double **blosum, **bl_ctrl;
double **qBlosum, *p_blsm;

double smatrix[21][21];

/* JP */
int debug;

/* The freqs in other (Dayhoff) order


A  		 R  	 N  	 D  	 C  		 Q
0.0866279 0.043972  0.0390894 0.0570451 0.0193078 0.0367281

E
0.0580589

G   		 H   I   	L  	 K   		M
0.0832518 0.0244313 0.048466  0.086209  0.0620286 0.0195027

F
0.0384319


P  		 S   T   	W   	Y  		 V
0.0457631 0.0695179 0.0610127 0.0143859 0.0352742 0.0708956

A   R   N   D   C   Q
E
G   H   I   L   K   M
F
P   S   T   W   Y   V  */


char **aname, **aname1, **aname2, **aseq, **aseq1, **aseq2;
int nal, nal1, nal2, nalmerge, alilen, alilen1, alilen2, maxalilen,
*astart, *astart1, *astart2, *alen;

int **align_mat1, **align_mat2;
int n_lowgaps, alilen_mat1, alilen_mat2;
char **aseq_mat1, **aseq_mat2;
char **aseqGapTr1, **aseqGapTr2;
int **tracebackDir;
int **flagNewGapQuery, **flagNewGapDb;
int *positive, **col_score;
int posGp, segment_len;

int *apos1, *apos2, *aposref1, *aposref2;
int start_ref1, start_ref2, end_ref1, end_ref2, reflen_nogp;
double coverage1, coverage2, falsecov, accuracy1, accuracy2;

static int **alignment1, **alignment2;
double **u_oaf,**h_oaf;
char *am="-WFYMLIVACGPTSNQDEHRKBZX*.wfymlivacgptsnqdehrkbzx";
char *am3[]={
"---",
"TRP",
"PHE",
"TYR",
"MET",
"LEU",
"ILE",
"VAL",
"ALA",
"CYS",
"GLY",
"PRO",
"THR",
"SER",
"ASN",
"GLN",
"ASP",
"GLU",
"HIS",
"ARG",
"LYS",
"ASX",
"GLX",
"UNK",
"***"
"...",
};

double  **read_aa_dmatrix(FILE *fmat);
int  **identity_imat(long n);
void argument();
void print_parameters(FILE *outfile,char *argi,char *argo,int nt,char *argt,int argb,char *args,int argm,int argf,int argc, int argw, char *argn,char *arga,double argg, char *argp,char *argd);

double gammln(double xx);
float ran1(long *idum);
double ScoreForRow(int *cntRow, double *pseudoCntRow);

double funcAl(double x,int len1,int len2, double *score_matrix_srt);

double funcAaFreq(double x, double *p_comp1, double *p_comp2);

double func_precise1(double x,int ntot, double *fixedPseudoCountRow);

double lambdaAl1(int len1, int len2, double *score_matrix_srt);

void **pseudoCounts(double **matrix, double n_eff, int len, double **pseudoCnt);
void **neffsForEachCol_maskGapReg(int **ali, int n, int len, double effgapmax, double effgapRegionMin, double **n_effAa, double *sum_eff_let, int *maskgapRegion, int *apos_filtr, int *len_lowgaps, double *n_eff);

void currdirlimits_walk (int nlow, int nhigh, void (*fcn)(char *));
void get_Neff_andQ(char *filename);

void sort(int n, double arr[]);

/*double *score_matrix_srt; */
double **score_matrix, *score_matrix_srt;
int **newScore_mat;
int **matrix1, **matrix2;

double *ident1, *ident2;
int **count;
int setsize = 1000;
int *maskgaps, *maskgaps1, *maskgaps2, *maskgapRegion, *maskgapRegion1, *maskgapRegion2;
double **pseudoCnt1, **pseudoCnt2;
double *p_comp1, *p_comp2;
double n_c;
double n_eff1, n_eff2;
double **n_effAa1, **n_effAa2;
double *sum_eff_let1, *sum_eff_let2;
int *apos_filtr1, *apos_filtr2;

/* int *ali1seq1, *ali1seq2; */

int scoreGivenEnd, score_final;

double Evalue;
double lambda_len[] = {0.277, 2.25}; /* coefficients for linear approximation of lambda(len) and K(len) */
double K_len[] = {0.044, 7.4};

double lambda_est, K_est; /* estimated lambda_g and K_g values */

double scorebylambda;

double lambda_u;
double b = 1.0;
double f= 32.0;

int **fV_RepeatOpenGapQuery, **fV_RepeatOpenGapDb;
int **fV_DbInClosestNewGapDb, **fV_QInClosestNewGapQuery;

int flag_errread = 0;

char **aseq_out1, **aseq_out2;
int *positive_out;
int inputpos1, inputpos2, outputpos, pos_beforegap1, pos_beforegap2, step1, step2;


static char str[MAXSTR+1];

	char ARG_I1[200],ARG_I2[200],ARG_O[100],ARG_P[100],ARG_D[50],ARG_S[100],ARG_Q[100],ARG_N[50],ARG_A[50];
	int ARG_F=1,ARG_C=0,ARG_V=0,ARG_M=0,ARG_B=60;
	int ARG_GO = 10, ARG_GE = 1; /* penalties for gap opening and extension */
	int ARG_E = 0; /* switch for adjustment of gap extension penalty according to gap content in opposite column */
	int ARG_R = 0; /* switch for reduction of column-column scores according to gap content */
	double ARG_G=1.0, ARG_T=1.0; /* thresholds of gap content for column excision and for "gapped regions" with waiving of 1st gap__open penalty */
	double ARG_L = LAMB_UNG; /* Ungapped lambda */

void prfprfdiag2(int fraglen, double *score) {

	double s=0;
	int i,j;

	for(i=1;i<=fraglen;i++) s = s + ScoreForTwoRows_smat3_21(i, i);

	*score = s;
}

void prfprfdiag1(int fraglen, double *score) {

        double s=0;
        int i,j;

        for(i=1;i<=fraglen;i++) {
		if(debug>1)fprintf(stdout, "%8.3f ", ScoreForTwoRows_smat1q_2d(i, i) );
		s = s + ScoreForTwoRows_smat1q_2d(i, i);
	}

        *score = s;
}





/*int main(int argc, char *argv[]) */
void prfprfmatrix(int **align1, int **align2, int alnlength1, int alnlength2, int nali1, int nali2, double **prfprfmat)
{
	FILE *smatrixfile,*qmatrixfile, *bl62file, *qbl62file;
	FILE *fout, *fpdb,*matrixfile,*fpdbout,*fp,*ft;
	int i,j,k,l,nt=0;
	int jposnogp, jmat;
	int fcount=0, fi=0;

	double av_sumlet1, av_sumlet2;

	char *sqbuf;
	int jj;
	int *Sequence2;

	char *bl62qijLocation = "blosum62.qij";
	char *bl62Location = "blosum62.sij";

	double (*scorefcn)(int pos1, int pos2);

	int nentry1, nentry2;

	int len_out;

	/* JP */
	double **score_matrix;
	double miniscore;
	static int int_scale=100;

	/*read input arguments */
	strcpy(ARG_S, bl62Location);
	strcpy(ARG_Q, bl62qijLocation);


        /* if(argc<=2) { argument(); exit(0);}
	for(i=1;i<argc;i++) {
	    if(strcmp(argv[i],"-i")==0) {strcpy(ARG_I1,argv[i+1]);i++;continue;}
	    if(strcmp(argv[i],"-j")==0) {strcpy(ARG_I2,argv[i+1]);i++;continue;}
	    if(strcmp(argv[i],"-o")==0) {strcpy(ARG_O,argv[i+1]);i++;continue;}
	    if(strcmp(argv[i],"-p")==0) {strcpy(ARG_P,argv[i+1]);i++;continue;}
	    if(strcmp(argv[i],"-d")==0) {strcpy(ARG_D,argv[i+1]);i++;continue;}
	    if(strcmp(argv[i],"-s")==0) {strcpy(ARG_S,argv[i+1]);i++;continue;}
	    if(strcmp(argv[i],"-q")==0) {strcpy(ARG_Q,argv[i+1]);i++;continue;}
	    if(strcmp(argv[i],"-n")==0) {strcpy(ARG_N,argv[i+1]);i++;continue;}
	    if(strcmp(argv[i],"-e")==0) {sscanf(argv[i+1],"%d",&ARG_E);i++;continue;}
	    if(strcmp(argv[i],"-r")==0) {sscanf(argv[i+1],"%d",&ARG_R);i++;continue;}
	    if(strcmp(argv[i],"-f")==0) {sscanf(argv[i+1],"%d",&ARG_F);i++;continue;}
    	    if(strcmp(argv[i],"-O")==0) {sscanf(argv[i+1],"%d",&ARG_GO);i++;continue;}
	    if(strcmp(argv[i],"-E")==0) {sscanf(argv[i+1],"%d",&ARG_GE);i++;continue;}
	    if(strcmp(argv[i],"-L")==0) {sscanf(argv[i+1],"%lf",&ARG_L);i++;continue;}
	    if(strcmp(argv[i],"-b")==0) {sscanf(argv[i+1],"%d",&ARG_B);i++;continue;}
	    if(strcmp(argv[i],"-c")==0) {sscanf(argv[i+1],"%d",&ARG_C);i++;continue;}
	    if(strcmp(argv[i],"-v")==0) {sscanf(argv[i+1],"%d",&ARG_V);i++;continue;}
	    if(strcmp(argv[i],"-m")==0) {sscanf(argv[i+1],"%d",&ARG_M);i++;continue;}
	    if(strcmp(argv[i],"-g")==0) {sscanf(argv[i+1],"%lf",&ARG_G);i++;continue;}
	    if(strcmp(argv[i],"-t")==0) {sscanf(argv[i+1],"%lf",&ARG_T);i++;continue;}

	} */

        if((ARG_F>8)||(ARG_F<0)){fprintf(stderr,"column-column score calculation method(-f): \n1, 3_21; 2, 3_22 ...., 8, 3_28\n");
                    exit(0);}
        if((ARG_E>1)||(ARG_E<0)){fprintf(stderr,"adjustment of gap extension penalty depending on gap content in opposite column (-e): \n0, no adjustment; 1, adjust gap extension\n");
                exit(0);}
        if((ARG_R>1)||(ARG_R<0)){fprintf(stderr,"reduction of column-column scores according to gap content (-r): \n0, no reduction; 1, adjust scores\n");
                exit(0);}
	if((ARG_G>1.0)||(ARG_G<=0)){fprintf(stderr,"gap content(-g) to eliminate a column must be no more than 1 and more than 0 \n");
		    exit(0);}

	gap__open = f*ARG_GO;
	gap__extend = f*ARG_GE;
	lambda_u = ARG_L/f;
	/* substitution matrix: should be bits (BLOSUM62 if reading from input failed) */

	for(i=0;i<=20;i++) for(j=0;j<=20;j++) smatrix[i][j] = blosum62_smatrix[i][j];


	for(i=0;i<=20;i++) {
		for(j=0;j<=20;j++) {
			 smatrix[i][j] *= LN2;
		}
	}


	/* q_ij matrix:  qij for BLOSUM62 if reading from input failed */

	alignment1 = align1;
	alignment2 = align2;
	alilen1 = alnlength1;
	alilen2 = alnlength2;
	nal1 = nali1;
	nal2 = nali2;
	if(debug>0) fprintf(stdout, "alcomp2: %d %d %d %d \n", alilen1, alilen2, nal1, nal2);

/* check switch for optional modes of comparison:
 1st seq VS full ali; full ali VS 1st seq; 1st seq VS 1st seq
*/
	switch (ARG_V) {
			case 1:
				nal1 = 1; break;
			case 2:
				nal2 = 1; break;
			case 3:
				nal1 =1; nal2= 1; break;
			}

/* filter "gapped" columns, get total n_eff, n_eff for each aa in each column and mask moderately "gapped regions" */
	apos_filtr1 = ivector(1,alilen1);
	n_effAa1 = dmatrix(1,alilen1,0,20);
	sum_eff_let1 = dvector(1,alilen1);
	maskgapRegion1 = ivector(0,alilen1+1);
	/*fprintf(stdout, "***********\n");
	for(i=1;i<=nal1;i++) {
	  for(j=1;j<=alilen1;j++) {
		fprintf(stdout, "%d ", alignment1[i][j]);
	  }
	  fprintf(stdout, "\n");
	}*/
	neffsForEachCol_maskGapReg(alignment1, nal1, alilen1, ARG_G, ARG_T, n_effAa1, sum_eff_let1, maskgapRegion1, apos_filtr1, &alilen_mat1, &n_eff1);
	//apos_filtr2 = ivector(1,alilen2);
	//n_effAa2 = dmatrix(1,alilen2,0,20);
	//sum_eff_let2 = dvector(1,alilen2);
	//maskgapRegion2 = ivector(0,alilen2+1);

	//neffsForEachCol_maskGapReg(alignment2, nal2, alilen2, ARG_G, ARG_T, n_effAa2, sum_eff_let2, maskgapRegion2, apos_filtr2, &alilen_mat2, &n_eff2);

/* Calculate target frequencies */
	pseudoCnt1 = dmatrix(1,alilen_mat1, 0, 20);
	//pseudoCnt2 = dmatrix(1,alilen_mat2, 0, 20);
	pseudoCounts(n_effAa1, n_eff1, alilen_mat1, pseudoCnt1);
	//pseudoCounts(n_effAa2, n_eff2, alilen_mat2, pseudoCnt2);

	/* JP */
  if(debug>1) {
	for(j=1;j<=20;j++) {
	  for(i=1;i<=alilen_mat1;i++) {
	 	fprintf(stdout, "%11.5g", n_effAa1[i][j]);
	  }
	  fprintf(stdout, "\n");
	}
	fprintf(stdout, "\n");
        for(j=1;j<=20;j++) {
	  for(i=1;i<=alilen_mat1;i++) {
                fprintf(stdout, "%11.5g", pseudoCnt1[i][j]);
	  }
	  fprintf(stdout, "\n");
	}

	/*for(i=1;i<=alilen_mat1;i++) {
	   fprintf(stdout, "%11.5g", sum_eff_let1[i]);
	}*/
	fprintf(stdout, "\n");
}

	return ;

	/*
	fprintf(stdout, "%8.5f\n", n_eff1);

	fprintf(stdout,".");*/

/* Choose column-to-column score formula to use */

		switch (ARG_F) {
			case 1:
				scorefcn = ScoreForTwoRows_smat3_21; break;
			case 2:
				scorefcn = ScoreForTwoRows_smat3_22; break;
			case 3:
				scorefcn = ScoreForTwoRows_smat3_23; break;
			case 7:
				scorefcn = ScoreForTwoRows_smat3_27; break;
			case 8:
				scorefcn = ScoreForTwoRows_smat3_28; break;
			default:
				scorefcn = ScoreForTwoRows_smat3_23;
			}


/* compute matrix of scores for column pairs, correct them with lambda */

	/*score_matrix = dmatrix(1,alilen_mat1,1,alilen_mat2); */
	score_matrix = prfprfmat;
	score_matrix_srt = dvector(1,alilen_mat1*alilen_mat2);



	k=1;
	for (i=1; i<=alilen_mat1; i++){
		for (j=1;j<=alilen_mat2;j++){
			score_matrix[i][j] = score_matrix_srt[k] = scorefcn(i, j);
			 /*fprintf(stdout, "%6.1f", score_matrix[i][j]);   */
			k++;
		}
		/*fprintf(stdout, "*\n");  */
	}

	free_dmatrix(n_effAa1,1,alilen1,0,20);
	free_dvector(sum_eff_let1,1,alilen1);
	free_dmatrix(pseudoCnt1,1,alilen1,0,20);
	free_dmatrix(n_effAa2,1,alilen2,0,20);
	free_dvector(sum_eff_let2,1,alilen2);
	free_dmatrix(pseudoCnt2,1,alilen2,0,20);

	lambda_al = lambdaAl1(alilen_mat1,alilen_mat2,score_matrix_srt);
	score_scale = lambda_al/ lambda_u;

	miniscore = 100000;
	for(i=1; i<=alilen_mat1; i++){
                for (j=1;j<=alilen_mat2;j++){
                        score_matrix[i][j] = score_matrix[i][j]*score_scale;
			if(miniscore > score_matrix[i][j]) miniscore = score_matrix[i][j];
			/*fprintf(stdout, "%5d%c%c", (int)(score_matrix[i][j]/30), am[alignment1[1][i]], am[alignment2[1][j]]); */
                }
		/*fprintf(stdout, "*\n"); */
        }
	if(debug>0) fprintf(stdout, "score_scale: %8.5f; lambda_al: %8.5f; lambda_u: %8.5f\n", score_scale, lambda_al, lambda_u);

	free_dvector(score_matrix_srt,1,alilen_mat1*alilen_mat2);
	free_ivector(apos_filtr1, 1,alilen1);
	free_ivector(apos_filtr2, 1,alilen2);
	free_ivector(maskgapRegion1, 0,alilen1+1);
	free_ivector(maskgapRegion2, 0,alilen2+1);

	/*exit(0);*/
	return;

}


/* from given alignment with aa as numbers, computes effective aa counts (PSIC->our formula)
and marks the columns with EFFECTIVE content of gaps > threshold (effgapmax) */

void **neffsForEachCol_maskGapReg(int **ali, int n, int len, double effgapmax, double effgapRegionMin, double **n_effAa, double *sum_eff_let, int *maskgapRegion, int *apos_filtr, int *len_lowgaps, double *nef)
{
	int i,j,k,l;
	int alilen_mat, nsymbols_col, nsymbols;
	double nef_loc;
	int ele;
        double *effnu;
        double sum_let;
        int *mark;
	int flagmark;

	if(debug>1) fprintf(stderr,".");

	effnu = dvector(0,20);
        mark = ivector(0,n+10);
	alilen_mat = 0;
	nsymbols = 0;
        for(j=1;j<=len;j++) {

                nsymbols_col = 0;
		sum_let=0;

                for(k=0;k<=20;++k){
/* Mark sequences that have amino acid  k (or gap, k=0) in this jth position */
                	 flagmark =0;
                        for(i=1;i<=n;++i){
                                mark[i]=0;

                                ele=ali[i][j];
                                if(ele==k){mark[i]=1; flagmark =1;}
                                ele=ali[i][j]-25;
                                if(ele==k) {mark[i]=1; flagmark =1;}
                        }

/* If aa k (or gap) is present in this position call compute k-th effective count */
                      if (flagmark == 1) {

				effnu[k]=effective_number_nogaps(ali,mark,n,1,len);
				nsymbols_col++;

			} else { effnu[k] = 0.0; }

                       if (k>0) sum_let += effnu[k];
                }


		if ( sum_let > 0 && 1.0*effnu[0]/(sum_let + effnu[0]) < effgapmax ) {
			alilen_mat++;
			for (k=0; k<=20; k++) {
				n_effAa[alilen_mat][k] = effnu[k];
			}
			sum_eff_let[alilen_mat] = sum_let;
			apos_filtr[alilen_mat] = j;
			nsymbols += nsymbols_col;

			if(1.0*effnu[0]/(sum_let + effnu[0]) < effgapRegionMin) {
				 maskgapRegion[alilen_mat] = 0;
			} else {
				maskgapRegion[alilen_mat] = 1;
			}

		}


	}


	nef_loc = 1.0*nsymbols/alilen_mat;
	*nef = nef_loc;
	*len_lowgaps = alilen_mat;

	maskgapRegion[0] = maskgapRegion[alilen_mat+1] = 0;

	free_dvector(effnu,0,20);
        free_ivector(mark,0,n+10);


}


void **pseudoCounts(double **matrix, double n_eff, int len, double **pseudoCnt)
{
	int i,j,k;
	double *f, *g;
	double sumN;
	double alpha, beta;

	alpha = n_eff-1;
	beta = 10.0;

	f = dvector(0,20);
	g = dvector(0,20);
	for (i=1;i<=len;i++) {
		sumN = 0;
		for (j=1;j<=20;j++) sumN += matrix[i][j];
		for (j=1;j<=20;j++) {
			f[j] = 1.0*matrix[i][j]/sumN;
		}
		for (j=1;j<=20;j++) {
			g[j] = 0;

	/***		for (k=1;k<=20;k++) g[j]+= qmatrix[j][k]*f[k]/p_dayhoff[k-1]; ***/

			for (k=1;k<=20;k++) g[j]+= qmatrix[j][k]*f[k]/p_rbnsn[k-1];

			pseudoCnt[i][j]= (alpha*f[j] + beta*g[j])/(alpha+beta);

	/**		fprintf (stderr, "f[%d] = %e  g[%d] = %e  pseudoCount = %e", j, f[j], j, g[j], pseudoCnt[i][j]); **/

		}
	}

}

/** Lambda finding - in version3_21, from all-to-all column comparisons (using funcAl()) ***/
double lambdaAl1(int len1, int len2, double *score_matrix_srt)
{
	int i,j,k;
	double dx,f,fmid,xmid,rtb;
	double x1, x2, xacc, scale;

/*** Sorting score_matrix_srt[] for the further summation of exp  ***/
	sort(len1*len2, score_matrix_srt);

	x1=1e-6;
	x2=1.0;
	xacc = 1e-10;

        f=funcAl(x1,len1,len2,score_matrix_srt);
        fmid=funcAl(x2,len1,len2,score_matrix_srt);

/***
       f=funcAaFreq(x1,p_comp1, p_comp2);
        fmid=funcAaFreq(x2,p_comp1, p_comp2);
***/

        if (f*fmid >= 0.0) nrerror("Root must be bracketed for bisection in rtbis");
        rtb = f < 0.0 ? (dx=x2-x1,x1) : (dx=x1-x2,x2);
        for (j=1;j<=JMAX;j++) {

        fmid=funcAl(xmid=rtb+(dx *= 0.5),len1,len2,score_matrix_srt);

/***                fmid=funcAaFreq(xmid=rtb+(dx *= 0.5),p_comp1, p_comp2); ***/

/***                fmid=func_precise1(xmid=rtb+(dx *= 0.5),n,qFixed); ***/


                if (fmid <= 0.0) rtb=xmid;
                if (fabs(dx) < xacc || fmid == 0.0) return rtb;
        }
        nrerror("Too many bisections in rtbis");
        return 0.0;

}

/** Calculation of the function: funcAl(lambda)=0  **/

/** Calculation of func based on all column combinations fom two alignments **/

double funcAl(double x, int len1, int len2, double *score_matrix_srt) {
	int i,j,k;
/***	double probln, prob, sumProb; ***/
	double f, ffin;
/***	int *countRow1, *countRow2;
	double *psCntRow1, *psCntRow2;
***/

/*    fprintf (stderr, "Entered func...\n");  */

	f=0.0;
	for (i=1; i<=len1*len2; i++) {

/*    fprintf (stderr, "i=%d\n",i);  */

/***		f+= prob*exp(x*ScoreForRow(countRow, fixedPseudoCountRow));  ***/

		f += exp(x*score_matrix_srt[i]);

/***		fprintf (stderr, "%d:%d->%e ", i,j,ScoreForTwoRows_Model6(i, j, 1.0, 0.0)); ***/

/***		if (f/f!=1.0) {
			fprintf (stderr, "pos1=%d pos2=%d::%e -> %e \n",i, j, ScoreForTwoRows_Model6(i,j, 1.0, 0.0), f);
			exit(0);
		}
***/
	}

/*		fprintf (stderr, "f = %e\n", f);  */

	ffin = f/(len1*len2) - 1.0;
	return ffin;
}


double funcAaFreq(double x,double *p_comp1, double *p_comp2)
{
 int i,j,k, sumPrevN;
        double prob1, prob2;
        double f, ffin;
        int *countRow;

/*    fprintf (stderr, "Entered func...\n");  */
        f=0.0;

        for (i=1; i<=20; i++) {
	        for (j=1; j<=20; j++) {

              	prob1 = p_comp1[i];
		prob2 = p_comp2[j];

/*        	  f+= prob1*prob2*exp(x*smatrix[i][j]);   */

             f+= prob1*prob2*exp(x*smatrix[i][j]/lambda_u);

/*              fprintf (stderr, "f = %e\n", f);  */

		 }
        }
        ffin = f - 1;
        return ffin;
}






double ScoreForRow(int *cntRow, double *pseudoCntRow) {
	int i,j;
	double s1, s2, s3;
	double sc, score;

	sc=0;
	for (i=1;i<=20;i++) {
/*		fprintf (stderr, "cntRow[%d]=%d  pseudoCntRow[%d]=%e  p_dayhoff[%d] = %e\n",i, cntRow[i], i,
pseudoCntRow[i], i, p_dayhoff[i]);
*/
		s1 = pseudoCntRow[i]/p_dayhoff[i];
		s2 = log(s1);
		s3=cntRow[i]*log(pseudoCntRow[i]/p_dayhoff[i]);
/***		 sc+=cntRow[i]*log(pseudoCntRow[i]/p_dayhoff[i-1]);  ***/

		 sc+=cntRow[i]*log(pseudoCntRow[i]/p_rbnsn[i-1]);

/**		fprintf (stderr, "s1= %e  s2= %e  s3= %e sc_%d = %e", s1, s2, s3, i, sc);  **/


	}

/*	fprintf (stderr, " %e ", s); */
	score = sc;
	return score;
}




float ran1(long *idum)
{
        int j;
        long k;
        static long iy=0;
        static long iv[NTAB];
        float temp;

        if (*idum <= 0 || !iy) {
                if (-(*idum) < 1) *idum=1;
                else *idum = -(*idum);
                for (j=NTAB+7;j>=0;j--) {
                        k=(*idum)/IQ;
                        *idum=IA*(*idum-k*IQ)-IR*k;
                        if (*idum < 0) *idum += IM;
                        if (j < NTAB) iv[j] = *idum;
                }
                iy=iv[0];
        }
        k=(*idum)/IQ;
        *idum=IA*(*idum-k*IQ)-IR*k;
        if (*idum < 0) *idum += IM;
        j=iy/NDIV;
        iy=iv[j];
        iv[j] = *idum;
        if ((temp=AM*iy) > RNMX) return RNMX;
        else return temp;
}

double gammln(double xx)
{
        double x,y,tmp,ser;
        static double cof[6]={76.18009172947146,-86.50532032941677,
                24.01409824083091,-1.231739572450155,
                0.1208650973866179e-2,-0.5395239384953e-5};
        int j;

        y=x=xx;
        tmp=x+5.5;
        tmp -= (x+0.5)*log(tmp);
        ser=1.000000000190015;
        for (j=0;j<=5;j++) ser += cof[j]/++y;
        return -tmp+log(2.5066282746310005*ser/x);
}

double  **read_aa_dmatrix(FILE *fmat){

/* read matrix from file *fmat */

int i,ncol,ri,rj,c,flag,j;
int col[31],row[31];
char stri[31];
double t;
double **mat;

mat=dmatrix(0,25,0,25);
for(i=0;i<=25;++i)for(j=0;j<=25;++j)mat[i][j]=0.0;

ncol=0;
i=0;
ri=0;
rj=0;
flag=0;

while( (c=getc(fmat)) != EOF){

if(flag==0 && c=='#'){flag=-1;continue;}
else if(flag==-1 && c=='\n'){flag=0;continue;}
else if(flag==-1){continue;}
else if(flag==0 && c==' '){flag=1;continue;}
else if(flag==1 && c=='\n'){flag=0;continue;}
else if(flag==1 && c==' '){continue;}
else if(flag==1){
                ++ncol;
                if(ncol>=25){nrerror("matrix has more than 24 columns: FATAL");exit(0);}
                col[ncol]=am2numBZX(c);
                continue;
                }
else if(flag==0 && c!=' ' && c!='#'){
                ri=0;
                ++rj;
                if(rj>=25){nrerror("matrix has more than 24 rows: FATAL");exit(0);}
                row[rj]=col[rj];
                for(i=0;i<=30;++i){stri[i]=' ';}
                stri[0]=c;
                j=0;
                flag=3;
                continue;
                }
else if (flag==2 && c==' '){for(i=0;i<=30;++i){stri[i]=' ';}j=0;continue;}
else if (flag==2 && c=='\n'){flag=0;continue;}
else if (flag==2){flag=3;stri[j]=c;if(j>30){nrerror("string too long:FATAL");exit(0);}continue;}
else if (flag==3 && c==' ' || flag==3 && c=='\n'){
                        j=0;
                        ++ri;
                        t=atof(stri);
                        mat[row[rj]][col[ri]]=t;
                        if (c=='\n')flag=0;else flag=2;
                        continue;
                        }
else if (flag==3){stri[++j]=c;continue;}

}

for(i=1;i<=ncol;i++) {
		for(j=i+1;j<=ncol;j++) {
/*			fprintf (stderr, "col[%d]=%d col[%d]=%d    ", i,col[i],j,col[j]);*/
			mat[col[i]][col[j]]=mat[col[j]][col[i]];
		}
	}

for(i=1;i<=ncol;i++) {
	for(j=1;j<=ncol;j++) {
		fprintf(stdout, "%7.4f,", mat[i][j]);
	}
	fprintf(stdout, "\n");
}
fprintf(stdout, "\n");
return mat;
}

int    **identity_imat(long n){
/* allocates square integer identity matrix of length n+1: m[0.n][0.n] */

int i,j;
int **m;
m=imatrix(0,n,0,n);
for(i=0;i<=n;++i)for(j=0;j<=n;++j){if(i==j)m[i][j]=1;else m[i][j]=0;}
return m;
}

void argument()
{
fprintf(stderr,"      al2co   arguments:\n");
fprintf(stderr,"\n");
fprintf(stderr,"  -i    Input alignment file [File in]\n");
fprintf(stderr,"        Format: ClustalW or simple alignment format\n");
fprintf(stderr,"        The title (first line) should begin with \"CLUSTAL W\", or\n");
fprintf(stderr,"        the title line should be deleted.\n");
fprintf(stderr,"\n");
fprintf(stderr,"  -o    Output file with conservation index for each position in the\n");
fprintf(stderr,"        alignment [File out] Optional\n");
fprintf(stderr,"        Default = STDOUT\n");
fprintf(stderr,"\n");
fprintf(stderr,"  -t    Output file with conservation index mapped to the alignment\n");
fprintf(stderr,"        [File out] Optional\n");
fprintf(stderr,"        Conservation indices are linearly rescaled to be from 0\n");
fprintf(stderr,"        to 9.99. C'=9.99*(C-MIN)/(MAX-MIN), where C and C' are the\n");
fprintf(stderr,"        the indices before and after rescaling respectively, MAX and\n");
fprintf(stderr,"        MIN are the highest index and lowest index before rescaling\n");
fprintf(stderr,"        respectively. The integer part of each rescaled index is\n");
fprintf(stderr,"        written out along with the sequence alignment.\n");
fprintf(stderr,"        Default = no output\n");
fprintf(stderr,"\n");
fprintf(stderr,"  -b    Block size of the output alignment file with conservation\n");
fprintf(stderr,"        [Integer] Optional\n");
fprintf(stderr,"        Default = 60\n");
fprintf(stderr,"\n");
fprintf(stderr,"  -s    Input file with the scoring matrix [File in] Optional\n");
fprintf(stderr,"        Format: NCBI\n");
fprintf(stderr,"        Notice: Scoring matrix is only used for sum-of-pairs measure\n");
fprintf(stderr,"        with option -c  2.\n");
fprintf(stderr,"        Default = identity matrix\n");
fprintf(stderr,"\n");
fprintf(stderr,"  -m    Scoring matrix transformation [Integer] Optional\n");
fprintf(stderr,"        Options:\n");
fprintf(stderr,"        0=no transformation,\n");
fprintf(stderr,"        1=normalization S'(a,b)=S(a,b)/sqrt[S(a,a)*S(b,b)],\n");
fprintf(stderr,"        2=adjustment S\"(a,b)=2*S(a,b)-(S(a,a)+S(b,b))/2\n");
fprintf(stderr,"        Default = 0\n");
fprintf(stderr,"\n");
fprintf(stderr,"  -f    Weighting scheme for amino acid frequency estimation [Integer] Optional\n");
fprintf(stderr,"        Options:\n");
fprintf(stderr,"        0=unweighted,\n");
fprintf(stderr,"        1=weighted by the modified method of Henikoff & Henikoff (2,3),\n");
fprintf(stderr,"        2=independent-count based (1)\n");
fprintf(stderr,"        Default = 2\n");
fprintf(stderr,"\n");
fprintf(stderr,"  -c    Conservation calculation method [Integer] Optional\n");
fprintf(stderr,"        Options:\n");
fprintf(stderr,"        0=entropy-based    C(i)=sum_{a=1}^{20}f_a(i)*ln[f_a(i)], where f_a(i)\n");
fprintf(stderr,"          is the frequency of amino acid a at position i,\n");
fprintf(stderr,"        1=variance-based   C(i)=sqrt[sum_{a=1}^{20}(f_a(i)-f_a)^2], where f_a\n");
fprintf(stderr,"          is the overall frequency of amino acid a,\n");
fprintf(stderr,"        2=sum-of-pairs measure   C(i)=sum_{a=1}^{20}sum_{b=1}^{20}f_a(i)*f_b(i)*S_{ab},\n");
fprintf(stderr,"          where S_{ab} is the element of a scoring matrix for amino acids a and b\n");
fprintf(stderr,"        Default = 0\n");
fprintf(stderr,"\n");
fprintf(stderr,"  -w    Window size used for averaging [Integer] Optional\n");
fprintf(stderr,"        Default = 1\n");
fprintf(stderr,"        Recommended value for motif analysis: 3\n");
fprintf(stderr,"\n");
fprintf(stderr,"  -n    Normalization option [T/F] Optional\n");
fprintf(stderr,"        Subtract the mean from each conservation index and divide by the\n");
fprintf(stderr,"        standard deviation.\n");
fprintf(stderr,"        Default = T\n");
fprintf(stderr,"\n");
fprintf(stderr,"  -a    All methods option [T/F] Optional\n");
fprintf(stderr,"        If set to true, the results of all 9 methods will be output.\n");
fprintf(stderr,"        1. unweighted entropy measure; 2. Henikoff entropy measure;\n");
fprintf(stderr,"        3. independent count entropy measure;\n");
fprintf(stderr,"        4. unweighted variance measure; 5. Henikoff variance measure;\n");
fprintf(stderr,"        6. independent count variance measure;\n");
fprintf(stderr,"        7. unweighted matrix-based sum-of-pairs measure;\n");
fprintf(stderr,"        8. Henikoff matrix-based sum-of-pairs measure;\n");
fprintf(stderr,"        9. independent count matrix-based sum-of-pairs measure;\n");
fprintf(stderr,"        Default = F\n");
fprintf(stderr,"\n");
fprintf(stderr,"  -g    Gap fraction to suppress conservation calculation [Real] Optional\n");
fprintf(stderr,"        The value should be more than 0 and no more than 1. Conservation\n");
fprintf(stderr,"        indices are calculated only for positions with gap fraction less\n");
fprintf(stderr,"        than the specified value. Otherwise, conservation indices will\n");
fprintf(stderr,"        be set to M-S, where M is the mean conservation value and S is\n");
fprintf(stderr,"        the standard deviation.\n");
fprintf(stderr,"        Default = 0.5\n");
fprintf(stderr,"\n");
fprintf(stderr,"  -p    Input pdb file [File in] Optional\n");
fprintf(stderr,"        The sequence in the pdb file should match exactly the first sequence of\n");
fprintf(stderr,"        the alignment.\n");
fprintf(stderr,"\n");
fprintf(stderr,"  -d    Output pdb file [File Out] Optional\n");
fprintf(stderr,"        The B-factors are replaced by the conservation indices.\n");
fprintf(stderr,"        Default = STDOUT\n");
fprintf(stderr,"\n");


}


void nrerror(char error_text[]){
fprintf(stderr,"%s\n",error_text);
fprintf(stderr,"FATAL - execution terminated\n");
exit(1);
}


char *cvector(long nl, long nh){
char *v;
v=(char *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(int)));
if (!v) nrerror("allocation failure in ivector()");
return v-nl+NR_END;
}


int *ivector(long nl, long nh){
int *v;
v=(int *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(int)));
if (!v) nrerror("allocation failure in ivector()");
return v-nl+NR_END;
}

/**** DUMP IN FAVOR OF NRUTIL.H ****/

long *lvector(long nl, long nh){
long int *v;
v=(long int *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(long int)));
if (!v) nrerror("allocation failure in lvector()");
return v-nl+NR_END;
}

double *dvector(long nl, long nh){
double *v;
v=(double *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(double)));
if (!v) nrerror("allocation failure in dvector()");
return v-nl+NR_END;
}

char **cmatrix(long nrl, long nrh, long ncl, long nch){
long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
char **m;
m=(char **)malloc((size_t)((nrow+NR_END)*sizeof(char*)));
if (!m) nrerror("allocation failure 1 in cmatrix()");
m += NR_END;
m -= nrl;

m[nrl]=(char *)malloc((size_t)((nrow*ncol+NR_END)*sizeof(char)));
if (!m[nrl]) nrerror("allocation failure 2 in cmatrix()");
m[nrl] += NR_END;
m[nrl] -= ncl;

for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;

return m;

}

int **imatrix(long nrl, long nrh, long ncl, long nch){
long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
int **m;
m=(int **)malloc((size_t)((nrow+NR_END)*sizeof(int*)));
if (!m) nrerror("allocation failure 1 in imatrix()");
m += NR_END;
m -= nrl;

m[nrl]=(int *)malloc((size_t)((nrow*ncol+NR_END)*sizeof(int)));
if (!m[nrl]) nrerror("allocation failure 2 in imatrix()");
m[nrl] += NR_END;
m[nrl] -= ncl;

for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;

return m;

}

double **dmatrix(long nrl, long nrh, long ncl, long nch){
long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
double **m;
m=(double **)malloc((size_t)((nrow+NR_END)*sizeof(double*)));
if (!m) nrerror("allocation failure 1 in dmatrix()");
m += NR_END;
m -= nrl;

m[nrl]=(double *)malloc((size_t)((nrow*ncol+NR_END)*sizeof(double)));
if (!m[nrl]) nrerror("allocation failure 2 in dmatrix()");
m[nrl] += NR_END;
m[nrl] -= ncl;

for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;

return m;
}


double ***d3tensor(long nrl,long nrh,long ncl,long nch,long ndl,long ndh){
long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1,ndep=ndh-ndl+1;
double ***t;

t=(double ***) malloc((size_t)((nrow+NR_END)*sizeof(double**)));
if(!t)nrerror("allocation failure 1 in d3tensor()");
t += NR_END;
t -= nrl;

t[nrl]=(double **) malloc((size_t)((nrow*ncol+NR_END)*sizeof(double*)));
if(!t[nrl])nrerror("allocation failure 2 in d3tensor()");
t[nrl] += NR_END;
t[nrl] -= ncl;

t[nrl][ncl]=(double *) malloc((size_t)((nrow*ncol*ndep+NR_END)*sizeof(double)));
if(!t[nrl][ncl])nrerror("allocation failure 3 in d3tensor()");
t[nrl][ncl] += NR_END;
t[nrl][ncl] -= ndl;

for(j=ncl+1;j<=nch;j++) t[nrl][j]=t[nrl][j-1]+ndep;
for(i=nrl+1;i<=nrh;i++){
	t[i]=t[i-1]+ncol;
	t[i][ncl]=t[i-1][ncl]+ncol*ndep;
	for(j=ncl+1;j<=nch;j++)t[i][j]=t[i][j-1]+ndep;
	}
return t;
}

void free_ivector(int *v, long nl, long nh)
/* free an int vector allocated with ivector() */
{
	free((FREE_ARG) (v+nl-NR_END));
}

void free_cvector(char *v, long nl, long nh)
/* free an unsigned char vector allocated with cvector() */
{
	free((FREE_ARG) (v+nl-NR_END));
}

void free_dvector(double *v, long nl, long nh)
/* free a double vector allocated with dvector() */
{
	free((FREE_ARG) (v+nl-NR_END));
}



void free_dmatrix(double **m, long nrl, long nrh, long ncl, long nch)
/* free a double matrix allocated by dmatrix() */
{
	free((FREE_ARG) (m[nrl]+ncl-NR_END));
	free((FREE_ARG) (m+nrl-NR_END));
}

void free_imatrix(int **m, long nrl, long nrh, long ncl, long nch)
/* free an int matrix allocated by imatrix() */
{
	free((FREE_ARG) (m[nrl]+ncl-NR_END));
	free((FREE_ARG) (m+nrl-NR_END));
}

void free_cmatrix(char **m, long nrl, long nrh, long ncl, long nch)
/* free a double matrix allocated by dmatrix() */
{
	free((FREE_ARG) (m[nrl]+ncl-NR_END));
	free((FREE_ARG) (m+nrl-NR_END));
}


int am2num(int c)
{
switch (c) {
           	 case 'W': case 'w':
                	c=1; break;
           	 case 'F': case 'f':
                	c=2; break;
           	 case 'Y': case 'y':
                	c=3; break;
           	 case 'M': case 'm':
                	c=4; break;
           	 case 'L': case 'l':
                	c=5; break;
           	 case 'I': case 'i':
          		c=6; break;
           	 case 'V': case 'v':
           		c=7; break;
          	 case 'A': case 'a':
			c=8; break;
           	 case 'C': case 'c':
                	c=9; break;
		 case 'G': case 'g':
			c=10; break;
           	 case 'P': case 'p':
             	 	c=11; break;
       		 case 'T': case 't':
			c=12; break;
	         case 'S': case 's':
			c=13; break;
           	 case 'N': case 'n':
                	c=14; break;
           	 case 'Q': case 'q':
                	c=15; break;
           	 case 'D': case 'd':
                	c=16; break;
           	 case 'E': case 'e':
                	c=17; break;
           	 case 'H': case 'h':
                	c=18; break;
           	 case 'R': case 'r':
                	c=19; break;
           	 case 'K': case 'k':
                	c=20; break;
           	 default :
			c=0;
		}
return (c);
}


int am2numBZX(c)
{
switch (c) {
                 case 'W': case 'w':
                        c=1; break;
                 case 'F': case 'f':
                        c=2; break;
                 case 'Y': case 'y':
                        c=3; break;
                 case 'M': case 'm':
                        c=4; break;
                 case 'L': case 'l':
                        c=5; break;
                 case 'I': case 'i':
                        c=6; break;
                 case 'V': case 'v':
                        c=7; break;
                 case 'A': case 'a':
                        c=8; break;
                 case 'C': case 'c':
                        c=9; break;
                 case 'G': case 'g':
                        c=10; break;
                 case 'P': case 'p':
                        c=11; break;
                 case 'T': case 't':
                        c=12; break;
                 case 'S': case 's':
                        c=13; break;
                 case 'N': case 'n':
                        c=14; break;
                 case 'Q': case 'q':
                        c=15; break;
                 case 'D': case 'd':
                        c=16; break;
                 case 'E': case 'e':
                        c=17; break;
                 case 'H': case 'h':
                        c=18; break;
                 case 'R': case 'r':
                        c=19; break;
                 case 'K': case 'k':
                        c=20; break;
                 case 'B': case 'b':
                        c=21; break;
                 case 'Z': case 'z':
                        c=22; break;
                 case 'X': case 'x':
                        c=23; break;
                 case '*':
                        c=24; break;
                 default :
                        c=0;
                }
return (c);
}


char am2lower(char inchr)
{
char c;
switch (inchr) {
              	 case '-':
                 	c='.'; break;
           	 case 'W':
                	c='w'; break;
           	 case 'F':
                	c='f'; break;
           	 case 'Y':
                	c='y'; break;
           	 case 'M':
                	c='m'; break;
           	 case 'L':
                	c='l'; break;
           	 case 'I':
          		c='i'; break;
           	 case 'V':
           		c='v'; break;
          	 case 'A':
			c='a'; break;
           	 case 'C':
                	c='c'; break;
		 case 'G':
			c='g'; break;
           	 case 'P':
             	 	c='p'; break;
       		 case 'T':
			c='t'; break;
	         case 'S':
			c='s'; break;
           	 case 'N':
                	c='n'; break;
           	 case 'Q':
                	c='q'; break;
           	 case 'D':
                	c='d'; break;
           	 case 'E':
                	c='e'; break;
           	 case 'H':
                	c='h'; break;
           	 case 'R':
                	c='r'; break;
           	 case 'K':
                	c='k'; break;
                 case 'B':
                 	c='b'; break;
                 case 'Z':
                 	c='z'; break;
                 case 'X':
                 	c='x'; break;
		 default :
                        c=inchr;
		}
return (c);
}

char am2upper(char inchr)
{
char c;
switch (inchr) {
              	 case '.':
                 	c='-'; break;
           	 case 'w':
                	c='W'; break;
           	 case 'f':
                	c='F'; break;
           	 case 'y':
                	c='Y'; break;
           	 case 'm':
                	c='M'; break;
           	 case 'l':
                	c='L'; break;
           	 case 'i':
          		c='I'; break;
           	 case 'v':
           		c='V'; break;
          	 case 'a':
			c='A'; break;
           	 case 'c':
                	c='C'; break;
		 case 'g':
			c='G'; break;
           	 case 'p':
             	 	c='P'; break;
       		 case 't':
			c='T'; break;
	         case 's':
			c='S'; break;
           	 case 'n':
                	c='N'; break;
           	 case 'q':
                	c='Q'; break;
           	 case 'd':
                	c='D'; break;
           	 case 'e':
                	c='E'; break;
           	 case 'h':
                	c='H'; break;
           	 case 'r':
                	c='R'; break;
           	 case 'k':
                	c='K'; break;
                 case 'b':
                 	c='B'; break;
                 case 'z':
                 	c='Z'; break;
                 case 'x':
                 	c='X'; break;
		 default :
                        c=inchr;
		}
return (c);
}





int **alignment;



static void *mymalloc(int size);
char *strsave(char *str);
char *strnsave(char *str, int l);
static char **incbuf(int n, char **was);
static int *incibuf(int n, int *was);

void readali(char *filename);
int **ali_char2int(char **aseq,int start_num, int start_seq);
int **read_alignment2int(char *filename,int start_num,int start_seq);

void counter(int b);
double effective_number(int **ali, int *marks, int n, int start, int end);
double effective_number_nogaps(int **ali, int *marks, int n, int start, int end);
double effective_number_nogaps_expos(int **ali, int *marks, int n, int start, int end, int pos);



static void *mymalloc(size)
int size;
{
	void *buf;

	if ((buf = malloc(size)) == NULL) {
		fprintf(stderr, "Not enough memory: %d\n", size);
		exit(1);
	}
	return buf;
}

char *strsave(str)
char *str;
{
	char *buf;
	int l;

	l = strlen(str);
	buf = mymalloc(l + 1);
	strcpy(buf, str);
	return buf;
}

char *strnsave(str, l)
char *str;
int l;
{
	char *buf;

	buf = mymalloc(l + 1);
	memcpy(buf, str, l);
	buf[l] = '\0';
	return buf;
}

static char **incbuf(n, was)
int n;
char **was;
{
	char **buf;
	char *aaa;

	buf = mymalloc((n+1) * sizeof(buf[0]));
	if (n > 0) {
		memcpy(buf, was, n * sizeof(was[0]));
		free(was);
	}
	buf[n] = NULL;
	return buf;
}

static int *incibuf(n, was)
int n, *was;
{
	int *ibuf;

	ibuf = mymalloc((n+1) * sizeof(ibuf[0]));

	if (n > 0) {
		memcpy(ibuf, was, n * sizeof(was[0]));
		free(was);
	}
	ibuf[n] = 0;
	return ibuf;
}
void err_readali(int err_num)
{
	fprintf(stderr,"Error with reading alignment: %d\n",err_num);
}

void readali(char *filename)

{
	FILE *fp;
	char *s, *ss, *seqbuf;
	int n, l, len, len0;
	int ii,mark=1;

	if ((fp = fopen(filename, "r")) == NULL) {
		fprintf(stderr, "No such file: \"%s\"\n", filename);
		err_readali(1);
/*		;exit(1); */
		flag_errread=1;
		return;
	}

	alilen = 0;
	nal = 0;
	n = 0;

	if(fgets(str, MAXSTR, fp) != NULL) {
		if(strncmp(str,"CLUSTAL ",8)!=0){rewind(fp);}
					}

	while (fgets(str, MAXSTR, fp) != NULL) {

/*	fprintf(stderr,"OK ");*/

		if (*str=='#' || strncmp(str,"//",2) == 0) {continue;}
		for (ss = str; isspace(*ss); ss++) ;
		if ((ii<=ss-str)&&(mark==0)) {continue;}

/*		fprintf(stderr, "n=%d ",n); */

		if (*ss == '\0') {
			if (n == 0) {
				continue;
			}
			if (nal == 0) {
				if (n == 0) {
					fprintf(stderr, "No alignments read\n");
					err_readali(2);
/*					exit(1); */
					flag_errread=1;
					return;
				}
				nal = n;
			} else if (n != nal) {
				fprintf(stderr, "Wrong nal, was: %d, now: %d\n", nal, n);
				err_readali(3);
/*				exit(1); */
				flag_errread=1;
				return;
			}
			n = 0;
			continue;
		}
		for (s = ss; *s != '\0' && !isspace(*s); s++) ;
		*s++ = '\0';

		if (nal == 0) {

			astart = incibuf(n, astart);
			alen = incibuf(n, alen);
			aseq = incbuf(n, aseq);
			aname = incbuf(n, aname);
			aname[n] = strsave(ss);

		} else {
			if (n < 0 || n >= nal) {
				fprintf(stderr, "Bad sequence number: %d of %d\n", n, nal);
				err_readali(4);
/*				exit(1); */
				flag_errread=1;
				return;
			}
			if (strcmp(ss, aname[n]) != 0) {
				fprintf(stderr, "Names do not match");
				fprintf(stderr, ", was: %s, now: %s\n", aname[n], ss);
				err_readali(5);
/*				exit(1); */
				flag_errread=1;
				return;
			}
		}
		for (ss = s; isspace(*ss); ss++);
		if(mark==1){
		ii = ss-str;
		mark=0;}

		for (s = ss; isdigit(*s); s++) ;
		if (isspace(*s)) {
			if (nal == 0) {
				astart[n] = atoi(ss);
			}
			for (ss = s; isspace(*ss); ss++);
		}
		for (s = ss, len=0, l = 0; *s != '\0' && !isspace(*s); s++) {
			if (isalpha(*s)) {
				l++;
			}

/*** Calculate len -- the full number of aa and gaps, excluding position numbers in the end ***/

			if (isalpha(*s) || *s == '-' || *s == '.') {
				len++;
			}


		}

/****		len = s - ss;  *************/



		if (n == 0) {
			len0 = len;
			alilen += len;
		} else if (len != len0) {
			fprintf(stderr, "wrong len for %s", aname[n]);
			fprintf(stderr, ", was: %d, now: %d\n", len0, len);
			err_readali(6);
/*			 exit(1); */
			flag_errread=1;
			return;
		}

		alen[n] += l;
		if (aseq[n] == NULL) {
			aseq[n] = strnsave(ss, len);
		} else {
			seqbuf = mymalloc(alilen+1);
			memcpy(seqbuf, aseq[n], alilen-len);
			free(aseq[n]);
			aseq[n] = seqbuf;
			memcpy(seqbuf+alilen-len, ss, len);
			seqbuf[alilen] = '\0';
		}
		n++;
	}
	if (nal == 0) {
		if (n == 0) {
			fprintf(stderr, "No alignments read\n");
			err_readali(7);
/*			exit(1); */
			flag_errread=1;
			return;
		}
		nal = n;
	} else if (n != 0 && n != nal) {
		fprintf(stderr, "Wrong nal, was: %d, now: %d\n", nal, n);
		err_readali(8);
/*		exit(1); */
		flag_errread=1;
		return;
	}

	fclose(fp);
}

/*** Print ali to stderr ****/
static void printali_ali(char *argo, int chunk, int n1, int n2, int len, char **aname1, char **aname2,
char **aseqGap1, char **aseqGap2, int *start1, int *start2, int *positive, int **col_score, int score)
{
        int i, j, k, jj, mlen, str_len, len_start;
	int ratio;
	int scoreFin;
	char arg_o[100], namebuf[100];
        char *sq;
	int *isq;
	char *sqn;
	FILE *fpp;
	strcpy(arg_o,argo);
	fpp=fopen(arg_o,"a");

	for (i=0; i<n1; i++) {
		start1[i] += apos1[1]-1;
		start2[i] += apos2[1]-1;
	}

        for (i=1, mlen=strlen(aname1[0]); i < n1; i++) {
                if (mlen < strlen(aname1[i])) {
                        mlen = strlen(aname1[i]);
                }
        }
        for (i=0 ; i < n2; i++) {
                if (mlen < strlen(aname2[i])) {
                        mlen = strlen(aname2[i]);
                }
        }

        jj = 0;

        do {

/* Print the chunk of the first alignment */
                if (jj > 0) {
                        fprintf(fpp, "\n");
		}

                for (i=0; i < n1; i++) {
     			strcpy(namebuf,aname1[i]);
     			fprintf(fpp, namebuf);
     			str_len = strlen(aname1[i]);
     			for(k=str_len;k<mlen+3;k++) fprintf(fpp," ");

    			if (jj==0) {

			for(len_start=0, ratio=start1[i]; ratio>0; ratio /= 10, len_start++);
     			fprintf(fpp, "%d", start1[i]);
     			} else {len_start=0;}

 			for(k=len_start;k<7;k++) fprintf(fpp," ");

                        sq = aseqGap1[i] + jj;

                        for (j=1; j+jj <=len && j <= chunk; j++) {

				 fprintf(fpp, "%c", sq[j]);
                        }

                       fprintf(fpp, "\n");
                }

			for(k=0;k<mlen+10;k++) fprintf(fpp," ");
			isq = positive + jj ;
			for (j=1; j+jj <= len && j <= chunk; j++) {
				if (isq[j]) {
					fprintf(fpp,"+");
				} else {
					fprintf(fpp," ");
				}
			}
			fprintf(fpp, "\n");


/*Print the chunk of the second alignment*/

                for (i=0; i < n2; i++) {
                        sqn = aname2[i];

                	for(k=0;k<mlen+3;k++){
				if(k<strlen(sqn)){fprintf(fpp,"%c",sqn[k]);}
				else fprintf(fpp," ");
			}

			if (jj==0) {

			for(len_start=0, ratio=start2[i]; ratio>0; ratio /= 10, len_start++);
/*     			len_start = log(start2[i])/log(10.0)+1; */
     			fprintf(fpp, "%d", start2[i]);
     			} else {len_start=0;}

 			for(k=len_start;k<7;k++) fprintf(fpp," ");


                        sq = aseqGap2[i] + jj;
                        for (j=1; j+jj <= len && j <= chunk; j++) {
                              fprintf(fpp, "%c", sq[j]);
                        }
                        fprintf(fpp, "\n");
                }
                fprintf(fpp, "\n");

                jj += chunk;

        } while (jj < len);

	fclose(fpp);
}

int **ali_char2int(char **aseq, int start_num, int start_seq){
/* fills the alignment ali[start_num..start_num+nal-1][start_seq..start_seq+alilen-1]
convetring charater to integer from aseq[0..nal-1][0..alilen-1]
*/

int i,j,end_num,end_seq;
int **ali;
end_num=start_num+nal-1;
end_seq=start_seq+alilen-1;
ali=imatrix(start_num,end_num,start_seq,end_seq);
for(i=start_num;i<=end_num;++i)for(j=start_seq;j<=end_seq;++j)ali[i][j]=am2num(aseq[i-start_num][j-start_seq]);
return ali;
}

int **read_alignment2int(char *filename,int start_num,int start_seq){
int **ali;
readali(filename);

if (flag_errread==1) return;

ali=ali_char2int(aseq,start_num,start_seq);
return ali;
}

double effective_number_nogaps(int **ali, int *marks, int n, int start, int end){

/* from the alignment of n sequences ali[1..n][1..l]
calculates effective number of sequences that are marked by 1 in mark[1..n]
for the segment of positions ali[][start..end]
Neff=ln(1-0.05*N-of-different-letters-per-site)/ln(0.95)
*/

int i,k,a,flag;
int *amco,lettercount=0,sitecount=0;
double letpersite=0,neff;

 amco=ivector(0,20);
for(k=start;k<=end;++k){
/******************DUMP the condition "consider only positions without gaps in the marked seqs" ***********/
/*****	flag=0;for(i=1;i<=n;++i)if(marks[i]==1 && ali[i][k]==0)flag=1;
	if(flag==1)continue;
*****/
	for(a=0;a<=20;++a)amco[a]=0;
	for(i=1;i<=n;++i)if(marks[i]==1)amco[ali[i][k]]++;
	flag=0;for(a=1;a<=20;++a)if(amco[a]>0){flag=1;lettercount++;}
	if(flag==1)sitecount++;

/*	else fprintf (stderr, "%d ", k); */
		       }
if(sitecount==0)letpersite=0;
else letpersite=1.0*lettercount/sitecount;

/*** neff = letpersite; ***/

 neff=-log(1.0-0.05*letpersite)/0.05129329438755;

 free_ivector(amco,0,20);
return neff;
}

int *letters;

void **freqInt(int **ali,int nal, int alilen, int **f,int *num_gaps,int
*effindiarr,double gap_threshold, double *p_comp)
{
	int i,j,k,effnumind, sumNC, fullCountNogaps;
	int count[21], sum_comp[21];

	fprintf (stderr, "freqInt started...\n");

	letters = ivector(0, alilen+1);
	letters[0]=1;

	/* find the number of frequences at each position */
	effnumind=0;
	sumNC = 0;
	for(i=0;i<=20;i++) sum_comp[i]=0;
	for(j=1;j<=alilen;j++){
		for(i=0;i<=20;i++) count[i]=0;
		for(i=1;i<=nal;i++) {
			if(ali[i][j]<=20) {
				count[ali[i][j]]++;
			}
			else {if(ali[i][j]>=25&&ali[i][j]<=45)
				{count[ali[i][j]-25]++;}
			      else {
				fprintf(stderr,"not good number for AA\n");
				fprintf(stderr,"%d", i);
				fprintf(stderr,"\n");
				fprintf(stderr,"%d", j);
				exit(0);
				   }
		}
	}
		/*** Adding to the sum of different symbols in the columns over the alignment, to derive N_C , and to the overall aa counts ***/

 		for (i=0; i<=20; i++) {
			 if(count[i]>0) {
				sumNC++;
				sum_comp[i]+=count[i];
			}

/*			fprintf  (stderr, "%d_%d ",count[i], sumNC);	*/
		}
/*		fprintf (stderr, "\n\n");	*/

		num_gaps[j] = count[0];
/***	   	 f[0][j] = count[0]*1.0/nal;    ***/
	   	f[j][0] = count[0];
		if(f[j][0]>nal) {
			fprintf(stderr,"gap number>total number\n");
			exit(0);
			      }

/* Eliminate the condition for small enough number of gaps */

/*		if(f[0][j]>=gap_threshold) {   /* ignore the case where gaps occur >= gap_threshold(percentage of gaps)
			f[0][j]=INDI;
			continue;
				}   */

		effnumind++;
		effindiarr[effnumind]=j;
		count[0]=nal-count[0];
		letters[j] = count[0];
/***		if(count[0]<=0){
			fprintf(stderr, "count[0] less than 0: %d  column = %d\n",count[0],j);
			exit(0);
			       }
***/
		for(k=1;k<=20;k++){
/* Eliminate the division by count[0] - not freqs but counts !!!   */
/*		f[k][j]=count[k]*1.0/count[0];			*/
		f[j][k]=count[k];
				  }
	}


	n_c = sumNC;
	n_c = n_c/alilen;

	fullCountNogaps = alilen*nal-sum_comp[0];
	for (i=1;i<=20;i++) p_comp[i] = 1.0*sum_comp[i]/fullCountNogaps;

	fprintf (stderr, "\n sumNC = %d  alilen = %d \n",sumNC, alilen);
	effindiarr[effnumind+1]=INDI;/*set the last element negative*/
	effindiarr[0]=effnumind;
}

/* Version for alscr_wwm1.c : no deleting columns wth gaps in 1st seq */
void **freqIntMaskGaps(int **ali,int nal, int alilen, int **f, double gapmax, double gapRegionMin, double *p_comp)
{
	int i,j,k, jnew, fullCountNogaps;
	int sumNC; /*sum of different symbols in columns over alignment, used to derived effective size n_c*/
	int count[21];
	int sum_comp[21]; /*overall aa counts*/
	letters = ivector(0, alilen+1);
	letters[0]=1;

	/* find the number of frequences at each position */
	sumNC = 0;
	for(i=0;i<=20;i++) sum_comp[i]=0;
	jnew = 0;
	for(j=1;j<=alilen;j++){
		for(i=0;i<=20;i++) count[i]=0;
		for(i=1;i<=nal;i++) {
			if(ali[i][j]<=20) {count[ali[i][j]]++;}
			else {if(ali[i][j]>=25&&ali[i][j]<=45) {count[ali[i][j]-25]++;}
			      else {
				fprintf(stderr,"not good number for AA\n");
				fprintf(stderr,"%d", i);
				fprintf(stderr,"\n");
				fprintf(stderr,"%d", j);
				exit(0);
				   }
			}
		}

 		for (i=0; i<=20; i++) {
			 if(count[i]>0) {
				sumNC++;
				sum_comp[i]+=count[i];
			}
		}

		if(count[0]>nal) {
			fprintf(stderr,"gap number>total number\n");
			exit(0);
			}

/* Eliminate the condition for small enough number of gaps */

/*		if(f[0][j]>=gap_threshold) {   /* ignore the case where gaps occur >= gap_threshold(percentage of gaps)
			f[0][j]=INDI;
			continue;
				}   */

		letters[j]=nal-count[0];

/* Raise flags for further deletion over higly gapped columns; mark "moderately gapped" regions */
		if(1.0*count[0]/nal >= gapmax /* || ali[1][j] == 0 || ali[1][j]==25 */ ) maskgaps[j] = 1;
		else {
			maskgaps[j] = 0;
			if(1.0*count[0]/nal >= gapRegionMin) maskgapRegion[j] = 1;
			else maskgapRegion[j] = 0;
			jnew++;
			for(k=0;k<=20;k++){f[jnew][k]=count[k];}
		}
	}
	n_lowgaps = jnew;
	n_c = sumNC;
	n_c = n_c/alilen;

	fullCountNogaps = alilen*nal-sum_comp[0];
	for (i=1;i<=20;i++) p_comp[i] = 1.0*sum_comp[i]/fullCountNogaps;


}

void **neffsForEachCol(int **ali, int n, int len, double **n_effAa, double *sum_eff_let)
{
	int i,j,k;
        int ele;
        double *effnu;
        double sum_let;
        int *mark;
	int flagmark;

	effnu = dvector(0,20);
        mark = ivector(0,n+10);
        for(j=1;j<=len;j++) {

        	sum_eff_let[j] = 0;
        	for(k=0;k<=20;++k){
        		 n_effAa[j][k]=0;
        		 }
        }

        for(j=1;j<=len;j++) {
                sum_let=0;

                for(k=0;k<=20;++k){
/* Mark sequences that have amino acid  k (or gap, k=0) in this jth position */
                	 flagmark =0;
                        for(i=1;i<=n;++i){
                                mark[i]=0;

                                ele=ali[i][j];
                                if(ele==k){mark[i]=1; flagmark =1;}
                                ele=ali[i][j]-25;
                                if(ele==k) {mark[i]=1; flagmark =1;}
                        }

/* If aa k (or gap) is present in this position call compute k-th effective count */
                        if (flagmark == 1) effnu[k]=effective_number_nogaps(ali,mark,n,1,len);
                        else effnu[k] = 0.0;

                       if (k>0) sum_let += effnu[k];
                }

      		for (k=0; k<=20; k++) {
			n_effAa[j][k] = effnu[k];
		}
		sum_eff_let[j] = sum_let;

	}
}


/*computes Smith-Waterman local alignment score and returns the
  evalue
  query is the query sequence
  queryLength is the length of query in amino acids
  dbSequence is the sequence corresponding to some matrix profile
  dbLength is the length of dbSequnece
  matrix is the position-specific matrix associated with dbSequence
  gapOpen is the cost of opening a gap
  gapExtend is the cost of extending an exisiting gap by 1 position
  queryEnd returns the final position in the query of an optimal
   local alignment
  dbEnd returns the final position in dbSequence of an optimal
   local alignment
  queryEnd and dbEnd can be used to run the local alignment in reverse
   to find optimal starting positions
  score is used to pass back the optimal score
  kbp holds the Karlin-Altschul paramters
  L holds an intermediate term for E-value computation
  adjustedDbLength is the adjusted database length used for e-value computation
  minGappedK holds the minimum gapped K for all matrices in the
  database, and is used for e-value computation */


 static int SmithWatermanScore(double **score_matrix,  int queryLength, int dbLength, int gapOpen, int gapextend, int queryEnd, int dbEnd, int **tracebackDir,
int **flagNewGapQuery, int **flagNewGapDb)
{
   int bestScore; /*best score seen so far*/
   int newScore;  /* score of next entry*/
   int bestQueryPos, bestDbPos; /*position ending best score in
                           query and database sequences*/
   int newGapCost; /*cost to have a gap of one character*/
   int gapExtend;
   int prevScoreNoGapQuery; /*score one row and column up
                               with no gaps*/
   int prevScoreGapQuery;   /*score if a gap already started in query*/
   int continueGapScore; /*score for continuing a gap in dbSequence*/
   int queryPos, dbPos; /*positions in query and dbSequence*/
/*   Nlm_FloatHi returnEvalue; /*e-value to return*/
   score_Vector scoreVector; /*keeps one row of the Smith-Waterman matrix
                           overwrite old row with new row*/
	int RowScore; /*score for match of two positions*/
	int gapDb2NoGap, gapQuery2NoGap, noGap2NoGap, score2NoGap;



/********************** Introduce arrays and variables:

int *fV_RepeatOpenGapQuery[queryPos] -- the current row of flagRepeatOpenGapQuery;

int *fV_RepeatOpenGapDb[queryPos] -- the current row of flagRepeatOpenGapDb;

int *fV_RepeatOpenGapDb1[queryPos] -- previous row of flagRepeatOpenGapDb, overwrite with the new row after
the SW matrix row is passed;

int *fV_RepeatOpenGapDb2[queryPos] -- the row of flagRepeatOpenGapDb 2 positions higher, overwrite with the
fV_RepeatOpenGapQuery1 after the SW matrix row is passed;

int *fV_RepeatOpenGapQuery1[queryPos] -- previous row of flagRepeatOpenGapQuery, overwrite with the new row after
the SW matrix row is passed;

int *fV_RepeatOpenGapQuery2[queryPos] -- the row of flagRepeatOpenGapQuery 2 positions higher, overwrite with the
fV_RepeatOpenGapQuery1 after the SW matrix row is passed;

int  *fV_QInClosestNewGapDb, *fV_DbInClosestNewGapQuery - current rows of the flags equal to the flags in the closest starting points of new gaps in Db (closest at the vertical queryPos = const) and Query (closest at the horizontal dbPos = const). They are used to calculate flagRepeatOpenGapQuery and flagRepeatOpenGapDb, respectively.

int  *fV_QInClosestNewGapDb1 - previous row of flagInClosestNewGapDb, overwrite with the new flagInClosestNewGapDb after the SW matrix row is passed;

int *fV_DbInClosestNewGapQuery1 -- previous row of flagInClosestNewGapQuery, overwrite with the new flagInClosestNewGapDb after the SW matrix row is passed;

??? Some of arrays fV_InClosestNewGap... can be probably replaced by some variables, since we use only positions (-1,-1), (0,-1), (-1,0) and (0,0), or we can use the scheme analogous to scoreVector.nogap and scoreVector.gapExist.

*****************************/
/***** int **fV_RepeatOpenGapQuery, **fV_RepeatOpenGapDb;  *****/
int **fV_QInClosestNewGapDb, **fV_DbInClosestNewGapQuery;

/*****int **fV_DbInClosestNewGapDb, **fV_QInClosestNewGapQuery; *****/


int flagRepeatOpenGapDb, flagRepeatOpenGapQuery;

   scoreVector.noGap = ivector (1,queryLength);
   scoreVector.gapExists = ivector (1,queryLength);

   fV_RepeatOpenGapQuery = imatrix (0,queryLength,0, dbLength);
   fV_RepeatOpenGapDb = imatrix (0,queryLength,0, dbLength);
   fV_QInClosestNewGapDb = imatrix (0,queryLength,0, dbLength);
   fV_DbInClosestNewGapQuery = imatrix (0,queryLength,0, dbLength);

   fV_DbInClosestNewGapDb = imatrix (0,queryLength,0, dbLength);
   fV_QInClosestNewGapQuery = imatrix (0,queryLength,0, dbLength);

   bestQueryPos = 0;
   bestDbPos = 0;
   bestScore = 0;
/***   newGapCost = gapOpen + gapExtend;  ***/
   for (queryPos = 1; queryPos <= queryLength; queryPos++) {
     scoreVector.noGap[queryPos] = 0;
     scoreVector.gapExists[queryPos] = -(gapOpen);
	}

   for (queryPos = 0; queryPos <= queryLength; queryPos++) {

	fV_RepeatOpenGapDb[queryPos][0] = 0;
	fV_RepeatOpenGapQuery[queryPos][0] = 0;

	fV_QInClosestNewGapDb[queryPos][0] = 0;
	fV_DbInClosestNewGapQuery[queryPos][0] = 0;
	fV_DbInClosestNewGapDb[queryPos][0] = 0;
	fV_QInClosestNewGapQuery[queryPos][0] = 0;

	fV_QInClosestNewGapDb[queryPos][1] = 0;
	fV_DbInClosestNewGapQuery[queryPos][1] = 0;
	fV_DbInClosestNewGapDb[queryPos][1] = 0;
	fV_QInClosestNewGapQuery[queryPos][1] = 0;
   }

   for(dbPos = 1; dbPos <= dbLength; dbPos++) {

     newScore = 0;
     noGap2NoGap = 0;
     prevScoreGapQuery = -(gapOpen);

        fV_RepeatOpenGapDb[0][dbPos]= 0;
	fV_RepeatOpenGapQuery[0][dbPos]= 0;

	fV_QInClosestNewGapDb[0][dbPos]= 0;
	fV_DbInClosestNewGapQuery[0][dbPos]= 0;
	fV_DbInClosestNewGapDb[0][dbPos]= 0;
	fV_QInClosestNewGapQuery[0][dbPos]= 0;

     for(queryPos = 1; queryPos <= queryLength; queryPos++) {

		flagNewGapQuery[queryPos][dbPos] = 0;
		flagNewGapDb[queryPos][dbPos] = 0;

/*** Check if we are in the gapped region of query and no gaps in db were opened against this region before;
	if TRUE eliminate gapOpen penalty; if this is the first position of the gapped region, reward the extending
	previous gap by compensating gapOpen in gapExtend ***/

	gapExtend = rint(g_e2(queryPos,b));

	flagRepeatOpenGapDb = fV_RepeatOpenGapDb[queryPos-1][dbPos];

	if (maskgapRegion1[queryPos]==1 && flagRepeatOpenGapDb==0) {newGapCost = gapExtend;}
	 else {	newGapCost = gapOpen + gapExtend; }

	if (maskgapRegion1[queryPos-1]==0 && maskgapRegion1[queryPos]==1 && fV_DbInClosestNewGapDb[queryPos-1][dbPos]==0) {
			gapExtend -= gapOpen;
	}

       /*testing scores with a gap in DB, either starting a new
         gap or extending an existing gap*/

       if ((newScore = newScore - newGapCost) >
	   (prevScoreGapQuery = prevScoreGapQuery - gapExtend)) {
         	prevScoreGapQuery = newScore;
         	flagNewGapQuery[queryPos][dbPos] = 1;
         }


/*** Check if we are in the gapped region of Db and no gaps were opened in query against this region before;
	if TRUE, eliminate gapOpen penalty ***/
         gapExtend = rint(g_e1(dbPos,b));

	flagRepeatOpenGapQuery = fV_RepeatOpenGapQuery[queryPos][dbPos-1];

	if (maskgapRegion2[dbPos]==1 && flagRepeatOpenGapQuery==0) {newGapCost = gapExtend;}
	else {newGapCost = gapOpen + gapExtend;}

/***	if this is the first position of the gapped region
 reward the extending previous gap by compensating gapOpen in gapExtend ***/
	if (maskgapRegion2[dbPos]==1 && maskgapRegion2[dbPos-1]==0 && fV_QInClosestNewGapQuery[queryPos][dbPos-1]==0) {
			gapExtend -= gapOpen;
		}

       /*testing scores with a gap in Query, either starting a new
         gap or extending an existing gap*/

       if ((newScore = scoreVector.noGap[queryPos] - newGapCost) >
           (continueGapScore = scoreVector.gapExists[queryPos] - gapExtend)) {
         continueGapScore = newScore;
         flagNewGapDb[queryPos][dbPos] = 1;
        }

       /*compute new score extending one position in query and db*/

       RowScore = rint(score_matrix[queryPos][dbPos]*score_scale - sgapfcn(queryPos,dbPos,b));

       newScore = noGap2NoGap + RowScore;

       if (newScore < 0)
       newScore = 0; /*Smith-Waterman locality condition*/

/*** TraceBackDir:
***/

	if (RowScore>0) {tracebackDir[queryPos][dbPos] = 6;}
	 else {	tracebackDir[queryPos][dbPos] = 5;}

	if (maskgapRegion1[queryPos] == 0) fV_RepeatOpenGapDb[queryPos][dbPos] = 0;
	else fV_RepeatOpenGapDb[queryPos][dbPos] = fV_RepeatOpenGapDb[queryPos-1][dbPos-1];

	if (maskgapRegion2[dbPos] == 0) fV_RepeatOpenGapQuery[queryPos][dbPos] = 0;
	else fV_RepeatOpenGapQuery[queryPos][dbPos] = fV_RepeatOpenGapQuery[queryPos-1][dbPos-1];

/**** Assign the flags coming from the closest NewGap in Query ******/

	if (flagNewGapDb[queryPos][dbPos] == 1) {
		 fV_DbInClosestNewGapQuery[queryPos][dbPos] = fV_RepeatOpenGapDb[queryPos][dbPos-1];
		 fV_QInClosestNewGapQuery[queryPos][dbPos] = fV_RepeatOpenGapQuery[queryPos][dbPos-1];
	} else {
	        fV_DbInClosestNewGapQuery[queryPos][dbPos] = fV_DbInClosestNewGapQuery[queryPos][dbPos-1];
        	fV_QInClosestNewGapQuery[queryPos][dbPos] = fV_QInClosestNewGapQuery[queryPos][dbPos-1];
        }

        if (maskgapRegion2[dbPos]==1) {fV_QInClosestNewGapQuery[queryPos][dbPos] = 1;}
        if (maskgapRegion1[queryPos]==1) {fV_DbInClosestNewGapQuery[queryPos][dbPos] = 1;}

/**** Assign the flags coming from the closest NewGap in Db ******/

	if (flagNewGapQuery[queryPos][dbPos] == 1) {
		fV_QInClosestNewGapDb[queryPos][dbPos] = fV_RepeatOpenGapQuery[queryPos-1][dbPos];
		fV_DbInClosestNewGapDb[queryPos][dbPos] = fV_RepeatOpenGapDb[queryPos-1][dbPos];
          } else {
        	fV_QInClosestNewGapDb[queryPos][dbPos] = fV_QInClosestNewGapDb[queryPos-1][dbPos];
        	fV_DbInClosestNewGapDb[queryPos][dbPos] = fV_DbInClosestNewGapDb[queryPos-1][dbPos];
        }

	if (maskgapRegion1[queryPos]==1) {fV_DbInClosestNewGapDb[queryPos][dbPos] = 1;}
	if (maskgapRegion2[dbPos]==1) {fV_QInClosestNewGapDb[queryPos][dbPos] = 1;}


       /*test two alternatives*/

/*** Gap in Db ***/

   if (newScore < prevScoreGapQuery) {
         newScore = prevScoreGapQuery;

/**** Determine tracebackDir pointer and the flags fV_RepeatOpenGapQuery and fV_RepeatOpenGapDb ******/

         if (flagNewGapQuery[queryPos][dbPos] == 1) { tracebackDir[queryPos][dbPos] = 1;}
          else {tracebackDir[queryPos][dbPos] = 2;}

	fV_RepeatOpenGapQuery[queryPos][dbPos] = fV_QInClosestNewGapDb[queryPos][dbPos];
	fV_RepeatOpenGapDb[queryPos][dbPos] = fV_DbInClosestNewGapDb[queryPos][dbPos];

   }

/*** Gap in Query ***/

   if (newScore < continueGapScore) {
         newScore = continueGapScore;

/**** Determine tracebackDir pointer and the flags fV_RepeatOpenGapQuery and fV_RepeatOpenGapDb ******/

         if (flagNewGapDb[queryPos][dbPos] == 1) {tracebackDir[queryPos][dbPos] = 3;}
         else {tracebackDir[queryPos][dbPos] = 4;}

	fV_RepeatOpenGapDb[queryPos][dbPos] = fV_DbInClosestNewGapQuery[queryPos][dbPos];
	fV_RepeatOpenGapQuery[queryPos][dbPos] = fV_QInClosestNewGapQuery[queryPos][dbPos];

       }

       noGap2NoGap = scoreVector.noGap[queryPos];
       scoreVector.noGap[queryPos] = newScore;
       scoreVector.gapExists[queryPos] = continueGapScore;

       if (newScore > bestScore) {
         bestScore = newScore;
         bestDbPos = dbPos;
         bestQueryPos = queryPos;
       }

	fprintf(stderr,"");
    }

 }

	fprintf(stderr,".");

  /* Closed assignments of arrays *score and *...End */

   /* MemFree(scoreVector);  */
   if (bestScore < 0)
     bestScore = 0;
   End1 = bestQueryPos;
   End2 = bestDbPos;

   return (bestScore);

/*   returnEvalue = scoreToEvalue(effSearchSpace, bestScore,  kbp); */
}


/*computes where optimal Smith-Waterman local alignment starts given the
  ending positions
  query is the query sequence
  queryLength is the length of query in amino acids
  dbSequence is the sequence corresponding to some matrix profile
  dbLength is the length of dbSequnece
  matrix is the position-specific matrix associated with dbSequence
  gapOpen is the cost of opening a gap
  gapExtend is the cost of extending an exisiting gap by 1 position
  queryEnd is the final position in the query of an optimal
   local alignment
  dbEnd is the final position in dbSequence of an optimal
   local alignment
  queryEnd and dbEnd can be used to run the local alignment in reverse
   to find optimal starting positions
  these are passed back in queryStart and dbStart
  the optimal score is passed in to check when it has
   been reached going backwards
  the score is also returned
  */

	static int SmithWatermanFindStart( double **score_matrix, int
queryLength, int dbLength, int gapOpen, int gapExtend,int queryEnd, int dbEnd, int score, int queryStart, int dbStart)
{

   int bestScore; /*best score seen so far*/
   int newScore;  /* score of next entry*/
   int bestQueryPos, bestDbPos; /*position starting best score in
                           query and database sequences*/
   int newGapCost; /*cost to have a gap of one character*/
   int prevScoreNoGapQuery; /*score one row and column up
                               with no gaps*/
   int prevScoreGapQuery;   /*score if a gap already started in query*/
   int continueGapScore; /*score for continuing a gap in dbSequence*/
   int queryPos, dbPos; /*positions in query and dbSequence*/

   score_Vector scoreVector; /*keeps one row of the Smith-Waterman matrix
                           overwrite old row with new row*/

   int flagNewGapQuery_Rev, flagNewGapDb_Rev;



/*  scoreVector = (SWpairs *) MemNew(queryLength * sizeof(SWpairs)); */
/********************** Introduce arrays and variables:

int *fV_RepeatOpenGapQuery[queryPos] -- the current row of flagRepeatOpenGapQuery;

int *fV_RepeatOpenGapDb[queryPos] -- the current row of flagRepeatOpenGapDb;

int *fV_RepeatOpenGapDb1[queryPos] -- previous row of flagRepeatOpenGapDb, overwrite with the new row after
the SW matrix row is passed;

int *fV_RepeatOpenGapDb2[queryPos] -- the row of flagRepeatOpenGapDb 2 positions higher, overwrite with the
fV_RepeatOpenGapQuery1 after the SW matrix row is passed;

int *fV_RepeatOpenGapQuery1[queryPos] -- previous row of flagRepeatOpenGapQuery, overwrite with the new row after
the SW matrix row is passed;

int *fV_RepeatOpenGapQuery2[queryPos] -- the row of flagRepeatOpenGapQuery 2 positions higher, overwrite with the
fV_RepeatOpenGapQuery1 after the SW matrix row is passed;

int  *fV_QInClosestNewGapDb, *fV_DbInClosestNewGapQuery - current rows of the flags equal to the flags in the closest starting points of new gaps in Db (closest at the vertical queryPos = const) and Query (closest at the horizontal dbPos = const). They are used to calculate flagRepeatOpenGapQuery and flagRepeatOpenGapDb, respectively.

int  *fV_QInClosestNewGapDb1 - previous row of flagInClosestNewGapDb, overwrite with the new flagInClosestNewGapDb after the SW matrix row is passed;

int *fV_DbInClosestNewGapQuery1 -- previous row of flagInClosestNewGapQuery, overwrite with the new flagInClosestNewGapDb after the SW matrix row is passed;

??? Some of arrays fV_InClosestNewGap... can be probably replaced by some variables, since we use only positions (-1,-1), (0,-1), (-1,0) and (0,0), or we can use the scheme analogous to scoreVector.nogap and scoreVector.gapExist.

*****************************/
int *fV_RepeatOpenGapQuery, *fV_RepeatOpenGapDb, *fV_RepeatOpenGapDb1, *fV_RepeatOpenGapDb2, *fV_RepeatOpenGapQuery1, *fV_RepeatOpenGapQuery2;
int *fV_QInClosestNewGapDb, *fV_DbInClosestNewGapQuery, *fV_QInClosestNewGapDb1 , *fV_DbInClosestNewGapQuery1;

int *fV_DbInClosestNewGapDb, *fV_QInClosestNewGapQuery, *fV_DbInClosestNewGapDb1 , *fV_QInClosestNewGapQuery1;

int flagRepeatOpenGapDb, flagRepeatOpenGapQuery;


   scoreVector.noGap = ivector (1,queryLength);
   scoreVector.gapExists = ivector (1,queryLength);

   fV_RepeatOpenGapQuery = ivector (1,queryLength+1);
   fV_RepeatOpenGapDb = ivector (1,queryLength+1);
   fV_RepeatOpenGapDb1 = ivector (1,queryLength+1);
   fV_RepeatOpenGapDb2 = ivector (1,queryLength+1);
   fV_RepeatOpenGapQuery1 = ivector (1,queryLength+1);
   fV_RepeatOpenGapQuery2 = ivector (1,queryLength+1);
   fV_QInClosestNewGapDb = ivector (1,queryLength+1);
   fV_QInClosestNewGapDb1 = ivector (1,queryLength+1);
   fV_DbInClosestNewGapQuery = ivector (1,queryLength+1);
   fV_DbInClosestNewGapQuery1 = ivector (1,queryLength+1);

   fV_DbInClosestNewGapDb = ivector (1,queryLength+1);
   fV_DbInClosestNewGapDb1 = ivector (1,queryLength+1);
   fV_QInClosestNewGapQuery = ivector (1,queryLength+1);
   fV_QInClosestNewGapQuery1 = ivector (1,queryLength+1);


   bestQueryPos = 0;
   bestDbPos = 0;
   bestScore = 0;

   for (queryPos = 1; queryPos <= queryLength; queryPos++) {
     scoreVector.noGap[queryPos] = 0;
     scoreVector.gapExists[queryPos] = -(gapOpen);
   }

      for (queryPos = 1; queryPos <= queryLength+1; queryPos++) {
	fV_RepeatOpenGapDb[queryPos] = 0;
	fV_RepeatOpenGapQuery[queryPos] = 0;

	fV_RepeatOpenGapDb1[queryPos] = 0;
   	fV_RepeatOpenGapDb2[queryPos] = 0;
	fV_RepeatOpenGapQuery1[queryPos] = 0;
	fV_RepeatOpenGapQuery2[queryPos] = 0;

	fV_QInClosestNewGapDb1[queryPos] = 0;
	fV_DbInClosestNewGapQuery1[queryPos] = 0;
	fV_DbInClosestNewGapDb1[queryPos] = 0;
	fV_QInClosestNewGapQuery1[queryPos] = 0;

	fV_QInClosestNewGapDb[queryPos] = 0;
	fV_DbInClosestNewGapQuery[queryPos] = 0;
	fV_DbInClosestNewGapDb[queryPos] = 0;
	fV_QInClosestNewGapQuery[queryPos] = 0;

   }

   for(dbPos = dbEnd; dbPos >= 1; dbPos--) {

     newScore = 0;
     prevScoreNoGapQuery = 0;
     prevScoreGapQuery = -(gapOpen);

     for(queryPos = queryEnd; queryPos >= 1; queryPos--) {
	flagNewGapQuery_Rev = 0;
	flagNewGapDb_Rev = 0;

/*** Check if we are in the gapped region of query and no gaps were opened against this region before;
	if TRUE, eliminate gapOpen penalty;
	if this is the first position of the gapped region, reward the extending
	previous gap by compensating gapOpen in gapExtend ***/

	gapExtend = rint(g_e2(queryPos,b));

	flagRepeatOpenGapDb = fV_RepeatOpenGapDb[queryPos+1];

	if (maskgapRegion1[queryPos]==1 && flagRepeatOpenGapDb==0) {newGapCost = gapExtend;}
	 else {newGapCost = gapOpen + gapExtend;}

	if (maskgapRegion1[queryPos+1]==0 && maskgapRegion1[queryPos]==1 && fV_DbInClosestNewGapDb[queryPos+1]==0) {
			gapExtend -= gapOpen;
	}

       /*testing scores with a gap in DB, either starting a new
         gap or extending an existing gap*/

      if ((newScore = newScore - newGapCost) >
	   (prevScoreGapQuery = prevScoreGapQuery - gapExtend)) {
         prevScoreGapQuery = newScore;
         flagNewGapQuery_Rev = 1;
       }


/*** Check if we are in the gapped region of Db and no gaps were opened against this region before;
	if TRUE, eliminate gapOpen penalty;
	if this is the first position of the gapped region, reward the extending
	previous gap by compensating gapOpen in gapExtend ***/
        gapExtend = rint(g_e1(dbPos,b));

	flagRepeatOpenGapQuery = fV_RepeatOpenGapQuery1[queryPos];

	if (maskgapRegion2[dbPos]==1 && flagRepeatOpenGapQuery==0) {newGapCost = gapExtend;}
	 else {	newGapCost = gapOpen + gapExtend;}

	if (maskgapRegion2[dbPos]==1 && maskgapRegion2[dbPos+1]==0 && fV_QInClosestNewGapQuery1[queryPos]==0) {
			gapExtend -= gapOpen;
		}

       /*testing scores with a gap in Query, either starting a new
         gap or extending an existing gap*/

       if ((newScore = scoreVector.noGap[queryPos] - newGapCost) >
           (continueGapScore = scoreVector.gapExists[queryPos] -gapExtend)) {
         continueGapScore = newScore;
         flagNewGapDb_Rev = 1;
        }

       /*compute new score extending one position in query and dbSequence*/

       newScore = prevScoreNoGapQuery + rint(score_matrix[queryPos][dbPos]*score_scale - sgapfcn(queryPos, dbPos, b));

       if (newScore < 0)
       newScore = 0; /*Smith-Waterman locality condition*/

	if (maskgapRegion1[queryPos] == 0) fV_RepeatOpenGapDb[queryPos] = 0;
	else fV_RepeatOpenGapDb[queryPos] = fV_RepeatOpenGapDb1[queryPos+1];

	if (maskgapRegion2[dbPos] == 0) fV_RepeatOpenGapQuery[queryPos] = 0;
	else fV_RepeatOpenGapQuery[queryPos] = fV_RepeatOpenGapQuery1[queryPos+1];

/**** Assign the flags coming from the closest NewGap in Query ******/

	if (flagNewGapDb_Rev == 1) {
		fV_DbInClosestNewGapQuery[queryPos] = fV_RepeatOpenGapDb1[queryPos];
		fV_QInClosestNewGapQuery[queryPos] = fV_RepeatOpenGapQuery1[queryPos];
	} else {
		fV_DbInClosestNewGapQuery[queryPos] = fV_DbInClosestNewGapQuery1[queryPos];
		fV_QInClosestNewGapQuery[queryPos] = fV_QInClosestNewGapQuery1[queryPos];
	}

	if (maskgapRegion2[dbPos]==1) {fV_QInClosestNewGapQuery[queryPos] = 1;}
        if (maskgapRegion1[queryPos]==1) {fV_DbInClosestNewGapQuery[queryPos] = 1;}


/**** Assign the flags coming from the closest NewGap in Db ******/

	if (flagNewGapQuery_Rev == 1) {
		fV_QInClosestNewGapDb[queryPos] = fV_RepeatOpenGapQuery[queryPos+1];
		fV_DbInClosestNewGapDb[queryPos] = fV_RepeatOpenGapDb[queryPos+1];
	} else {
		fV_QInClosestNewGapDb[queryPos] = fV_QInClosestNewGapDb[queryPos+1];
		fV_DbInClosestNewGapDb[queryPos] = fV_DbInClosestNewGapDb[queryPos+1];
	}

	if (maskgapRegion1[queryPos]==1) {fV_DbInClosestNewGapDb[queryPos] = 1;}
	if (maskgapRegion2[dbPos]==1) {fV_QInClosestNewGapDb[queryPos] = 1;}

       /*test two alternatives*/
       if (newScore < prevScoreGapQuery) {
         newScore = prevScoreGapQuery;
	fV_RepeatOpenGapQuery[queryPos] = fV_QInClosestNewGapDb[queryPos];
	fV_RepeatOpenGapDb[queryPos] = fV_DbInClosestNewGapDb[queryPos];
       }


       if (newScore < continueGapScore) {
         newScore = continueGapScore;
	fV_RepeatOpenGapDb[queryPos] = fV_DbInClosestNewGapQuery[queryPos];
	fV_RepeatOpenGapQuery[queryPos] = fV_QInClosestNewGapQuery[queryPos];
       }

       prevScoreNoGapQuery = scoreVector.noGap[queryPos];
       scoreVector.noGap[queryPos]= newScore;
       scoreVector.gapExists[queryPos] = continueGapScore;

       if (newScore > bestScore) {
         bestScore = newScore;
         bestDbPos = dbPos;
         bestQueryPos = queryPos;
       }

       if (bestScore >= score) break;

     }
     if (bestScore >= score)  break;

                /***** Reassignments in "the array stack" of flags for previous rows of SW matrix *************/

	for (queryPos=1; queryPos<=queryLength; queryPos++) {
		fV_RepeatOpenGapQuery2[queryPos] = fV_RepeatOpenGapQuery1[queryPos];
		fV_RepeatOpenGapQuery1[queryPos] = fV_RepeatOpenGapQuery[queryPos];
		fV_RepeatOpenGapDb2[queryPos] = fV_RepeatOpenGapDb1[queryPos];
		fV_RepeatOpenGapDb1[queryPos] = fV_RepeatOpenGapDb[queryPos];

		fV_QInClosestNewGapDb1[queryPos] = fV_QInClosestNewGapDb[queryPos];
		fV_DbInClosestNewGapQuery1[queryPos] = fV_DbInClosestNewGapQuery[queryPos];

		fV_DbInClosestNewGapDb1[queryPos] = fV_DbInClosestNewGapDb[queryPos];
		fV_QInClosestNewGapQuery1[queryPos] = fV_QInClosestNewGapQuery[queryPos];
	}

   }

   free(scoreVector.noGap);
   free(scoreVector.gapExists);

   if (bestScore < 0)
     bestScore = 0;

   Start1 = bestQueryPos;
   Start2 = bestDbPos;

   scoreGivenEnd = bestScore;

   return(bestScore);
}

/* Traces back the best alignment path using tracebackDir[][], flagNewGapDb[][] and flagNewGapQuery[][];
output is the set of arrays: aligned portions of the aseq... arrays, with gaps inserted,
scores for each position in alignment,
flags for positive matches,
positions in the initial alignments that are aligned
*/

	void **traceback_outputPos(int start_ali1, int start_ali2, int end_ali1, int end_ali2, int **tracebackDir, int **flagNewGapQuery, int **flagNewGapDb, int *apos1, int *apos2)
{
	int pos1, pos2, posGapped, dir, i, j;
	char **aseqGapTrInt1, **aseqGapTrInt2;
	int *positiveInt, *apos1Int, *apos2Int;
	int **col_scoreInt;
	int gapOpen, gapExtend, newGapCost, colScore, d0, d1, d2, d3;
	int sctrl;
	int ascore[10];
	int jnogp1, jnogp2;

	int flagRepeatOpenGap1, flagRepeatOpenGap2;

	positiveInt = ivector(0,alilen_mat1+alilen_mat2);

	apos1Int = ivector(0,alilen_mat1+alilen_mat2);
	apos2Int = ivector(0,alilen_mat1+alilen_mat2);

	col_scoreInt = imatrix(0,alilen_mat1+alilen_mat2, 0,9);

	gapOpen = gap__open;

	sctrl = 0;
	segment_len = 0;
	flagRepeatOpenGap1 = 0;
	flagRepeatOpenGap2 = 0;
	pos1 = end_ali1;
	pos2 = end_ali2;
	posGapped = alilen_mat1+alilen_mat2;


/*** TraceBackDir: 6 - positive match; 5 - non-positive match;
1 - previousScoreGapQuery wins, the gap is new; 2- previousScoreGapQuery wins, the gap is extended from existing;
3 - continueGapScore wins, the gap is new; 4- continueGapScore wins, the gap is extended from existing;
***/

	do {
		dir = tracebackDir[pos1][pos2];
		if (dir==3) {

			positiveInt[posGapped]=0;

			apos1Int[posGapped] = 0;
			apos2Int[posGapped] = apos_filtr2[pos2];

			gapExtend = rint(g_e1(pos2,b));
			newGapCost = gapOpen + gapExtend;
			colScore = -newGapCost;
		        if (maskgapRegion2[pos2]==1 && flagRepeatOpenGap1==0) {
				colScore += gapOpen;
				flagRepeatOpenGap1 = 1;
		        }
			sctrl += colScore;
			if (sctrl < 0) sctrl = 0;

			ScoreOverColumn (colScore, fV_QInClosestNewGapQuery[pos1][pos2], fV_DbInClosestNewGapDb[pos1][pos2], fV_RepeatOpenGapDb[pos1][pos2], fV_RepeatOpenGapQuery[pos1][pos2], flagRepeatOpenGap1, flagRepeatOpenGap2, ascore);
			for (i=0; i<=9; i++) 	col_scoreInt[posGapped][i] = ascore[i];

			pos2--;
			posGapped--;
			segment_len ++;

		}

		if (dir==4) {
			do {
			positiveInt[posGapped]=0;

			apos1Int[posGapped] = 0;
			apos2Int[posGapped] = apos_filtr2[pos2];

		        gapExtend = rint(g_e1(pos2,b));
		        colScore = -gapExtend;
		        if (maskgapRegion2[pos2]==1 && maskgapRegion2[pos2-1]==0 && flagRepeatOpenGap1 == 0) {
		                colScore += gapOpen;
		                flagRepeatOpenGap1 = 1;
		        }

			sctrl += colScore;
			if (sctrl < 0) sctrl = 0;

			ScoreOverColumn (colScore, fV_QInClosestNewGapQuery[pos1][pos2], fV_DbInClosestNewGapDb[pos1][pos2], fV_RepeatOpenGapDb[pos1][pos2], fV_RepeatOpenGapQuery[pos1][pos2], flagRepeatOpenGap1, flagRepeatOpenGap2, ascore);
			for (i=0; i<=9; i++) 	col_scoreInt[posGapped][i] = ascore[i];

			pos2--;
			posGapped--;
			segment_len ++;

			} while (flagNewGapDb[pos1][pos2]!= 1);

			positiveInt[posGapped]=0;

			apos1Int[posGapped] = 0;
			apos2Int[posGapped] = apos_filtr2[pos2];

		        gapExtend = rint(g_e1(pos2,b));
		        newGapCost = gapOpen + gapExtend;
		        colScore = -newGapCost;
		        if (maskgapRegion2[pos2]==1 && flagRepeatOpenGap1==0) {
		                colScore += gapOpen;
		                flagRepeatOpenGap1 = 1;
		        }

			sctrl += colScore;
			if (sctrl < 0) sctrl = 0;

			ScoreOverColumn (colScore, fV_QInClosestNewGapQuery[pos1][pos2], fV_DbInClosestNewGapDb[pos1][pos2], fV_RepeatOpenGapDb[pos1][pos2], fV_RepeatOpenGapQuery[pos1][pos2], flagRepeatOpenGap1, flagRepeatOpenGap2, ascore);
			for (i=0; i<=9; i++) 	col_scoreInt[posGapped][i] = ascore[i];

			pos2--;
			posGapped--;
			segment_len ++;

		}

		if (dir==1) {

			positiveInt[posGapped]=0;

			apos1Int[posGapped] = apos_filtr1[pos1];
			apos2Int[posGapped] = 0;

		        gapExtend = rint(g_e2(pos1,b));
		        newGapCost = gapOpen + gapExtend;
		        colScore = -newGapCost;
		        if (maskgapRegion1[pos1]==1 && flagRepeatOpenGap2==0) {
		                colScore += gapOpen;
		                flagRepeatOpenGap2 = 1;
		        }

			sctrl += colScore;
			if (sctrl < 0) sctrl = 0;

			ScoreOverColumn (colScore, fV_QInClosestNewGapQuery[pos1][pos2], fV_DbInClosestNewGapDb[pos1][pos2], fV_RepeatOpenGapDb[pos1][pos2], fV_RepeatOpenGapQuery[pos1][pos2], flagRepeatOpenGap1, flagRepeatOpenGap2, ascore);
			for (i=0; i<=9; i++) 	col_scoreInt[posGapped][i] = ascore[i];

			pos1--;
			posGapped--;
			segment_len ++;
		}

		if (dir==2) {
		do {

			positiveInt[posGapped]=0;

			apos1Int[posGapped] = apos_filtr1[pos1];
			apos2Int[posGapped] = 0;

		        gapExtend = rint(g_e2(pos1,b));
		        colScore = -gapExtend;
		        if (maskgapRegion1[pos1]==1 && maskgapRegion1[pos1-1]==0 && flagRepeatOpenGap2 == 0) {
		                colScore += gapOpen;
		                flagRepeatOpenGap2 = 1;
		        }

			sctrl += colScore;
			if (sctrl < 0) sctrl = 0;

			ScoreOverColumn (colScore, fV_QInClosestNewGapQuery[pos1][pos2], fV_DbInClosestNewGapDb[pos1][pos2], fV_RepeatOpenGapDb[pos1][pos2], fV_RepeatOpenGapQuery[pos1][pos2], flagRepeatOpenGap1, flagRepeatOpenGap2, ascore);
			for (i=0; i<=9; i++) 	col_scoreInt[posGapped][i] = ascore[i];

			pos1--;
			posGapped--;
			segment_len ++;

		} while (flagNewGapQuery[pos1][pos2] != 1 && sctrl<score);

			positiveInt[posGapped]=0;

			apos1Int[posGapped] = apos_filtr1[pos1];
			apos2Int[posGapped] = 0;

		        gapExtend = rint(g_e2(pos1,b));
		        newGapCost = gapOpen + gapExtend;
		        colScore = -newGapCost;
		        if (maskgapRegion1[pos1]==1 && flagRepeatOpenGap2==0) {
		                colScore += gapOpen;
		                flagRepeatOpenGap2 = 1;
		        }

			sctrl += colScore;
			if (sctrl < 0) sctrl = 0;

			ScoreOverColumn (colScore, fV_QInClosestNewGapQuery[pos1][pos2], fV_DbInClosestNewGapDb[pos1][pos2], fV_RepeatOpenGapDb[pos1][pos2], fV_RepeatOpenGapQuery[pos1][pos2], flagRepeatOpenGap1, flagRepeatOpenGap2, ascore);
			for (i=0; i<=9; i++) 	col_scoreInt[posGapped][i] = ascore[i];

			pos1--;
			posGapped--;
			segment_len ++;

		}

		if (dir==5) {

			positiveInt[posGapped]=0;

			apos1Int[posGapped] = apos_filtr1[pos1];
			apos2Int[posGapped] = apos_filtr2[pos2];

			colScore = rint(score_matrix[pos1][pos2]*score_scale - sgapfcn(pos1,pos2,b));

			sctrl += colScore;
			if (sctrl < 0) sctrl = 0;

		        if (maskgapRegion2[pos2] == 0) flagRepeatOpenGap1 = 0;
		        if (maskgapRegion1[pos1] == 0) flagRepeatOpenGap2 = 0;

			pos2--;
			pos1--;
			posGapped--;

			segment_len++;
		}

		if (dir==6) {

			positiveInt[posGapped]=1;

			apos1Int[posGapped] = apos_filtr1[pos1];
			apos2Int[posGapped] = apos_filtr2[pos2];

			colScore = rint(score_matrix[pos1][pos2]*score_scale - sgapfcn(pos1,pos2,b));

			sctrl += colScore;
			if (sctrl < 0) sctrl = 0;

			ScoreOverColumn (colScore, fV_QInClosestNewGapQuery[pos1][pos2], fV_DbInClosestNewGapDb[pos1][pos2], fV_RepeatOpenGapDb[pos1][pos2], fV_RepeatOpenGapQuery[pos1][pos2], flagRepeatOpenGap1, flagRepeatOpenGap2, ascore);
			for (i=0; i<=9; i++) 	col_scoreInt[posGapped][i] = ascore[i];

		        if (maskgapRegion2[pos2] == 0) flagRepeatOpenGap1 = 0;
		        if (maskgapRegion1[pos1] == 0) flagRepeatOpenGap2 = 0;

			pos2--;
			pos1--;
			posGapped--;

			segment_len++;
		}


/***	} while ((pos1>=start_ali1) && (pos2>=start_ali2)); ***/
	} while (sctrl<score && pos1>0 && pos2>0);

	posGp = posGapped+1;
	jnogp1 = jnogp2 = 1;
	for (j=posGp;j<posGp+segment_len;j++) {
		positive[j-posGp+1] = positiveInt[j];

		apos1[j-posGp+1] = apos1Int[j];
		apos2[j-posGp+1] = apos2Int[j];

		for (i=0;i<=9; i++) col_score[i][j-posGp+1] = col_scoreInt[j][i];
	}

	free_ivector(positiveInt, 0,alilen_mat1+alilen_mat2);
	free_imatrix(col_scoreInt, 0,alilen_mat1+alilen_mat2, 0,9);
	free_ivector(apos1Int, 0,alilen_mat1+alilen_mat2);
	free_ivector(apos2Int, 0,alilen_mat1+alilen_mat2);

}

/* Computes the score for column match using "beta-function" formula */
	int ScoreForTwoRows(double *subjectRow, double *queryRow)
	{
	double effnum_subj, effnum_query; /* Total counts in each row (N1 and N2) */
	double sum_mu; /* Variable used in the calculation of P */
	double p, pnorm, p_back, lgodd; /* Logarithms of : the probability P, normalized P and the background probability */
/*	double *p_dayhoff; /* background probabilities of amino acid occurrence (from Whelan & Goldman 2001) */
	int i_lgodd; /* Returned logg-odd score*/
	int i;
	double gamma1, gamma2, gamma3;

	effnum_subj = effnum_query = 0;

	for(i=1;i<=20;i++) {
		effnum_subj += subjectRow[i];
		effnum_query += queryRow[i];
	}

	sum_mu = effnum_subj+effnum_query+21;
/*	p = lgamma(effnum_subj+1) + lgamma(effnum_query+1) - lgamma(sum_mu);
	p_back = lgamma(effnum_subj+1) + lgamma(effnum_query+1);  */
	lgodd = -lgamma(sum_mu);

	for (i=1; i<=20; i++) {
		gamma1 = lgamma(subjectRow[i]+queryRow[i]+1.0) ;
/******		gamma2 = lgamma(subjectRow[i]+1.0);
		gamma3 = lgamma(queryRow[i]+1.0);
		p += gamma1 - gamma2 - gamma3;
		p_back += (subjectRow[i]+queryRow[i])*log(p_dayhoff[i-1]) - gamma2 - gamma3;  *****/
		lgodd += -(subjectRow[i]+queryRow[i])*log(p_dayhoff[i-1]) + gamma1;

/*		fprintf(stderr, "ln p_%d = %e    ln p_back_%d = %e   lgodd_%d = %e \n", i, p, i, p_back, i, lgodd);
		fprintf(stderr, "subject_row[%d] = %e, query_row[%d] = %e", i, subjectRow[i], i, queryRow[i]);
		fprintf(stderr, "gamma1_%d = %e ; gamma2_%d = %e ; gamma3_%d = %e \n", i, gamma1, i, gamma2, i, gamma3);  */

	}

/*	pnorm = p + lgamma(21);  */
/*	fprintf(stderr, "p= %e    pnorm= %e \n", p, pnorm);
	fprintf(stderr, "p_back= %e \n", p_back); */

	lgodd += lgamma(21);
/*	fprintf(stderr, "lgodd = %e \n", lgodd);  */

	i_lgodd = rint(lgodd);
	return(i_lgodd);
}


double ScoreForTwoRows_Model6(int pos1, int pos2, double score_scale, double b)
{
	int i, k1, k2;
	double s;
	double ngap1, ngap2, g1, g2;
	s=0.0;

	if ((sum_eff_let1[pos1]==1.0) && (sum_eff_let2[pos2] == 1.0)) {
		k1=k2=1;
		for (i=1;i<=20;i++) {
			if (matrix1[pos1][i]!= 0.0) k1=i;
			if (matrix2[pos2][i]!= 0.0) k2=i;
		}
		s = smatrix[k1][k2];
	} else {
		for (i=1;i<=20;i++) {
		s+= n_effAa1[pos1][i]*(sum_eff_let2[pos2]-1)/sum_eff_let1[pos1]*log(pseudoCnt2[pos2][i]/p_rbnsn[i-1]) + n_effAa2[pos2][i]*(sum_eff_let1[pos1]-1)/sum_eff_let2[pos2]*log(pseudoCnt1[pos1][i]/p_rbnsn[i-1]);
		}
/***		s= s/(sum_eff_let1[pos1]+sum_eff_let2[pos2]-2);   ***/
	}
	s = s*score_scale/lambda_u;

	g1 = n_effAa1[pos1][0]/(n_effAa1[pos1][0]+sum_eff_let1[pos1]);
	g2 = n_effAa2[pos2][0]/(n_effAa2[pos2][0]+sum_eff_let2[pos2]);
	s -= b*((1-g1)*g2 + (1-g2)*g1);

	return s;
}

double ScoreForTwoRows_smat1q_2d(int pos1, int pos2) /* 1 serves as query; 2 serves as database */
{
        int i, k1, k2;
        double s;
        double ngap1, ngap2, g1, g2;
        double comp1, comp2;
        double s1,s2;
        s1=s2=0.0;

        comp1 = sum_eff_let1[pos1]-1.0;
        comp2 = sum_eff_let2[pos2]-1.0;

	//fprintf(stdout, "\n# ");
        for (i=1;i<=20;i++) {

		//fprintf(stdout, "%9.3f%9.3f", n_effAa1[pos1][i], pseudoCnt2[pos2][i]);
                s1 +=  n_effAa1[pos1][i]*log(pseudoCnt2[pos2][i]/p_rbnsn[i-1]);
        }
	//fprintf(stdout, "\n");

	s1 = s1/sum_eff_let1[pos1];

        s = s1/lambda_u;

        return s;
}

double ScoreForTwoRows_smat3_21(int pos1, int pos2)
{
	int i, k1, k2;
	double s;
	double ngap1, ngap2, g1, g2;
	double comp1, comp2;
	double s1,s2;
	s1=s2=0.0;

	comp1 = sum_eff_let1[pos1]-1.0;
	comp2 = sum_eff_let2[pos2]-1.0;

	if (abs(comp1)<1e-5 && abs(comp2)<1e-5) {
		k1=k2=1;
		for (i=1;i<=20;i++) {
			if (n_effAa1[pos1][i]!= 0.0) k1=i;
			if (n_effAa2[pos2][i]!= 0.0) k2=i;
		}

		s = smatrix[k1][k2];
	} else {
		for (i=1;i<=20;i++) {

		s1 +=  n_effAa1[pos1][i]*log(pseudoCnt2[pos2][i]/p_rbnsn[i-1]);
		s2 +=  n_effAa2[pos2][i]*log(pseudoCnt1[pos1][i]/p_rbnsn[i-1]);

/*		s+= n_effAa1[pos1][i]*(sum_eff_let2[pos2]-1)/sum_eff_let1[pos1]*log(pseudoCnt2[pos2][i]/p_rbnsn[i-1]) +
n_effAa2[pos2][i]*(sum_eff_let1[pos1]-1)/sum_eff_let2[pos2]*log(pseudoCnt1[pos1][i]/p_rbnsn[i-1]);
*/
		}

/*** Do normalization of s: ***/
		s = s1*(sum_eff_let2[pos2]-1)/sum_eff_let1[pos1] + s2*(sum_eff_let1[pos1]-1)/sum_eff_let2[pos2];
		s= s/(sum_eff_let1[pos1]+sum_eff_let2[pos2]-2);

	}
	s = s/lambda_u;

	return s;
}

double ScoreForTwoRows_smat3_22(int pos1, int pos2)
{
	int i, k1, k2;
	double s;
	double ngap1, ngap2, g1, g2;
	double comp1, comp2;
	s=0.0;

	comp1 = sum_eff_let1[pos1]-1.0;
	comp2 = sum_eff_let2[pos2]-1.0;
	if (abs(comp1)<1e-5 && abs(comp2)<1e-5) {
		k1=k2=1;
		for (i=1;i<=20;i++) {
			if (matrix1[pos1][i]!= 0.0) k1=i;
			if (matrix2[pos2][i]!= 0.0) k2=i;
		}
		s = smatrix[k1][k2];
	} else {
		for (i=1;i<=20;i++) {
		s+= n_effAa1[pos1][i]*(sum_eff_let2[pos2]-1)/sum_eff_let1[pos1]*log(pseudoCnt2[pos2][i]/p_rbnsn[i-1]) + n_effAa2[pos2][i]*(sum_eff_let1[pos1]-1)/sum_eff_let2[pos2]*log(pseudoCnt1[pos1][i]/p_rbnsn[i-1]);
		}

/*** No normalization of s ***/
	}
	s = s/lambda_u;

	return s;
}

double ScoreForTwoRows_smat3_23(int pos1, int pos2)
{
	int i, k1, k2;
	double s;
	double ngap1, ngap2, g1, g2;
	double comp1, comp2;
	s=0.0;

	comp1 = sum_eff_let1[pos1]-1.0;
	comp2 = sum_eff_let2[pos2]-1.0;
	if (abs(comp1)<1e-5 && abs(comp2)<1e-5) {
		k1=k2=1;
		for (i=1;i<=20;i++) {
			if (matrix1[pos1][i]!= 0.0) k1=i;
			if (matrix2[pos2][i]!= 0.0) k2=i;
		}
		s = smatrix[k1][k2];
	} else {
		for (i=1;i<=20;i++) {
/* No division by the opposite sum_eff_let in each of two terms (formula 3_23): */
		s+= n_effAa1[pos1][i]*(sum_eff_let2[pos2]-1)*log(pseudoCnt2[pos2][i]/p_rbnsn[i-1]) + n_effAa2[pos2][i]*(sum_eff_let1[pos1]-1)*log(pseudoCnt1[pos1][i]/p_rbnsn[i-1]);
		}

/*** Do normalization of s: ***/
		s= s/(sum_eff_let1[pos1]+sum_eff_let2[pos2]-2);
	}
	s = s/lambda_u;

	return s;
}


double ScoreForTwoRows_smat3_27(int pos1, int pos2)
{
	int i, k1, k2;
	double s;
	double ngap1, ngap2, g1, g2;
	double comp1, comp2;

	s=0.0;
	comp1 = sum_eff_let1[pos1]-1.0;
	comp2 = sum_eff_let2[pos2]-1.0;
	if (abs(comp1)<1e-5 && abs(comp2)<1e-5) {
		k1=k2=1;
		for (i=1;i<=20;i++) {
			if (matrix1[pos1][i]!= 0.0) k1=i;
			if (matrix2[pos2][i]!= 0.0) k2=i;
		}
		s = smatrix[k1][k2];
	} else {
		for (i=1;i<=20;i++) {

/* Simplest scoring formula 3_27 */
		s+= n_effAa1[pos1][i]*log(pseudoCnt2[pos2][i]/p_rbnsn[i-1]) + n_effAa2[pos2][i]*log(pseudoCnt1[pos1][i]/p_rbnsn[i-1]);
		}

/*** No normalization of s ***/
	}
	s = s/lambda_u;

	return s;
}

double ScoreForTwoRows_smat3_28(int pos1, int pos2)
{
	int i, k1, k2;
	double s;
	double ngap1, ngap2, g1, g2;
	double comp1, comp2;
	s=0.0;

	comp1 = sum_eff_let1[pos1]-1.0;
	comp2 = sum_eff_let2[pos2]-1.0;
	if (abs(comp1)<1e-5 && abs(comp2)<1e-5) {
		k1=k2=1;
		for (i=1;i<=20;i++) {
			if (matrix1[pos1][i]!= 0.0) k1=i;
			if (matrix2[pos2][i]!= 0.0) k2=i;
		}
		s = smatrix[k1][k2];
	} else {
		for (i=1;i<=20;i++) {
/* NO normalization and no division by the opposite sum_eff_let in each of two terms (formula 3_28): */
		s+= n_effAa1[pos1][i]*(sum_eff_let2[pos2]-1)*log(pseudoCnt2[pos2][i]/p_rbnsn[i-1]) + n_effAa2[pos2][i]*(sum_eff_let1[pos1]-1)*log(pseudoCnt1[pos1][i]/p_rbnsn[i-1]);
		}
/*** No normalization of s ***/
	}
	s = s/lambda_u;

	return s;
}

/* version of S_g = 0 */
double Sgap6_smat_off(int pos1, int pos2, double b)
{
	double sg;
	sg = 0.0;
	return sg;
}

/* computes S_g - reduction of col-col score due to gap content */
double Sgap6_smat(int pos1, int pos2, double b)
{
	double g1, g2, sg;
	g1 = n_effAa1[pos1][0]/(n_effAa1[pos1][0]+sum_eff_let1[pos1]);
	g2 = n_effAa2[pos2][0]/(n_effAa2[pos2][0]+sum_eff_let2[pos2]);
	sg = f*b*((1-g1)*g2 + (1-g2)*g1);
	return sg;
}


/* computes gap extension penalty in 1 depending on gap content in 2 */
double GapExtend1(int pos2, double b)
{
	double ge;
	ge = gap__extend*b*(sum_eff_let2[pos2]/(n_effAa2[pos2][0]+sum_eff_let2[pos2]));
	return ge;
}

/* version not depending on gap content */
double GapExtend1_off(int pos2, double b)
{
	double ge;
	ge = gap__extend;
	return ge;
}

/* computes gap extension penalty in 2 depending on gap content in 1 */
double GapExtend2(int pos1, double b)
{
	double ge;
	ge = gap__extend*b*(sum_eff_let1[pos1]/(n_effAa1[pos1][0]+sum_eff_let1[pos1]));
	return ge;
}

/* version not depending on gap content */
double GapExtend2_off(int pos1, double b)
{
	double ge;
	ge = gap__extend;
	return ge;
}

void sort(int n, double arr[])
{
	unsigned long i,ir=n,j,k,l=1;
	int jstack=0,*istack;
	double a,temp;

	istack=ivector(1,NSTACK);
	for (;;) {
		if (ir-l < M) {
			for (j=l+1;j<=ir;j++) {
				a=arr[j];
				for (i=j-1;i>=1;i--) {
					if (arr[i] <= a) break;
					arr[i+1]=arr[i];
				}
				arr[i+1]=a;
			}
			if (jstack == 0) break;
			ir=istack[jstack--];
			l=istack[jstack--];
		} else {
			k=(l+ir) >> 1;
			SWAP(arr[k],arr[l+1])
			if (arr[l+1] > arr[ir]) {
				SWAP(arr[l+1],arr[ir])
			}
			if (arr[l] > arr[ir]) {
				SWAP(arr[l],arr[ir])
			}
			if (arr[l+1] > arr[l]) {
				SWAP(arr[l+1],arr[l])
			}
			i=l+1;
			j=ir;
			a=arr[l];
			for (;;) {
				do i++; while (arr[i] < a);
				do j--; while (arr[j] > a);
				if (j < i) break;
				SWAP(arr[i],arr[j]);
			}
			arr[l]=arr[j];
			arr[j]=a;
			jstack += 2;
			if (jstack > NSTACK) nrerror("NSTACK too small in sort.");
			if (ir-i+1 >= j-l) {
				istack[jstack]=ir;
				istack[jstack-1]=i;
				ir=j-1;
			} else {
				istack[jstack]=j-1;
				istack[jstack-1]=l;
				l=i;
			}
		}
	}
	free (istack);
}


int *ScoreOverColumn (int colScore, int flag1, int flag2, int flag3, int flag4, int flag5, int flag6, int *column_score)
{
	int d0, d1, d2, d3;
			if (colScore>=0) d0 = 1;
			else d0 = 0;
			d1 = abs(colScore)/100;
			d2 = abs(colScore)/10 - 10*d1;
			d3 = abs(colScore) - 100*d1 - 10*d2;
			column_score[0] = d0;
			column_score[1] = d1;
			column_score[2] = d2;
			column_score[3] = d3;

                        column_score[4] = flag1;
                        column_score[5] = flag2;
                        column_score[7] = flag3;
                        column_score[6] = flag4;

                        column_score[8] = flag5;
                        column_score[9] = flag6;

}

/* Reads reference FSSP alignment, creates the array of positions of capital letter pairs */
void *ReadRef (char *inputfile)
{
/*	char *path = "/home/sadreyev/fssp_test/db_id/pairaln/"; */
	int i, pos, pos1, pos2;

/*	strcat (path, inputfile);
	readali(path);
*/
	readali (inputfile);

	aposref1 = ivector (0,alilen);
	aposref2 = ivector (0,alilen);
	pos=pos1=pos2=0;
	for (i=0; i<alilen; i++) {
			if (isalpha(aseq[0][i])) pos1++;
			if (isalpha(aseq[1][i])) pos2++;
			if (isupper(aseq[0][i]) && isupper(aseq[1][i])) {
				if (pos==0) {start_ref1 = pos1; start_ref2 = pos2;}
				pos++;
				aposref1[pos]=pos1;
				aposref2[pos]=pos2;
				end_ref1 = pos1; end_ref2 = pos2;
			}
		}

	reflen_nogp = pos;
}

int CompareAlnVsReferenceAln (int *apos1, int *apos2, int *aposref1, int *aposref2, int start_ref1, int start_ref2, int end_ref1, int end_ref2 /*, int coverage1, int coverage2, int accuracy1, int accuracy2 */)
{
	int i,j,k;
	int start_ali1, start_ali2, end_ali1, end_ali2;
	int len_common1, len_common2;
/*** Starting/ending positions of intersections	between ali1,2 and ref1,2 : ****/
	int start_common1, end_common1, start_common2, end_common2;
/*** Indexes of elements closest to the starting/ending points of intersections between ali1,2 and ref1,2 : ****/
	int ind_ali1_startc1, ind_ali1_endc1, ind_ali2_startc2, ind_ali2_endc2;
	int ind_ref1_startc1, ind_ref1_endc1, ind_ref2_startc2, ind_ref2_endc2;
	int ind_ref_startc, ind_ref_endc;
	double factor = 0.6931472;

/*** Calculation of coverage1 and coverage2.
In parallel to finding common points of starts and ends corresponding to the intersection,
find the indexes in arrays (apos[] and aposref[]) of the elements closest to these points  *****/

	for (i=1; apos1[i]==0; i++);
	start_ali1 = apos1[i];

	for (i=1; apos2[i]==0; i++);
	start_ali2 = apos2[i];

	for (i=segment_len; apos1[i]==0; i--);
	end_ali1 = apos1[i];

	for (i=segment_len; apos2[i]==0; i--);
	end_ali2 = apos2[i];

	if (end_ali1<end_ref1) {
		end_common1 = end_ali1;
		ind_ali1_endc1 = segment_len;
		for (i=reflen_nogp; aposref1[i]>end_common1 && i>1; i--);
		ind_ref1_endc1 = i;
	} else {
		end_common1 = end_ref1;
		ind_ref1_endc1 = reflen_nogp;
		for (i=segment_len; (apos1[i]>end_common1 || apos1[i]==0) && i>1; i--);
		ind_ali1_endc1 = i;
	}

	if (end_ali2<end_ref2) {
		end_common2 = end_ali2;
		ind_ali2_endc2 = segment_len;
		for (i=reflen_nogp; aposref2[i]>end_common2 && i>1; i--);
		ind_ref2_endc2 = i;
	} else {
		end_common2 = end_ref2;
		ind_ref2_endc2 = reflen_nogp;
		for (i=segment_len; (apos2[i]>end_common2 || apos2[i]==0) && i>1 ; i--);
		ind_ali2_endc2 = i;
	}

	if (start_ali1>start_ref1) {
		start_common1 = start_ali1;
		ind_ali1_startc1 = 1;
		for (i=1; aposref1[i]<start_common1 && i<reflen_nogp; i++);
		ind_ref1_startc1 = i;
	} else {
		start_common1 = start_ref1;
		ind_ref1_startc1 = 1;
		for (i=1; (apos1[i]<start_common1 || apos1[i]==0) && i<segment_len ; i++);
		ind_ali1_startc1 = i;
	}

	if (start_ali2>start_ref2) {
		start_common2 = start_ali2;
		ind_ali2_startc2 = 1;
		for (i=1; aposref2[i]<start_common2 && i<reflen_nogp; i++);
		ind_ref2_startc2 = i;
	} else {
		start_common2 = start_ref2;
		ind_ref2_startc2 = 1;
		for (i=1; (apos2[i]<start_common2 || apos2[i]==0) && i<segment_len ; i++);
		ind_ali2_startc2 = i;
	}


	len_common1 = end_common1 - start_common1 + 1;
	len_common2 = end_common2 - start_common2 + 1;

	if (len_common1<=0 || len_common2<=0) {coverage1 = coverage2 = accuracy1 = accuracy2 = 0.0; return;}



	coverage1 = 1.0*(len_common1 + len_common2)/((end_ref1 - start_ref1 + 1) + (end_ref2 - start_ref2 + 1));
	coverage2 = 0.5*len_common1/(end_ref1 - start_ref1 + 1) + 0.5*len_common2/(end_ref2 - start_ref2 + 1);
	falsecov = 1.0*((end_ali1-start_ali1+1) - len_common1 + (end_ali2-start_ali2+1) - len_common2)/ (len_common1+len_common2);


/*** Calculation of accuracy1 and 2 ***/

	if (ind_ref1_startc1 < ind_ref2_startc2) ind_ref_startc = ind_ref1_startc1;
	else ind_ref_startc = ind_ref2_startc2;

	if (ind_ref1_endc1 > ind_ref2_endc2) ind_ref_endc = ind_ref1_endc1;
	else ind_ref_endc = ind_ref2_endc2;


	accuracy1 = accuracy2 = 0.0;
/*	Let us do it simpler and search not throught the intersection regions but through all aposref1,2 and
apos1,2. So let us temporarily change the cycle limits.

	for (i=ind_ref_startc; i<=ind_ref_endc; i++) {
		for (j=ind_ali1_startc1; j<=ind_ali1_endc1; j++) {
*/
	for (i=1; i<=reflen_nogp; i++) {
		for (j=1; j<=segment_len; j++) {

			if (apos1[j]>aposref1[i]) break;
			if (apos1[j] == aposref1[i]) {
/*				for (k=ind_ali2_startc2; k<=ind_ali2_endc2; k++) {
*/
				for (k=1; k<=segment_len; k++) {
					if (apos2[k]>aposref2[i]) break;
					if (apos2[k] == aposref2[i]) {
						if (k==j) {
							accuracy1 += 1.0;
							 accuracy2 +=1.0;
/**							 fprintf (stderr, "\n match: posref1=%d\n",aposref1[i]); **/
							 }
						else { accuracy2 += exp(-factor*(abs(k-j))); }
					}
				}
			}
		}
	}

	accuracy1 /= ind_ref_endc - ind_ref_startc + 1.0;
	accuracy2 /= ind_ref_endc - ind_ref_startc + 1.0;

}

