%% General pipeline for assmebly
% Here we load data, clean it up, apply masks and try several variations of
% pairwise comparison experiments, we can vary the window "w", stretch
% factors, or just the subset of data we do this on...
% Main output from here is the final barcodeGen you are actually using and
% the oS (overlap Structure).

% 1. Load raw kymographs
% 2. Merge neighboring barcodes
% 3. Masking of Peaks (extreme values) and Repetitive Finder
% 4. Pairwise comparison
% 5. Parametrization.m (separate script)
% 5. Bargrouping.m (separate script)
% 6. Compare bargroups to reference (separate script)
%% 1. Load Raw Kymos using HCA importing stuff
        % import settings
        import CBT.Hca.Import.import_hca_settings;
        [sets] = import_hca_settings('hca_settings.txt');

          % load all user selected settings
        import CBT.Hca.Settings.get_user_settings;
        sets = get_user_settings(sets);
         
        % timestamp for the results
        sets.timestamp = datestr(clock(), 'yyyy-mm-dd_HH_MM_SS');
 
        % align
        sets.minOverlap =  120;
        sets.maxShift = 20;
        sets.skipPreAlign = 0;
        sets.detPeaks = 1;
    
        % generate barcodes
        sets.maxLen=Inf;
        sets.skipEdgeDetection = 0;
        sets.bitmasking.untrustedPx = 6; % depending on nm/bp
        sets.minLen = max(sets.minOverlap); % dependent on nm/px ratio
        sets.minLen=120;
        % add kymographs

        import CBT.Hca.Import.add_kymographs_fun;
        [kymoStructs] = add_kymographs_fun(sets);

        %  put the kymographs into the structure
        import CBT.Hca.Core.edit_kymographs_fun;
        kymoStructs = edit_kymographs_fun(kymoStructs,sets.timeFramesNr);

        names=(sets.kymosets.filenames).';

        % cut out bitmasks
        for i=1:length(kymoStructs)
            kymoStructs{i}.unalignedBitmask=kymoStructs{i}.unalignedBitmask(1:sets.timeFramesNr,:) ;
        end

        % Kill hca from path here please
        % align kymos
        [kymoStructs, badMasks] = DBM4.align_kymos(kymoStructs, sets, names);


    %     % sets.minLen
    import DBM4.gen_barcodes_from_kymo;
    barcodeGen =  gen_barcodes_from_kymo(kymoStructs, sets,sets.maxLen);

%% Found two duplicates, gonna fish them out
l=0;
for i=1:length(barcodeGen)
names_fish{i,1}=barcodeGen{i}.name;
names_fish{i,2}=length(barcodeGen{i}.rawBarcode)==377 || length(barcodeGen{i}.rawBarcode)==433 ;
if names_fish{i,2}==1 
l=l+1;
names_sus{l,1}=names_fish{i,1};
names_sus{l,2}=i;
end
end
barcodeGen{575}=[];
barcodeGen{574}=[];


index = cellfun(@isempty, barcodeGen) == logical(false);
barcodeGen = barcodeGen(index);

%% 2. Merge neighbor barcodes
% This uses the same system to merge as local alignment
% but here is based in the naming as well, so we know they do belong
% together, we can allow minimal stretching to stitch them

    import DBM4.Bargrouping.merge_neighbor_barcodes;
    sets.minOverlap=100;
    %sF=1;
    sF=0.97:0.01:1.03;
    [barGenMerged, posMulti, cnt_unique] = merge_neighbor_barcodes(barcodeGen,sets.minOverlap,sF);
    disp('Done stitching neighbor barcodes');

 % to check some...
 for i=1:length(barGenMerged)
 
     mergId=i;
 mergPair = 1;
 figure;   
 barmerged_plot(barGenMerged,barcodeGen,mergId,true)
  
   % import Core.plot_match_simple;
   %[f] = plot_match_simple(barGenMerged{mergId}.pairBars{mergPair}, barGenMerged{mergId}.overlapStruct{mergPair},2,1); 

    % Plot of Partial and total overlap vs Overlap.size vs sF 
end

% Ideally barGenMerged becomes your new base barcodeGen
barcodeGen_allmerged=barGenMerged;
        
%% 3. Masking of peaks and Repetitive finder
% This is super important
% Depending on the data one might want to be more or less conservative
% peak and rep masks are created but not forced
% one can also separate into barcodeGens of rep / partiallyrep/ norep
barcodeGen_allmerged = mask_peaks(barcodeGen_allmerged);


[barcodeGen_allmerged_periodic,barcodeGen_allmerged_notperiodic,barcodeGen_allmerged_partiallyperiodic,barcodeGen_allmerged] = periodic_finder(barcodeGen_allmerged);


%% 4. Pairwise comparison - Overlap Structure Experiments
% Experiment 1 OS Experiment overlap 150 nonrepetitives only perioic finder used neighbors merged masked POMBE
% Os settings
timestamp=datestr(clock(), 'yyyy-mm-dd_HH_MM_SS');
sF = 0.85:0.01:1.15;
minOverlap = 150; 
bitmask_buff=5;
masked =0;

% barcodeGen selection many different subdatasets were tested, here we choose one 
 % barcodeGen selection
%  barcodeGen=barcodeGen_periodic_unmerged
%  barcodeGen=barcodeGen_notperiodic_unmerged
%  barcodeGen=barcodeGen_edgy_unmerged
%  barcodeGen=barcodeGen_allunmerged
%  barcodeGen=barcodeGen_periodic_merged
  barcodeGen=barcodeGen_notperiodic_merged;
%  barcodeGen=barcodeGen_edgy_merged
%  barcodeGen=barcodeGen_allmerged

for k=1:length(barcodeGen)
    if isfield(barcodeGen{1},'repmask')
    barcodeGen{k}.rawBitmask= barcodeGen{k}.rawBitmask & ~barcodeGen{k}.repmask;
    end
    if masked ==1 
    barcodeGen{k}.rawBitmask= barcodeGen{k}.rawBitmask & barcodeGen{k}.peakBitmask;
    end
end

barcodeGen = expand_mask(barcodeGen, bitmask_buff);

% kill intermediate unmasked stuff
barcodeGen = keeplongestbitmask(barcodeGen);

% size check
barcodeGen_sizefiltered={};
l=0;
for i=1:length(barcodeGen)

    if sum(diff(find(barcodeGen{i}.rawBitmask==0))>=(minOverlap*max(sF)))>=1  % This weird thing here makes sure all barcodes have a continuous bitmask thats longer that the minOverlap
    l=l+1;
    barcodeGen_sizefiltered{l}=barcodeGen{i};
    end

end
barcodeGen=barcodeGen_sizefiltered;

% OS run
% Big exp
bars = barcodeGen;
tic
[oS] = calc_overlap_mp(bars,sF, minOverlap,timestamp);
toc

% Save oS and bars
save('E:\Spombe\reprocess after bargi\oS_results_pombe_2025\1_oS_norepallmergedpeakmasked_150_extrabuff.mat', ...
    "barcodeGen","oS")


