%% Fig 4a, b? Suppl fig 1 suppl Fig 2 , Suppl fig5 , Suppl Table 3
% Run repetitive bcodes into HCA to see where they land
% Try with all reps but only reps

c=0;
barcodeGen_allreps_pombe={};
for i = 1:length(barcodeGen_allmerged_fixed)
if barcodeGen_allmerged_fixed{i}.rep_check==1  
c=c+1;
barcodeGen_allreps_pombe{c}=barcodeGen_allmerged_fixed{i};
barcodeGen_allreps_pombe{c}.original='pombe';
end
end

c=0;
barcodeGen_allreps_scer={};
for i = 1:length(barcodeGen_allmerged_scer)
if barcodeGen_allmerged_scer{i}.rep_check==1  
c=c+1;
barcodeGen_allreps_scer{c}=barcodeGen_allmerged_scer{i};
barcodeGen_allreps_scer{c}.original='scer';
end
end

%barcodeGen_allreps_mix=[barcodeGen_allreps_pombe barcodeGen_allreps_scer];

% now let's global HCA
%theoryStruct=[theoryStruct_newpsf_scer
%   theoryStruct_newpsf];

tokeep=[2 11 20 21]; %chose only these 4 to show, to many gets cluttery
barcodeGen_allreps_mix=barcodeGen_allreps_pombe(tokeep);

% now let's global HCA
theoryStruct= theoryStruct_newpsf;

        import CBT.Hca.Import.import_hca_settings;
       
        [hcaSets] = import_hca_settings('hca_settings.txt');
    % change names to support previous hca
     hcaSets.bitmaskSettings=0;
     hcaSets.timestamp = datestr(clock(), 'yyyy-mm-dd_HH_MM_SS');
     hcaSets.theory.stretchFactors=0.85:0.01:1.15;

    % compare theory to experiment
    barcodeGenC=barcodeGen_allreps_mix;
    import CBT.Hca.Core.Comparison.compare_distance;
    [rezMax,bestBarStretch,bestLength] = compare_distance(barcodeGenC,theoryStruct, hcaSets );
    % Extract data
    comparisonStructAll = rezMax;
    for i=1:length(comparisonStructAll)
        for j=1:length(bestBarStretch{i})
            comparisonStructAll{i}{j}.bestBarStretch = bestBarStretch{i}(j);
            comparisonStructAll{i}{j}.length = bestLength{i}(j);
        end
    end
    import CBT.Hca.Core.Comparison.combine_theory_results;
    [comparisonStruct] = combine_theory_results(theoryStruct, rezMax,bestBarStretch,bestLength);
   
% Best chromosome extract
for n = 1:numel(theoryStruct)
 names{n}=theoryStruct(n).name;
 theoryStruct(n).rawBitmask=ones(length(theoryStruct(n).rawBarcode),1)';
 theoryStruct(n).rawBitmask(1:6)=0;
 theoryStruct(n).rawBitmask(end-6:end)=0;
 theoryStruct(n).rawBitmask=logical(theoryStruct(n).rawBitmask);
end
best_chromosomes=[];
for i = 1:numel(comparisonStruct)
      found = 0; % default if nothing matches
    for j = 1:numel(names)
        if strcmp(comparisonStruct{i}.name, names{j})
            found = j;
            break; % stop after first match
        end
    end   
    best_chromosomes(i) = found; 
end
% hotizontal plot
    font_size=8;
cross_hca=0; %0 no local, % 1local and global % 2 only local
chr_to_plot=0; %(0=all)
%chr_to_plot=17:19;
makethemlines=0; %
layerwidth=4;
threshold=-Inf;%4.5;% -7.88
unmask=0;
coverageline=0;
%whole_genome_coverage(barcodeGen, comparisonStruct, best_chromosomes, theoryStruct,layerwidth, chromosome,cross_hca,font_size,threshold,consensus_line,unmask,hcaSets)

%[f,chrconsensus,b,y]
 [f] = whole_genome_coverage(barcodeGenC, comparisonStruct,best_chromosomes, theoryStruct, layerwidth,chr_to_plot,cross_hca,font_size,threshold,makethemlines,unmask,coverageline);% other arguments are offset, chromosome (0=all), highlight (to highlight local or not ) 1/0
 set(gca,'Color','white')

 %% now  run them masked?
    for k=1:length(barcodeGen_allreps_mix)
     if isfield(barcodeGen_mix{1},'repmask')
     barcodeGen_allreps_mix{k}.rawBitmask= barcodeGen_allreps_mix{k}.rawBitmask & ~barcodeGen_allreps_mix{k}.repmask;
     barcodeGen_allreps_mix{k}.rawBitmask= barcodeGen_allreps_mix{k}.rawBitmask & barcodeGen_allreps_mix{k}.peakBitmask;
     end
    end

        % compare theory to experiment
    barcodeGenC=barcodeGen_allreps_mix;
    import CBT.Hca.Core.Comparison.compare_distance;
    [rezMax,bestBarStretch,bestLength] = compare_distance(barcodeGenC,theoryStruct, hcaSets );
    % Extract data
    comparisonStructAll = rezMax;
    for i=1:length(comparisonStructAll)
        for j=1:length(bestBarStretch{i})
            comparisonStructAll{i}{j}.bestBarStretch = bestBarStretch{i}(j);
            comparisonStructAll{i}{j}.length = bestLength{i}(j);
        end
    end
    import CBT.Hca.Core.Comparison.combine_theory_results;
    [comparisonStruct] = combine_theory_results(theoryStruct, rezMax,bestBarStretch,bestLength);
   
