%% Script used in DOGMA for parametrization of bargrouping (null-model extreme values distribution fitting)
% Input for this is barcodeGen directly after assembly_pipeline_yeast
% Output of this is 'alphaN','alphav','betaN','betanu'
%% Parameters 
minlen=250; 
masked=1;
  barcodeGen=barcodeGen_allmerged_fixed;

for k=1:length(barcodeGen)
    if masked ==1 
    barcodeGen{k}.rawBitmask= barcodeGen{k}.rawBitmask & barcodeGen{k}.peakBitmask;
    end
end

for i = 1:length(barcodeGen)
    old_mask = barcodeGen{i}.rawBitmask;
    new_mask = expand_mask(old_mask, bitmask_buff);
    barcodeGen{i}.rawBitmask = new_mask;
end

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

    if sum(diff(find(barcodeGen{i}.rawBitmask==0))>=(minlen*max(1.1)))>=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;

bGrand = barcodeGen;


lengths = cellfun(@(x) sum(x.rawBitmask),bGrand);
[maxL,idx] = max(lengths);
% as bGrand2, the longest barcode
bGrand2 = bGrand(idx);
bGrand(idx) = [];
% TO choose some random ones not all
try
bGrand = bGrand(1: 100);
catch
try
    bGrand=bGrand(1:50);
catch
end
end

rng("default")

% Parameters to generate random data:
N = 100;% number of random fragments % default 200
pxpsf = 2.72; % psf
mL = 1250; % albertas has it as this
snr = 0.5:0.5:5;
stdE = 0.05; %0.02; % fragment length-rescale factor std 0.02
isC = 0; % whether circular barcode
tL = 200000; % rand length in px % default 10000
stdM = 20; % additive noise mean
sF = 0.9:0.025:1.1; % length re-scaling factor
minL = 90; % minimum length / same as minoverlap?
stdL = 100; %std of length of fragment

% LA <= LB!!!
lA = 100:10:160; % 3 
lB = 250:50:maxL; % will keep % 6
% w < lA2
w = 50:20:150; % My try optimum 150? so 5 or 8 

%lA2 = 500:50:800; % his try
lA2= 170:10:230; % 5, 210

%% Global pcc distribution
% Fix for linearity longVec, shortVec, shortVecBit,longVecBit,k,isLinearTF,numPixelsAroundBestTheoryMask, numBestCoef
comparisonFun = @(x,y,z,w,u) unmasked_MASS_PCC(y,x,z,w,2^(4+nextpow2(length(x))),1,50);

