% This Function applies the block code model on a given genome for a
% specific n amd k values. It plots 4 curves: A- the Translated Curve, B-
% the Hypothetical GeneMark Curve, C- The Hypothetical Glimmer Curve, and
% D- the Nontranslated Curve according to May's paper.

function [A,B,C,D] = BlockCodeModel(n,k,Genome,DistMetric)

% n = 5; 
% k = 2;
% Genome = 1;
% DistMetric = 1;

% Load your genome:
switch Genome
    case 1, load 'MG1655.mat';
    case 2, load 'O157H7.mat';
    case 3, load 'Bacillus_Subtilis.mat';
    case 4, load 'Staphylococcus_Aureus_Mu50.mat';
    case 5, load 'Salmonella_Typhimurium_LT2.mat';
end
% Initializations
InitPos = -30; 
% -30 means that we consider the 30 bases upstream the initiation codon. 
% [b_-30, b_-29, ..., b_-1, b_0, b_1, b_2, ... b_28, b_29].

clc, close all;
[InfoSeq,ParitySeq,codebook] = BlockCode(n,k);      
for i = 1:length(AGL)   
    AnalysisSeq = dna2rna(DNA{1}(AGL{i,2}+InitPos:AGL{i,2}-InitPos-1));
    SeqLength = length(AnalysisSeq);  
    NumValid = SeqLength - n + 1;
    % NumValid = The number of valid comparisons between the analysis 
    % sequence and the codebook. i.e. number of n-mers in the analysis 
    % sequence (For n = 5, InitPos = -30,NumValid = 56).
    PosValid = InitPos + NumValid - 1; 
    % The corresponding number of NumValid considering the negative values.
    % We do not need it here. We just found it as described in the paper.
    % For (5,2) block code, PosValid = 25.
    for j = 1:NumValid
        r{i,j} = AnalysisSeq(j:j+n-1); 
        % r{i,j} = the jth n-bases-long string of the analysis sequence.
        switch DistMetric
            case 1, HD = HamDist(r{i,j},codebook); 
%                   HD = Hamming Distance between r{i,j} and the codebook
            case 2, HD = FreeEnergyDist(r{i,j},codebook);
            case 3, HD = FreeEnergyDistNew(r{i,j},codebook,2);
        end
        MDI(i,j) = min(HD);  
        % MDI = Minmum Distance Information matrix
    end
    {'AGL', n, k, Genome, DistMetric, length(AGL),i}
end
% plot([InitPos:PosValid],sum(MDI)/length(MDI),'*-')
% TranlatedCurve = sum(MDI)/length(MDI);
A = sum(MDI)/length(MDI);



clear MDI;
for i = 2:length(HGMGL)   
    AnalysisSeq = dna2rna(DNA{1}(HGMGL{i,2}+InitPos:HGMGL{i,2}-InitPos-1));
    SeqLength = length(AnalysisSeq);
    NumValid = SeqLength - n + 1;
    PosValid = InitPos + NumValid - 1;
    for j = 1:NumValid
        r{i,j} = AnalysisSeq(j:j+n-1);
        switch DistMetric
            case 1, HD = HamDist(r{i,j},codebook); 
%                   HD = Hamming Distance between r{i,j} and the codebook
            case 2, HD = FreeEnergyDist(r{i,j},codebook);
            case 3, HD = FreeEnergyDistNew(r{i,j},codebook,2);
        end
        MDI(i,j) = min(HD);
    end
    {'HGMGL',  n, k, Genome, DistMetric, length(HGMGL),i}
end
% hold on, plot([InitPos:PosValid],sum(MDI)/length(MDI),'ro-')
% HypotheticalGeneMark = sum(MDI)/length(MDI);
B = sum(MDI)/length(MDI);

clear MDI;
for i = 2:length(HGGL)   
    AnalysisSeq = dna2rna(DNA{1}(HGGL{i,2}+InitPos:HGGL{i,2}-InitPos-1));
    SeqLength = length(AnalysisSeq);
    NumValid = SeqLength - n + 1;
    PosValid = InitPos + NumValid - 1;
    for j = 1:NumValid
        r{i,j} = AnalysisSeq(j:j+n-1);
        switch DistMetric
            case 1, HD = HamDist(r{i,j},codebook); 
%                   HD = Hamming Distance between r{i,j} and the codebook
            case 2, HD = FreeEnergyDist(r{i,j},codebook);
            case 3, HD = FreeEnergyDistNew(r{i,j},codebook,2);
        end
        MDI(i,j) = min(HD);
    end
    {'HGGL',  n, k, Genome, DistMetric,length(HGGL),i}
end
% hold on, plot([InitPos:PosValid],sum(MDI)/length(MDI),'kd-')
% HypotheticalGlimmer = sum(MDI)/length(MDI);
C = sum(MDI)/length(MDI);


clear MDI;
for i = 2:length(NTGL)   
    AnalysisSeq = dna2rna(DNA{1}(NTGL{i,2}+InitPos:NTGL{i,2}-InitPos-1));
    SeqLength = length(AnalysisSeq);
    NumValid = SeqLength - n + 1;
    PosValid = InitPos + NumValid - 1;
    for j = 1:NumValid
        r{i,j} = AnalysisSeq(j:j+n-1);
        switch DistMetric
            case 1, HD = HamDist(r{i,j},codebook); 
%                   HD = Hamming Distance between r{i,j} and the codebook
            case 2, HD = FreeEnergyDist(r{i,j},codebook);
            case 3, HD = FreeEnergyDistNew(r{i,j},codebook,2);
        end
        MDI(i,j) = min(HD);
    end
    {'NTGL',  n, k, Genome, DistMetric,length(NTGL), i}
end
% plot([InitPos:PosValid],sum(MDI)/length(MDI),'gx-')
% NonTranlatedCurve = sum(MDI)/length(MDI);
D = sum(MDI)/length(MDI);



PlotFigure(A,B,C,D,n,k);



% legend('Translated Curve','Hypothetical GeneMark Curve', 'Hypothetical Glimmer Curve','Non-Translated Curve');
% grid
% xlabel('Base Position');
% ylabel('Minimum Hamming Distance');
% title(['Minimum Hamming Distance (dmin) for (',num2str(n),',',num2str(k),') Block Code Model']);
% axis([InitPos,PosValid,ylim])

% figure,
% plot([InitPos:PosValid],A,'*-');
% hold on, plot([InitPos:PosValid],B,'ro-');
% plot([InitPos:PosValid],C,'kd-');
% plot([InitPos:PosValid],D,'gx-');
% legend('Translated Curve','Hypothetical GeneMark Curve', 'Hypothetical Glimmer Curve','Non-Translated Curve');
% grid
% xlabel('Base Position');
% ylabel('Minimum Hamming Distance');
% title(['Minimum Hamming Distance (dmin) for (',num2str(n),',',num2str(k),') Block Code Model']);
% axis([InitPos,PosValid,ylim]);