function [result_SEK elapsedtime_READ elapsedtime_REF elapsedtime_algo] = SEK_implementation(k_SEK,READ_seq,REF_seq, NoOfSpecies)



k = k_SEK;  % Note that we used k= 4
I=100; % Maximum allowable iteration in OMP^{+,1}    
Nu=0.00001; % allowable tolerance in l1 norm computation for OMP^{+,1} 


%% Generation of k-mers from reads (measurement); 
% the output kmers are stored in kmer_vector_reads_SEK.mat (the main variable is kmer_vector_reads)
tic;
tstart_READ = tic;
generate_read_kmer_vector_SEK(READ_seq, k);
elapsedtime_READ = toc(tstart_READ);

disp('Elapsed time to generate k-mers for reads'); elapsedtime_READ


%% Generation of k-mers from reference; 
% the output k-mers are stored in kmer_vector_allref_SEK.mat (the main variable is all_species_kmer_vector)
tic
tstart_REF = tic;
generate_ref_kmer_vector_SEK(REF_seq, k);
elapsedtime_REF = toc(tstart_REF);

disp('Elapsed time to generate k-mers for reference'); elapsedtime_REF


%% Compilation of reference k-mers; this provides correct indexing 
% the output is stored in kmer_vector_allref_trans_SEK.mat (the main variables are all_species_kmer_trans and fragment2species)
allref_vector_compile_SEK('kmer_vector_allref_SEK.mat');


%% SEK (using OMP^{+,1}) execution

% Loading the relevant mat files
load('kmer_vector_reads_SEK.mat'); load('kmer_vector_allref_trans_SEK.mat'); 
X = all_species_kmer_trans;   
Mu = mean(kmer_vector_reads,1)';

tstart_algo=tic;
[gamma ~]= OMP_plus_1(X,Mu,Nu,I); gamma = gamma / sum(gamma);

result_SEK = zeros(1,NoOfSpecies);

for i =1:length(gamma)
    if gamma(i) ~= 0, result_SEK(fragment2species(i)) = result_SEK(fragment2species(i)) + gamma(i); end
end
elapsedtime_algo=toc(tstart_algo);

disp('Elapsed time to execute OMP^{+,1} algorithm:'); elapsedtime_algo

end