%%
msafile = 'SAPAN_MSA.fasta';
msafile_type = 'faln';

switch msafile_type 
    case 'faln'
    [~,nmsa] = faln_to_nmsa(msafile);
    case 'aln'        
    [~,nmsa] = aln_to_nmsa(msafile);
end

fcov = 50;   % percentage of total positions in the msa based on which a 
            % corresponding number of top coevolution scores is reported.
npos = size(nmsa,2);
ncov = round(npos*fcov/100);


%%
% nbZPX2
        [nbZPX2] = NMSA_to_fnb4_ZPX2(nmsa);

% PSICOV        
        [PSICOV] = NMSA_to_PSICOV(nmsa,...
            'NOGAPS',0.8,1.0,'ZERO','FRO',...
            0.006,0.0,0.0001,'QUIC','RHO',20);

%DCA        
        nmsa_to_faln(nmsa,'msa.faln');
        % Here we perform DCA using the original code from Andrea Pagnani.
        dca('msa.faln','dca_mat.matrix');
        % Here we read in the dca matrix
        import_dca_matrix('dca_mat.matrix');
        % Here we convert the dca_mat.matrix into our internal matrix
        [ ~,DCA,~,~ ] = dca_to_wMI_DCA( dca_mat,fcov );
        % Here we remove the external files generated at each cycle.
        clear dca_mat
        !rm msa.faln dca_mat.matrix


%%
% Here we view the coevolution matrices before any correction for gaps.
figure;imagesc(nbZPX2);
figure;imagesc(PSICOV);
figure;imagesc(DCA);

%%
% Here we calculate a matrix of weights to correct for the presence of
% gaps.

gapW0 = ones(npos,npos);
gapW1 = correct_coevmat_forgaps(nmsa);
gapW2 = gapW1.^2;
gapW3 = gapW1.^3;

% Here we choose how to corrects the effect of gaps on the coevolution 
% matrices calculated by various methods. Three progressively stronger
% corrections (gapW1,gapW2,gapW3) can be applied.
gW = gapW3;

%%
    
 nbZPX2_orig_1 = nbZPX2;
 PSICOV_orig_1 = PSICOV;
 DCA_orig_1 = DCA;
 
%%
% Here we make sure all the matrices are positive to prevent odd phenomena
% in the sorting.

 nbZPX2_orig_2 = nbZPX2_orig_1 - min(nbZPX2_orig_1(:));
 PSICOV_orig_2 = PSICOV_orig_1 - min(PSICOV_orig_1(:));
 DCA_orig_2 = DCA_orig_1 - min(DCA_orig_1(:));
 
%%

 nbZPX2 = nbZPX2_orig_2.*gW;
 PSICOV = PSICOV_orig_2.*gW;
 DCA = DCA_orig_2.*gW;

%%
figure;imagesc(nbZPX2);
figure;imagesc(PSICOV);
figure;imagesc(DCA);

%% Coevolving positions and scores
        sorted_nbZPX2 = sort_matrix_descend(nbZPX2);
        nbZPX2_covar_vec = sorted_nbZPX2(1:ncov,1:3);            
	%fp = fopen('outputfile', 'w');
	%fprintf(fp,' %g', nbZPX2_covar_vec);
	csvwrite('nbZPX2covarvec.csv',nbZPX2_covar_vec);
        sorted_PSICOV = sort_matrix_descend(PSICOV);
        PSICOV_covar_vec = sorted_PSICOV(1:ncov,1:3);            
	csvwrite('PSICOV_covar_vec.csv',PSICOV_covar_vec);
	%fp1 = fopen(PSICOV_covar_vec_outputfile, 'w');
	%fprintf(fp1, PSICOV_covar_vec);
        sorted_DCA = sort_matrix_descend(DCA);
        DCA_covar_vec = sorted_DCA(1:ncov,1:3);            
	csvwrite('DCA_covar_vec.csv',DCA_covar_vec);
        %fp2 = fopen(DCA_covar_vec_outputfile, 'w');
        %fprintf(fp2, DCA_covar_vec);