maxPCC = cell(length(lA),length(lB));
osOutput = cell(length(lA),length(lB));
pccs =  cell(length(lA),length(lB)); 
for j = 1:length(lB)
    for k=1:length(lA)
    j
    k
    curLen = lA(k);
    curLenB = lB(j);
    %%

    % cut bGrand to wanted length:
    barcodes = cellfun(@(x) x.rawBarcode(x.rawBitmask),bGrand','un',false);
    barcodes = cellfun(@(x) x(randi(length(x)-curLen)+(1:curLen)),barcodes,'un',false); % random cut-out of length curLen
    bitmasks = cellfun(@(x) true(size(x)), barcodes, 'un', 0);
    barSynth = cell2struct([barcodes, bitmasks],{'rawBarcode','rawBitmask'},2);

        % cut bGrand to wanted length:
    barcodes2 = cellfun(@(x) x.rawBarcode(x.rawBitmask),bGrand2,'un',false);
    barcodes2 = cellfun(@(x) x(randi(length(x)-curLenB)+(1:curLenB)),barcodes2,'un',false); % random cut-out of length curLen
    bitmasks2 = cellfun(@(x) true(size(x)), barcodes2, 'un', 0);
    barSynth2 = cell2struct([barcodes2, bitmasks2],{'rawBarcode','rawBitmask'},2);


    for l=1:length(barSynth)
        A = arrayfun(@(y) imresize(barSynth(l).rawBarcode(barSynth(l).rawBitmask),'Scale' ,[1 y]),sF,'un',false);
        pcCur = 0;
        for d=1:length(A)
                pccThis = comparisonFun(A{d},barcodes2{1},ones(1,length(A{d})),ones(1,length(barcodes2{1})),2^16);
                pcCur = max([pcCur,pccThis]);
        end
        pccs{j,k}(l) = pcCur;
    end
    end
end

%% plot of PCC score distribution
TT = length(lB);
NN = length(lA);
maxPCCAll = pccs;

import Zeromodel.beta_ev_params;

intV = 5; % degrees of freedom for chi2 test
a_fit2 = zeros(TT,NN);n_fit2 = zeros(TT,NN);chi2Score = zeros(TT,NN);
for j=1:TT
    for k=1:NN
        [parameters] = beta_ev_params(maxPCCAll{j,k}, lA(k)/2);

        a_fit2(j,k) = parameters(1);
        n_fit2(j,k) = parameters(2);
    end
end

% PLOT FIGURE P-VAL fit

par1 = mean(a_fit2);
par1CI =  1.98*[std(a_fit2,[],1)./sqrt(NN)];
par3 = mean(n_fit2,2);
par3CI =  1.98*[std(n_fit2,[],2)./sqrt(TT)];

iy = 3; % La
j = 2;% Lb
%alphav=mean(par1./lA);
alphav=par1(iy)/lA(iy);
%alphav=a_fit2(j,iy)/lA(iy); %Exact value

%alphaN=mean(par3'./(2*lB));
alphaN=par3(j)/(2*lB(j)); 


f=figure('Position',[1 1 800 800])
tiledlayout(3,2,'TileSpacing','tight')

nexttile
hold on
errorbar(lA,par1./lA,par1CI./lA)
xlabel('$L_A$','Interpreter','latex')
ylabel('$\nu_{eff}/L_A$','Interpreter','latex')
title('(A) $\alpha_{\nu}$ estimation','Interpreter','latex')
plot(lA,alphav*ones(1,length(lA)),'--')
xlim([lA(1) lA(end)])

nexttile
hold on

errorbar(lB,par3'./(2*lB),par3CI'./(2*lB))
xlabel('$L_B$','Interpreter','latex')
ylabel('$N_{eff}/(2\cdot L_B)$','Interpreter','latex')
title('(B) $\alpha_{N}$ estimation','Interpreter','latex')
plot(lB,alphaN*ones(1,length(lB)),'--')
xlim([lB(1) lB(end)])

nexttile

imagesc(lA,lB,a_fit2);
title('(C) $\nu_{eff}$ heatmap','Interpreter','latex')
ylabel('$L_B$','Interpreter','latex')
xlabel('$L_A$','Interpreter','latex')
colorbar;colormap(gray) 
nexttile
imagesc(lB,lA,n_fit2');
colorbar;colormap(gray) 
title('(D) $N_{eff}$ heatmap','Interpreter','latex')
ylabel('$L_A$','Interpreter','latex')
xlabel('$L_B$','Interpreter','latex')

nexttile

xx=0.1:0.001:0.9;
import Zeromodel.beta_ev_pdf;
[p] = beta_ev_pdf(xx,a_fit2(j,iy), 1, n_fit2(j,iy));
[p2] = beta_ev_pdf(xx, alphav*lA(iy), 1, alphaN*2*lB(j));
% f=figure,
plot(xx,p)
hold on
histogram(pccs{j,iy},'Normalization','pdf')
title(['(E) EVD fit histogram,', ' $L_A=$ ',num2str(lA(iy)),' $L_B=$ ',num2str(lB(j))  ],'Interpreter','latex')
xlim([min(pccs{j,iy})-0.1 max(pccs{j,iy})+0.1])
plot(xx,p2,'LineWidth',2)

lgd = legend({'Fitted EVD','PCC histogram','Pre-calculated parameters EVD'},'Interpreter','latex','Location','southoutside')

%% local score distribution

bGrandN = bGrand;

maxPCC = cell(length(lA2),length(w));
osOutput = cell(length(lA2),length(w));
for j = 1:length(w)
    for k=1:length(lA2)
    j
    k
    curLen = lA2(k);
    minOverlap = w(j);
    
    % cut bGrand to wanted length:
    barcodes = cellfun(@(x) x.rawBarcode(x.rawBitmask),bGrandN','un',false);
    barcodes = cellfun(@(x) x(randi(length(x)-curLen)+(1:curLen)),barcodes,'un',false); % random cut-out of length curLen
    bitmasks = cellfun(@(x) true(size(x)), barcodes, 'un', 0);
    barSynth = cell2struct([barcodes, bitmasks],{'rawBarcode','rawBitmask'},2);
    tic
    % bars = barGenMerged(cellfun(@(x) sum(x.rawBitmask),barGenMerged)>minOverlap); % only for those larger than min overlap
    [oS] = calc_overlap_mp(barSynth,sF, minOverlap, [],3);
    toc
    localScore = [oS(:).score]; % local
    lenA = [oS(:).lenA]; % lenA
    lenB = [oS(:).lenB]; % lenB\
    partialScore = [oS(:).partialScore];
    partialLength = [oS(:).partialLength];
    overlapLength = [oS(:).overlaplen];

    idx = reshape(1:size(oS,1)*size(oS,2), size(oS,1),size(oS,2)); % from bg_test_1
    idx = tril(idx);

    localScore(idx(idx~=0)) = nan;
    partialScore(idx(idx~=0)) = nan;

    maxPCC{k,j} = localScore(~isnan(localScore));
    osOutput{k,j} = oS;
    end
end


%% PLOT FOR LOCAL DIST

import Zeromodel.beta_ev_params;

casetest=1;
a_fit2 = zeros(size(maxPCC,1),size(maxPCC,2)); % first dim lA2, second dim w
n_fit2 = zeros(size(maxPCC,1),size(maxPCC,2));
for j = 1:length(w)
    for k=1:length(lA2)
            [parameters] = beta_ev_params(maxPCC{k,j},w(j)/2);%(randLenMin+(k-1)*gap)/2);
            a_fit2(k,j) = parameters(1);
            n_fit2(k,j) = parameters(2);
    end
end

a_fit2 = a_fit2';
n_fit2 = n_fit2';

% Similar plot for local pcc
barL2 = lA2;
barL1 = lA2;

% figure,plot(mean(n_fit2')./(shortL-minOverlapLen))
lW = length(w);
par1 = zeros(1,lW);
par2 = zeros(1,lW);
par1CI = zeros(1,lW);
% choose 1
for j = 1:lW
    nPar = n_fit2(j,:);
    parEq = (2*(barL1));
    
    % fnpar = polyfit(parEq,nPar,1) % fit a line.    
    mdl = fitlm(parEq,nPar,'intercept',false);
    fnpar(1) = mdl.Coefficients.Estimate(1);
    fnpar(2) = 0;

    par2(j) = [mdl.Coefficients.Estimate(1)];
    par1(j) = [mean(a_fit2(j,:))]; % par1 mean
    par1CI(j) =  1.98*[std(a_fit2(j,:))/sqrt(length(a_fit2(j,:)))];
end


ix = 3; % this is w = 120?

par3 = n_fit2(ix,:); % can only do par3 for fixed w or fixed len2, since it depends on both. Write this in the text!
overlapL = w(ix);
% minL = randLenMin+gap*([1:TT]-1);

iy =5; %LA2 3
j=3; % w
xx=0.1:0.001:0.9;

%betanu = mean(par1./w);
%betanu=a_fit2(j,iy)/w(j);
betanu=par1(j)/w(j);

%betaN = mean(par3./((2.*((barL1-w(ix)).*((barL2-w(ix)))))));
%betaN=n_fit2(j,iy)/(2.*(lA2(iy)-w(j)).*(lA2(iy)-w(j)));
betaN=par3(iy)/(2.*(lA2(iy)-w(j)).*(lA2(iy)-w(j)));

f=figure('Position',[1 1 800 800])
tiledlayout(3,2,'TileSpacing','tight')

nexttile
hold on

errorbar(w,par1./w,par1CI./w)
xlabel('$w$','Interpreter','latex')
ylabel('$\nu_{eff}/w$','Interpreter','latex')
title('(A) $\beta_{\nu}$ estimation','Interpreter','latex')
plot(w,betanu*ones(1,length(w)),'--')

nexttile
hold on
% plot(barL1,par3./(2.*(max(barL1-overlapL(ix),barL2-overlapL(ix)))))
plot(barL1,par3./((2.*((barL1-w(ix)).*((barL2-w(ix)))))))
% plot(barL1,par3./((2.*((barL1).*((barL2))))))

xlabel('$L_A$','Interpreter','latex')
ylabel('$N_{eff}/(2(\cdot L_B-w)(L_A-w)$','Interpreter','latex') %(2(\cdot L_B-w)(L_A-w)
title('(B) $\beta_{N}$ estimation','Interpreter','latex')
% 
plot(barL1,betaN*ones(1,length(barL1)),'--')

nexttile
imagesc(barL1,w,a_fit2);
title('(C) $\nu_{eff}$ heatmap','Interpreter','latex')
colorbar;colormap(gray) 
xlabel('$L_A$','Interpreter','latex')
ylabel('$w$','Interpreter','latex')
nexttile
imagesc(barL1,w,n_fit2);
colorbar;colormap(gray) 
title('(D) $N_{eff}$ heatmap','Interpreter','latex')
xlabel('$L_A$','Interpreter','latex')
ylabel('$w$','Interpreter','latex')

nexttile
import Zeromodel.beta_ev_pdf;
[p] = beta_ev_pdf(xx,a_fit2(j,iy), 1, n_fit2(j,iy));
[p2] = beta_ev_pdf(xx,betanu*w(j), 1, betaN*(2.*(lA2(iy)-w(j)).*(lA2(iy)-w(j))));
plot(xx,p)
hold on
histogram(maxPCC{iy,j},'Normalization','pdf')
xlabel('Max local PCC','Interpreter','latex')
title(['E) EVD fit histogram,', ' $L_A=$ ',num2str(lA2(iy)),', $w=$ ',num2str(w(j)) ],'Interpreter','latex')
xlim([min(maxPCC{iy,j})-0.1 max(maxPCC{iy,j})+0.1])
plot(xx,p2,'LineWidth',2)
lgd = legend({'Fitted EVD','PCC histogram','Pre-calculated parameters EVD'},'Interpreter','latex','Location','southoutside');



save('E:/Spombe/reprocess after bargi/oS_results_pombe/parameters/5_mini_fixed2.mat','alphaN','alphav','betaN','betanu');
