close all;
clear all;
clc;
%*************************************************************************
%*  Peter Miščík
%*  Program KompresDNA - kompresia sekvencie
%*  vstup - DNA sekvencia (FASTA, EMBL, GENBANK)
%*  výstup - *.komdan, *.komdnal
%*************************************************************************
reply1 = input('Zvoľte typ formátu DNA sekvencie (FASTA,EMBL,GENBANK):', 's');
switch reply1
    case 'FASTA'
        sequence = input('\nZadajte názov sekvencie s priponou fasta: ','s');
        [header seq] = fastaread (sequence);   % Načítanie DNA sekvencie fasta formatu
        N = length(seq); % Dĺžka sekvencie
        alpha = 'ACGT'; % abeceda sekvencie
    case 'EMBL'
        sequence = input('\nZadajte názov sekvencie s priponou: ','s');
        seq=emblread (sequence, 'SequenceOnly',true); 
        N = length(seq); % načítanie z EMBL
        alpha = 'acgt';
    case 'GENBANK'
        sequence = input('\nZadajte názov sekvencie s priponou: ','s');
        seq=genbankread(sequence); % načítanie z genbank
        seq=seq.Sequence;
        N = length(seq);
        alpha = 'acgt';
end

reply2 = input('\nChcete zadať úsek sekvencie?(A/N): ', 's');
if reply2=='A'
    str1=sprintf('\nZvoľte interval sekvencie od 1 do %d. Od: ',N);
    a=input(str1);
    str2=sprintf('\nDo: ');
    b=input(str2);
    seq=seq(a:b);
    N=length(seq);
end

str=sprintf('\nPrebieha kompresia...');
disp(str)

k=1;    % inicializácia premennych
j=1;
compress='';
compress{1}='';
N=length(seq);
v=seq(1:15);    % vloženie pre prehľadové okno
i=16;          
dlzka=[];
h=1;

while i<=N-14   % algoritmus vyhľadávania
    q=v(1:i-1);
    if i>2^14       % podmienka pre prehľadové okno max. 2^14
        q=v(i-2^14:i-1);
    end
    p=seq(i:i+14);  % vloženie do náhľadového okna
    for t=(1:length(q)-14)
        [R1,R2,R]=findR(q(t:t+14),p);   % funkcia vyhľadávania zhody medzi
        if R>2                          % nahľadového a prehľadového okna
            [R1,R2,R]=findR(fliplr(q(t:t+14)),p); % otočenie sekvencie
            flip='1';
        else
            flip='0';
        end
        switch R
            case 0  % 0 editačných operácií
                if t==1
                    compress{k}=(v(j:i-1));
                end
                if cellfun('length',compress(k))~=0
                    k=k+1;      
                end
                compress{k}=strcat(flip,dec2bin(t,14),'000000000000');  % kódovanie zhody na 27 bitový reťazec
                i=i+14;
                j=i+1;
                break
            case 1  % 1 editačných operácií
                if t==1
                    compress{k}=(v(j:i-1));
                end
                if cellfun('length',compress(k))~=0
                    k=k+1;
                end
                code=encode(p(R1),alpha); % kódovanie zhody na 27 bitový reťazec
                compress{k}=strcat(flip,dec2bin(t,14),dec2bin(R1,4),code,'000000');
                i=i+14;
                j=i+1;
                break
            case 2  % 2 editačných operácií
                if t==1
                    compress{k}=(v(j:i-1));
                end
                if cellfun('length',compress(k))~=0
                    k=k+1;
                end
                code1=encode(p(R1),alpha);
                code2=encode(p(R2),alpha);% kódovanie zhody na 27 bitový reťazec
                compress{k}=strcat(flip,dec2bin(t,14),dec2bin(R1,4),code1,dec2bin(R2,4),code2);
                i=i+14;
                j=i+1;
                break
            otherwise % nenájdenie zhody
                if cellfun('length',compress(k))~=0
                    if compress{k}(1)=='1' || compress{k}(1)=='0'
                        k=k+1;
                    end
                end
                compress{k}=(v(j:i-1));
        end
    end
    v=seq(1:i);
    i=i+1;
    if i>=7315
        ahoj=1;
    end
end

code=encode(seq(i-1:end),alpha); % doplnenie zyšku sekvencie
compress{k+1}=code;

for a=1:length(compress)    % úprava sekvencie na bitový reťazec
    if cellfun('length',compress(a))~=0
        if compress{a}(1)~='1'&&compress{a}(1)~='0'
            code=encode(compress{a},alpha);
            compress{a}=code;
        end
    end
end

CompressionRatio=sum(cellfun('length',compress))/N % kompresný pomer

dlzka = cellfun('length',compress); % príprava súboru na uloženie
i=1;
j=1;
k=0;

dlzky = [];

while i<=length(dlzka)
    popka=dlzka(i);
    if popka==27
        k=k+1;
    else
        if k~=0
            dlzky(j)=(k*10)+27;
            k=0;
            j=j+1;
        end
        dlzky(j)=dlzka(i);
        j=j+1;
    end
i=i+1;    
end

commat=cell2mat(compress);
Z=mod(length(commat),8);
commat=commat(1:end-Z);
if Z==0
    dlzky(end+1)=0;
    dlzky(end+1)=0;
else
    zvysok=commat(end-Z+1:end);
    dlzky(end+1)=bin2dec(zvysok);
    dlzky(end+1)=Z;
end
    commat=reshape(commat,8,length(commat)/8)';
    commat=bin2dec(commat);
    
uloz=input('Zadajte názov súboru pre uloženie: ','s'); % uloženie výst. režazca
uloz1=strcat(uloz,'.komdna');
uloz2=strcat(uloz,'.komdnal');

fid = fopen(uloz1,'wb'); % uloženie súboru dĺžok
fwrite(fid,commat,'uint8');
fclose(fid);

save(uloz2,'dlzky');
str=sprintf('\nUložené súbory sú: %s, %s',uloz1,uloz2);
disp(str)