text:
function [post] = spm_mci_post (mcmc,M,U,Y,true_P)
% Estimate posterior density
% FORMAT [post] = spm_mci_post (mcmc,M,U,Y,true_P)
%
% mcmc          .inference = 'amc','ais','vl' or 'langevin'
%               .verbose = 0 or 1 to plot progress (default 0)
%               .maxits = max number of iterations for sampling
%               .init = init parameter values (default is prior mean)
% M             model structure
% U             inputs (shouldn't be empty)
% Y             data
% true_P        true parameters (if known)
%
% post          structure containing posterior (mean, samples etc)
%__________________________________________________________________________
% Copyright (C) 2014 Wellcome Trust Centre for Neuroimaging

% Will Penny
% $Id: spm_mci_post.m 6697 2016-01-27 14:57:28Z spm $

try, verbose=mcmc.verbose; catch, verbose=0; end
if nargin < 5 || isempty(true_P)
    tp=0;
else
    tp=1;
end

tic;
switch mcmc.inference,

    case 'ais',
        disp('Annealed Importance Sampling');
        disp(' ');

        mcmc.nsamp=mcmc.maxits;

        tic;
        post = spm_mci_ais (mcmc,M,U,Y);
        toc

        Nsamp=size(post.P,2);
        post.ind=[1:Nsamp];
        post.Ep=mean(post.P(:,post.ind)')';
        post.mcmc=mcmc;

        % Eigenrep needed for spm_mci_joint later
        M = spm_mci_minit (M);

    case 'multi-amc',
        disp('Multiple chains of Adaptive Monte Carlo');
        disp(' ');
        Nsamp=ceil(0.5*mcmc.J);
        Nscale=0;
        Ntune=Nsamp;
        mc = spm_mci_popdef (Nscale,Ntune,Nsamp);
        mc.verbose=verbose;

        MM{1}=M;
        UU{1}=U;

        for it=1:mcmc.maxits,
            % Loop over chains
            mcmc.init{1}=spm_normrnd(M.pE,M.pC,1);
            Psamp = spm_mci_pop (mc,MM,UU,Y);
            if it ==1
                P=Psamp{1}.theta;
            else
                P=[P,Psamp{1}.theta];
            end
        end
        post.ind=[1:size(P,2)];
        post.Ep=mean(P(:,post.ind)')';
        post.P=P;
        post.mcmc=mcmc;

        % Eigenrep needed for spm_mci_joint later
        M = spm_mci_minit (M);

    case 'amc',
        disp('Adaptive Monte Carlo');
        disp(' ');

        % MH defaults
        tune=1;
        if tune
            Nscale=ceil(0.25*mcmc.maxits);
            Ntune=Nscale;
            Nsamp=ceil(0.5*mcmc.maxits);
        else
            Nscale=ceil(0.5*mcmc.maxits);
            Ntune=0;
            Nsamp=Nscale;
        end
        mc = spm_mci_popdef (Nscale,Ntune,Nsamp);
        mc.verbose=verbose;

        % Draw initialisation point from prior ?
        %mcmc.init{1}=spm_normrnd(M.pE,M.pC,1);
        try, mc.init{1}=mcmc.init; catch, mc.init{1}=spm_vec(M.pE); end

        MM{1}=M;
        UU{1}=U;
        tic;
        [Psamp,logev,D,MM] = spm_mci_pop (mc,MM,UU,Y);
        toc
        M=MM{1};

        if tp
            [post.Ep,post.SDp]=spm_mci_report (Psamp,mc,true_P);
        else
            disp('Initial params:');
            disp(mc.init{1});
            [post.Ep,post.SDp]=spm_mci_report (Psamp,mc);
        end

        post.Ep=post.Ep';
        post.P=Psamp{1}.theta;
        post.E=-Psamp{1}.logq;
        post.dL=Psamp{1}.dL;
        post.bayes_fb=zeros(1,length(post.E));
        post.acc=Psamp{1}.acc;
        post.logev=logev;
        post.D=D;
        post.mcmc=mcmc;
        post.ind=mc.ind_samp;

    case 'vl',
        disp('Variational Laplace');
        disp(' ');

        D.y=Y;
        if ~verbose
            M.nograph=1;
        end

        if isstruct(U)
            UI=U;
        else
            UI.u=U';
            UI.dt=M.T/M.N;
        end

        % Change VL defaults
        if isfield(mcmc,'maxits'), M.Nmax=mcmc.maxits; end
        if isfield(mcmc,'init'), M.P=mcmc.init; end
        if isfield(mcmc,'hE'), M.hE=mcmc.hE; end
        if isfield(mcmc,'hC'), M.hC=mcmc.hC; end

        [Ep,Cp,Eh,F,tmp1,tmp2,tmp3,k] = spm_nlsi_GN (M,UI,D);

        post.Ep=spm_vec(Ep);
        post.Cp=Cp;
        post.Eh=Eh;
        post.Ce=diag(1./exp(Eh));
        post.logev=F;
        post.its=k;
        if isfield(M,'P')
            post.init=M.P;
        end

        % Eigenrep needed for spm_mci_joint later
        M = spm_mci_minit (M);

    case 'langevin',
        disp('Langevin Monte Carlo');
        disp(' ');

        [M,stats]=spm_mci_lgv(mcmc,M,U,Y);
        Psamp=stats.P';
        Nsamp=size(Psamp,1);
        burn_in=round(0.3*size(Psamp,1));
        post=stats;
        post.ind=[burn_in+1:Nsamp];
        post.targ=Psamp(post.ind,:);
        post.Ep=mean(post.targ)';

        % 2.5%, 50% and 97.5% quantiles
        q = [.025 .5 .975];
        for j=1:M.Np,
            sx = post.P(j,post.ind);
            post.quantiles(j,:) = quantile(sx,q);
        end

    otherwise
        disp('Unknown inference method');
end

if strcmp(mcmc.inference,'VL')
    Up=UI;
else
    Up=U;
end

% Generate data fit from posterior mean
if isfield(M,'IS')
    if strcmp(M.IS,'spm_gen_erp')
        Ps=spm_unvec(post.Ep,M.pE);
        post.Yhat=feval('spm_gen_erp',Ps,M,Up);
    else
        post.Yhat = feval(M.IS,post.Ep,M,Up);
    end
else
    post.Yhat = spm_mci_fwd (post.Ep,M,Up);
end

post.els=toc;
disp(sprintf('Optimisation time = %1.2f seconds',post.els));

lw=2;
if ~isfield(M,'IS')
    % Plot time series
    figure
    rm=ceil(sqrt(M.l));
    for i=1:M.l,
        if M.l>3
            subplot(rm,rm,i);
        else
            subplot(M.l,1,i);
        end
        plot(M.t,Y(:,i),'LineWidth',lw);
        hold on
        plot(M.t,post.Yhat(:,i),'r','LineWidth',lw);
        grid on
        set(gca,'FontSize',16);
        legend('Data','Fit');
        xlabel('Time');
        ylabel(sprintf('y(%d)',i));
    end
end

% get parameters in reduced space
Pr=M.V'*(post.Ep-M.vpE);
post.L_est = spm_mci_joint (Pr,M,U,Y);
disp(sprintf('Estimated Log Joint=%1.2f',post.L_est));

if tp
    % get parameters in reduced space
    Pr=M.V'*(spm_vec(true_P)-M.vpE);
    post.L_true = spm_mci_joint (Pr,M,U,Y);
    disp(sprintf('True Log Joint=%1.2f',post.L_true));
    disp(' ');
end

pt=4;
if M.Np > pt
    hp=figure;
    set(hp,'Name','Parameters');
    plot(post.Ep,'r','LineWidth',lw);
    xlabel('Parameter');
    set(gca,'FontSize',16);
    grid on
else
    disp('Estimated (latent) params:');
    disp(post.Ep);
end

if tp
    if M.Np > pt
        hold on
        plot(spm_vec(true_P),'LineWidth',lw);
        legend('Estimated','True');
    else
        disp('True (latent) params:');
        disp(spm_vec(true_P));
    end
end

switch mcmc.inference,
    case {'amc','langevin'},
        post.type='sample';
    otherwise
        post.type='gaussian';
end

author: spm


text:
function [cls,M1] = spm_preproc_write8(res,tc,bf,df,mrf,cleanup,bb,vx,odir)
% Write out VBM preprocessed data
% FORMAT [cls,M1] = spm_preproc_write8(res,tc,bf,df,mrf,cleanup,bb,vx,odir)
%__________________________________________________________________________
% Copyright (C) 2008-2016 Wellcome Trust Centre for Neuroimaging

% John Ashburner
% $Id: spm_preproc_write8.m 7415 2018-09-10 18:24:16Z john $

% Prior adjustment factor.
% This is a fudge factor to weaken the effects of the tissue priors.  The
% idea here is that the bias from the tissue priors probably needs to be
% reduced because of the spatial smoothing typically used in VBM studies.
% Having the optimal bias/variance tradeoff for each voxel is not the same
% as having the optimal tradeoff for weighted averages over several voxels.

if isfield(res,'mg')
    lkp = res.lkp;
    Kb  = max(lkp);
else
    Kb  = size(res.intensity(1).lik,2);
end
N       = numel(res.image);

if nargin<2, tc = true(Kb,4); end % native, import, warped, warped-mod
if nargin<3, bf = false(N,2); end % field, corrected
if nargin<4, df = false(1,2); end % inverse, forward
if nargin<5, mrf = 1;         end % MRF parameter
if nargin<6, cleanup = 1;     end % Run the ad hoc cleanup
if nargin<7, bb = NaN(2,3);   end % Default to TPM bounding box
if nargin<8, vx = NaN;        end % Default to TPM voxel size
if nargin<9, odir = [];       end % Output directory

% Read essentials from tpm (it will be cleared later)
tpm = res.tpm;
if ~isstruct(tpm) || ~isfield(tpm, 'bg1')
    tpm = spm_load_priors8(tpm);
end
d1      = size(tpm.dat{1});
d1      = d1(1:3);
M1      = tpm.M;

% Define orientation and field of view of any "normalised" space
% data that may be generated (wc*.nii, mwc*.nii, rc*.nii & y_*.nii).
if any(isfinite(bb(:))) || any(isfinite(vx))
    % If a bounding box is supplied, combine this with the closest
    % bounding box derived from the dimensions and orientations of
    % the tissue priors.
    [bb1,vx1] = spm_get_bbox(tpm.V(1), 'old');
    bb(~isfinite(bb)) = bb1(~isfinite(bb));
    if ~isfinite(vx), vx = abs(prod(vx1))^(1/3); end
    bb(1,:) = vx*round(bb(1,:)/vx);
    bb(2,:) = vx*round(bb(2,:)/vx);
    odim    = abs(round((bb(2,1:3)-bb(1,1:3))/vx))+1;

    mm  = [[bb(1,1) bb(1,2) bb(1,3)
            bb(2,1) bb(1,2) bb(1,3)
            bb(1,1) bb(2,2) bb(1,3)
            bb(2,1) bb(2,2) bb(1,3)
            bb(1,1) bb(1,2) bb(2,3)
            bb(2,1) bb(1,2) bb(2,3)
            bb(1,1) bb(2,2) bb(2,3)
            bb(2,1) bb(2,2) bb(2,3)]'; ones(1,8)];
    vx3 = [[1       1       1
            odim(1) 1       1
            1       odim(2) 1
            odim(1) odim(2) 1
            1       1       odim(3)
            odim(1) 1       odim(3)
            1       odim(2) odim(3)
            odim(1) odim(2) odim(3)]'; ones(1,8)];
    mat    = mm/vx3;
else
    % Use the actual dimensions and orientations of
    % the tissue priors.
    odim = tpm.V(1).dim;
    mat  = tpm.V(1).mat;
end

[pth,nam] = fileparts(res.image(1).fname);
if ~isempty(odir) && ischar(odir), pth = odir; end
ind  = res.image(1).n;
d    = res.image(1).dim(1:3);

[x1,x2,o] = ndgrid(1:d(1),1:d(2),1);
x3 = 1:d(3);

chan(N) = struct('B1',[],'B2',[],'B3',[],'T',[],'Nc',[],'Nf',[],'ind',[]);
for n=1:N
    d3         = [size(res.Tbias{n}) 1];
    chan(n).B3 = spm_dctmtx(d(3),d3(3),x3);
    chan(n).B2 = spm_dctmtx(d(2),d3(2),x2(1,:)');
    chan(n).B1 = spm_dctmtx(d(1),d3(1),x1(:,1));
    chan(n).T  = res.Tbias{n};

    % Need to fix writing of bias fields or bias corrected images, when the data used are 4D.
    [pth1,nam1] = fileparts(res.image(n).fname);
    if ~isempty(odir) && ischar(odir), pth1 = odir; end
    chan(n).ind = res.image(n).n;

    if bf(n,2)
        chan(n).Nc      = nifti;
        chan(n).Nc.dat  = file_array(fullfile(pth1,['m', nam1, '.nii']),...
                                     res.image(n).dim(1:3),...
                                     [spm_type('float32') spm_platform('bigend')],...
                                     0,1,0);
        chan(n).Nc.mat  = res.image(n).mat;
        chan(n).Nc.mat0 = res.image(n).mat;
        chan(n).Nc.descrip = 'Bias corrected';
        create(chan(n).Nc);
    end

    if bf(n,1),
        chan(n).Nf      = nifti;
        chan(n).Nf.dat  = file_array(fullfile(pth1,['BiasField_', nam1, '.nii']),...
                                     res.image(n).dim(1:3),...
                                     [spm_type('float32') spm_platform('bigend')],...
                                     0,1,0);
        chan(n).Nf.mat  = res.image(n).mat;
        chan(n).Nf.mat0 = res.image(n).mat;
        chan(n).Nf.descrip = 'Estimated Bias Field';
        create(chan(n).Nf);
    end
end

do_cls   = any(tc(:)) || nargout>=1;
tiss(Kb) = struct('Nt',[]);
for k1=1:Kb
    if tc(k1,4) || any(tc(:,3)) || tc(k1,2) || nargout>=1,
        do_cls  = true;
    end
    if tc(k1,1),
        tiss(k1).Nt      = nifti;
        tiss(k1).Nt.dat  = file_array(fullfile(pth,['c', num2str(k1), nam, '.nii']),...
                                      res.image(n).dim(1:3),...
                                      [spm_type('uint8') spm_platform('bigend')],...
                                      0,1/255,0);
        tiss(k1).Nt.mat  = res.image(n).mat;
        tiss(k1).Nt.mat0 = res.image(n).mat;
        tiss(k1).Nt.descrip = ['Tissue class ' num2str(k1)];
        create(tiss(k1).Nt);
        do_cls = true;
    end
end

prm     = [3 3 3 0 0 0];
Coef    = cell(1,3);
Coef{1} = spm_bsplinc(res.Twarp(:,:,:,1),prm);
Coef{2} = spm_bsplinc(res.Twarp(:,:,:,2),prm);
Coef{3} = spm_bsplinc(res.Twarp(:,:,:,3),prm);

do_defs = any(df);
do_defs = do_cls | do_defs;
if do_defs
    if df(1)
        Ndef      = nifti;
        Ndef.dat  = file_array(fullfile(pth,['iy_', nam, '.nii']),...
                               [res.image(1).dim(1:3),1,3],...
                               [spm_type('float32') spm_platform('bigend')],...
                               0,1,0);
        Ndef.mat  = res.image(1).mat;
        Ndef.mat0 = res.image(1).mat;
        Ndef.descrip = 'Inverse Deformation';
        create(Ndef);
    end
    if df(2) || any(any(tc(:,[2,3,4]))) || nargout>=1,
        y = zeros([res.image(1).dim(1:3),3],'single');
    end
end

spm_progress_bar('init',length(x3),['Working on ' nam],'Planes completed');
M = M1\res.Affine*res.image(1).mat;

if do_cls
    Q = zeros([d(1:3),Kb],'single');
end

for z=1:length(x3)

    % Bias corrected image
    cr = cell(1,N);
    for n=1:N
        f          = spm_sample_vol(res.image(n),x1,x2,o*x3(z),0);
        bf         = exp(transf(chan(n).B1,chan(n).B2,chan(n).B3(z,:),chan(n).T));
        cr{n}      = bf.*f;
        if ~isempty(chan(n).Nc),
            % Write a plane of bias corrected data
            chan(n).Nc.dat(:,:,z,chan(n).ind(1),chan(n).ind(2)) = cr{n};
        end
        if ~isempty(chan(n).Nf),
            % Write a plane of bias field
            chan(n).Nf.dat(:,:,z,chan(n).ind(1),chan(n).ind(2)) = bf;
        end
    end

    if do_defs
        % Compute the deformation (mapping voxels in image to voxels in TPM)
        [t1,t2,t3] = defs(Coef,z,res.MT,prm,x1,x2,x3,M);

        if exist('Ndef','var')
            % Write out the deformation to file, adjusting it so mapping is
            % to voxels (voxels in image to mm in TPM)
            Ndef.dat(:,:,z,1,1) = M1(1,1)*t1 + M1(1,2)*t2 + M1(1,3)*t3 + M1(1,4);
            Ndef.dat(:,:,z,1,2) = M1(2,1)*t1 + M1(2,2)*t2 + M1(2,3)*t3 + M1(2,4);
            Ndef.dat(:,:,z,1,3) = M1(3,1)*t1 + M1(3,2)*t2 + M1(3,3)*t3 + M1(3,4);
        end

        if exist('y','var')
            % If needed later, save in variable y
            y(:,:,z,1) = t1;
            y(:,:,z,2) = t2;
            y(:,:,z,3) = t3;
        end

        if do_cls
            % Generate variable Q if tissue classes are needed
            msk = any((f==0) | ~isfinite(f),3);

            if isfield(res,'mg')
                % Parametric representation of intensity distributions
                q   = zeros([d(1:2) Kb]);
                q1  = likelihoods(cr,[],res.mg,res.mn,res.vr);
                q1  = reshape(q1,[d(1:2),numel(res.mg)]);
                b   = spm_sample_priors8(tpm,t1,t2,t3);
                wp  = res.wp;
                s   = zeros(size(b{1}));
                for k1 = 1:Kb
                    b{k1} = wp(k1)*b{k1};
                    s     = s + b{k1};
                end
                for k1=1:Kb
                    tmp       = sum(q1(:,:,lkp==k1),3);
                    tmp(msk)  = 1e-3;
                    q(:,:,k1) = tmp.*(b{k1}./s);
                end
            else
                % Nonparametric representation of intensity distributions
                q   = spm_sample_priors8(tpm,t1,t2,t3);
                wp  = res.wp;
                s   = zeros(size(q{1}));
                for k1 = 1:Kb
                    q{k1} = wp(k1)*q{k1};
                    s     = s + q{k1};
                end
                for k1 = 1:Kb
                    q{k1} = q{k1}./s;
                end
                q   = cat(3,q{:});

                for n=1:N
                    tmp = round(cr{n}*res.intensity(n).interscal(2) + res.intensity(n).interscal(1));
                    tmp = min(max(tmp,1),size(res.intensity(n).lik,1));
                    for k1=1:Kb
                        likelihood = res.intensity(n).lik(:,k1);
                        q(:,:,k1)  = q(:,:,k1).*likelihood(tmp);
                    end
                end
            end
            Q(:,:,z,:) = reshape(q,[d(1:2),1,Kb]);
        end
    end
    spm_progress_bar('set',z);
end
spm_progress_bar('clear');

cls   = cell(1,Kb);
if do_cls
    P = zeros([d(1:3),Kb],'uint8');
    if mrf==0
        % Normalise to sum to 1
        sQ = (sum(Q,4)+eps)/255;
        for k1=1:size(Q,4)
            P(:,:,:,k1) = uint8(round(Q(:,:,:,k1)./sQ));
        end
        clear sQ
    else
        % Use a MRF cleanup procedure
        nmrf_its = 10;
        spm_progress_bar('init',nmrf_its,['MRF: Working on ' nam],'Iterations completed');
        G   = ones([Kb,1],'single')*mrf;
        vx2 = 1./single(sum(res.image(1).mat(1:3,1:3).^2));
        for iter=1:nmrf_its
            spm_mrf(P,Q,G,vx2);
            spm_progress_bar('set',iter);
        end
    end
    clear Q

    if cleanup
        % Use an ad hoc brain cleanup procedure
        if size(P,4)>3
            P = clean_gwc(P,cleanup);
        else
            warning('Cleanup not done.');
        end
    end

    % Write tissues if necessary
    for k1=1:Kb
        if ~isempty(tiss(k1).Nt)
            for z=1:length(x3)
                tmp = double(P(:,:,z,k1))/255;
                tiss(k1).Nt.dat(:,:,z,ind(1),ind(2)) = tmp;
            end
        end
    end
    spm_progress_bar('clear');

    % Put tissue classes into a cell array...
    for k1=1:Kb
        if tc(k1,4) || any(tc(:,3)) || tc(k1,2) || nargout>=1
            cls{k1} = P(:,:,:,k1);
        end
    end
    clear P % ...and remove the original 4D array
end

clear tpm
M0  = res.image(1).mat;

if any(tc(:,2))
    % "Imported" tissue class images

    % Generate mm coordinates of where deformations map from
    x      = affind(rgrid(d),M0);

    % Generate mm coordinates of where deformation maps to
    y1     = affind(y,M1);

    % Procrustes analysis to compute the closest rigid-body
    % transformation to the deformation, weighted by the
    % interesting tissue classes.
    ind    = find(tc(:,2)); % Saved tissue classes
    [dummy,R] = spm_get_closest_affine(x,y1,single(cls{ind(1)})/255);
    clear x y1

    mat0   = R\mat; % Voxel-to-world of original image space

    fwhm   = max(vx./sqrt(sum(res.image(1).mat(1:3,1:3).^2))-1,0.01);
    for k1=1:size(tc,1)
        if tc(k1,2)

            % Low pass filtering to reduce aliasing effects in downsampled images,
            % then reslice and write to disk
            tmp1    = decimate(single(cls{k1}),fwhm);
            Ni      = nifti;
            Ni.dat  = file_array(fullfile(pth,['rc', num2str(k1), nam, '.nii']),...
                                 odim,...
                                 [spm_type('float32') spm_platform('bigend')],...
                                 0,1,0);
            Ni.mat         = mat;
            Ni.mat_intent  = 'Aligned';
            Ni.mat0        = mat0;
            Ni.mat0_intent = 'Aligned';
            Ni.descrip     = ['Imported Tissue ' num2str(k1)];
            create(Ni);

            for i=1:odim(3)
                tmp           = spm_slice_vol(tmp1,M0\mat0*spm_matrix([0 0 i]),odim(1:2),[1,NaN])/255;
                Ni.dat(:,:,i) = tmp;
            end
            clear tmp1
        end
    end
end

if any(tc(:,3)) || any(tc(:,4)) || nargout>=1 || df(2)
    % Adjust stuff so that warped data (and deformations) have the
    % desired bounding box and voxel sizes, instead of being the same
    % as those of the tissue probability maps.
    M = mat\M1;
    for i=1:size(y,3)
        t1         = y(:,:,i,1);
        t2         = y(:,:,i,2);
        t3         = y(:,:,i,3);
        y(:,:,i,1) = M(1,1)*t1 + M(1,2)*t2 + M(1,3)*t3 + M(1,4);
        y(:,:,i,2) = M(2,1)*t1 + M(2,2)*t2 + M(2,3)*t3 + M(2,4);
        y(:,:,i,3) = M(3,1)*t1 + M(3,2)*t2 + M(3,3)*t3 + M(3,4);
    end
    M1 = mat;
    d1 = odim;
end

if any(tc(:,3)) || any(tc(:,4)) || nargout>=1

    if any(tc(:,3))
        C = zeros([d1,Kb],'single');
    end

    spm_progress_bar('init',Kb,'Warped Tissue Classes','Classes completed');
    for k1 = 1:Kb
        if ~isempty(cls{k1})
            c = single(cls{k1})/255;
            if any(tc(:,3))
                [c,w]       = spm_diffeo('push',c,y,d1(1:3));
                vx          = sqrt(sum(M1(1:3,1:3).^2));
                spm_field('boundary',1);
                C(:,:,:,k1) = spm_field(w,c,[vx  1e-6 1e-4 0  3 2]);
                clear w
            else
                c      = spm_diffeo('push',c,y,d1(1:3));
            end
            if nargout>=1
                cls{k1} = c;
            end
            if tc(k1,4)
                N      = nifti;
                N.dat  = file_array(fullfile(pth,['mwc', num2str(k1), nam, '.nii']),...
                                    d1,...
                                    [spm_type('float32') spm_platform('bigend')],...
                                    0,1,0);
                N.mat  = M1;
                N.mat0 = M1;
                N.descrip = ['Jac. sc. warped tissue class ' num2str(k1)];
                create(N);
                N.dat(:,:,:) = c*abs(det(M0(1:3,1:3))/det(M1(1:3,1:3)));
            end
            spm_progress_bar('set',k1);
        end
    end
    spm_progress_bar('Clear');

    if any(tc(:,3))
        spm_progress_bar('init',Kb,'Writing Warped Tis Cls','Classes completed');
        C = max(C,eps);
        s = sum(C,4);
        for k1=1:Kb
            if tc(k1,3)
                N      = nifti;
                N.dat  = file_array(fullfile(pth,['wc', num2str(k1), nam, '.nii']),...
                                    d1,'uint8',0,1/255,0);
                N.mat  = M1;
                N.mat0 = M1;
                N.descrip = ['Warped tissue class ' num2str(k1)];
                create(N);
                N.dat(:,:,:) = C(:,:,:,k1)./s;
            end
            spm_progress_bar('set',k1);
        end
        spm_progress_bar('Clear');
        clear C s
    end
end

if df(2)
    y         = spm_diffeo('invdef',y,d1,eye(4),M0);
    y         = spm_extrapolate_def(y,M1);
    N         = nifti;
    N.dat     = file_array(fullfile(pth,['y_', nam, '.nii']),...
                           [d1,1,3],'float32',0,1,0);
    N.mat     = M1;
    N.mat0    = M1;
    N.descrip = 'Deformation';
    create(N);
    N.dat(:,:,:,:,:) = reshape(y,[d1,1,3]);
end

return;

%==========================================================================
% function [x1,y1,z1] = defs(sol,z,MT,prm,x0,y0,z0,M)
%==========================================================================
function [x1,y1,z1] = defs(sol,z,MT,prm,x0,y0,z0,M)
iMT = inv(MT);
x1  = x0*iMT(1,1)+iMT(1,4);
y1  = y0*iMT(2,2)+iMT(2,4);
z1  = (z0(z)*iMT(3,3)+iMT(3,4))*ones(size(x1));
x1a = x0    + spm_bsplins(sol{1},x1,y1,z1,prm);
y1a = y0    + spm_bsplins(sol{2},x1,y1,z1,prm);
z1a = z0(z) + spm_bsplins(sol{3},x1,y1,z1,prm);
x1  = M(1,1)*x1a + M(1,2)*y1a + M(1,3)*z1a + M(1,4);
y1  = M(2,1)*x1a + M(2,2)*y1a + M(2,3)*z1a + M(2,4);
z1  = M(3,1)*x1a + M(3,2)*y1a + M(3,3)*z1a + M(3,4);
return;

%==========================================================================
% function t = transf(B1,B2,B3,T)
%==========================================================================
function t = transf(B1,B2,B3,T)
if ~isempty(T)
    d2 = [size(T) 1];
    t1 = reshape(reshape(T, d2(1)*d2(2),d2(3))*B3', d2(1), d2(2));
    t  = B1*t1*B2';
else
    t = zeros(size(B1,1),size(B2,1),size(B3,1));
end
return;

%==========================================================================
% function p = likelihoods(f,bf,mg,mn,vr)
%==========================================================================
function p = likelihoods(f,bf,mg,mn,vr)
K  = numel(mg);
N  = numel(f);
M  = numel(f{1});
cr = zeros(M,N);
for n=1:N
    if isempty(bf)
        cr(:,n) = double(f{n}(:));
    else
        cr(:,n) = double(f{n}(:).*bf{n}(:));
    end
end
p  = ones(numel(f{1}),K);
for k=1:K
    amp    = mg(k)/sqrt((2*pi)^N * det(vr(:,:,k)));
    d      = bsxfun(@minus,cr,mn(:,k)')/chol(vr(:,:,k));
    p(:,k) = amp*exp(-0.5*sum(d.*d,2)) + eps;
end
return;

%==========================================================================
% function dat = decimate(dat,fwhm)
%==========================================================================
function dat = decimate(dat,fwhm)
% Convolve the volume in memory (fwhm in voxels).
lim = ceil(2*fwhm);
x  = -lim(1):lim(1); x = spm_smoothkern(fwhm(1),x); x  = x/sum(x);
y  = -lim(2):lim(2); y = spm_smoothkern(fwhm(2),y); y  = y/sum(y);
z  = -lim(3):lim(3); z = spm_smoothkern(fwhm(3),z); z  = z/sum(z);
i  = (length(x) - 1)/2;
j  = (length(y) - 1)/2;
k  = (length(z) - 1)/2;
spm_conv_vol(dat,dat,x,y,z,-[i j k]);
return;

%==========================================================================
% function y1 = affind(y0,M)
%==========================================================================
function y1 = affind(y0,M)
y1 = zeros(size(y0),'single');
for d=1:3
    y1(:,:,:,d) = y0(:,:,:,1)*M(d,1) + y0(:,:,:,2)*M(d,2) + y0(:,:,:,3)*M(d,3) + M(d,4);
end
return;

%==========================================================================
% function x = rgrid(d)
%==========================================================================
function x = rgrid(d)
x = zeros([d(1:3) 3],'single');
[x1,x2] = ndgrid(single(1:d(1)),single(1:d(2)));
for i=1:d(3)
    x(:,:,i,1) = x1;
    x(:,:,i,2) = x2;
    x(:,:,i,3) = single(i);
end
return;

%==========================================================================
% function [P] = clean_gwc(P,level)
%==========================================================================
function [P] = clean_gwc(P,level)
if nargin<2, level = 1; end

b    = P(:,:,:,2);

% Build a 3x3x3 seperable smoothing kernel
%--------------------------------------------------------------------------
kx=[0.75 1 0.75];
ky=[0.75 1 0.75];
kz=[0.75 1 0.75];
sm=sum(kron(kron(kz,ky),kx))^(1/3);
kx=kx/sm; ky=ky/sm; kz=kz/sm;

th1 = 0.15;
if level==2, th1 = 0.2; end
% Erosions and conditional dilations
%--------------------------------------------------------------------------
niter  = 32;
niter2 = 32;
spm_progress_bar('Init',niter+niter2,'Extracting Brain','Iterations completed');
for j=1:niter
    if j>2, th=th1; else th=0.6; end  % Dilate after two its of erosion
    for i=1:size(b,3)
        gp       = double(P(:,:,i,1));
        wp       = double(P(:,:,i,2));
        bp       = double(b(:,:,i))/255;
        bp       = (bp>th).*(wp+gp);
        b(:,:,i) = uint8(round(bp));
    end
    spm_conv_vol(b,b,kx,ky,kz,-[1 1 1]);
    spm_progress_bar('Set',j);
end

% Also clean up the CSF.
if niter2 > 0
    c = b;
    for j=1:niter2
        for i=1:size(b,3)
            gp       = double(P(:,:,i,1));
            wp       = double(P(:,:,i,2));
            cp       = double(P(:,:,i,3));
            bp       = double(c(:,:,i))/255;
            bp       = (bp>th).*(wp+gp+cp);
            c(:,:,i) = uint8(round(bp));
        end
        spm_conv_vol(c,c,kx,ky,kz,-[1 1 1]);
        spm_progress_bar('Set',j+niter);
    end
end

th = 0.05;
for i=1:size(b,3)
    slices = cell(1,size(P,4));
    for k1=1:size(P,4)
        slices{k1} = double(P(:,:,i,k1))/255;
    end
    bp        = double(b(:,:,i))/255;
    bp        = ((bp>th).*(slices{1}+slices{2}))>th;
    slices{1} = slices{1}.*bp;
    slices{2} = slices{2}.*bp;

    if niter2>0
        cp        = double(c(:,:,i))/255;
        cp        = ((cp>th).*(slices{1}+slices{2}+slices{3}))>th;
        slices{3} = slices{3}.*cp;
    end
    tot       = zeros(size(bp))+eps;
    for k1=1:size(P,4)
        tot   = tot + slices{k1};
    end
    for k1=1:size(P,4)
        P(:,:,i,k1) = uint8(round(slices{k1}./tot*255));
    end
end
spm_progress_bar('Clear');


text:
classdef FDV < ALGORITHM
% <multi/many> <real/integer> <large/none>
% Fuzzy decision variable framework with various internal optimizers
% Rate      --- 0.8 --- Fuzzy evolution rate. Default = 0.8
% Acc       --- 0.4 --- Step acceleration. Default = 0.4
% optimizer ---   5 --- Internal optimisation algorithm. 1 = NSGA-II, 2 = NSGA-III, 3 = MOEA/D, 4 = CMOPSO, 5 = LMOCSO
% type      ---   1 --- The type of aggregation function for MOEA/D

%------------------------------- Reference --------------------------------
% X. Yang, J. Zou, S. Yang, J. Zheng, and Y. Liu, A fuzzy decision
% variables framework for large-scale multiobjective optimization, IEEE
% Transactions on Evolutionary Computation, 2021.
%--------------------------------------------------------------------------

%  Copyright (C) 2021 Xu Yang
%  Xu Yang <xuyang.busyxu@qq.com> or <xuyang369369@gmail.com>

    methods
        function main(Algorithm,Problem)
            %% Set the default parameters
            [Rate,Acc,optimizer,type] = Algorithm.ParameterSet(0.8,0.4,5,1);

            %% NSGAII
            if optimizer==1
                % Generate random population
                Population = Problem.Initialization();
                [~,FrontNo,CrowdDis] = EnvironmentalSelection_NSGAII(Population,Problem.N);

                % Optimization
                while Algorithm.NotTerminated(Population)
                    MatingPool = TournamentSelection(2,Problem.N,FrontNo,-CrowdDis);
                    OffDec     = OperatorGA(Problem,Population(MatingPool).decs);
                    %% FDV
                    if Problem.FE/Problem.maxFE <= Rate
                        Offspring = FDVOperator(Problem,Rate,Acc,OffDec);
                    else
                        Offspring = Problem.Evaluation(OffDec);
                    end
                    %%
                    [Population,FrontNo,CrowdDis] = EnvironmentalSelection_NSGAII([Population,Offspring],Problem.N);
                end
            end

            %% NSGAIII
            if optimizer==2
                % Generate the reference points and random population
                [Z,Problem.N] = UniformPoint(Problem.N,Problem.M);
                Population    = Problem.Initialization();
                Zmin          = min(Population(all(Population.cons<=0,2)).objs,[],1);

                % Optimization
                while Algorithm.NotTerminated(Population)
                    MatingPool = TournamentSelection(2,Problem.N,sum(max(0,Population.cons),2));
                    OffDec     = OperatorGA(Problem,Population(MatingPool).decs);
                    %% FDV
                    if Problem.FE/Problem.maxFE <= Rate
                        Offspring = FDVOperator(Problem,Rate,Acc,OffDec);
                    else
                        Offspring = Problem.Evaluation(OffDec);
                    end
                    Zmin       = min([Zmin;Offspring(all(Offspring.cons<=0,2)).objs],[],1);
                    Population = EnvironmentalSelection_NSGAIII([Population,Offspring],Problem.N,Z,Zmin);
                end
            end

            %% MOEA/D
            if optimizer==3
                % Generate the weight vectors
                [W,Problem.N] = UniformPoint(Problem.N,Problem.M);
                T = ceil(Problem.N/10);

                % Detect the neighbours of each solution
                B = pdist2(W,W);
                [~,B] = sort(B,2);
                B = B(:,1:T);

                % Generate random population
                Population = Problem.Initialization();
                Z = min(Population.objs,[],1);

                % Optimization
                while Algorithm.NotTerminated(Population)
                    % For each solution
                    for i = 1 : Problem.N
                        % Choose the parents
                        P = B(i,randperm(size(B,2)));

                        % Generate an offspring
                        OffDec = OperatorGAhalf(Problem,Population(P(1:2)).decs);
                        %% FDV
                        if Problem.FE/Problem.maxFE <= Rate
                            Offspring = FDVOperator(Problem,Rate,Acc,OffDec);
                        else
                            Offspring = Problem.Evaluation(OffDec);
                        end

                        % Update the ideal point
                        Z = min(Z,Offspring.obj);

                        % Update the neighbours
                        switch type
                            case 1
                                % PBI approach
                                normW   = sqrt(sum(W(P,:).^2,2));
                                normP   = sqrt(sum((Population(P).objs-repmat(Z,T,1)).^2,2));
                                normO   = sqrt(sum((Offspring.obj-Z).^2,2));
                                CosineP = sum((Population(P).objs-repmat(Z,T,1)).*W(P,:),2)./normW./normP;
                                CosineO = sum(repmat(Offspring.obj-Z,T,1).*W(P,:),2)./normW./normO;
                                g_old   = normP.*CosineP + 5*normP.*sqrt(1-CosineP.^2);
                                g_new   = normO.*CosineO + 5*normO.*sqrt(1-CosineO.^2);
                            case 2
                                % Tchebycheff approach
                                g_old = max(abs(Population(P).objs-repmat(Z,T,1)).*W(P,:),[],2);
                                g_new = max(repmat(abs(Offspring.obj-Z),T,1).*W(P,:),[],2);
                            case 3
                                % Tchebycheff approach with normalization
                                Zmax  = max(Population.objs,[],1);
                                g_old = max(abs(Population(P).objs-repmat(Z,T,1))./repmat(Zmax-Z,T,1).*W(P,:),[],2);
                                g_new = max(repmat(abs(Offspring.obj-Z)./(Zmax-Z),T,1).*W(P,:),[],2);
                            case 4
                                % Modified Tchebycheff approach
                                g_old = max(abs(Population(P).objs-repmat(Z,T,1))./W(P,:),[],2);
                                g_new = max(repmat(abs(Offspring.obj-Z),T,1)./W(P,:),[],2);
                        end
                        Population(P(g_old>=g_new)) = Offspring;
                    end
                end
            end

            %% CMOPSO
            if optimizer == 4
                % Generate random population
                Population = Problem.Initialization();

                % Optimization
                while Algorithm.NotTerminated(Population)
                    [OffDec,OffVel] = Operator_CMOPSO(Problem,Population);
                    %% FDV
                    if Problem.FE/Problem.maxFE <= Rate
                        Offspring = FDVOperator(Rate,Acc,OffDec,OffVel);
                    else
                        Offspring = Problem.Evaluation(OffDec,OffVel);
                    end
                    Population = EnvironmentalSelection_CMOPSO([Population,Offspring],Problem.N);
                end
            end

            %% LMOCSO
            if optimizer == 5
                 % Generate random population
                [V,Problem.N] = UniformPoint(Problem.N,Problem.M);
                Population    = Problem.Initialization();
                Population    = EnvironmentalSelection_LMOCSO(Population,V,(Problem.FE/Problem.maxFE)^2);

                % Optimization
                while Algorithm.NotTerminated(Population)
                    % Calculate the fitness by shift-based density   SDE (the shift-based density estimation strategy)
                    PopObj = Population.objs;
                    N      = size(PopObj,1);
                    fmax   = max(PopObj,[],1);
                    fmin   = min(PopObj,[],1);
                    PopObj = (PopObj-repmat(fmin,N,1))./repmat(fmax-fmin,N,1);
                    Dis    = inf(N);
                    for i = 1 : N
                        SPopObj = max(PopObj,repmat(PopObj(i,:),N,1));
                        for j = [1:i-1,i+1:N]
                            Dis(i,j) = norm(PopObj(i,:)-SPopObj(j,:));
                        end
                    end
                    Fitness = min(Dis,[],2);

                    if length(Population) >= 2
                        Rank = randperm(length(Population),floor(length(Population)/2)*2);
                    else
                        Rank = [1,1];
                    end
                    Loser  = Rank(1:end/2);
                    Winner = Rank(end/2+1:end);
                    Change = Fitness(Loser) >= Fitness(Winner);
                    Temp   = Winner(Change);
                    Winner(Change) = Loser(Change);
                    Loser(Change)  = Temp;

                    [OffDec,OffVel] = Operator_LMOCSO(Problem,Population(Loser),Population(Winner),Rate);
                    %% FDV
                    iter = Problem.FE/Problem.maxFE;
                    if iter <= Rate
                        Offspring = FDVOperator(Problem,Rate,Acc,OffDec,OffVel);
                    else
                        Offspring = Problem.Evaluation(OffDec,OffVel);
                    end
                    Population = EnvironmentalSelection_LMOCSO([Population,Offspring],V,(Problem.FE/Problem.maxFE)^2);
                end
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function Choose = LastSelection(PopObj1,PopObj2,K,Z,Zmin)
% Select part of the solutions in the last front

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

% This function is written by Zhenshou Song
% Email: zssong@stu.xidian.edu.cn

    PopObj = [PopObj1;PopObj2] - repmat(Zmin,size(PopObj1,1)+size(PopObj2,1),1);
    [N,M]  = size(PopObj);
    N1     = size(PopObj1,1);
    N2     = size(PopObj2,1);
    NZ     = size(Z,1);

    %% Normalization
    % Detect the extreme points
    Extreme = zeros(1,M);
    w       = zeros(M)+1e-6+eye(M);
    for i = 1 : M
        [~,Extreme(i)] = min(max(PopObj./repmat(w(i,:),N,1),[],2));
    end
    % Calculate the intercepts of the hyperplane constructed by the extreme
    % points and the axes
    Hyperplane = PopObj(Extreme,:)\ones(M,1);
    a = 1./Hyperplane;
    if any(isnan(a))
        a = max(PopObj,[],1)';
    end
    % Normalization
    PopObj = PopObj./repmat(a',N,1);

    %% Associate each solution with one reference point
    % Calculate the distance of each solution to each reference vector
    Cosine   = 1 - pdist2(PopObj,Z,'cosine');
    Distance = repmat(sqrt(sum(PopObj.^2,2)),1,NZ).*sqrt(1-Cosine.^2);
    % Associate each solution with its nearest reference point
    [d,pi] = min(Distance',[],1);

    %% Calculate the number of associated solutions except for the last front of each reference point
    rho = hist(pi(1:N1),1:NZ);

    %% Environmental selection
    Choose  = false(1,N2);
    Zchoose = true(1,NZ);
    % Select K solutions one by one
    while sum(Choose) < K
        % Select the least crowded reference point
        Temp = find(Zchoose);
        Jmin = find(rho(Temp)==min(rho(Temp)));
        j    = Temp(Jmin(randi(length(Jmin))));
        I    = find(Choose==0 & pi(N1+1:end)==j);
        % Then select one solution associated with this reference point
        if ~isempty(I)
            if rho(j) == 0
                [~,s] = min(d(N1+I));
            else
                s = randi(length(I));
            end
            Choose(I(s)) = true;
            rho(j) = rho(j) + 1;
        else
            Zchoose(j) = false;
        end
    end
end


text:
function [model] = BCPF_IC(Y, varargin)
% Bayesian CP Factorization for Image Completion
% Author : Qibin Zhao  2014
%
% -----------------------------------------------------------------------
%  [model] = BCPF_IC(Y, 'PARAM1', val1, 'PARAM2', val2, ...)
%
%  INPUTS
%     Y              - Input tensor
%     'obs'          - Binary (0-1) tensor indicating missing entries
%                      (0: missing; 1: observed)
%     'init'         - Initialization method
%                     - 'ml'  : SVD initilization (default)
%                     - 'rand': Random matrices
%     'maxRank'      - The initialization of rank (larger than true rank)
%     'dimRed'       - 1: Remove unnecessary components automaticly (default)
%                    - 0: Not remove
%     'maxiters'     - max number of iterations (default: 100)
%     'tol'          - lower band change tolerance for convergence dection
%                      (default: 1e-5)
%     'noise'        - whether noise is updated
%                        - 'on': update noise parameter (default)
%                        - 'off': fixed noise parameter (1e-5)
%     'predVar'      - Predictive distribution
%                         - 1:  compute and output
%                         - 0:  doesnot compute  (default)
%     'verbose'      - visualization of results
%                       - 0: no
%                       - 1: text (default)
%                       - 2: online display image
%                       - 3: show factors by image
%                       - 4: show factors by hinton plot (very slow)
%   OUTPUTS
%      model         - Model parameters and hyperparameters
% -----------------------------------------------------------------------
%
%   Example:
%
%     [model] = BCPF_IC(Y, 'obs', O, 'init', 'rand', 'maxRank', 10, 'dimRed', 1, 'maxiters', 100, ...
%                                'tol', 1e-6, 'verbose', 3);
%
% < Bayesian CP Factorization of Incomplete Image >
% Copyright (C) 2014  Qibin Zhao
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%
warning off; %#ok<WNOFF>
randn('state',1); rand('state',1); %#ok<RAND>
dimY = size(Y);
N = ndims(Y);

%% Set parameters from input or by using defaults
ip = inputParser;
ip.addParamValue('obs', ones(dimY), @(x) (isnumeric(x) || islogical(x)) );
ip.addParamValue('init', 'rand', @(x) (ismember(x,{'ml','rand'})));
ip.addParamValue('maxRank', max(dimY), @isscalar);
ip.addParamValue('maxiters', 100, @isscalar);
ip.addParamValue('tol', 1e-5, @isscalar);
ip.addParamValue('verbose', 1, @isscalar);
ip.addParamValue('noise', 'on', @(x)ismember(x,{'on','off'}));
ip.addParamValue('dimRed', 1, @isscalar);
ip.addParamValue('predVar', 0, @isscalar);
ip.parse(varargin{:});

O     = ip.Results.obs;
init  = ip.Results.init;
R   = ip.Results.maxRank;
maxiters  = ip.Results.maxiters;
tol   = ip.Results.tol;
verbose  = ip.Results.verbose;
DIMRED   = ip.Results.dimRed;
noise = ip.Results.noise;
predVar = ip.Results.predVar;

%% Initialization
Y = tensor(Y.*O);
O = tensor(O);
nObs = sum(O(:));

a_gamma0     = 1e-6;
b_gamma0     = 1e-6;
if  strcmp(noise,'on')
    a_beta0      = 1e-6;
    b_beta0      = 1e-6;
else
    a_beta0      = 1e-1;
    b_beta0      = 1e-6;
end
gammas = ones(R,1);
beta = 1e4;
dscale = 1;

switch init,
    case 'ml'    % Maximum likelihood
        Z = cell(N,1);
        ZSigma = cell(N,1);
        if ~isempty(find(O==0))
            Y(find(O==0)) = sum(Y(:))/nObs;
        end
        for n = 1:N
            ZSigma{n} = (repmat(eye(R), [1 1 dimY(n)]));
            [U, S, V] = svd(double(tenmat(Y,n)), 'econ');
            if R <= size(U,2)
                Z{n} = U(:,1:R)*(S(1:R,1:R)).^(0.5);
            else
                Z{n} = [U*(S.^(0.5)) randn(dimY(n), R-size(U,2))];
            end
        end
        Y = Y.*O;
    case 'rand'   % Random initialization
        Z = cell(N,1);
        ZSigma = cell(N,1);
        for n = 1:N
            Z{n} = randn(dimY(n),R);
            ZSigma{n} = repmat(eye(R), [1 1 dimY(n)]);
        end
end
% --------- E(aa') = cov(a,a) + E(a)E(a')----------------
EZZT = cell(N,1);
for n=1:N
    EZZT{n} = (reshape(ZSigma{n}, [R*R, dimY(n)]))';
end

Fit =0;
LB = 0;
X = double(ktensor(Z));

%% Create figures
if verbose >2,
    scrsz = get(0,'ScreenSize');
    h1 = figure('Position',[scrsz(3)*0.2 scrsz(4)*0.3 scrsz(3)*0.6 scrsz(4)*0.4]);
    figure(h1);
    switch verbose,
        case 4,
            subplot(2,3,1); hintonDiagram(Z{1}); title('Mode-1'); ylabel('Length of #-mode');
            subplot(2,3,2); hintonDiagram(Z{2}); title('Mode-2'); xlabel('Latent dimensions');
            if N>=3, subplot(2,3,3); hintonDiagram(Z{3}); title('Mode-3'); end
        case 3,
            subplot(2,3,1); imagesc(Z{1}); title('Mode-1'); ylabel('Length of #-mode');
            subplot(2,3,2); imagesc(Z{2}); title('Mode-2'); xlabel('Latent dimensions');
            if N>=3, subplot(2,3,3); imagesc(Z{3}); title('Mode-3');end
    end
    subplot(2,3,4); bar(gammas); title('Posterior mean of \lambda'); xlabel('Latent components'); ylabel(''); axis tight;
    subplot(2,3,5); plot(LB, '-r.','LineWidth',1.5,'MarkerSize',10 ); title('Lower bound'); xlabel('Iteration');  grid on;
    subplot(2,3,6); plotGamma(a_beta0, a_beta0); title('Posterior pdf'); xlabel('Noise precision \tau');grid on;
    set(findall(h1,'type','text'),'fontSize',12);
    drawnow;
end
if verbose ==2;
    h3 = figure;
    temp = 255.*(X-min(X(:)))/(max(X(:))-min(X(:)));
    imshow(uint8(temp));
    title(['Iter.= '  num2str(0),',  Rank = ' num2str(R)],'FontSize', 13, 'color','b');
    tic;
    xlabel(['(BCPF)  Time: ' num2str(round(toc)), ' seconds'],'FontSize', 13, 'color','b');
    drawnow;
end

%% Model learning
for it=1:maxiters,
    %% Update factor matrices
    Aw = diag(gammas);
    for n=1:N
        ENZZT = reshape(khatrirao_fast(EZZT{[1:n-1, n+1:N]},'r')' * double(tenmat(O,n)'), [R,R,dimY(n)]);
        FslashY = khatrirao_fast(Z{[1:n-1, n+1:N]},'r')' * tenmat(Y.*O, n)';
        for i=1:dimY(n)
            ZSigma{n}(:,:,i) = (beta * ENZZT(:,:,i) + Aw )^(-1);
            Z{n}(i,:) = (beta * ZSigma{n}(:,:,i) * FslashY(:,i))';
        end
        EZZT{n} = (reshape(ZSigma{n}, [R*R, dimY(n)]) + khatrirao_fast(Z{n}',Z{n}'))';
    end

    %% Update latent tensor X
    X = double(ktensor(Z));

    %% Update hyperparameters gamma
    a_gammaN = (0.5*sum(dimY) + a_gamma0)*ones(R,1);
    b_gammaN = 0;
    for n=1:N
        b_gammaN = b_gammaN + diag(Z{n}'*Z{n}) + diag(sum(ZSigma{n},3));
    end
    b_gammaN = b_gamma0 + 0.5.* b_gammaN;
    gammas = a_gammaN./b_gammaN;

    %% update noise beta
    %  The most time and space consuming part
    if 0 % save time but large space needed
        EX2 =  O(:)' * khatrirao_fast(EZZT,'r') * ones(R*R,1);
    else  % save space but slow
        temp1 = cell(N,1);
        EX2 =0;
        for i =1:R
            for n=1:N
                temp1{n} = EZZT{n}(:,(i-1)*R+1: i*R);
            end
            EX2 = EX2 + O(:)' * khatrirao_fast(temp1,'r')* ones(R,1);
        end
    end
    err = Y(:)'*Y(:) - 2*Y(:)'*X(:) + EX2;
    if  strcmp(noise,'on')
        a_betaN = a_beta0 + 0.5*nObs;
        b_betaN = b_beta0 + 0.5*err;
    else
        a_betaN = a_beta0;
        b_betaN = b_beta0;
    end
    beta = a_betaN/b_betaN;
    Fit = 1 - sqrt(sum(err(:)))/norm(Y(:));

    %% Lower bound
    temp1 = -0.5*nObs*safelog(2*pi) + 0.5*nObs*(psi(a_betaN)-safelog(b_betaN)) - 0.5*(a_betaN/b_betaN)*err;
    temp22 =0;
    for n=1:N
        temp22= temp22 + Z{n}'*Z{n} + sum(ZSigma{n},3);
    end
    temp2 = -0.5*R*sum(dimY)*safelog(2*pi) + 0.5*sum(dimY)*sum(psi(a_gammaN)-safelog(b_gammaN)) -0.5*trace(diag(gammas)* temp22);
    temp3 = sum(-safelog(gamma(a_gamma0)) + a_gamma0*safelog(b_gamma0) -  b_gamma0.*(a_gammaN./b_gammaN) + (a_gamma0-1).*(psi(a_gammaN)-safelog(b_gammaN)));
    temp4 = -safelog(gamma(a_beta0)) + a_beta0*safelog(b_beta0) + (a_beta0-1)*(psi(a_betaN)-safelog(b_betaN)) - b_beta0*(a_betaN/b_betaN);
    temp5=0;
    for n=1:N
        for i=1:size(ZSigma{n},3)
            temp5 = temp5 + 0.5*safelog(det(ZSigma{n}(:,:,i))) + 0.5*R*(1+safelog(2*pi));
        end
    end
    temp6 = sum(safelog(gamma(a_gammaN)) - (a_gammaN-1).*psi(a_gammaN) -safelog(b_gammaN) + a_gammaN);
    temp7 = safelog(gamma(a_betaN)) - (a_betaN-1)*psi(a_betaN) -safelog(b_betaN) + a_betaN;
    LB(it) = temp1 + temp2 + temp3 + temp4 + temp5 + temp6 + temp7;

    %% Prune irrelevant dimensions?
    Zall = cell2mat(Z);
    comPower = diag(Zall' * Zall);
    comTol = sum(dimY)*eps(norm(Zall,'fro'));
    rankest = sum(comPower> comTol );
    if max(rankest)==0
        disp('Rank becomes 0 !!!');
        break;
    end
    if DIMRED==1  && it >=2,
        if R~= max(rankest)
            indices = comPower > comTol;
            gammas = gammas(indices);
            temp = ones(R,R);
            temp(indices,indices) = 0;
            temp = temp(:);
            for n=1:N
                Z{n} = Z{n}(:,indices);
                ZSigma{n} = ZSigma{n}(indices,indices,:);
                EZZT{n} = EZZT{n}(:, temp == 0);
            end
            R = max(rankest);
        end
    end

    %% Display progress
    if it>2
        LBRelChan = abs(LB(it) - 2*LB(it-1) + LB(it-2))/-LB(2);
    else
        LBRelChan = NaN;
    end
    if verbose,
        fprintf('Iter. %d: RelChan = %g, Fit = %g, R = %d \n', it, LBRelChan, Fit, rankest);
    end

    %% visualize online results
    if verbose >2 ,
        switch verbose,
            case 4,
                set(0,'CurrentFigure',h1);
                subplot(2,3,1); hintonDiagram(Z{1}); title('Mode-1'); ylabel('Length of #-mode');
                subplot(2,3,2); hintonDiagram(Z{2}); title('Mode-2'); xlabel('Latent dimensions');
                if N>=3, subplot(2,3,3); hintonDiagram(Z{3}); title('Mode-3'); end
            case 3,
                set(0,'CurrentFigure',h1);
                subplot(2,3,1); imagesc(Z{1}); title('Mode-1'); ylabel('Length of #-mode');
                subplot(2,3,2); imagesc(Z{2}); title('Mode-2'); xlabel('Latent dimensions');
                if N>=3, subplot(2,3,3); imagesc(Z{3}); title('Mode-3'); end
        end
        subplot(2,3,4); bar(gammas); title('Posterior mean of \lambda'); xlabel('Latent components'); ylabel(''); axis tight;
        subplot(2,3,5); plot(LB, '-r.','LineWidth',1.5,'MarkerSize',10 ); title('Lower bound'); xlabel('Iteration');  grid on;
        subplot(2,3,6); plotGamma(a_betaN, b_betaN); title('Posterior pdf'); xlabel('Noise precision \tau');grid on;
        set(findall(h1,'type','text'),'fontSize',12);
        drawnow;
    end
    if verbose==2
        set(0,'CurrentFigure',h3);
        figure(h3);
        %        temp = (X-min(X(:)))/(max(X(:))-min(X(:)));
        %        image(temp);
        %        axis off;
        imshow(uint8(X));
        title(['Iter.= '  num2str(it),',  Rank = ' num2str(max(rankest))],'FontSize', 13, 'color','b');
        xlabel(['(BCPF)  Time: ' num2str(round(toc)), ' seconds'],'FontSize', 13, 'color','b');
        drawnow;
    end

    %% Convergence check
    if it>5 && abs(LBRelChan) < tol
        disp('\\\======= Converged===========\\\');
        break;
    end
end

%% Predictive distribution
if predVar==1
    Xvar =  tenzeros(size(Y));
    for n=1:N
        Xvar = tenmat(Xvar,n);
        Fslash = khatrirao_fast(Z{[1:n-1, n+1:N]},'r');
        if 1
            temp1 = double(tenmat(tensor(ZSigma{n}),3));
            temp2 = khatrirao_fast(Fslash', Fslash');
            Xvar(:,:) = Xvar(:,:) + temp1*temp2;
        else
            % ---  slow computation ------
            for i=1:size(Xvar,1)     %#ok
                Xvar(i,:) = Xvar(i,:) + diag(Fslash * ZSigma{n}(:,:,i) *Fslash')';
            end
            % ---  slow computation ------
        end
        Xvar = tensor(Xvar);
    end
    Xvar = Xvar + beta^(-1);
    Xvar = Xvar.*(2*a_betaN)/(2*a_betaN-2);
    Xvar = Xvar.*(dscale^2);
else
    Xvar =[];
end

%% Prepare the results
SNR = 10*log10(var(X(:))*beta);
X = ktensor(Z)*dscale;
X = arrange(X);

%% Output
model.X = X;
model.ZSigma = ZSigma;
model.gammas = gammas;
model.Fit = Fit;
model.SNR = SNR;
model.Xvar = double(Xvar);
model.TrueRank = rankest;
model.LowBound = max(LB);

function y = safelog(x)
x(x<1e-300)=1e-200;
x(x>1e300)=1e300;
y=log(x);

author: andrewssobral


text:
function showHOG(w)

hogim = HOGpicture(w);
imagesc(hogim)
axis image
axis off
grid on
drawnow

author: zouchuhang

repo: LayoutNet


text:
function p08_header ( )

%*****************************************************************************80
%
%% P08_HEADER prints some information about problem 08.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Per-Olof Persson and Gilbert Strang,
%    A Simple Mesh Generator in MATLAB,
%    SIAM Review,
%    Volume 46, Number 2, June 2004, pages 329-345.
%
%  Parameters:
%
%    None.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'P08:\n' );
  fprintf ( 1, '  Strang and Persson example #7\n' );
  fprintf ( 1, '  Pie slice with notch and hole.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  A uniform mesh density is requested.\n' );

  fprintf ( 1, '\n' );
  boundary_segment_num = p08_boundary_segment_num ( );
  fprintf ( 1, '  Number of boundary segments = %d\n', boundary_segment_num );
  fixed_num = p08_fixed_num ( );
  fprintf ( 1, '  Number of fixed points = %d\n', fixed_num );
  hole_num = p08_hole_num ( );
  fprintf ( 1, '  Number of holes = %d\n', hole_num );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [s, err_mse, iter_time]=greed_omp_linsolve(x,A,m,varargin)
% greed_omp_linsolve: Orthogonal Matching Pursuit algorithm based on matlab
% linsolve solution. For reference only!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Usage
% [s, err_mse, iter_time]=greed_omp_linsolve(x,P,m,'option_name','option_value')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input
%   Mandatory:
%               x   Observation vector to be decomposed
%               P   Either:
%                       1) An nxm matrix (n must be dimension of x)
%                       2) A function handle (type "help function_format"
%                          for more information)
%                          Also requires specification of P_trans option.
%                       3) An object handle (type "help object_format" for
%                          more information)
%               m   length of s
%
%   Possible additional options:
%   (specify as many as you want using 'option_name','option_value' pairs)
%   See below for explanation of options:
%__________________________________________________________________________
%   option_name    |     available option_values                | default
%--------------------------------------------------------------------------
%   stopCrit       | M, corr, mse, mse_change                   | M
%   stopTol        | number (see below)                         | n/4
%   P_trans        | function_handle (see below)                |
%   maxIter        | positive integer (see below)               | n
%   verbose        | true, false                                | false
%   start_val      | vector of length m                         | zeros
%
%   Available stopping criteria :
%               M           -   Extracts exactly M = stopTol elements.
%               corr        -   Stops when maximum correlation between
%                               residual and atoms is below stopTol value.
%               mse         -   Stops when mean squared error of residual
%                               is below stopTol value.
%               mse_change  -   Stops when the change in the mean squared
%                               error falls below stopTol value.
%
%   stopTol: Value for stopping criterion.
%
%   P_trans: If P is a function handle, then P_trans has to be specified and
%            must be a function handle.
%
%   maxIter: Maximum number of allowed iterations.
%
%   verbose: Logical value to allow algorithm progress to be displayed.
%
%   start_val: Allows algorithms to start from partial solution.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Outputs
%	 s              Solution vector
%    err_mse        Vector containing mse of approximation error for each
%                   iteration
%    iter_time      Vector containing computation times for each iteration
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Description
%   greed_omp_linsolve performs a greedy signal decomposition.
%   In each iteration a new element is selected depending on the inner
%   product between the current residual and columns in P.
%   The non-zero elements of s are approximated by orthogonally projecting
%   x onto the selected elements in each iteration.
%   The inverse problem is solved using matlabs linsolve command.
%   This is slow and is only provided for reference.
%
% See Also
%   greed_qr, greed_omp_chol, greed_omp_cg, greed_omp_cgp, greed_omp_pinv,
%   greed_gp, greed_nomp
%
% Copyright (c) 2007 Thomas Blumensath
%
% The University of Edinburgh
% Email: thomas.blumensath@ed.ac.uk
% Comments and bug reports welcome
%
% This file is part of sparsity Version 0.1
% Created: April 2007
%
% Part of this toolbox was developed with the support of EPSRC Grant
% D000246/1
%
% Please read COPYRIGHT.m for terms and conditions.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                    Default values and initialisation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[n1 n2]=size(x);
if n2 == 1
    n=n1;
elseif n1 == 1
    x=x';
    n=n2;
else
   display('x must be a vector.');
   return
end

sigsize     = x'*x/n;
initial_given=0;
err_mse     = [];
iter_time   = [];
STOPCRIT    = 'M';
STOPTOL     = ceil(n/4);
MAXITER     = n;
verbose     = false;
s_initial   = zeros(m,1);
vectnfact   = ones(m,1);

if verbose
   display('Initialising...')
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                           Output variables
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

switch nargout
    case 3
        comp_err=true;
        comp_time=true;
    case 2
        comp_err=true;
        comp_time=false;
    case 1
        comp_err=false;
        comp_time=false;
    case 0
        error('Please assign output variable.')
    otherwise
        error('Too many output arguments specified')
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                       Look through options
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Put option into nice format
Options={};
OS=nargin-3;
c=1;
for i=1:OS
    if isa(varargin{i},'cell')
        CellSize=length(varargin{i});
        ThisCell=varargin{i};
        for j=1:CellSize
            Options{c}=ThisCell{j};
            c=c+1;
        end
    else
        Options{c}=varargin{i};
        c=c+1;
    end
end
OS=length(Options);
if rem(OS,2)
   error('Something is wrong with argument name and argument value pairs.')
end

for i=1:2:OS
   switch Options{i}
        case {'stopCrit'}
            if (strmatch(Options{i+1},{'M'; 'corr'; 'mse'; 'mse_change'},'exact'));
                STOPCRIT    = Options{i+1};
            else error('stopCrit must be char string [M, corr, mse, mse_change]. Exiting.'); end
        case {'stopTol'}
            if isa(Options{i+1},'numeric') ; STOPTOL     = Options{i+1};
            else error('stopTol must be number. Exiting.'); end
        case {'P_trans'}
            if isa(Options{i+1},'function_handle'); Pt = Options{i+1};
            else error('P_trans must be function _handle. Exiting.'); end
        case {'maxIter'}
            if isa(Options{i+1},'numeric'); MAXITER     = Options{i+1};
            else error('maxIter must be a number. Exiting.'); end
        case {'verbose'}
            if isa(Options{i+1},'logical'); verbose     = Options{i+1};
            else error('verbose must be a logical. Exiting.'); end
        case {'vecNormFac'}
            if isa(Options{i+1},'numeric')& length(Options{i+1}) == m , vectnfact = Options{i+1};
            else error('verbose must be a logical. Exiting.'); end
        case {'start_val'}
            if isa(Options{i+1},'numeric') & length(Options{i+1}) == m ;
                s_initial     = Options{i+1};
                initial_given=1;
            else error('start_val must be a vector of length m. Exiting.'); end
        otherwise
            error('Unrecognised option. Exiting.')
   end
end

if strcmp(STOPCRIT,'M')
    maxM=STOPTOL;
else
    maxM=MAXITER;
end

if nargout >=2
    err_mse = zeros(maxM,1);
end
if nargout ==3
    iter_time = zeros(maxM,1);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                        Make P and Pt functions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if          isa(A,'float')      P =@(z) A*z;  Pt =@(z) A'*z;
elseif      isobject(A)         P =@(z) A*z;  Pt =@(z) A'*z;
elseif      isa(A,'function_handle')
    try
        if          isa(Pt,'function_handle'); P=A;
        else        error('If P is a function handle, Pt also needs to be a function handle. Exiting.'); end
    catch error('If P is a function handle, Pt needs to be specified. Exiting.'); end
else        error('P is of unsupported type. Use matrix, function_handle or object. Exiting.'); end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                        Do we start from zero or not?
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if initial_given ==1;
    IN          = find(s_initial);
    s=zeros(m,1);
    if isa(A,'function_handle') || isobject(A)
         Pmat=zeros(n,length(IN));
         for i=1:length(IN)
             mask=zeros(m,1);
             mask(IN(i))=1;
             Pmat(:,i)=P(mask);
         end
         s(IN)=linsolve(Pmat,x);
     else
         s(IN)=linsolve(A(:,IN),x);
     end
    Residual    = x-P(s);
    oldERR      = Residual'*Residual/n;

else
    IN          = [];
    Residual    = x;
    s           = s_initial;
    sigsize     = x'*x/n;
    oldERR      = sigsize;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                 Random Check to see if dictionary is normalised
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%         mask=zeros(m,1);
%         mask(ceil(rand*m))=1;
%         nP=norm(P(mask));
%         if abs(1-nP)>1e-3;
%             display('Dictionary appears not to have unit norm columns.')
%         end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                        Main algorithm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if verbose
   display('Main iterations...')
end
tic
t=0;
DR=Pt(Residual).*vectnfact;
done = 0;
iter=1;
while ~done
     DR(IN)=0;
     [v I]=max(abs(DR));
     IN=[IN I];
     if isa(A,'function_handle') || isobject(A)
         Pmat=zeros(n,length(IN));
         for i=1:length(IN)
             mask=zeros(m,1);
             mask(IN(i))=1;
             Pmat(:,i)=P(mask);
         end
         s(IN)=linsolve(Pmat,x);
     else
         s(IN)=linsolve(A(:,IN),x);
     end
     Residual=x-P(s);
     DR=Pt(Residual).*vectnfact;

     ERR=Residual'*Residual/n;
     if comp_err
         err_mse(iter)=ERR;
     end

     if comp_time
         iter_time(iter)=toc;
     end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                        Are we done yet?
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

     if strcmp(STOPCRIT,'M')
         if iter >= STOPTOL
             done =1;
         elseif verbose && toc-t>10
            display(sprintf('Iteration %i. --- %i iterations to go',iter ,STOPTOL-iter))
            t=toc;
         end
    elseif strcmp(STOPCRIT,'mse')
         if comp_err
            if err_mse(iter)<STOPTOL;
                done = 1;
            elseif verbose && toc-t>10
                display(sprintf('Iteration %i. --- %i mse',iter ,err_mse(iter)))
                t=toc;
            end
         else
             if ERR<STOPTOL;
                done = 1;
             elseif verbose && toc-t>10
                display(sprintf('Iteration %i. --- %i mse',iter ,ERR))
                t=toc;
             end
         end
     elseif strcmp(STOPCRIT,'mse_change') && iter >=2
         if comp_err && iter >=2
              if ((err_mse(iter-1)-err_mse(iter))/sigsize <STOPTOL);
                done = 1;
             elseif verbose && toc-t>10
                display(sprintf('Iteration %i. --- %i mse change',iter ,(err_mse(iter-1)-err_mse(iter))/sigsize ))
                t=toc;
             end
         else
             if ((oldERR - ERR)/sigsize < STOPTOL);
                done = 1;
             elseif verbose && toc-t>10
                display(sprintf('Iteration %i. --- %i mse change',iter ,(oldERR - ERR)/sigsize))
                t=toc;
             end
         end
     elseif strcmp(STOPCRIT,'corr')
          if max(abs(DR)) < STOPTOL;
             done = 1;
          elseif verbose && toc-t>10
                display(sprintf('Iteration %i. --- %i corr',iter ,max(abs(DR))))
                t=toc;
          end
     end

    % Also stop if residual gets too small or maxIter reached
     if comp_err
         if err_mse(iter)<1e-16
             display('Stopping. Exact signal representation found!')
             done=1;
         end
     else

         if iter>1
             if ERR<1e-16
                 display('Stopping. Exact signal representation found!')
                 done=1;
             end
         end
     end

     if iter >= MAXITER
         display('Stopping. Maximum number of iterations reached!')
         done = 1;
     end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                    If not done, take another round
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

     if ~done
        iter=iter+1;
        oldERR=ERR;
     end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                  Only return as many elements as iterations
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargout >=2
    err_mse = err_mse(1:iter);
end
if nargout ==3
    iter_time = iter_time(1:iter);
end

if verbose
   display('Done')
end

% Change history
%
% 8 of Februray: Algo does no longer stop if dictionary is not normaliesd.

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
%% Example social DB

clear all;
close all;
clear classes
clc;

graph = prtGraphDataGenKarate;
plot(graph);

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function test_plu ( )

%*****************************************************************************80
%
%% TEST_PLU tests the PLU factors.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 April 2015
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST_PLU\n' );
  fprintf ( 1, '  A = a test matrix of order M by N\n' );
  fprintf ( 1, '  P, L, U are the PLU factors.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  ||A|| = Frobenius norm of A.\n' );
  fprintf ( 1, '  ||A-PLU|| = Frobenius norm of A-P*L*U.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Title                    M     N      ' );
  fprintf ( 1, '||A||            ||A-PLU||\n' );
  fprintf ( 1, '\n' );
%
%  A123
%
  title = 'A123';
  m = 3;
  n = 3;
  a = a123 ( );
  [ p, l, u ] = a123_plu ( );
  error_frobenius = r8mat_is_plu ( m, n, a, p, l, u );
  norm_a_frobenius = r8mat_norm_fro ( m, n, a );
  fprintf ( 1, '  %-20s  %4d  %4d  %14g  %14g\n', ...
    title, m, n, norm_a_frobenius, error_frobenius );
%
%  BODEWIG
%
  title = 'BODEWIG';
  m = 4;
  n = 4;
  a = bodewig ( );
  [ p, l, u ] = bodewig_plu ( );
  error_frobenius = r8mat_is_plu ( m, n, a, p, l, u );
  norm_a_frobenius = r8mat_norm_fro ( m, n, a );
  fprintf ( 1, '  %-20s  %4d  %4d  %14g  %14g\n', ...
    title, m, n, norm_a_frobenius, error_frobenius );
%
%  BORDERBAND
%
  title = 'BORDERBAND';
  m = 5;
  n = 5;
  a = borderband ( n );
  [ p, l, u ] = borderband_plu ( n );
  error_frobenius = r8mat_is_plu ( m, n, a, p, l, u );
  norm_a_frobenius = r8mat_norm_fro ( m, n, a );
  fprintf ( 1, '  %-20s  %4d  %4d  %14g  %14g\n', ...
    title, m, n, norm_a_frobenius, error_frobenius );
%
%  DIF2
%
  title = 'DIF2';
  m = 5;
  n = 5;
  a = dif2 ( m, n );
  [ p, l, u ] = dif2_plu ( n );
  error_frobenius = r8mat_is_plu ( m, n, a, p, l, u );
  norm_a_frobenius = r8mat_norm_fro ( m, n, a );
  fprintf ( 1, '  %-20s  %4d  %4d  %14g  %14g\n', ...
    title, m, n, norm_a_frobenius, error_frobenius );
%
%  GFPP
%
  title = 'GFPP';
  m = 5;
  n = 5;
  r8_lo = -5.0;
  r8_hi = +5.0;
  seed = 123456789;
  [ alpha, seed ] = r8_uniform_ab ( r8_lo, r8_hi, seed );
  a = gfpp ( n, alpha );
  [ p, l, u ] = gfpp_plu ( n, alpha );
  error_frobenius = r8mat_is_plu ( m, n, a, p, l, u );
  norm_a_frobenius = r8mat_norm_fro ( m, n, a );
  fprintf ( 1, '  %-20s  %4d  %4d  %14g  %14g\n', ...
    title, m, n, norm_a_frobenius, error_frobenius );
%
%  GIVENS
%
  title = 'GIVENS';
  m = 5;
  n = 5;
  a = givens ( m, n );
  [ p, l, u ] = givens_plu ( n );
  error_frobenius = r8mat_is_plu ( m, n, a, p, l, u );
  norm_a_frobenius = r8mat_norm_fro ( m, n, a );
  fprintf ( 1, '  %-20s  %4d  %4d  %14g  %14g\n', ...
    title, m, n, norm_a_frobenius, error_frobenius );
%
%  KMS
%
  title = 'KMS';
  m = 5;
  n = 5;
  r8_lo = -5.0;
  r8_hi = +5.0;
  seed = 123456789;
  [ alpha, seed ] = r8_uniform_ab ( r8_lo, r8_hi, seed );
  a = kms ( alpha, m, n );
  [ p, l, u ] = kms_plu ( alpha, n );
  error_frobenius = r8mat_is_plu ( m, n, a, p, l, u );
  norm_a_frobenius = r8mat_norm_fro ( m, n, a );
  fprintf ( 1, '  %-20s  %4d  %4d  %14g  %14g\n', ...
    title, m, n, norm_a_frobenius, error_frobenius );
%
%  LEHMER
%
  title = 'LEHMER';
  m = 5;
  n = 5;
  a = lehmer ( m, n );
  [ p, l, u ] = lehmer_plu ( n );
  error_frobenius = r8mat_is_plu ( m, n, a, p, l, u );
  norm_a_frobenius = r8mat_norm_fro ( m, n, a );
  fprintf ( 1, '  %-20s  %4d  %4d  %14g  %14g\n', ...
    title, m, n, norm_a_frobenius, error_frobenius );
%
%  MAXIJ
%
  title = 'MAXIJ';
  m = 5;
  n = 5;
  a = maxij ( m, n );
  [ p, l, u ] = maxij_plu ( n );
  error_frobenius = r8mat_is_plu ( m, n, a, p, l, u );
  norm_a_frobenius = r8mat_norm_fro ( m, n, a );
  fprintf ( 1, '  %-20s  %4d  %4d  %14g  %14g\n', ...
    title, m, n, norm_a_frobenius, error_frobenius );
%
%  MINIJ
%
  title = 'MINIJ';
  m = 5;
  n = 5;
  a = minij ( n, n );
  [ p, l, u ] = minij_plu ( n );
  error_frobenius = r8mat_is_plu ( m, n, a, p, l, u );
  norm_a_frobenius = r8mat_norm_fro ( m, n, a );
  fprintf ( 1, '  %-20s  %4d  %4d  %14g  %14g\n', ...
    title, m, n, norm_a_frobenius, error_frobenius );
%
%  MOLER1
%
  title = 'MOLER1';
  m = 5;
  n = 5;
  r8_lo = -5.0;
  r8_hi = +5.0;
  seed = 123456789;
  [ alpha, seed ] = r8_uniform_ab ( r8_lo, r8_hi, seed );
  a = moler1 ( alpha, n, n );
  [ p, l, u ] = moler1_plu ( alpha, n );
  error_frobenius = r8mat_is_plu ( m, n, a, p, l, u );
  norm_a_frobenius = r8mat_norm_fro ( m, n, a );
  fprintf ( 1, '  %-20s  %4d  %4d  %14g  %14g\n', ...
    title, m, n, norm_a_frobenius, error_frobenius );
%
%  MOLER3
%
  title = 'MOLER3';
  m = 5;
  n = 5;
  a = moler3 ( m, n );
  [ p, l, u ] = moler3_plu ( n );
  error_frobenius = r8mat_is_plu ( m, n, a, p, l, u );
  norm_a_frobenius = r8mat_norm_fro ( m, n, a );
  fprintf ( 1, '  %-20s  %4d  %4d  %14g  %14g\n', ...
    title, m, n, norm_a_frobenius, error_frobenius );
%
%  OTO
%
  title = 'OTO';
  m = 5;
  n = 5;
  a = oto ( m, n );
  [ p, l, u ] = oto_plu ( n );
  error_frobenius = r8mat_is_plu ( m, n, a, p, l, u );
  norm_a_frobenius = r8mat_norm_fro ( m, n, a );
  fprintf ( 1, '  %-20s  %4d  %4d  %14g  %14g\n', ...
    title, m, n, norm_a_frobenius, error_frobenius );
%
%  PASCAL2
%
  title = 'PASCAL2';
  m = 5;
  n = 5;
  a = pascal2 ( n );
  [ p, l, u ] = pascal2_plu ( n );
  error_frobenius = r8mat_is_plu ( m, n, a, p, l, u );
  norm_a_frobenius = r8mat_norm_fro ( m, n, a );
  fprintf ( 1, '  %-20s  %4d  %4d  %14g  %14g\n', ...
    title, m, n, norm_a_frobenius, error_frobenius );
%
%  PLU
%
  title = 'PLU';
  n = 5;
  seed = 123456789;
  for i = 1 : n
    i4_lo = i;
    i4_hi = n;
    [ j, seed ] = i4_uniform_ab ( i4_lo, i4_hi, seed );
    pivot(i) = j;
  end
  a = plu ( n, pivot );
  [ p, l, u ] = plu_plu ( n, pivot );
  error_frobenius = r8mat_is_plu ( m, n, a, p, l, u );
  norm_a_frobenius = r8mat_norm_fro ( m, n, a );
  fprintf ( 1, '  %-20s  %4d  %4d  %14g  %14g\n', ...
    title, m, n, norm_a_frobenius, error_frobenius );
%
%  VAND2
%
  title = 'VAND2';
  m = 4;
  n = 4;
  r8_lo = -5.0;
  r8_hi = +5.0;
  seed = 123456789;
  [ x, seed ] = r8vec_uniform_ab ( m, r8_lo, r8_hi, seed );
  a = vand2 ( m, x );
  [ p, l, u ] = vand2_plu ( m, x );
  error_frobenius = r8mat_is_plu ( m, n, a, p, l, u );
  norm_a_frobenius = r8mat_norm_fro ( m, n, a );
  fprintf ( 1, '  %-20s  %4d  %4d  %14g  %14g\n', ...
    title, m, n, norm_a_frobenius, error_frobenius );
%
%  WILSON
%
  title = 'WILSON';
  m = 4;
  n = 4;
  a = wilson ( );
  [ p, l, u ] = wilson_plu ( );
  error_frobenius = r8mat_is_plu ( m, n, a, p, l, u );
  norm_a_frobenius = r8mat_norm_fro ( m, n, a );
  fprintf ( 1, '  %-20s  %4d  %4d  %14g  %14g\n', ...
    title, m, n, norm_a_frobenius, error_frobenius );

  return
end


text:
function g = sigmoidGradient(z)
%SIGMOIDGRADIENT returns the gradient of the sigmoid function
%evaluated at z
%   g = SIGMOIDGRADIENT(z) computes the gradient of the sigmoid function
%   evaluated at z. This should work regardless if z is a matrix or a
%   vector. In particular, if z is a vector or matrix, you should return
%   the gradient for each element.

g = zeros(size(z));
g = sigmoid(z) .* (1 - sigmoid(z));
% ====================== YOUR CODE HERE ======================
% Instructions: Compute the gradient of the sigmoid function evaluated at
%               each value of z (z can be a matrix, vector or scalar).

% =============================================================

end

author: Ayatans

repo: Machine-Learning-homework


text:
function [p,feasible,vol_reduction,seen_x] = propagate_bounds_lp(p,upper,lower,lpsolver,xmin);
LU = [p.lb p.ub];
seen_x = {};
if ~p.options.bmibnb.lpreduce | ((size(p.lpcuts,1)==0) & (any(p.lb(p.linears)<-1e8) & any(p.ub(p.linears)>1e8)))
    vol_reduction = 1;
    p.feasible = 1;
    p.lb(p.integer_variables) = ceil(p.lb(p.integer_variables)-1e-7);
    p.ub(p.integer_variables) = floor(p.ub(p.integer_variables)+1e-7);
    p.lb(p.binary_variables) = ceil(p.lb(p.binary_variables)-1e-7);
    p.ub(p.binary_variables) = floor(p.ub(p.binary_variables)+1e-7);
else
    try
        [p,p.feasible,seen_x] =  boxreduce(p,upper,lower,lpsolver,p.options,xmin);
    catch
        % LINPROG easiy crashes in poorly constrained cases
    end
end
if ~isequal(LU,[p.lb p.ub])
    p.changedbounds = 1;
end
feasible = p.feasible;
vol_reduction = 0;

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
%SerialLink.gravload Gravity load on joints
%
% TAUG = R.gravload(Q) is the joint gravity loading (1xN) for the robot R
% in the joint configuration Q (1xN), where N is the number of robot
% joints.  Gravitational acceleration is a property of the robot object.
%
% If Q is a matrix (MxN) each row is interpreted as a joint configuration
% vector, and the result is a matrix (MxN) each row being the corresponding
% joint torques.
%
% TAUG = R.gravload(Q, GRAV) as above but the gravitational
% acceleration vector GRAV is given explicitly.
%
% See also SerialLink.gravjac, SerialLink.rne, SerialLink.itorque, SerialLink.coriolis.

% Copyright (C) 1993-2017, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

function tg = gravload(robot, q, grav)
	assert(numcols(q) == robot.n, 'RTB:SerialLink:gravload:badarg', 'Insufficient columns in q');

	if nargin == 2
		tg = rne(robot, q, zeros(size(q)), zeros(size(q)));
	elseif nargin == 3
		tg = rne(robot, q, zeros(size(q)), zeros(size(q)), grav);
    end

author: petercorke

repo: robotics-toolbox-matlab


text:
Network UNet {
	Layer CONV1_1 {
		Type: CONV
		Dimensions { K: 64, C:1, R: 3, S: 3, Y: 572, X: 572 }
		Dataflow {
			SpatialMap(1,1) K;
			TemporalMap(64,64) C;
			TemporalMap(7,7) R;
			TemporalMap(7,7) S;
			TemporalMap(7,1) Y;
			TemporalMap(7,1) X;
			Cluster(64, P);
			SpatialMap(1,1) C;
		}
	}

	Layer CONV1_2 {
		Type: CONV
		Dimensions { K: 64, C: 64, R: 3, S: 3, Y: 570, X: 570 }
		Dataflow {
			SpatialMap(1,1) K;
			TemporalMap(64,64) C;
			TemporalMap(3,3) R;
			TemporalMap(3,3) S;
			TemporalMap(3,1) Y;
			TemporalMap(3,1) X;
			Cluster(64, P);
			SpatialMap(1,1) C;
		}
	}

	Layer CONV2_1 {
		Type: CONV
		Dimensions { K: 128, C: 64, R: 3, S: 3, Y: 284, X: 284 }
		Dataflow {
			SpatialMap(1,1) K;
			TemporalMap(64,64) C;
			TemporalMap(3,3) R;
			TemporalMap(3,3) S;
			TemporalMap(3,1) Y;
			TemporalMap(3,1) X;
			Cluster(64, P);
			SpatialMap(1,1) C;
		}
	}

	Layer CONV2_2 {
		Type: CONV
		Dimensions { K: 128, C: 128, R: 3, S: 3, Y: 282, X: 282 }
		Dataflow {
			SpatialMap(1,1) K;
			TemporalMap(64,64) C;
			TemporalMap(3,3) R;
			TemporalMap(3,3) S;
			TemporalMap(3,1) Y;
			TemporalMap(3,1) X;
			Cluster(64, P);
			SpatialMap(1,1) C;
		}
	}

	Layer CONV3_1 {
		Type: CONV
		Dimensions { K: 256, C: 128, R: 3, S: 3, Y: 140, X: 140 }
		Dataflow {
			SpatialMap(1,1) K;
			TemporalMap(64,64) C;
			TemporalMap(3,3) R;
			TemporalMap(3,3) S;
			TemporalMap(3,1) Y;
			TemporalMap(3,1) X;
			Cluster(64, P);
			SpatialMap(1,1) C;
		}
	}

	Layer CONV3_2 {
		Type: CONV
		Dimensions { K: 256, C: 256, R: 3, S: 3, Y: 138, X: 138 }
		Dataflow {
			SpatialMap(1,1) K;
			TemporalMap(64,64) C;
			TemporalMap(3,3) R;
			TemporalMap(3,3) S;
			TemporalMap(3,1) Y;
			TemporalMap(3,1) X;
			Cluster(64, P);
			SpatialMap(1,1) C;
		}
	}

	Layer CONV4_1 {
		Type: CONV
		Dimensions { K: 512, C: 256, R: 3, S: 3, Y: 68, X: 68 }
		Dataflow {
			SpatialMap(1,1) K;
			TemporalMap(64,64) C;
			TemporalMap(3,3) R;
			TemporalMap(3,3) S;
			TemporalMap(3,1) Y;
			TemporalMap(3,1) X;
			Cluster(64, P);
			SpatialMap(1,1) C;
		}
	}

	Layer CONV4_2 {
		Type: CONV
		Dimensions { K: 512, C: 512, R: 3, S: 3, Y: 66, X: 66 }
		Dataflow {
			SpatialMap(1,1) K;
			TemporalMap(64,64) C;
			TemporalMap(3,3) R;
			TemporalMap(3,3) S;
			TemporalMap(3,1) Y;
			TemporalMap(3,1) X;
			Cluster(64, P);
			SpatialMap(1,1) C;
		}
	}

	Layer CONV5_1 {
		Type: CONV
		Dimensions { K: 1024, C: 512, R: 3, S: 3, Y: 32, X: 32 }
		Dataflow {
			SpatialMap(1,1) K;
			TemporalMap(64,64) C;
			TemporalMap(3,3) R;
			TemporalMap(3,3) S;
			TemporalMap(3,1) Y;
			TemporalMap(3,1) X;
			Cluster(64, P);
			SpatialMap(1,1) C;
		}
	}

	Layer CONV5_2 {
		Type: CONV
		Dimensions { K: 1024, C: 1024, R: 3, S: 3, Y: 30, X: 30 }
		Dataflow {
			SpatialMap(1,1) K;
			TemporalMap(64,64) C;
			TemporalMap(3,3) R;
			TemporalMap(3,3) S;
			TemporalMap(3,1) Y;
			TemporalMap(3,1) X;
			Cluster(64, P);
			SpatialMap(1,1) C;
		}
	}

	Layer TRCONV1 {
		Type: TRCONV
		Dimensions { K: 512, C: 1024, R: 2, S: 2, Y: 28, X: 28 }
		Dataflow {
			SpatialMap(1,1) K;
			TemporalMap(64,64) C;
			TemporalMap(2,2) R;
			TemporalMap(2,2) S;
			TemporalMap(2,1) Y;
			TemporalMap(2,1) X;
			Cluster(64, P);
			SpatialMap(1,1) C;
		}
	}

	Layer CONV4_3 {
		Type: CONV
		Dimensions { K: 512, C: 1024, R: 3, S: 3, Y: 56, X: 56 }
		Dataflow {
			SpatialMap(1,1) K;
			TemporalMap(64,64) C;
			TemporalMap(3,3) R;
			TemporalMap(3,3) S;
			TemporalMap(3,1) Y;
			TemporalMap(3,1) X;
			Cluster(64, P);
			SpatialMap(1,1) C;
		}
	}

	Layer CONV4_4 {
		Type: CONV
		Dimensions { K: 512, C: 512, R: 3, S: 3, Y: 54, X: 54 }
		Dataflow {
			SpatialMap(1,1) K;
			TemporalMap(64,64) C;
			TemporalMap(3,3) R;
			TemporalMap(3,3) S;
			TemporalMap(3,1) Y;
			TemporalMap(3,1) X;
			Cluster(64, P);
			SpatialMap(1,1) C;
		}
	}

	Layer TRCONV2 {
		Type: TRCONV
		Dimensions { K: 512, C: 512, R: 2, S: 2, Y: 52, X: 52 }
		Dataflow {
			SpatialMap(1,1) K;
			TemporalMap(64,64) C;
			TemporalMap(2,2) R;
			TemporalMap(2,2) S;
			TemporalMap(2,1) Y;
			TemporalMap(2,1) X;
			Cluster(64, P);
			SpatialMap(1,1) C;
		}
	}

	Layer CONV3_3 {
		Type: CONV
		Dimensions { K: 256, C: 512, R: 3, S: 3, Y: 104, X: 104 }
		Dataflow {
			SpatialMap(1,1) K;
			TemporalMap(64,64) C;
			TemporalMap(3,3) R;
			TemporalMap(3,3) S;
			TemporalMap(3,1) Y;
			TemporalMap(3,1) X;
			Cluster(64, P);
			SpatialMap(1,1) C;
		}
	}

	Layer CONV3_4 {
		Type: CONV
		Dimensions { K: 256, C: 256, R: 3, S: 3, Y: 102, X: 102 }
		Dataflow {
			SpatialMap(1,1) K;
			TemporalMap(64,64) C;
			TemporalMap(3,3) R;
			TemporalMap(3,3) S;
			TemporalMap(3,1) Y;
			TemporalMap(3,1) X;
			Cluster(64, P);
			SpatialMap(1,1) C;
		}
	}

	Layer TRCONV3 {
		Type: TRCONV
		Dimensions { K: 128, C: 256, R: 2, S: 2, Y: 100, X: 100 }
		Dataflow {
			SpatialMap(1,1) K;
			TemporalMap(64,64) C;
			TemporalMap(2,2) R;
			TemporalMap(2,2) S;
			TemporalMap(2,1) Y;
			TemporalMap(2,1) X;
			Cluster(64, P);
			SpatialMap(1,1) C;
		}
	}

	Layer CONV2_3 {
		Type: CONV
		Dimensions { K: 128, C: 256, R: 3, S: 3, Y: 200, X: 200 }
		Dataflow {
			SpatialMap(1,1) K;
			TemporalMap(64,64) C;
			TemporalMap(3,3) R;
			TemporalMap(3,3) S;
			TemporalMap(3,1) Y;
			TemporalMap(3,1) X;
			Cluster(64, P);
			SpatialMap(1,1) C;
		}
	}

	Layer CONV2_4 {
		Type: CONV
		Dimensions { K: 128, C: 128, R: 3, S: 3, Y: 198, X: 198 }
		Dataflow {
			SpatialMap(1,1) K;
			TemporalMap(64,64) C;
			TemporalMap(3,3) R;
			TemporalMap(3,3) S;
			TemporalMap(3,1) Y;
			TemporalMap(3,1) X;
			Cluster(64, P);
			SpatialMap(1,1) C;
		}
	}

	Layer TRCONV4 {
		Type: TRCONV
		Dimensions { K: 64, C: 128, R: 2, S: 2, Y: 196, X: 196 }
		Dataflow {
			SpatialMap(1,1) K;
			TemporalMap(64,64) C;
			TemporalMap(2,2) R;
			TemporalMap(2,2) S;
			TemporalMap(2,1) Y;
			TemporalMap(2,1) X;
			Cluster(64, P);
			SpatialMap(1,1) C;
		}
	}

	Layer CONV1_3 {
		Type: CONV
		Dimensions { K: 64, C: 128, R: 3, S: 3, Y: 392, X: 392 }
		Dataflow {
			SpatialMap(1,1) K;
			TemporalMap(64,64) C;
			TemporalMap(3,3) R;
			TemporalMap(3,3) S;
			TemporalMap(3,1) Y;
			TemporalMap(3,1) X;
			Cluster(64, P);
			SpatialMap(1,1) C;
		}
	}

	Layer CONV1_4 {
		Type: CONV
		Dimensions { K: 64, C: 64, R: 3, S: 3, Y: 390, X: 390 }
		Dataflow {
			SpatialMap(1,1) K;
			TemporalMap(64,64) C;
			TemporalMap(3,3) R;
			TemporalMap(3,3) S;
			TemporalMap(3,1) Y;
			TemporalMap(3,1) X;
			Cluster(64, P);
			SpatialMap(1,1) C;
		}
	}

	Layer CONV1_5 {
		Type: CONV
		Dimensions { K: 2, C: 64, R: 1, S: 1, Y: 388, X: 388 }
		Dataflow {
			SpatialMap(1,1) K;
			TemporalMap(64,64) C;
			TemporalMap(1,1) R;
			TemporalMap(1,1) S;
			TemporalMap(1,1) Y;
			TemporalMap(1,1) X;
			Cluster(64, P);
			SpatialMap(1,1) C;
		}
	}
}

// Accelerator {
//  PE { NumPEs: 128; VectorWidth: 4; MultPrecision: INT8, AddPrecision: INT16 }
//  Buffer { GlobalL2: 2048, LocalL1: 64}
//  NoC {Bandwidth: 64; AvgLatency: 2}
// }

author: maestro-project

repo: maestro


text:
function [accuracy,threshold] = Sys_accuracy(match, nonmatch)
    [TAR,FAR,nonmatch_score] = rocplot(match, nonmatch);
    %[length(TAR),length(FAR),length(match),length(nonmatch_score),length(nonmatch)]
    accuracies=(TAR*length(match)+(1-FAR)*length(nonmatch))/(length(nonmatch)+length(match));
    [~,idx]=sort(abs(TAR-(1-FAR)),1,'ascend');

    [accuracy,a_idx]=max(accuracies(idx(1:100)));
    threshold=nonmatch_score(idx(a_idx(1)));

author: happynear

repo: FaceVerification


text:
function [kWh] = eV2kWh(eV)
% Convert energy or work from electron volts to kilowatt-hours.
% Chad A. Greene 2012
kWh = eV*4.4504925e-26;

author: Sable

repo: mcbench-benchmarks


text:
function res = conditions(this, varargin)
% Method for getting condition labels, over trials
% FORMAT res = conditions(this, ind, conditionlabels)
% _______________________________________________________________________
% Copyright (C) 2008-2012 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak
% $Id: conditions.m 5025 2012-10-31 14:44:13Z vladimir $

res = getset(this, 'trials', 'label', varargin{:});

if nargin == 1 && ~iscell(res)
    res = {res};
end


text:
classdef StressMinimizationWithGivenTopology < handle

    properties (Access = public)

    end

    properties (Access = private)
        mesh
        dataRes
        designVariable
    end

    properties (Access = private)
       microCase
       path
       fileResName
       initialIter
       vademecumName
       topOptProblem
    end

    methods (Access = public)

        function obj = StressMinimizationWithGivenTopology()
            obj.init();
            obj.wrapResAndMesh();
            obj.createTopOpt();
            obj.createDesignVariable();
            obj.solveTopOpt();
        end

    end

    methods (Access = private)

        function init(obj)
            obj.microCase = 'Rectangle';
            obj.vademecumName = 'SuperEllipseQMax';
            meshType  = 'Medium';
            fCase = [obj.microCase,meshType];
            obj.path = ['/media/alex/My Passport/LatticeResults/StressNorm',fCase,'/'];
            obj.fileResName  = 'ExperimentingPlot';
            obj.initialIter = 3;
        end

        function createTopOpt(obj)
            s = SettingsTopOptProblem('ExperimentingPlotRectangle');
            s.designVarSettings.creatorSettings.m1 = obj.dataRes.DesignVar1;
            s.designVarSettings.creatorSettings.m2 = obj.dataRes.DesignVar2;
            s.designVarSettings.creatorSettings.alpha0 =  obj.dataRes.AlphaGauss';
            obj.topOptProblem = TopOpt_Problem(s);
        end

        function solveTopOpt(obj)
            obj.topOptProblem.computeVariables();
            obj.topOptProblem.postProcess();
        end

        function wrapResAndMesh(obj)
            s.folderPath = obj.path;
            s.fileName   = [obj.fileResName,num2str(obj.initialIter)];
            w = WrapperMshResFiles(s);
            w.compute();
            obj.mesh    = w.mesh;
            obj.dataRes = w.dataRes;
        end

        function createDesignVariable(obj)
            s = SettingsDesignVariable();
            s.mesh                  = Mesh_Total(obj.mesh);
            s.type                  = 'MicroParams';
            s.initialCase           = 'full';
            s.scalarProductSettings = obj.createScalarProductParams();
            s.creatorSettings       = obj.createCreatorSettings();
            desVar = DesignVariable.create(s);
            desVar.alpha = obj.dataRes.AlphaGauss';
            obj.designVariable = desVar;
        end

        function s = createScalarProductParams(obj)
            s.epsilon = [];
            s.mesh = obj.mesh;
        end

        function s = createCreatorSettings(obj)
            s.m1 = obj.dataRes.DesignVar1;
            s.m2 = obj.dataRes.DesignVar2;
            s.homogSettings.fileName = obj.vademecumName;
            s.homogSettings.type     = 'ByVademecum';
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [Btuph] = TW2Btuph(TW)
% Convert power from terawatts to British thermal units per hour.
% Chad A. Greene 2012
Btuph = TW*3.412141633e+12;


function minus_log_dnsty = blp_opt_hyperpara(hyperpara,hh,prior,olsreg_,F,G,Fo,positions_nylags,position_constant)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 'blp_opt_hyperpara' computes the marginal likelihood over the
% hyperparameters of for the LP

% Inputs:
% - hyperpara, prior shrinkage
% likelihood

% Output: marginal data density

% Filippo Ferroni, 3/21/2020
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

hyperpara       = exp(hyperpara);
log_dnsty       = blp_ml(hyperpara,hh,prior,olsreg_,F,G,Fo,positions_nylags,position_constant);
minus_log_dnsty = -log_dnsty;


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% chvtr.m %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function chvtr(f,vtr)
% checks whether a required value to reach has already been reached; in
% that case flag is set to 0, otherwise it is set to 1
% Input:
% f	function value to be checked
% vtr	value to reach
% Output:
% flag	= 0  vtr has been reached
%	= 1  otherwise

function flag = chvtr(f,vtr)
if f <= vtr
  flag = 0;
else
  flag = 1;
end

author: lacerbi

repo: optimviz

save_path: github-repos/MATLAB/lacerbi-optimviz


text:
function plotProjectedAmps(view, scanList,projectionPhase)
%
% plotProjectedAmplitudes(view, [scanList])
%
% Bar plot of the amplitudes for each scan, averaging across
% all pixels (in all slices) in the current ROI.  The bar heights
% and a coarse SEM can be obtained from get(gca,'UserData').
%
% gmb  5/25/98
% bw   2/19/99  Added seY field to the UserData field.
%	    seY is an estimate of the variability in the
%      amplitudes.  It is the SEM of the in the complex
%      (amp*exp(-i*ph)) representation.  The values are
%      computed in vectorMean.m
% rfd  7/18/02 Added optional scanList.

if(~exist('scanList','var'))
    scanList = [1:numScans(view)];
end
if (~exist('projectionPhase','var'))
    computeProjPhase=1;
else
    computeProjPhase=0;
end
% Compute means across scans, for all pixels in the
% currently selected ROI.  The seZ value is the mean
% distance from the mean.
[meanAmps,meanPhs,seZ] = vectorMeans(view);

%Reference scan is the current scan
refScan = getCurScan(view);

% Compute the amplitude projected onto the reference phase
if (computeProjPhase)
    projectionPhase=meanPhs(refScan);
end
meanProjectedAmps = meanAmps.*cos(meanPhs-projectionPhase);

selectGraphWin;

% Header
ROIname = view.ROIs(view.selectedROI).name;
headerStr = ['Mean of projected amplitudes, ROI ',ROIname];
set(gcf,'Name',headerStr);

%plot the bar graph
clf
fontSize = 14;
h = mybar(meanProjectedAmps(scanList), zeros(size(scanList)), num2str(scanList'));
xlabel('Scan','FontSize',fontSize);
ylabel('Mean Projected Amplitude','FontSize',fontSize);
ylim =get(gca,'YLim');
set(gca,'YLim',ylim*1.1);
set(gca,'FontSize',fontSize);

foo = cell2struct(h,'bar');
hbar = foo(find(scanList==refScan)).bar;
set(hbar,'FaceColor','r')

%Save the data in gca('UserData')
data.y = meanProjectedAmps;
data.refScan = refScan;
data.seY = seZ;
set(gca,'UserData',data);

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% Copyright 2011 Zdenek Kalal
%
% This file is part of TLD.
%
% TLD is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% TLD is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with TLD.  If not, see <http://www.gnu.org/licenses/>.

function [BB1, shift] = bb_predict(BB0,pt0,pt1)

of  = pt1 - pt0;
dx  = median(of(1,:));
dy  = median(of(2,:));

d1  = pdist(pt0','euclidean');
d2  = pdist(pt1','euclidean');
s   = median(d2./d1);

s1  = 0.5*(s-1)*bb_width(BB0);
s2  = 0.5*(s-1)*bb_height(BB0);

BB1  = [BB0(1)-s1; BB0(2)-s2; BB0(3)+s1; BB0(4)+s2] + [dx; dy; dx; dy];
shift = [s1; s2];

author: yuxng


text:
function [d,fp,dt,tc,t]=v_readhtk(file)
%V_READHTK  read an HTK parameter file [D,FP,DT,TC,T]=(FILE)
%
% Input:
%    FILE = name of HTX file
% Outputs:
%       D = data: column vector for waveforms, one row per frame for other types
%      FP = frame period in seconds
%      DT = data type (also includes Voicebox code for generating data)
%             0  WAVEFORM     Acoustic waveform
%             1  LPC          Linear prediction coefficients
%             2  LPREFC       LPC Reflection coefficients:  -v_lpcar2rf([1 LPC]);LPREFC(1)=[];
%             3  LPCEPSTRA    LPC Cepstral coefficients
%             4  LPDELCEP     LPC cepstral+delta coefficients (obsolete)
%             5  IREFC        LPC Reflection coefficients (16 bit fixed point)
%             6  MFCC         Mel frequency cepstral coefficients
%             7  FBANK        Log Fliter bank energies
%             8  MELSPEC      linear Mel-scaled spectrum
%             9  USER         User defined features
%            10  DISCRETE     Vector quantised codebook
%            11  PLP          Perceptual Linear prediction
%            12  ANON
%      TC = full type code = DT plus (optionally) one or more of the following modifiers
%               64  _E  Includes energy terms
%              128  _N  Suppress absolute energy
%              256  _D  Include delta coefs
%              512  _A  Include acceleration coefs
%             1024  _C  Compressed
%             2048  _Z  Zero mean static coefs
%             4096  _K  CRC checksum (not implemented yet)
%             8192  _0  Include 0'th cepstral coef
%            16384  _V  Attach VQ index
%            32768  _T  Attach delta-delta-delta index
%       T = text version of type code e.g. LPC_C_K

%   Thanks to Dan Ellis (ee.columbia.edu) for sorting out decompression.
%   Thanks to Stuart Anderson (whispersys.com) for making it work on 64 bit machines.

%      Copyright (C) Mike Brookes 2005
%      Version: $Id: v_readhtk.m 10865 2018-09-21 17:22:45Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fid=fopen(file,'r','b');
if fid < 0
    error(sprintf('Cannot read from file %s',file));
end
nf=fread(fid,1,'int32');             % number of frames
fp=fread(fid,1,'int32')*1.E-7;       % frame interval (converted to seconds)
by=fread(fid,1,'int16');            % bytes per frame
tc=fread(fid,1,'int16');            % type code (see comments above for interpretation)
tc=tc+65536*(tc<0);
cc='ENDACZK0VT';                    % list of suffix codes
nhb=length(cc);                     % number of suffix codes
ndt=6;                              % number of bits for base type
hb=floor(tc*pow2(-(ndt+nhb):-ndt));
hd=hb(nhb+1:-1:2)-2*hb(nhb:-1:1);   % extract bits from type code
dt=tc-pow2(hb(end),ndt);            % low six bits of tc represent data type

% hd(7)=1 CRC check
% hd(5)=1 compressed data
if (dt==5)  % hack to fix error in IREFC files which are sometimes stored as compressed LPREFC
    fseek(fid,0,'eof');
    flen=ftell(fid);        % find length of file
    fseek(fid,12,'bof');
    if flen>14+by*nf        % if file is too long (including possible CRCC) then assume compression constants exist
        dt=2;               % change type to LPREFC
        hd(5)=1;            % set compressed flag
        nf=nf+4;            % frame count doesn't include compression constants in this case
    end
end

if any(dt==[0,5,10])        % 16 bit data for waveforms, IREFC and DISCRETE
    d=fread(fid,[by/2,nf],'int16').';
    if ( dt == 5),
        d=d/32767;                    % scale IREFC
    end
else
    if hd(5)                            % compressed data - first read scales
        nf = nf - 4;                    % frame count includes compression constants
        ncol = by / 2;
        scales = fread(fid, ncol, 'float');
        biases = fread(fid, ncol, 'float');
        d = ((fread(fid,[ncol, nf], 'int16')+repmat(biases,1,nf)).*repmat(1./scales,1,nf)).';
    else                              % uncompressed data
        d=fread(fid,[by/4,nf],'float').';
    end
end;
fclose(fid);
if nargout > 4
    ns=sum(hd);                 % number of suffixes
    kinds={'WAVEFORM' 'LPC' 'LPREFC' 'LPCEPSTRA' 'LPDELCEP' 'IREFC' 'MFCC' 'FBANK' 'MELSPEC' 'USER' 'DISCRETE' 'PLP' 'ANON' '???'};
    t=[kinds{min(dt+1,length(kinds))} reshape(['_'*ones(1,ns);cc(hd>0)],1,2*ns)];
end

author: ImperialCollegeLondon

repo: sap-voicebox