% correlation to the best match within genome
nFields = 5; % each species, each hybrid allele, and the mean of both hybrid alleles
fieldNames = {'cer', 'par', 'hyc', 'hyp', 'hybMean'};
N = 5;
PCM2 = nan(length(intint), length(intint), nFields);
PCM = nan(length(intint), length(intint), nFields);
for i = 1:4; PCM(:, :, i) = PCMs.(F{i}).PCM; end; 
hybPCM = nan(length(intint), length(intint), 2);
hybPCM(:, :, 1) = PCMs.WCDhyc.PCM;
hybPCM(:, :, 2) = PCMs.WCDhyp.PCM;
hybPCMmean = nanmean(hybPCM, 3);
PCM(:, :, 5) = hybPCMmean;
nnWithinGenome = struct;
fieldLabels = {'firstOrder', 'secondOrder'};
disp(['pre-allocate nnWithinGenome']);
for i = 1:length(fieldLabels)
    nnWithinGenome.(fieldLabels{i}).corrs = nan(length(intint), nFields);
    nnWithinGenome.(fieldLabels{i}).geneIdx = nan(length(intint), nFields);
    nnWithinGenome.(fieldLabels{i}).corrsN = nan(length(intint), nFields);
    nnWithinGenome.(fieldLabels{i}).geneIdxN = cell(1, nFields);
end
%%
for i = 1:nFields
    disp(['generating first order corr for: ',fieldNames{i}]);
    currPCM = PCM(:, :, i);
    currPCM(1:(length(intint)+1):end) = nan;
    [currMax, maxIdx] = max(currPCM, [], 2);
    nnWithinGenome.firstOrder.corrs(:, i) = currMax;
    nnWithinGenome.firstOrder.geneIdx(:, i) = maxIdx;
    corrsN = nan(length(intint), 1);
    geneIdxN = nan(length(intint), N);
    for j = 1:length(intint)
        [sortd, sortdidx] = sort(currPCM(j, :), 'descend');
        corrsN(j) = nanmedian(sortd(2:N+1));
        geneIdxN(j, :) = sortdidx(2:N+1);   
    end
    nnWithinGenome.firstOrder.corrsN(:, i) = corrsN;
    nnWithinGenome.firstOrder.geneIdxN{i} = geneIdxN;
    
    disp(['generating second order corr for: ',fieldNames{i}]);
    currPCM2 = corr(currPCM, 'rows', 'pairwise');
    currPCM2(1:(length(intint)+1):end) = nan;
    [currMax2, maxIdx2] = max(currPCM2, [], 2);
    nnWithinGenome.secondOrder.corrs(:, i) = currMax2;
    nnWithinGenome.secondOrder.geneIdx(:, i) = maxIdx2;
    PCM2(:, :, i) = currPCM2;
    disp(['sortind and finding NNs for: ',fieldNames{i}]);
    for j = 1:length(intint)
        [sortd, sortdidx] = sort(currPCM2(j, :), 'descend');
        corrsN(j) = nanmedian(sortd(2:N+1));
        geneIdxN(j, :) = sortdidx(2:N+1);   
    end
    nnWithinGenome.secondOrder.corrsN(:, i) = corrsN;
    nnWithinGenome.secondOrder.geneIdxN{i} = geneIdxN;
end

mats = {PCM, PCM2};
for i = 1:2
    mat = mats{i};
    byCer = nnWithinGenome.(fieldLabels{i}).geneIdxN{1};
    parByCer = nan(length(intint), N);
    for j = 1:length(intint)
        currIdx = byCer(j, :);
        parByCer(j, :) = mat(j, currIdx, 2);
    end
    nnWithinGenome.(fieldLabels{i}).parByCerN = nanmedian(parByCer, 2);
    
    byPar = nnWithinGenome.(fieldLabels{i}).geneIdxN{2};
    cerByPar = nan(length(intint), N);
    for j = 1:length(intint)
        currIdx = byPar(j, :);
        cerByPar(j, :) = mat(j, currIdx, 1);
    end
    nnWithinGenome.(fieldLabels{i}).cerByParN = nanmedian(cerByPar, 2);
end
    
figure; 
subplot(1, 2, 1);
dscatter(nnWithinGenome.firstOrder.corrs(:, 1), nnWithinGenome.secondOrder.corrs(:, 1));
xlabel('1st order corr'); ylabel('2nd order corr');
subplot(1, 2, 2);
dscatter(nnWithinGenome.secondOrder.corrs(:, 1), nnWithinGenome.secondOrder.corrsN(:, 1));

name = datestr(datetime);
name = strrep(name(1:11), '-', '_');
save([homeDir, 'revision/PCM2_', name, '.mat'], 'PCM2');
save([homeDir, 'revision/nnWithinGenome_', name, '.mat'], 'nnWithinGenome');