function inv_patchdef = spm_cfg_eeg_inv_patchdef
% Configuration file for taking a number of previous inversion results
% (maybe based on different data), smoothing and creating an approximate posterior
%__________________________________________________________________________

% Gareth Barnes
% Copyright (C) 2010-2022 Wellcome Centre for Human Neuroimaging


inv_patchdef          = cfg_exbranch;
inv_patchdef.tag      = 'inv_patchdef';
inv_patchdef.name     = 'Set prior sampling';
inv_patchdef.val      = @inv_patchdef_cfg;
inv_patchdef.help     = {'Use to make up a patch definition file'};
inv_patchdef.prog     = @run_inv_patchdef;
inv_patchdef.vout     = @vout_inv_patchdef;
inv_patchdef.modality = {'MEG'};


%==========================================================================
function varargout = inv_patchdef_cfg

persistent cfg
if ~isempty(cfg), varargout = {cfg}; return; end

D = cfg_files;
D.tag = 'D';
D.name = 'Prior patch density file';
D.filter = 'gii';
D.num = [1 Inf];
D.help = {'Select the prior patch distribution file'};

npatch = cfg_entry;
npatch.tag = 'npatch';
npatch.name = 'Number of patches';
npatch.strtype = 'n';
npatch.help = {'Total number of patches'};
npatch.val = {512};

niter = cfg_entry;
niter.tag = 'niter';
niter.name = 'Number of iterations';
niter.strtype = 'n';
niter.help = {'Total number of random patch sets'};
niter.val = {32};

[cfg,varargout{1}] = deal({D,npatch,niter});


%==========================================================================
function  out = run_inv_patchdef(job)

inverse = [];
if numel(job.D)<1
    error('Need to add a number of files to combine');
end

prior_mesh=gifti(job.D);

pprob=prior_mesh.cdata(:);
[sprior,sortind]=sort(pprob);
cumprior=cumsum(sprior);

figure;
subplot(2,1,1);
plot(1:length(pprob),pprob,'o');

rng('shuffle');
Ip=zeros(job.niter,job.npatch);
dum=zeros(size(sprior));
for j=1:job.niter
    for k=1:job.npatch
         pval=(k-1)./job.npatch;
         sind=find(cumprior>=pval);
         randind=randperm(length(sind));
         Ip(j,k)=sortind(sind(randind(1)));
         dum(Ip(j,k))=dum(Ip(j,k))+1;
    end
    subplot(2,1,2); hold on;
    plot(1:length(pprob),dum,'.');
    drawnow;
    
end

[a1,b1,c1]=fileparts(job.D{1});
patchfilename=[a1 filesep b1 sprintf('N%d_Np%d_%d.mat',job.niter,job.npatch,randi(1e6))];
save(patchfilename,'Ip');

out.D = {patchfilename};


%==========================================================================
function dep = vout_inv_patchdef(job)
% Output is always in field "D", no matter how job is structured
dep = cfg_dep;
dep.sname = 'Smooth posterior from a number of imaging source reconstructions';
% reference field "D" from output
dep.src_output = substruct('.','D');
% this can be entered into any evaluated input
dep.tgt_spec   = cfg_findspec({{'filter','mat'}});
