/*** 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 "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


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 *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 **aname1, **aname2, **aseq1, **aseq2;
int nal1, nal2, nalmerge, alilen1, alilen2, maxalilen,
*astart1, *astart2;

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;

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 prfprfmatrix(int **align1, int alnlength1, int nali1, double **pseudoCnt1)
{
	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);}

	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;
	alilen1 = alnlength1;
	nal1 = nali1;
	if(debug>0) fprintf(stdout, "alcomp2: %d %d %d %d \n", alilen1, alilen2, nal1, nal2);


/* 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);
	neffsForEachCol_maskGapReg(alignment1, nal1, alilen1, ARG_G, ARG_T, n_effAa1, sum_eff_let1, maskgapRegion1, apos_filtr1, &alilen_mat1, &n_eff1);
/*	for(i=1;i<=alilen1;i++){
		printf("i=%3d ",i);
		for(j=0;j<=20;j++) printf("%8.5f ",n_effAa1[i][j]);
		printf("\n");
	}
*/
/* Calculate target frequencies */
//	pseudoCnt1 = dmatrix(1,alilen_mat1, 0, 20);
	pseudoCounts(n_effAa1, n_eff1, alilen_mat1, pseudoCnt1);


	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;

	/*for(i=1;i<=n;i++){
		printf("%3d ",i);
		for(j=1;j<=len;j++) printf("%d ",ali[i][j]);
		printf("\n");
	}*/


	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=1;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);
				//fprintf(stdout, "%d %f\n", k, effnu[k]);
				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]); **/

		}
	}

}

/** 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 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 **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;

	}
}
