%% Code from here was used for Fig 4a-c
% Fig where we take some bargroups and dissect them and look at the
% alignment

% We did many overlap experiments (different pretreatment to barcodes and different w's
% First we try with 1-5 oS
% Parameters for the run. 
%load('D:\Spombe\reprocess after bargi\oS_results_pombe_2025\parameters\1_mini')
%Pombe exp 1
alphaNu1=alphav
alphaNu2=betanu
alphaN1=alphaN
alphaN2=betaN

sets.minOverlap =150; %

% For consistency check
sets.scDiffSetting = 0.01; % local stretch when adding a barcode to a group
sets.pxDifSetting = 10; % dist of pixels from actual bargroup 
pthresh =  0.05; % false positive control 
showFig = 1; % whether to show output figure

timestamp = datestr(clock(), 'yyyy-mm-dd_HH_MM_SS'); % timestamp for results

[sortedVals, sortedIds,localScore,partialScore,lenA,lenB,partialLength,pvalLocal,pvalLeftOver, pvalCombined, sortedValsBad, sortedIdsBad] = ...
    calculate_sorted_pvals(oS,sets.minOverlap, alphaNu1, pthresh,[],alphaN1,alphaNu2,alphaN2); %localCCStruct

% bargroup plots
% number of minimum amount of bargroup members
nbar=3;
pxCov=1; % coverage per pixel

sortedValsGood = sortedVals(~isnan(sortedVals));
sortedIdsGood  = sortedIds(~isnan(sortedVals));
localScore = localScore(~isnan(sortedVals));

% Main barcode islands script:
import Core.barcode_island_output;
[barsetGen, outConsensus, coverage, consensus, islandElts, islandPx,cGenAll,barcodeIslandsData, barStruct,barIslands,bargroup_meancov,meanCC,varCC,SNR] =...
    barcode_island_output(sortedValsGood,sortedIdsGood, oS, barcodeGen,timestamp,sets.scDiffSetting,sets.pxDifSetting, [],showFig,nbar,pxCov);

accepted = cellfun(@(x) x.accepted, barsetGen.currentData); % all the

%% Choose which bargroups to check for


tokeep=[6,9,10,16,18];


meanCC(tokeep)          
varCC(tokeep)                    
SNR(tokeep)             
w =150;% minimum overlap length
per_stretch=5;
to_run="bargroups"; % nice bargroups only
%to_run="merged"; %barcodes and nice bargroups
%to_run="barcodes"; %barcodes only
bGeni=barcodeGen;
timestamp = datestr(clock(), 'yyyy-mm-dd_HH_MM_SS');
sF = (1-(per_stretch+8)/100):0.01:(1+per_stretch/100); 

bTi=bT_newpsf_pombe;
for i=1:length(bTi)
bTi{i}.rawBitmask = logical(bTi{i}.rawBitmask);
try
bTi{i}.name=theoryStruct_newpsf(i).name;
catch
end
bTi{i}.theory=i;
end
% Expand mask
bitmask_buff=0;
bTi = expand_mask(bTi, bitmask_buff,1); % last argument adds mask (6px) on the edges

for i=1:length(bGeni)
    bGeni{i}.index=i;
end

if (to_run == "bargroups" || to_run == "merged" )
for i=1:length(cGenAll)
    cGenAll{i}.bargroupid=i;
end
mini_cGenAll=cGenAll(tokeep);
if isempty(mini_cGenAll)
mini_cGenAll=cGenAll;
end

if (to_run == "merged")
 [bGeni_grouped,bGeni_trimmed,allIdx,bGeni_excluded]  = barcode_bargroup_cleaner(bGeni, mini_cGenAll);
bGeni=bGeni_grouped;
else
bGeni=mini_cGenAll;
end

end

for i = 1:numel(bGeni)
    bm = bGeni{i}.rawBitmask;
    firstOne = find(bm == 1, 1, 'first');
    lastOne  = find(bm == 1, 1, 'last');
    bm(firstOne:lastOne) = 1;
    bGeni{i}.rawBitmask = bm;
end

cGen =[bGeni bTi];

l=0;
cGen_sizefiltered={};
numTheoryBar=0;
for i=1:length(cGen)

    if max(diff(find([0 cGen{i}.rawBitmask 0]==0)))>=(w*max(sF)) % This weird thing here makes sure all barcodes have a continuous bitmask thats longer that w 
    l=l+1;
    cGen_sizefiltered{l}=cGen{1,i};
    if isfield(cGen_sizefiltered{l},'theory')
    numTheoryBar=numTheoryBar+1;
    end
    end
end

cGen=cGen_sizefiltered;
numExpBar=length(cGen)-numTheoryBar;


%% align to theory
%different w for each cGen element? 
% Problem here is that theories are short!
% maybe add a check on mp_to_struct and make stuff NA if shorter than W
minioS_pombe_variable=[];
numExpBar=numel(bGeni);
numTheoryBar=numel(bTi);
%parpool
tic
parfor c = 1:numel(bGeni)
tGen=[];
w=round(max(diff(find([0 bGeni{c}.rawBitmask 0]==0)))/max(sF)*0.8);
tGen{1}=bGeni{c};
tGen=[tGen bTi];
[tempoS] = calc_overlap_mp(tGen,sF, w,timestamp);
tempoS=tempoS(1,2:end);
minioS_pombe_variable=[minioS_pombe_variable; tempoS];
end
toc
clear tempoS

%% Best match for each chromosome
minioS=minioS_pombe_variable;
unmask=0;
% Let's take best chromosome for each bargroup
% To do this we need to cut down oS to (exp, theory)
% Then for each row, we take the highest collumn
% What do we take? pval? cclocal? ccreminder?
% stouffer score? (not weighted!!!)

%localscore
alllocalscores= reshape([minioS(:).score], size(minioS,1),size(minioS,2)); % from bg_test_1
alllocalscores=alllocalscores(:);
%globalscore
allglobalscores= reshape([minioS(:).fullscore], size(minioS,1),size(minioS,2)); % from bg_test_1
allglobalscores=allglobalscores(:);
w = minioS(1,2).h;
alllenA = reshape([minioS(:).lenA], size(minioS,1),size(minioS,2)); 
alllenA = alllenA(:);
alllenA(alllenA==0) = nan;
alllenB = reshape([minioS(:).lenB], size(minioS,1),size(minioS,2)); % from bg_test_1
alllenB(alllenB==0) = nan;
%partialscoresandlengths
allpartialScore= reshape([minioS(:).partialScore], size(minioS,1),size(minioS,2)); % from bg_test_1
allpartialScore=allpartialScore(:);
allpartialLength= reshape([minioS(:).partialLength], size(minioS,1),size(minioS,2)); % from bg_test_1
allpartialLength=allpartialLength(:);

% p-value for local
import Zeromodel.beta_ev_cdf;
pvalfun = @(x,l1,l2,nuF,w) 1-beta_ev_cdf(x,nuF*w,1,max(alphaN1*2*(max(l1,l2)-w),alphaN2*2*(l1-w+1)*(l2-w+1)),0);

k=length(alllocalscores);
minL = 50;
allpval = nan(1,k);

        for i=1:k
            %allpval(i) = pvalfun(alllocalscores(i),alllenA(i),alllenB(i),0.10); %whats this 0.07 doing?!
            allpval(i) = pvalfun(alllocalscores(i),alllenA(i),alllenB(i),alphaNu2,w);
        end
        %p-pvalue for single
        import Zeromodel.beta_cdf;
        pvalfunLeftover = @(x,l1,nuF) 1-beta_cdf(x,nuF*l1,1,0);
        allpvalLeftOver = nan(1,k);
        for i=1:k
            if allpartialLength(i) > minL
                allpvalLeftOver(i) = pvalfunLeftover(allpartialScore(i),allpartialLength(i),alphaNu1);
            end
        end
%         pval = [];
        allpvalCombined = (norminv(1-allpval)+norminv(1-allpvalLeftOver))/sqrt(2);
        allpvalCombined(isnan(allpvalCombined)) = (norminv(1-allpval(isnan(allpvalCombined))));

allpvalCombined(allpvalCombined==-Inf)=-8;
allpvalCombined(allpvalCombined==Inf)=8;
allscores=allpvalCombined;

%% Best is sorted according to local score (not stouffer)
% maybe ccdif no, but stscore or something else
ccdif_thres=0;
s=[];

stscores=[];
pvalLs=[];
pvalLOs=[];
Vals=[];
ccdifs=[];
best_chromosomes=[];
for r=1:numExpBar
 [sortedVals,sortedIds,pscores,fullScores,overlaplen,partialScore,partialLength,pvalL,pvalLO,stscore]=sorted_scores(minioS(r,:),alphaN1,alphaN2,alphaNu1,alphaNu2);
 [~,maxind]=max(stscore);
 
 %how does this look?
 stscores(r)=max(stscore);
 pvalLs(r)=min(pvalL); 
 pvalLOs(r)=min(pvalLO);
Vals(r)=max(sortedVals);
ccdifs(r)=sortedVals(1)-sortedVals(2);
% WHAT TO DO HERE???
if ccdifs(r)<=ccdif_thres
best_chromosomes(r)= 0 ;% if ccdif < 0 kick them out!
else
best_chromosomes(r)=sortedIds(1); % best local score
%best_chromosomes(r)=sortedIds(maxind); %best stouffer score
end
end



font_size=8;
see_local=0; %0 no local, % 1local and global % 2 only local
chr_to_plot=0 ; %(0=all)
consensus_line=1; % 1 for full , 2 for local
layerwidth=1.1;
threshold=0;%6.5;%4.5;% -7.88
horizontalgap=600;
unmask=0;

%% Example vars 

 cGen=cGen_example_1_5oS;
 minioS= minioS_1_5;
 best_chromosomes=best_chromosomes_1_5_examples;
%% horizontal plot

 [f,chromosomeConsensus_bargroupsonly] = bargroup_detail_vs_reference(cGen, minioS,best_chromosomes,layerwidth,chr_to_plot,see_local,font_size,threshold,consensus_line,unmask,horizontalgap,alphaN1,alphaN2,alphaNu1,alphaNu2);% other arguments are offset, chromosome (0=all), highlight (to highlight local or not ) 1/0