% Best chromosome extract
for n = 1:numel(theoryStruct)
 names{n}=theoryStruct(n).name;
 theoryStruct(n).rawBitmask=ones(length(theoryStruct(n).rawBarcode),1)';
 theoryStruct(n).rawBitmask(1:6)=0;
 theoryStruct(n).rawBitmask(end-6:end)=0;
 theoryStruct(n).rawBitmask=logical(theoryStruct(n).rawBitmask);
end
best_chromosomes=[];
for i = 1:numel(comparisonStruct)
      found = 0; % default if nothing matches
    for j = 1:numel(names)
        if strcmp(comparisonStruct{i}.name, names{j})
            found = j;
            break; % stop after first match
        end
    end   
    best_chromosomes(i) = found; 
end
% hotizontal plot
    font_size=8;
cross_hca=0; %0 no local, % 1local and global % 2 only local
chr_to_plot=0; %(0=all)
%chr_to_plot=17:19;
makethemlines=0; %
layerwidth=4;
threshold=-Inf;%4.5;% -7.88
unmask=0;
coverageline=0;
%whole_genome_coverage(barcodeGen, comparisonStruct, best_chromosomes, theoryStruct,layerwidth, chromosome,cross_hca,font_size,threshold,consensus_line,unmask,hcaSets)

%[f,chrconsensus,b,y]
 [f] = whole_genome_coverage(barcodeGenC, comparisonStruct,best_chromosomes, theoryStruct, layerwidth,chr_to_plot,cross_hca,font_size,threshold,makethemlines,unmask,coverageline);% other arguments are offset, chromosome (0=all), highlight (to highlight local or not ) 1/0
 set(gca,'Color','white')






%% supl fig of merge molecules
  import DBM4.Bargrouping.merge_neighbor_barcodes;
    sets.minOverlap=130;
    %sF=1;
    sF=0.90:0.01:1.10;
    [barGenMerged, posMulti, cnt_unique] = merge_neighbor_barcodes(barcodeGen_tomerge_new_scer,sets.minOverlap,sF);
    disp('Done stitching neighbor barcodes');

%for i=1:length(barGenMerged)

mergId=17;
mergPair = 1;
figure;   
barmerged_plot(barGenMerged,barcodeGen_tomerge_new_scer,mergId,true)
set(gcf,'Color','white')


%% Maybe better Suppl fig of the rep region detection algorithm
minreps = 5; %default, can be omitted,  minimum ammount of contiguous reps to mask
plots=1;%default=0, enabeled just for suppl fig
periodic_finder(barcodeGen_notperiodic_unmerged_scer{4},minreps,plots);
set(gcf,'Color','white')
periodic_finder(barcodeGen_edgyrep{3},minreps,plots);
set(gcf,'Color','white')
periodic_finder(barcodeGen_regrep{3},minreps,plots);
set(gcf,'Color','white')




%% Supl fig of gaps in chrI and chrII

% Go for the outConsensus and kill 1 last barcode for chr II and 2 last
% bcodes for chr I, those where the last added.... without them we get the
% original gaps

% Let's go the long way and extract all the members into individual
% barcodes
cGen{1}=cGen_final_pombe{1};
cGen{1}.outConsensus=cGen{1}.outConsensus(1:end-2,:);
bGen=[];
for i = 1:height(cGen{1}.outConsensus)
bGen{i}.rawBarcode=cGen{1}.outConsensus(i,:);
bGen{i}.rawBitmask=~isnan(bGen{i}.rawBarcode);
end

cGen{1}=cGen_final_pombe{2};
cGen{1}.outConsensus=cGen{1}.outConsensus(1:end-1,:);
bGen=[];
for i = 1:height(cGen{1}.outConsensus)
bGen{i}.rawBarcode=cGen{1}.outConsensus(i,:);
bGen{i}.rawBitmask=~isnan(bGen{i}.rawBarcode);
end







%% Mean cc std cc snr for ALL bargroups


cGenAll=cGen_final;

cGenAll=cGen_final_pombe;

meanCC=[];
varCC=[];
SNR=[];
for ix=1:length(cGenAll)
    nBars = size(cGenAll{ix}.outConsensus,1);
    ccVals = nan(nBars,1);

    for i = 1:nBars
    ccVals(i) = corr( ...
        cGenAll{ix}.rawBarcode(:), ...
        cGenAll{ix}.outConsensus(i,:)', ...
        'rows','complete' ...
    );
    end

    meanCC(ix) = nanmean(ccVals);
    varCC(ix)  = nanvar(ccVals);


    %Option 3 SNR approach?
    cov = sum(~isnan(cGenAll{ix}.outConsensus),1);
    valid = cov > 1;

    cons = nanmean(cGenAll{ix}.outConsensus(:,valid),1);
    signal = nanstd(cons);
    noise  = nanmean(nanstd(cGenAll{ix}.outConsensus(:,valid),0,1));

    SNR(ix) = signal / noise;
end

bargroups_table=table([1:length(meanCC)]',meanCC(:),varCC(:),SNR(:))

%% longest rep?
bgen=barcodeGen_edgyrep;
lens={};
for i =1:length(bgen)
    i
    lens{i}=bgen{i}.name;
    bgen{i}.name
end
 [maxl,which_max]=max(lens)
figure
 plot(bgen{36}.rawBarcode)
% 4556 annd 4557
 bgen_bigrep{2}=bgen{36}.rawBarcode;
