text:
function talCoords = mrAnatAcpc2Tal(talScale, acpcCoords)
%
% talCoords = mrAnatAcpc2Tal(talScale, acpcCoords)
%
% Converts acpc coords to Talairach coords.
%
% HISTORY:
% 2004.10.11 RFD: pulled code out of dtiFiberUI.
% 2005.04.29 RFD: added pcReference field to make this more general.

if(size(acpcCoords,1)~=3) acpcCoords = acpcCoords'; end
if(size(acpcCoords,1)~=3)
    error('The talCoords are the wrong size!');
end
talCoords = acpcCoords;
if(~isfield(talScale,'acpc') || isempty(talScale.acpc))
    return;
end
if(~isfield(talScale,'pcReference') || isempty(talScale.pcReference))
    tal = mrAnatGetTalairachDists;
    talScale.pcReference = tal.acpc;
end

% Do the easy ones first
% LEFT-RIGHT
inds = acpcCoords(1,:)>0;
talCoords(1,inds) = talCoords(1, inds).*talScale.rac;
inds = ~inds;
talCoords(1,inds) = talCoords(1, inds).*talScale.lac;
% SUP-INF
inds = acpcCoords(3,:)>0;
talCoords(3,inds) = talCoords(3,inds).*talScale.sac;
inds = ~inds;
talCoords(3,inds) = talCoords(3,inds).*talScale.iac;

pc = talScale.pcReference/talScale.acpc;
inds = acpcCoords(2,:)>0;
talCoords(2,inds) = talCoords(2,inds).*talScale.aac;
% between ac-pc
inds = acpcCoords(2,:)>=pc & acpcCoords(2,:)<0;
talCoords(2,inds) = talCoords(2,inds).*talScale.acpc;
inds = acpcCoords(2,:)<pc;
talCoords(2,inds) = talScale.pcReference + (talCoords(2,inds)-pc).*talScale.ppc;
return;


text:
function placetitlebar(fig)
%PLACETITLEBAR ensures that a figure's titlebar is on screen.
%
%   Copyright 1993-2004 The MathWorks, Inc.
%   $Revision: 1.1.4.1 $  $Date: 2004/01/24 09:36:24 $

oldRootUnits = get(0, 'Units');
oldFigUnits = get(fig, 'Units');

set(0, 'Units', 'pixels');
set(fig, 'Units', 'pixels');

screenpos = get(0, 'Screensize');
outerpos = get(fig, 'Outerposition');
if outerpos(2) + outerpos(4) > screenpos(4)
    outerpos(2) = screenpos(4) - outerpos(4);
    set(fig, 'Outerposition', outerpos);
end
%restore units
set(0, 'Units', oldRootUnits);
set(fig, 'Units', oldFigUnits);

author: zouchuhang

repo: LayoutNet


text:
function [ diff_median_mask]= diff_median_blk( In, win, thres)

p = 2*win + 1;

median_val = medfilt2(In, [p p]);

diff_median_mat = abs( median_val - In);

diff_median_mask= diff_median_mat >= ( thres* max( diff_median_mat(:)));

author: xingchenzhang

repo: VIFB

save_path: github-repos/MATLAB/xingchenzhang-VIFB


text:
function value = normal_ms_moment ( order, mu, sigma )

%*****************************************************************************80
%
%% NORMAL_MS_MOMENT evaluates moments of the Normal PDF.
%
%  Discussion:
%
%    The formula was posted by John D Cook.
%
%    Order  Moment
%    -----  ------
%      0    1
%      1    mu
%      2    mu^2 +         sigma^2
%      3    mu^3 +  3 mu   sigma^2
%      4    mu^4 +  6 mu^2 sigma^2 +   3      sigma^4
%      5    mu^5 + 10 mu^3 sigma^2 +  15 mu   sigma^4
%      6    mu^6 + 15 mu^4 sigma^2 +  45 mu^2 sigma^4 +  15      sigma^6
%      7    mu^7 + 21 mu^5 sigma^2 + 105 mu^3 sigma^4 + 105 mu   sigma^6
%      8    mu^8 + 28 mu^6 sigma^2 + 210 mu^4 sigma^4 + 420 mu^2 sigma^6 + 105 sigma^8
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 August 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer ORDER, the order of the moment.
%    0 <= ORDER.
%
%    Input, real MU, the mean of the distribution.
%
%    Input, real SIGMA, the standard deviation of the distribution.
%
%    Output, real VALUE, the value of the central moment.
%
  j_hi = floor ( order / 2 );

  value = 0.0;
  for j = 0 : j_hi
    value = value ...
      + r8_choose ( order, 2 * j ) ...
      * r8_factorial2 ( 2 * j - 1 ) ...
      * mu ^ ( order - 2 * j ) * sigma ^ ( 2 * j );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function removeNSxData(timestamps, filename)
%
%   removeNSxData(timestamps, filename)
%
%   This function removes chunks of data given in a 2-column timestamps
%   variable from a NSx file and then saves it as another NSx. This is
%   useful when removing useless data from the file to decrease the size of
%   the NSx file without modifying the file timestamps. The file size only
%   changes when the NSx file is compressed.
%
%   timestamps: This variable contains beginning and ending timestamps of
%               the data that need to be removed from the NSx file. The
%               format is:
%                          [BegTimeStamp1 EndTimeStamp1
%                           BegTimeStamp2 EndTimeStamp2
%                           BegTimeStamp3 EndTimeStamp3
%                                ...          ...
%                           BegTimeStampN EndTimeStampN]
%
%               The data in the NSx between the pairs of timestamps will be
%               removed and set to 0.
%               (REQUIRED)
%
%   filename:   The name of the NSx file to be used. This input is also
%               optional. In its absense, a dialog will open and will
%               prompt the user to select an NSx file.
%               (OPTIONAL)
%
%   Example:    removeNSxData([100,600; 1000,2000], 'c:\datafile\sampleNSx.ns5')
%
%               In the example above, the timestamps between 100:600 and
%               1000:2000 will be removed from the file sampleNSx.ns5 and
%               the resulting NSx will be saved in sampleNSx_chunked.ns5
%               file.
%
%   Kian Torab
%   ktorab@blackrockmicro.com
%   Blackrock Microsystems
%
%   Version 1.0.1.0

% Validating the timestamps variable
if ~exist('timestamps', 'var')
    disp('The variable timestamps is required.');
    return;
end

if size(timestamps, 2) ~= 2
    disp('The timestamps variable should be a Nx2 variable containing beginning timestamps in the first column and ending timestamps in the second column.');
    return;
end

if ~all(timestamps(:,1) < timestamps(:,2))
    disp('The timestamps in the first column (beginning) have to be smaller than the ones in the second column (ending).')
    return;
end

% Figuring out the filename, if not passed on to the function
if ~exist('filename', 'var')
	[fname, path] = getFile('*.*', 'Choose an NSx file...');
    filename = [path fname];
end

if exist(filename, 'file') ~= 2
    disp('The file does not exist.');
    return;
end

% Openning and reading the file
FID = fopen(filename, 'r', 'ieee-le');
readData = fread(FID, '*uint8');
FID = fclose(FID);

% Extracting some information from the NSx file
fileTypeID   = char(readData(1:8))';
headerBytes  = typecast(readData(11:14), 'uint32')+9;
samplingFreq = double(typecast(readData(291:294), 'uint32'))/...
               double(typecast(readData(287:290), 'uint32'));
channelCount = typecast(readData(311:314), 'uint32');

% Converting timestamps into uint16 (2 bytes per data point)
timestamps = timestamps*2;

% Removing given timestamps
for idx = 1:size(timestamps, 1)
    readData((timestamps(idx,1)-2)*channelCount+1+headerBytes:timestamps(idx,2)*channelCount+headerBytes) = NaN;
end

% Saving the chunked file
writeFilename = [filename(1:end-4) '_chunked', filename(end-3:end)];
FID = fopen(writeFilename, 'w+', 'ieee-le');
fwrite(FID, readData);
fclose(FID);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function rez = fullMPMU(rez, DATA)

ops = rez.ops;

Nfilt   = ops.Nfilt;
lam =  ones(Nfilt, 1, 'single');
lam(:)    = ops.lam(3);

[W, U, mu, UtU, nu] = decompose_dWU(ops, rez.dWU, ops.Nrank,rez.ops.kcoords);

pm = exp(-ops.momentum(2));
Params = double([ops.NT ops.Nfilt ops.Th(3) ops.maxFR 10 ops.Nchan ops.Nrank pm ops.epu ops.nt0]);

Params(3) = ops.Th(3);
Params(4) = 50000;
Params(5) = 50;

if ops.GPU
    U0 = gpuArray(U);
else
    U0 = U;
end
%%
nt0     = rez.ops.nt0;
Nrank   = ops.Nrank;
WtW     = zeros(Nfilt,Nfilt,2*nt0-1, 'single');
for i = 1:Nrank
    for j = 1:Nrank
        utu0 = U0(:,:,i)' * U0(:,:,j);
        if ops.GPU
            wtw0 =  gather_try(mexWtW2(Params, W(:,:,i), W(:,:,j), utu0));
        else
            wtw0 =  getWtW2(Params, W(:,:,i), W(:,:,j), utu0);
            wtw0 = permute(wtw0, [2 3 1]);
        end
        WtW = WtW + wtw0;
        clear wtw0 utu0
%         wtw0 = squeeze(wtw(:,i,:,j,:));

    end
end

mWtW = max(WtW, [], 3);
WtW = permute(WtW, [3 1 2]);

if ops.GPU
    WtW = gpuArray(WtW);
end
%%
Nbatch_buff = rez.temp.Nbatch_buff;
Nbatch      = rez.temp.Nbatch;
Nchan       = ops.Nchan;
if ~ops.GPU
   fW = rez.fW; % load fft-ed templates
end
% mWtW = mWtW - diag(diag(mWtW));

% rez.WtW = gather_try(WtW);
%
clear wtw0 utu0 U0
%
clear nspikes2
st3 = [];
rez.st3 = [];

if ops.verbose
   fprintf('Time %3.0fs. Running the final template matching pass...\n', toc)
end

if Nbatch_buff<Nbatch
    fid = fopen(ops.fproc, 'r');
end
msg = [];

if ~isempty(ops.nNeigh)
    nNeigh    = ops.nNeigh;

    rez.cProj = zeros(5e6, nNeigh, 'single');

    % sort pairwise templates
    nsp = sum(rez.nspikes,2);
    vld = single(nsp>100);
    cr    = mWtW .* (vld * vld');
    cr(isnan(cr)) = 0;
    [~, iNgsort] = sort(cr, 1, 'descend');

    % save full similarity score
    rez.simScore = cr;
    maskTT = zeros(Nfilt, 'single');
    rez.iNeigh = iNgsort(1:nNeigh, :);
    for i = 1:Nfilt
        maskTT(rez.iNeigh(:,i),i) = 1;
    end
end
if ~isempty(ops.nNeighPC)
    nNeighPC  = ops.nNeighPC;
    load PCspikes
    ixt = round(linspace(1, size(Wi,1), ops.nt0));
    Wi = Wi(ixt, 1:3);
    rez.cProjPC = zeros(5e6, 3*nNeighPC, 'single');

    % sort best channels
    [~, iNch]       = sort(abs(U(:,:,1)), 1, 'descend');
    maskPC          = zeros(Nchan, Nfilt, 'single');
    rez.iNeighPC    = iNch(1:nNeighPC, :);
    for i = 1:Nfilt
        maskPC(rez.iNeighPC(:,i),i) = 1;
    end
    maskPC = repmat(maskPC, 3, 1);
end

irun = 0;
i1nt0 = int32([1:nt0])';
%%
LAM = lam .* (20./mu).^2;

NT = ops.NT;
batchstart = 0:NT:NT*(Nbatch-Nbatch_buff);

for ibatch = 1:Nbatch
    if ibatch>Nbatch_buff
        offset = 2 * ops.Nchan*batchstart(ibatch-Nbatch_buff); % - ioffset;
        fseek(fid, offset, 'bof');
        dat = fread(fid, [NT ops.Nchan], '*int16');
    else
       dat = DATA(:,:,ibatch);
    end
    if ops.GPU
        dataRAW = gpuArray(dat);
    else
        dataRAW = dat;
    end
    dataRAW = single(dataRAW);
    dataRAW = dataRAW / ops.scaleproc;

    % project data in low-dim space
    if ops.GPU
        data    = gpuArray.zeros(NT, Nfilt, Nrank, 'single');
    else
        data   = zeros(NT, Nfilt, Nrank, 'single');
    end
    for irank = 1:Nrank
        data(:,:,irank) = dataRAW * U(:,:,irank);
    end
    data                = reshape(data, NT, Nfilt*Nrank);

    if ops.GPU
        [st, id, x, errC, PCproj] ...
                        = mexMPmuFEAT(Params,data,W,WtW, mu, lam .* (20./mu).^2, nu);
    else
         [st, id, x, errC, PCproj]= cpuMPmuFEAT(Params,data,fW,WtW, mu, lam .* (20./mu).^2, nu, ops);
    end

    if ~isempty(st)
        if ~isempty(ops.nNeighPC)
            % PCA coefficients
            inds            = repmat(st', nt0, 1) + repmat(i1nt0, 1, numel(st));
            try  datSp      = dataRAW(inds(:), :);
            catch
                datSp       = dataRAW(inds(:), :);
            end
            datSp           = reshape(datSp, [size(inds) Nchan]);
            coefs           = reshape(Wi' * reshape(datSp, nt0, []), size(Wi,2), numel(st), Nchan);
            coefs           = reshape(permute(coefs, [3 1 2]), [], numel(st));
            coefs           = coefs .* maskPC(:, id+1);
            iCoefs          = reshape(find(maskPC(:, id+1)>0), 3*nNeighPC, []);
            rez.cProjPC(irun + (1:numel(st)), :) = gather_try(coefs(iCoefs)');
        end
        if ~isempty(ops.nNeigh)
            % template coefficients
            % transform coefficients
            PCproj          = bsxfun(@rdivide, ...
                bsxfun(@plus, PCproj, LAM.*mu), sqrt(1+LAM));

            PCproj          = maskTT(:, id+1) .* PCproj;
            iPP             = reshape(find(maskTT(:, id+1)>0), nNeigh, []);
            rez.cProj(irun + (1:numel(st)), :) = PCproj(iPP)';
        end
        % increment number of spikes
        irun            = irun + numel(st);

        if ibatch==1;
            ioffset         = 0;
        else
            ioffset         = ops.ntbuff;
        end
        st                  = st - ioffset;

        %     nspikes2(1:size(W,2)+1, ibatch) = histc(id, 0:1:size(W,2));
        STT = cat(2, ops.nt0min + double(st) +(NT-ops.ntbuff)*(ibatch-1), ...
            double(id)+1, double(x), ibatch*ones(numel(x),1));
        st3             = cat(1, st3, STT);
    end
    if rem(ibatch,100)==1
%         nsort = sort(sum(nspikes2,2), 'descend');
        fprintf(repmat('\b', 1, numel(msg)));
        msg             = sprintf('Time %2.2f, batch %d/%d,  NTOT %d\n', ...
            toc, ibatch,Nbatch, size(st3,1));
        fprintf(msg);

    end
end
%%
[~, isort]      = sort(st3(:,1), 'ascend');
st3             = st3(isort,:);

rez.st3         = st3;
if ~isempty(ops.nNeighPC)
    % re-sort coefficients for projections
    rez.cProjPC(irun+1:end, :)  = [];
    rez.cProjPC                 = reshape(rez.cProjPC, size(rez.cProjPC,1), [], 3);
    rez.cProjPC                 = rez.cProjPC(isort, :,:);
    for ik = 1:Nfilt
        iSp                     = rez.st3(:,2)==ik;
        OneToN                  = 1:nNeighPC;
        [~, isortNeigh]         = sort(rez.iNeighPC(:,ik), 'ascend');
        OneToN(isortNeigh)      = OneToN;
        rez.cProjPC(iSp, :,:)   = rez.cProjPC(iSp, OneToN, :);
    end

    rez.cProjPC                 = permute(rez.cProjPC, [1 3 2]);
end
if ~isempty(ops.nNeigh)
    rez.cProj(irun+1:end, :)    = [];
    rez.cProj                   = rez.cProj(isort, :);

    % re-index the template coefficients
    for ik = 1:Nfilt
        iSp                     = rez.st3(:,2)==ik;
        OneToN                  = 1:nNeigh;
        [~, isortNeigh]         = sort(rez.iNeigh(:,ik), 'ascend');
        OneToN(isortNeigh)      = OneToN;
        rez.cProj(iSp, :)       = rez.cProj(iSp, OneToN);
    end
end

%%
% rez.ops             = ops;
rez.W               = W;
rez.U               = U;
rez.mu              = mu;

rez.t2p = [];
for i = 1:Nfilt
    wav0            = W(:,i,1);
    wav0            = my_conv(wav0', .5)';
   [~, itrough]     = min(wav0);
    [~, t2p]        = max(wav0(itrough:end));
    rez.t2p(i,1)    = t2p;
    rez.t2p(i,2)    = itrough;
end

rez.nbins           = histc(rez.st3(:,2), .5:1:Nfilt+1);

[~, rez.ypos]       = max(rez.U(:,:,1), [], 1);
if Nbatch_buff<Nbatch
    fclose(fid);
end

% center the templates
rez.W               = cat(1, zeros(nt0 - (ops.nt0-1-ops.nt0min), Nfilt, Nrank), rez.W);
rez.WrotInv         = (rez.Wrot/200)^-1;
%%
Urot = U;
for k = 1:size(U,3)
   Urot(:,:,k)  = rez.WrotInv' * Urot(:,:,k);
end
for n = 1:size(U,2)
    rez.Wraw(:,:,n) = mu(n) * sq(Urot(:,n,:)) * sq(rez.W(:,n,:))';
end
%

author: cortex-lab

repo: KiloSort

save_path: github-repos/MATLAB/cortex-lab-KiloSort


text:
% Optimal doping profile optimization with current gain constraint.
% Boyd, Kim, Vandenberghe, and Hassibi, "A tutorial on geometric programming"
% Joshi, Boyd, and Dutton, "Optimal doping profiles via geometric programming"
% Written for CVX by Almir Mutapcic 02/08/06
% (a figure is generated)
%
% Determines the optimal doping profile that minimizes base transit
% time subject to a lower bound constraint on the current gain (beta).
% This problem can be posed as a GP:
%
%   minimize   tau_B
%       s.t.   Nmin <= v <= Nmax
%              y_(i+1) + v_i^const1 <= y_i
%              w_(i+1) + v_i^const2 <= w_i, etc...
%              beta => beta_min
%
% where variables are v_i, y_i, and w_i.

% problem size
M = 20;

% problem constants
g1 = 0.42;
g2 = 0.69;
Nmax = 5*10^18;
Nmin = 5*10^16;
Nref = 10^17;
Dn0 = 20.72;
ni0 = 1.4*(10^10);
WB = 10^(-5);
C =  WB^2/((M^2)*(Nref^g1)*Dn0);

% minimum current gain values
beta_min_GE = [1 1.4 1.8 2.2 2.6 3.0 3.4 3.43]*(1e-11);

% exponent powers
pwi = g2 -1;
pwj = 1+g1-g2;

v_array = [];
for k = 1:length(beta_min_GE)
    fprintf( 'beta_min_GE = %g: ', beta_min_GE(k) );
    cvx_begin gp quiet
        % optimization variables
        variables v(M) y(M) w(M)

        % objective function is the base transmit time
        tau_B = C*w(1);

        minimize( tau_B )
        subject to
        % fixed problem constraints
        Nmin <= v <= Nmax;

        for i=1:M-1
            y(i+1) + v(i)^pwj <= y(i);
            w(i+1) + y(i)*v(i)^pwi <= w(i);
        end

        % equalities
        y(M) == v(M)^pwj;
        w(M) == y(M)*v(M)^pwi;

        % changing constraint
        (WB*beta_min_GE(k)/(M*Nref^(g1-g2)*Dn0))*y(1) <= 1;
    cvx_end
    fprintf( '%s\n', cvx_status );
    % keep the optimal solution
    v_array = [v_array v];
end

% plot the basic optimal doping profile
figure, clf
nbw = 0:1/M:1-1/M;
for k = 1:length(beta_min_GE)
  semilogy(nbw,v_array(:,k),'LineWidth',2); hold on;
end
axis([0 1 1e16 1e19]);
xlabel('base');
ylabel('doping');
text(0,Nmin,'Nmin ', 'HorizontalAlignment','right');
text(0,Nmax,'Nmax ', 'HorizontalAlignment','right');
hold off;

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function pred_label=SpatialJSRC(feat, train_data, train_label, sp_label, block_size, param)
% Classification of superpixels with joint sparse representation
% Input:
%    feat: extracted img features, nr*nc*nd where nd is the number of feature channel
%    dict: the dictionary
%    dict_label: the label of dictionary atoms
%    sp_label: the label of superpixels, nr*nc
%    block_size: the size for block-wise processing
%    param: parameter for sparse representation
%            param.L (optional, maximum number of elements in each decomposition,
%               min(m,p) by default)
%            param.eps (optional, threshold on the squared l2-norm of the residual,
%               0 by default
%            param.lambda (optional, penalty parameter, 0 by default
%            param.numThreads (optional, number of threads for exploiting
%            multi-core / multi-cpus. By default, it takes the value -1,
%            which automatically selects all the available CPUs/cores).
%   2016-10-20, jlfeng

[nr,nc,~]=size(feat);
sp_label=reshape(sp_label,[nr*nc 1]);
[group_label,sortidx]=sort(sp_label);
feat_sort=VectorIndexing3D(feat,sortidx);
idx_group_start=zeros(length(unique(group_label)),1);
idx_group_start(2:end)=int32(find(diff(group_label)~=0));

num_sp=length(idx_group_start);
num_block=ceil(num_sp/block_size);
disp('Block-wise classification strategy is used. ');
disp(['Number of Blocks: ', num2str(num_block)]);
pause(0.05)
pred_label_sort=zeros(nr*nc,1);
for nn=1:num_block
    disp(['Processing block ',num2str(nn)]);tic
    idxstart=(nn-1)*block_size;
    idxend=min(nn*block_size-1,num_sp-1);
    idx_block=(group_label>=idxstart)&(group_label<=idxend);
    test_data=feat_sort(:,idx_block);
    idx_group_start_block=idx_group_start(idxstart+1:idxend+1);
    idx_group_start_block=idx_group_start_block-idx_group_start_block(1);
    pred_label_block=JSRClassifier(train_data,train_label,test_data,idx_group_start_block,param);
    pred_label_sort(idx_block)=pred_label_block;
end
pred_label=zeros(nr*nc,1);
pred_label(sortidx)=pred_label_sort;
pred_label=reshape(pred_label,[Nx Ny]);

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function compareMapsByFields(map1, fields1, map2, fields2)
    % create features based on field properties, compare featrues

    numFields1 = length(fields1);
    numFields2 = length(fields2);
    distances = zeros(numFields2, 1);

    matchedPoints1 = zeros(numFields1, 2);
    matchedPoints2 = zeros(numFields1, 2);
    toRemove = false(numFields1, 1); % above matched data will be truncated based on toRemove variable

    matchedIndices = zeros(numFields1, 1);

    for i = 1:numFields1
        matchedPoints1(i, :) = [fields1(i).peakX fields1(i).peakY];
        features1 = fieldFeatures(fields1(i));

        for j = 1:numFields2
            features2 = fieldFeatures(fields2(j));
            distances(j) = hausdorffDistance(features1, features2);
%             peakPoint = [fields2(j).peakX fields2(j).peakY];
%             distances(j) = pdist2(matchedPoints1(i, :), peakPoint);
        end

        [minDist, minInd] = min(distances);
        if minDist < 21 && ~any(matchedIndices == minInd)
%         if ~any(matchedIndices == minInd)
            features2 = fieldFeatures(fields2(minInd));

            distance = pdist2(features1, features2);

            matchedPoints2(i, :) = [fields2(minInd).peakX fields2(minInd).peakY];
            matchedIndices(i) = minInd;
        else
            toRemove(i) = true;
        end
        fprintf('Field %u, distance %f, matched index %u\n', i, minDist, minInd);

    end
    matchedPoints1(toRemove, :) = [];
    matchedPoints2(toRemove, :) = [];

    if ~isempty(matchedPoints1)
        figure;
        plot.matchedPoints(map1, map2, matchedPoints1, matchedPoints2);
    else
        fprintf('No matches\n');
    end
end

function features = fieldFeatures(field)
%     features = zeros(5, 1);

    features(1) = field.x;
    features(2) = field.y;
%     features(3) = field.x;
%     features(4) = field.y;
%     features(3) = field.area;
%     features(1) = field.size;
    features(3) = field.Perimeter;
    features(4) = field.meanRate;
%     features(3) = field.Eccentricity;
%     features(4) = field.Extent;
%     features(5) = field.area;
    features(5) = field.Orientation;
end

author: kavli-ntnu


text:
function varargout = aibcutpush(varargin)
% VL_AIBCUTPUSH  Quantize based on VL_AIB cut
%  Y = VL_AIBCUTPUSH(MAP, X) maps the data X to elements of the AIB
%  cut specified by MAP.
%
%  The function is equivalent to Y = MAP(X).
%
%  See also: VL_HELP(), VL_AIB().
[varargout{1:nargout}] = vl_aibcutpush(varargin{:});

author: yihui-he

repo: panorama

save_path: github-repos/MATLAB/yihui-he-panorama


text:
function kern = sqexpKernExpandParam(kern, params)

% SQEXPKERNEXPANDPARAM Create kernel structure from SQEXP kernel's parameters.
% FORMAT
% DESC returns a pre-built compound squared exponential kernel structure filled with the
% parameters in the given vector. This is used as a helper function to
% enable parameters to be optimised in, for example, the NETLAB
% optimisation functions.
% ARG kern : the kernel structure in which the parameters are to be
% placed.
% ARG param : vector of parameters which are to be placed in the
% kernel structure.
% RETURN kern : kernel structure with the given parameters in the
% relevant locations.
%
% SEEALSO : sqexpKernParamInit, sqexpKernExtractParam, kernExpandParam
%
% COPYRIGHT : Neil D. Lawrence, 2004

% KERN

kern.inverseWidth = params(1);
kern.rbfVariance = params(2);
kern.biasVariance = params(3);
kern.whiteVariance = params(4);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function test_issue1410

% MEM 6gb
% WALLTIME 00:20:00
% DEPENDENCY

load(dccnpath('/home/common/matlab/fieldtrip/data/test/latest/raw/meg/preproc_ctf151.mat'), 'data');
grad = data.grad;
clear data

%%

vol = [];
vol.o = [0 0 4];
vol.r = 12;
vol.unit = 'cm';
vol.type = 'singlesphere';

%%

cfg = [];
cfg.dip.pos = [0 0 9];
cfg.dip.mom = [1 0 0];
cfg.dip.unit = 'cm';
cfg.headmodel = vol;
cfg.grad = grad;
data = ft_dipolesimulation(cfg);

cfg = [];
timelock = ft_timelockanalysis(cfg, data);

cfg = [];
cfg.method = 'mtmfft';
cfg.taper = 'hanning';
cfg.output = 'fourier';
freq = ft_freqanalysis(cfg, data);

%%

tfdata = timelock;
method = 'lcmv';

cfg = [];
cfg.xgrid = -8:2:8;
cfg.ygrid = -8:2:8;
cfg.zgrid = (-8:2:8) + 4;
cfg.unit = 'cm';
cfg.headmodel = vol;
cfg.method = method;
cfg.keepleadfield = 'yes';
cfg.(method).keepfilter = 'yes';
cfg.(method).lambda = '10%';
cfg.channel = 'MEG';
source1 = ft_sourceanalysis(cfg, tfdata);

%%

% in the next code the labels are missing for the precomputed filters
% this should assume that they were computed with the same channel selection

cfg = [];
cfg.method = method;
cfg.sourcemodel.pos = source1.pos;
cfg.sourcemodel.filter = source1.avg.filter;
cfg.channel = 'MEG';
source2 = ft_sourceanalysis(cfg, tfdata);

clear source1 source2

%%
% this is mostly from http://www.fieldtriptoolbox.org/tutorial/beamformingextended/

clear all
load(dccnpath('/home/common/matlab/fieldtrip/data/ftp/tutorial/sensor_analysis/subjectK.mat'));

data_combined = ft_appenddata([], data_left, data_right);

cfg           = [];
cfg.toilim    = [-0.8 1.1];
cfg.minlength = 'maxperlen'; % this ensures all resulting trials are equal length
data          = ft_redefinetrial(cfg, data_combined);

[ftver, ftdir] = ft_version;

%%
% the first analysis is for visual gamma

cfg        = [];
cfg.toilim = [-0.8 0];
data_bsl   = ft_redefinetrial(cfg, data);

cfg.toilim = [0.3 1.1];
data_exp   = ft_redefinetrial(cfg, data);

cfg      = [];
data_cmb = ft_appenddata(cfg, data_bsl, data_exp);

% give a number to each trial: 0 = baseline, 1 = experimental condition
data_cmb.trialinfo = [zeros(length(data_bsl.trial), 1); ones(length(data_exp.trial), 1)];

cfg            = [];
cfg.method     = 'mtmfft';
cfg.output     = 'fourier';
cfg.keeptrials = 'yes';
cfg.tapsmofrq  = 15;
cfg.foi        = 55;
freq_cmb       = ft_freqanalysis(cfg, data_cmb);

cfg                = [];
cfg.trials         = freq_cmb.trialinfo == 0;
freq_bsl           = ft_selectdata(cfg, freq_cmb);

cfg.trials         = freq_cmb.trialinfo == 1;
freq_exp           = ft_selectdata(cfg, freq_cmb);

%%
% this is used both for visual gamma and beta coherence

load(dccnpath('/home/common/matlab/fieldtrip/data/ftp/tutorial/beamformer_extended/hdm.mat'))
load(dccnpath('/home/common/matlab/fieldtrip/data/ftp/tutorial/beamformer_extended/sourcemodel.mat'))

cfg             = [];
cfg.sourcemodel = sourcemodel;
cfg.headmodel   = hdm;
cfg.channel     = {'MEG'};
cfg.grad        = data.grad;
cfg.singleshell.batchsize = 10;
sourcemodel_lf  = ft_prepare_leadfield(cfg);

templatedir = fullfile(ftdir, 'template', 'sourcemodel');
template = load(fullfile(templatedir, 'standard_sourcemodel3d8mm')); % 8mm spacing grid

templatedir = fullfile(ftdir, 'external', 'spm8', 'templates');
template_mri = ft_read_mri(fullfile(templatedir, 'T1.nii'));
template_mri.coordsys = 'mni'; % we know it's in MNI space

%%

cfg                   = [];
cfg.frequency         = freq_cmb.freq;
cfg.grad              = freq_cmb.grad;
cfg.method            = 'dics';
cfg.keeptrials        = 'yes';
cfg.channel           = 'MEG';
cfg.sourcemodel       = sourcemodel_lf;
cfg.keeptrials        = 'yes';
cfg.dics.lambda       = '5%';
cfg.dics.keepfilter   = 'yes';
cfg.dics.fixedori     = 'yes';
cfg.dics.realfilter   = 'yes';
source                = ft_sourceanalysis(cfg, freq_cmb);

% beam pre- and poststim by using the common filter
cfg.sourcemodel.filter  = source.avg.filter;
cfg.sourcemodel.label   = source.avg.label;
source_bsl       = ft_sourceanalysis(cfg, freq_bsl);
source_exp       = ft_sourceanalysis(cfg, freq_exp);

cfg = [];
cfg.parameter = 'avg.pow';
cfg.operation = '(x1 ./ x2) - 1';
source_diff = ft_math(cfg, source_exp, source_bsl);

source_diff.pos = template.sourcemodel.pos;
source_diff.dim = template.sourcemodel.dim;

cfg              = [];
cfg.parameter    = 'pow';
cfg.interpmethod = 'nearest';
source_diff_int  = ft_sourceinterpolate(cfg, source_diff, template_mri);

cfg               = [];
cfg.method        = 'slice';
cfg.funparameter  = 'pow';
cfg.maskparameter = cfg.funparameter;
cfg.funcolorlim   = [0.0 1.2];
cfg.opacitylim    = [0.0 1.2];
cfg.opacitymap    = 'rampup';
ft_sourceplot(cfg, source_diff_int);

clear source_bsl source_exp source_diff source_diff_int

%%
% now for the EMG coherence

cfg                 = [];
cfg.toilim          = [-1 -0.0025];
cfg.minlength       = 'maxperlen'; % this ensures all resulting trials are equal length
data_stim           = ft_redefinetrial(cfg, data);

cfg                 = [];
cfg.output          = 'powandcsd';
cfg.method          = 'mtmfft';
cfg.taper           = 'dpss';
cfg.tapsmofrq       = 5;
cfg.foi             = 20;
cfg.keeptrials      = 'yes';
cfg.channel         = {'MEG' 'EMGlft' 'EMGrgt'};
cfg.channelcmb      = {'MEG' 'MEG'; 'MEG' 'EMGlft'; 'MEG' 'EMGrgt'};
freq_csd            = ft_freqanalysis(cfg, data_stim);

%%

cfg                 = [];
cfg.method          = 'dics';
cfg.refchan         = 'EMGlft';
cfg.channel         = 'MEG';
cfg.frequency       = 20;
% with precomputed leadfields
cfg.sourcemodel     = sourcemodel_lf;
source_coh_lft      = ft_sourceanalysis(cfg, freq_csd);

% without precomputed leadfields
cfg.sourcemodel     = removefields(sourcemodel_lf, {'leadfield', 'leadfielddimord', 'label'});
cfg.headmodel       = hdm;
source_coh          = ft_sourceanalysis(cfg, freq_csd);

figure; hold on
plot(source_coh.avg.coh(:), '.')
plot(source_coh_lft.avg.coh(:), 'ro')

assert(isalmostequal(source_coh_lft.avg.coh, source_coh.avg.coh, 'reltol', 1e-6));

%%

cfg                  = [];
cfg.method           = 'dics';
cfg.refchan          = 'EMGlft';
cfg.channel          = 'MEG';
cfg.frequency        = 20;
cfg.headmodel        = hdm;
cfg.sourcemodel.pos  = [3.6987   -3.1192   11.2426];
cfg.sourcemodel.unit = 'cm';
cfg.reducerank       = 2;
cfg.keepleadfield    = 'yes';
cfg.dics.keepfilter  = 'yes';
cfg.backproject      = 'yes';
source3              = ft_sourceanalysis(cfg, freq_csd);
cfg.backproject      = 'no';
source2              = ft_sourceanalysis(cfg, freq_csd);

assert(isalmostequal(source3.avg.coh, source2.avg.coh, 'reltol', 1e-6));
assert(~isequal(source3.leadfield, source2.leadfield)); % should have different number of columns
assert(~isequal(source3.avg.filter, source2.avg.filter)); % should have different number of rows

clear source3 source2

%%

source_coh_lft.pos = template.sourcemodel.pos;
source_coh_lft.dim = template.sourcemodel.dim;

source_coh.pos = template.sourcemodel.pos;
source_coh.dim = template.sourcemodel.dim;

cfg              = [];
cfg.parameter    = 'coh';
cfg.interpmethod = 'nearest';
cfg.coordsys     = 'mni';
source_coh_int   = ft_sourceinterpolate(cfg, source_coh_lft, template_mri);

cfg               = [];
cfg.method        = 'ortho';
cfg.funparameter  = 'coh';
cfg.maskparameter = 'coh';
cfg.headmodel     = hdm;
ft_sourceplot(cfg, source_coh_int);

cfg               = [];
cfg.method        = 'surface';
cfg.funparameter  = 'coh'; % use it to represent color
cfg.maskparameter = 'coh'; % and use it for transparency
ft_sourceplot(cfg, source_coh_int);

%%
% let's also try dics_refdip

% these are approximately in the motor cortex
diplft = [4  3 10];
diprgt = [4 -3 10];

% take the nearest grid point
d = sqrt(sum((sourcemodel_lf.pos - repmat(diplft, size(sourcemodel_lf.pos,1), 1)).^2, 2));
[mind, indx] = min(d);
diplft = sourcemodel_lf.pos(indx,:);

cfg                 = [];
cfg.method          = 'dics';
cfg.refdip          = diplft;
cfg.channel         = 'MEG';
cfg.frequency       = 20;
% with precomputed leadfields
cfg.sourcemodel     = sourcemodel_lf; % the leadfields will be discarded
cfg.headmodel       = hdm;
source_coh_diplft   = ft_sourceanalysis(cfg, freq_csd);

cfg = [];
cfg.funparameter = 'coh'; % use it to represent color
ft_sourceplot(cfg, source_coh_diplft); % it should be 1 on the location of diplft


text:
classdef PTKDensityInterpolation < PTKPlugin
    % PTKDensityInterpolation. Plugin for interpolating density values to a
    % different voxel size
    %
    %     This is a plugin for the Pulmonary Toolkit. Plugins can be run using
    %     the gui, or through the interfaces provided by the Pulmonary Toolkit.
    %     See PTKPlugin.m for more information on how to run plugins.
    %
    %     Plugins should not be run directly from your code.
    %
    %
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2012.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    properties
        ButtonText = 'Density Interpolation'
        ToolTip = 'Recomputes image density with different sized voxels'
        Category = 'Lungs'

        AllowResultsToBeCached = true
        AlwaysRunPlugin = false
        PluginType = 'ReplaceOverlay'
        HidePluginInDisplay = false
        FlattenPreviewImage = false
        PTKVersion = '1'
        ButtonWidth = 6
        ButtonHeight = 1
        GeneratePreview = false
        Visibility = 'Developer'
    end

    methods (Static)
        function results = RunPlugin(application, reporting)

            % The size of the 'voxels' for calculating the density.
            % Note if you want to change this, you should change the above
            % property
            %     AlwaysRunPlugin = false
            % so that the plug is forced to re-run

            interp_voxel_size_mm = [5, 5, 5];

            % Fetch the intensity of just the lung regions
            roi = application.GetResult('PTKLungROI');
            left_and_right_lungs = application.GetResult('PTKLeftAndRightLungs');

            reporting.ShowProgress('Finding lung region for density');
            mask = left_and_right_lungs.RawImage > 0;
            roi_lung = int16(roi.RawImage).*int16(mask);

            interp_voxel_size_units = interp_voxel_size_mm./roi.VoxelSize;
            image_size = roi.ImageSize;
            i_span = 1 : interp_voxel_size_units(1) : image_size(1);
            j_span = 1 : interp_voxel_size_units(2) : image_size(2);
            k_span = 1 : interp_voxel_size_units(3) : image_size(3);
            [interp_i, interp_j, interp_k] = ndgrid(i_span, j_span, k_span);

            % Interpolate dentisty values
            reporting.ShowProgress('Interpolating to new voxel grid');
            roi_interp = interpn(single(roi_lung), interp_i, interp_j, interp_k, '*linear');

            % Interpolate mask (gives a measure of how much of each voxel is
            % within the lung vs outside the lung)
            reporting.ShowProgress('Interpolating to new voxel size');
            roi_mask = interpn(single(mask), interp_i, interp_j, interp_k, '*linear');

            % Rescale voxels that are partially outside the lung, and remove
            % those more than 50% outside of the lung
            roi_interp = roi_interp./roi_mask;
            roi_interp(roi_mask < 0.5) = 0;

            roi_interp = roi_interp/max(roi_interp(:));

            reporting.ShowProgress('Finding coordinates on original grid');
            i_span_r = single(1 + (1:image_size(1))/interp_voxel_size_units(1));
            j_span_r = single(1 + (1:image_size(2))/interp_voxel_size_units(2));
            k_span_r = single(1 + (1:image_size(3))/interp_voxel_size_units(3));

            [interp_i, interp_j, interp_k] = ndgrid(i_span_r, j_span_r, k_span_r);
            interp_i = min(size(roi_interp, 1), round(interp_i));
            interp_j = min(size(roi_interp, 2), round(interp_j));
            interp_k = min(size(roi_interp, 3), round(interp_k));

            reporting.ShowProgress('Interpolation to original grid');
            indices = sub2ind(size(roi_interp), (interp_i(:)), (interp_j(:)), (interp_k(:)));
            results_raw = zeros(image_size, 'single');
            results_raw(:) = roi_interp(indices);

            results = roi.BlankCopy;

            results.ChangeRawImage(single(results_raw));
            results.ImageType = PTKImageType.Scaled;

        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function y=polyValVec(P,x)
%%POLYVALVEC Evaluate a set of polynomials at a given scalar point or set
%            of points. Evaluating a collection of polynomials produces a
%            vector estimate. This could, for example, be used to evaluate
%            polynomials to interpolate each component of a vector quantity
%            to a desired point given the ointerpolating polynomials for
%            each dimension.
%
%INPUTS: P A numDimX(numDeg+1) matrix of numDim polynomial coefficients for
%          each of the numDim components. numDeg is the degree of the
%          polynomials. The coefficients for each row are arranged for
%          exponents in the order [x^numDeg, x^(numDeg-1),...,1,0].
%        x A scalar value of x or a linear vector of numX points at which
%          the polynomials should be evaluated.
%
%OUTPUTS: y The numDim X numX set of polynomials evaluated at each of the
%           numX points.
%
%The algorithm just calls Matlab's polyval function multiple times in a
%loop.
%
%February 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

numDim=size(P,1);
numX=length(x);

%Make sure that it is a row vector.
x=reshape(x,1,numX);

y=zeros(numDim,numX);
for curDim=1:numDim
    y(curDim,:)=polyval(P(curDim,:),x);
end
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
% Script to reproduce the experiments leading to the results provided in the
% Table 2 of the paper "Deep Scattering Spectrum" by J. Andén and S. Mallat.

% Delta-Delta-MFCCs for window size 370 ms

run_name = 'DSS_Table2_GTZAN_mfcc_370ms';

N = 5*2^17;

src = gtzan_src('/path/to/gtzan');

filt1_opt.wavelet_type = {'gabor','morlet'};
filt1_opt.Q = [8 2];
filt1_opt.J = T_to_J(8192,filt1_opt);
filt1_opt.boundary = 'symm';

sc1_opt.M = 1;

filters = filter_bank(N, filt1_opt);

features = {@(x)(feval(@(x)([x; circshift(x,[0 +1]); circshift(x,[0 -1])]), ...
	format_scat(log_scat(spec_freq_average(x,filters,sc1_opt)))))};

db = prepare_database(src,features);
db.features = single(db.features);
db = svm_calc_kernel(db,'gaussian','square',1:2:size(db.features,2));

rs = RandStream.create('mt19937ar','Seed',floor(pi*1e9));
RandStream.setGlobalStream(rs);
[train_set{1}, test_set{1}] = create_partition([src.objects.class], 0.9);
for k = 2:10
	[train_set{k}, test_set{k}] = ...
		next_fold([src.objects.class], train_set{k-1}, test_set{k-1});
end

optt.kernel_type = 'gaussian';
optt.C = 2.^[0:4:8];
optt.gamma = 2.^[-16:4:-8];
optt.search_depth = 3;
optt.full_test_kernel = 1;

for k = 1:10
	[dev_err_grid,C_grid,gamma_grid] = ...
		svm_adaptive_param_search(db,train_set{k},[],optt);

	[dev_err(k),ind] = min(mean(dev_err_grid{end},2));
	C(k) = C_grid{end}(ind);
	gamma(k) = gamma_grid{end}(ind);

	optt1 = optt;
	optt1.C = C(k);
	optt1.gamma = gamma(k);

	model = svm_train(db,train_set{k},optt1);
	labels(:,k) = svm_test(db,model,test_set{k});
	err(k) = classif_err(labels(:,k),test_set{k},db.src);

	fprintf('dev err = %f, test err = %f\n',dev_err(k),err(k));

	save([run_name '.mat'],'labels','dev_err','err','C','gamma');
end

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function [tri, x, map] = tri_squeeze(tri, x)
% TRI_SQUEEZE  Remove disconnected vertices from triangular mesh
%
% [TRI2, X2, I] = tri_squeeze(TRI, X)
%
%   TRI is a 3-column matrix. Each row contains the 3 nodes that form one
%   triangular facet in the mesh.
%
%   X is a 3-column matrix. X(i, :) contains the xyz-coordinates of the
%   i-th node in the mesh.
%
%   Assuming that TRI, X represents a triangular mesh where some of the
%   vertices in X are not part of any triangle in TRI, this function
%   returns TRI2, X2, that represent the same mesh, but with the
%   disconnected vertices removed.
%
%   I is the mapping of preserved vertices such that X2 = X(I), given as a
%   list of indices.

% Author: Ramon Casero <rcasero@gmail.com>
% Copyright © 2012-2014 University of Oxford
% Version: 0.2.0
%
% University of Oxford means the Chancellor, Masters and Scholars of
% the University of Oxford, having an administrative office at
% Wellington Square, Oxford OX1 2JD, UK.
%
% This file is part of Gerardus.
%
% 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. The offer of this
% program under the terms of the License is subject to the License
% being interpreted in accordance with English Law and subject to any
% action against the University of Oxford being under the jurisdiction
% of the English Courts.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see
% <http://www.gnu.org/licenses/>.

% check arguments
narginchk(2, 2);
nargoutchk(0, 3);

% number of vertices
N = size(x, 1);

% list of all vertices that are part of at least a triangle
idx = unique(tri(:));

% vector that will do the mapping
map = zeros(N, 1);
map(idx) = (1:length(idx))';

% rename vertices
tri = map(tri);

% delete disconnected vertices
x(~map, :) = [];

if (nargout > 2)
    % format the mapping of preserved vertices as indices instead of a
    % boolean vector
    map = find(map);
end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function [strategy, MEU, niter] = solve_limid(engine, varargin)
% SOLVE_LIMID Find the (locally) optimal strategy for a LIMID
% [strategy, MEU, niter] = solve_limid(inf_engine, ...)
%
% strategy{d} = stochastic policy for node d (a decision node)
% MEU = maximum expected utility
% niter = num iterations used
%
% The following optional arguments can be specified in the form of name/value pairs:
% [default in brackets]
%
% max_iter - max. num. iterations [ 1 ]
% tol - tolerance required of consecutive MEU values, used to assess convergence [1e-3]
% order - order in which decision nodes are optimized [ reverse numerical order ]
%
% e.g., solve_limid(engine, 'tol', 1e-2, 'max_iter', 10)

bnet = bnet_from_engine(engine);

% default values
max_iter = 1;
tol = 1e-3;
D = bnet.decision_nodes;
order = D(end:-1:1);

args = varargin;
nargs = length(args);
for i=1:2:nargs
  switch args{i},
   case 'max_iter', max_iter  = args{i+1};
   case 'tol',      tol = args{i+1};
   case 'order',    order = args{i+1};
   otherwise,
    error(['invalid argument name ' args{i}]);
  end
end

CPDs = bnet.CPD;
ns = bnet.node_sizes;
N = length(ns);
evidence = cell(1,N);
strategy = cell(1, N);

iter = 1;
converged = 0;
oldMEU = 0;
while ~converged & (iter <= max_iter)
  for d=order(:)'
    engine = enter_evidence(engine, evidence, 'exclude', d);
    [m, pot] = marginal_family(engine, d);
    %pot = marginal_family_pot(engine, d);
    [policy, score] = upot_to_opt_policy(pot);
    e = bnet.equiv_class(d);
    CPDs{e} = set_fields(CPDs{e}, 'policy', policy);
    engine = update_engine(engine, CPDs);
    strategy{d} = policy;
  end
  engine = enter_evidence(engine, evidence);
  [m, pot] = marginal_nodes(engine, []);
  %pot = marginal_family_pot(engine, []);
  [dummy, MEU] = upot_to_opt_policy(pot);
  if approxeq(MEU, oldMEU, tol)
    converged = 1;
  end
  oldMEU = MEU;
  iter = iter + 1;
end
niter = iter - 1;

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [rVectB, vVectB] = getPositOfBodyWRTSun(time, bodyInfo, celBodyData)
    try
        if(bodyInfo.propTypeEnum == BodyPropagationTypeEnum.TwoBody || (numel(time) == 1 && time == bodyInfo.epoch))
            chain = bodyInfo.getOrbitElemsChain();
            [rVectB, vVectB] = getPositOfBodyWRTSun_alg_fast_mex(time, chain{:});
        elseif(bodyInfo.propTypeEnum == BodyPropagationTypeEnum.Numerical)
            parentBodyInfo = bodyInfo.getParBodyInfo(celBodyData);

            [rVect, vVect] = getStateAtTime(bodyInfo, time, []);
            cartElem = CartesianElementSet(time, rVect, vVect, parentBodyInfo.getBodyCenteredInertialFrame(), true);
            cartElem = convertToFrame(cartElem, celBodyData.getTopLevelBody().getBodyCenteredInertialFrame());

            rVectB = [cartElem.rVect];
            vVectB = [cartElem.vVect];

        else
            error('Unknown celestial body prop sim type.');
        end
    catch ME %#ok<NASGU> %if bad then use old way
        numTimes = length(time);
        loop = true;
        rVectB = zeros(3,numTimes);
        vVectB = zeros(3,numTimes);
        while(loop)
            try
                parentBodyInfo = bodyInfo.getParBodyInfo(celBodyData);
            catch
                parentBodyInfo = getParentBodyInfo(bodyInfo, celBodyData);
            end

            if(isempty(parentBodyInfo))
                break;
            end

            try
                parentGM = bodyInfo.getParentGmuFromCache();
            catch
                parentGM = getParentGM(bodyInfo,celBodyData);
            end

            [rVect, vVect] = getStateAtTime(bodyInfo, time, parentGM); %getParentGM(bodyInfo, celBodyData)
            rVectB = rVectB + rVect;
            vVectB = vVectB + vVect;

            bodyInfo = parentBodyInfo;
        end
    end

%     if(numel(time) > 1 && numel(time) == numel(bodyInfo))
%         [uniBodies,~,ic] = unique(bodyInfo,'stable');
%
%         rVectB = NaN(3, numel(time));
%         vVectB = rVectB;
%         for(i=1:length(uniBodies))
%             subBodyInfo = uniBodies(i);
%             bool = i == ic;
%             subTimes = time(bool);
%
%             [rVectB(:,bool), vVectB(:,bool)] = getPositOfBodyWRTSun(subTimes, subBodyInfo, celBodyData);
%         end
%     else
%         try %try new way
%             if(bodyInfo.propTypeEnum == BodyPropagationTypeEnum.TwoBody || (numel(time) == 1 && time == bodyInfo.epoch))
% %                 if(numel(time) == 1 && numel(bodyInfo.lastComputedTime) > 0 && bodyInfo.lastComputedTime == time)
% %                     rVectB = bodyInfo.lastComputedRVect;
% %                     vVectB = bodyInfo.lastComputedVVect;
% %                 else
% %                     chain = bodyInfo.getOrbitElemsChain();
% %                     [rVectB, vVectB] = getPositOfBodyWRTSun_alg(time, chain{:});
% %
% %                     if(numel(time) == 1)
% %                         bodyInfo.lastComputedTime = time;
% %                         bodyInfo.lastComputedRVect = rVectB;
% %                         bodyInfo.lastComputedVVect = vVectB;
% %                     end
% %                 end
%
%                 if(numel(time) == 1)
%                     bool = time == bodyInfo.lastComputedTime;
%
%                     if(any(bool))
%                         rVectVVectB = bodyInfo.lastComputedRVectVVect(:,bool);
%                         rVectB = rVectVVectB(1:3,1);
%                         vVectB = rVectVVectB(4:6,1);
%
%                     else
%                         chain = bodyInfo.getOrbitElemsChain();
%                         try
%                             [rVectB, vVectB] = getPositOfBodyWRTSun_alg_fast_mex(time, chain{:});
%                         catch ME %#ok<NASGU>
%                             [rVectB, vVectB] = getPositOfBodyWRTSun_alg(time, chain{:});
%                         end
%
%                         bodyInfo.lastComputedTime = horzcat(bodyInfo.lastComputedTime, time);
%                         bodyInfo.lastComputedRVectVVect = horzcat(bodyInfo.lastComputedRVectVVect, [rVectB;vVectB]);
%
%                         [bodyInfo.lastComputedTime,ia,~] = unique(bodyInfo.lastComputedTime,'stable');
%                         bodyInfo.lastComputedRVectVVect = bodyInfo.lastComputedRVectVVect(:, ia);
%                     end
%                 else
%                     chain = bodyInfo.getOrbitElemsChain();
%                     try
%                         [rVectB, vVectB] = getPositOfBodyWRTSun_alg_fast_mex(time, chain{:});
%                     catch ME %#ok<NASGU>
%                         [rVectB, vVectB] = getPositOfBodyWRTSun_alg(time, chain{:});
%                     end
%
%                     bodyInfo.lastComputedTime = horzcat(bodyInfo.lastComputedTime, time);
%                     bodyInfo.lastComputedRVectVVect = horzcat(bodyInfo.lastComputedRVectVVect, [rVectB;vVectB]);
%
%                     [bodyInfo.lastComputedTime,ia,~] = unique(bodyInfo.lastComputedTime,'stable');
%                     bodyInfo.lastComputedRVectVVect = bodyInfo.lastComputedRVectVVect(:, ia);
%                 end
%
%             elseif(bodyInfo.propTypeEnum == BodyPropagationTypeEnum.Numerical)
%                 parentBodyInfo = bodyInfo.getParBodyInfo(celBodyData);
%
%                 [rVect, vVect] = getStateAtTime(bodyInfo, time, []);
%                 cartElem = CartesianElementSet(time, rVect, vVect, parentBodyInfo.getBodyCenteredInertialFrame(), true);
%                 cartElem = convertToFrame(cartElem, celBodyData.getTopLevelBody().getBodyCenteredInertialFrame());
%
%                 rVectB = [cartElem.rVect];
%                 vVectB = [cartElem.vVect];
%
%             else
%                 error('Unknown celestial body prop sim type.');
%             end
%         catch ME %if bad then use old way
%             numTimes = length(time);
%             loop = true;
%             rVectB = zeros(3,numTimes);
%             vVectB = zeros(3,numTimes);
%             while(loop)
%                 try
%                     parentBodyInfo = bodyInfo.getParBodyInfo(celBodyData);
%                 catch
%                     parentBodyInfo = getParentBodyInfo(bodyInfo, celBodyData);
%                 end
%
%                 if(isempty(parentBodyInfo))
%                     break;
%                 end
%
%                 try
%                     parentGM = bodyInfo.getParentGmuFromCache();
%                 catch
%                     parentGM = getParentGM(bodyInfo,celBodyData);
%                 end
%
%                 [rVect, vVect] = getStateAtTime(bodyInfo, time, parentGM); %getParentGM(bodyInfo, celBodyData)
%                 rVectB = rVectB + rVect;
%                 vVectB = vVectB + vVect;
%
%                 bodyInfo = parentBodyInfo;
%             end
%         end
%     end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function truncated_normal_test22 ( )

%*****************************************************************************80
%
%% TRUNCATED_NORMAL_B_PDF_TEST tests TRUNCATED_NORMAL_B_PDF;
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 September 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TRUNCATED_NORMAL_B_PDF_TEST\n' );
  fprintf ( 1, '  TRUNCATED_NORMAL_B_PDF evaluates the Upper Truncated Normal PDF.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The "parent" normal distribution has\n' );
  fprintf ( 1, '    mean = mu\n' );
  fprintf ( 1, '    standard deviation = sigma\n' );
  fprintf ( 1, '  The parent distribution is truncated to\n' );
  fprintf ( 1, '  the interval (-oo,b]\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '                                                  Stored         Computed\n' );
  fprintf ( 1, '       X        Mu         S         B             PDF             PDF\n' );
  fprintf ( 1, '\n');

  n_data = 0;

  while ( 1 )

    [ n_data, mu, sigma, b, x, pdf1 ] ...
      = truncated_normal_b_pdf_values ( n_data );

    if ( n_data == 0 )
      break
    end

    pdf2 = truncated_normal_b_pdf ( x, mu, sigma, b );

    fprintf( 1, '  %8.1f  %8.1f  %8.1f  %8.1f  %14g  %14g\n', x, mu, sigma, b, pdf1, pdf2 );

  end

  return
end


text:
function [PSNR, SSIM, IFC] = evaluate_SR(img_GT, img_HR, scale, compute_ifc)
% -------------------------------------------------------------------------
%   Description:
%       Compute PSNR, SSIM and IFC for SR
%       We convert RGB image to grayscale and crop boundaries for 'scale'
%       pixels
%
%   Input:
%       - img_GT        : Ground truth image
%       - img_HR        : predicted HR image
%       - scale         : upsampling scale
%       - compute_ifc   : evaluate IFC [default = 0 since it's slow]
%
%   Citation:
%       Fast and Accurate Image Super-Resolution with Deep Laplacian Pyramid Networks
%       Wei-Sheng Lai, Jia-Bin Huang, Narendra Ahuja, and Ming-Hsuan Yang
%       arXiv, 2017
%
%   Contact:
%       Wei-Sheng Lai
%       wlai24@ucmerced.edu
%       University of California, Merced
% -------------------------------------------------------------------------

    if ~exist('compute_ifc', 'var')
        compute_ifc = 0;
    end

    %% quantize pixel values
    img_GT = im2double(im2uint8(img_GT));
    img_HR = im2double(im2uint8(img_HR));

    %% convert to gray scale
    if( size(img_GT, 3) > 1 )
        img_GT = rgb2ycbcr(img_GT); img_GT = img_GT(:, :, 1);
        img_HR = rgb2ycbcr(img_HR); img_HR = img_HR(:, :, 1);
    end

    %% crop boundary
    img_GT = shave_bd(img_GT, scale);
    img_HR = shave_bd(img_HR, scale);

    % evaluate
    PSNR = psnr(img_GT, img_HR);
    SSIM = ssim(img_GT, img_HR);

    % comment IFC to speed up testing
    IFC = 0;
    if compute_ifc
        IFC = ifcvec(img_GT, img_HR);
        if( ~isreal(IFC) )
            IFC = 0;
        end
    end

end

repo: LapSRN


text:
% The COBRAToolbox: testEfmBackboneExtraction.m
%
% Purpose:
%     - test extraction of backbone from a given set of EFMs
%
% Authors:
%     - Created initial test script: Chaitra Sarathy 2 Dec 19
%     - Updates to header docs: Chaitra Sarathy 19 Dec 19

solvers = prepareTest('requiredToolboxes', {'statistics_toolbox'});

global CBTDIR

% save the current path and initialize the test
currentDir = cd(fileparts(which(mfilename)));

% determine the test path for references
testPath = pwd;

% Load reference data
load([testPath filesep 'testData_efmBackboneExtraction.mat']);

% This is a function with two required inputs:
% Case 1: Test whether efmBackboneExtraction gives the expected output when using 2 input and 2 output arguments
[selectedRxns,  rxnDistribution] = efmBackboneExtraction(testEFMRxns, percentage);
assert(isequal(selectedRxns, selectedRxns_ref));
assert(isequal(rxnDistribution,rxnDistribution_ref));

% Case 2: Test whether efmBackboneExtraction throws an error when using 1 input argument only
assert(verifyCobraFunctionError('efmBackboneExtraction', 'inputs', {testEFMRxns}, 'outputArgCount', 1));

% Case 3: Test whether efmBackboneExtraction throws an error when using NO input arguments
assert(verifyCobraFunctionError('efmBackboneExtraction', 'inputs', {}, 'outputArgCount', 1));

author: opencobra

repo: cobratoolbox


text: %% S2FunHarmonicSym

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
clear;
netStruct = load('../data/vgg16_cuhk_batch32_Rankloss_2:1:0.5_margin1_re/net-epoch-60.mat');
net = dagnn.DagNN.loadobj(netStruct.net);
clear netStruct;
net.mode = 'test' ;
net.move('gpu') ;
net.removeLayer('RankLoss');
net.conserveMemory = true;
im_mean = reshape(net.meta.normalization.averageImage,1,1,3);

load('../dataset/CUHK-PEDES-prepare/url_data.mat');
p = imdb.images.data(imdb.images.set==3);
ff = [];
%%------------------------------

for i = 1:numel(p)
    disp(i);
    str = p{i};
    im = imread(str);
    oim = im;
    f = getFeature2(net,oim,im_mean,'x0','fc1_1bn');
    f = sum(sum(f,1),2);
    f2 = getFeature2(net,fliplr(oim),im_mean,'x0','fc1_1bn');
    f2 = sum(sum(f2,1),2);
    f = f+f2;
    size4 = size(f,4);
    f = reshape(f,[],size4)';
    f = norm_zzd(f);
    ff = cat(1,ff,f);
end
save('./resnet_cuhk_img.mat','ff','-v7.3');
%}

ff = [];
load('../dataset/CUHK-PEDES-prepare/cuhk_word2.mat');
test_word = wordcnn(:,end-6155:end);

for i = 1:6156
    disp(i);
    content = test_word(:,i);
    len = numel(find(content>0));
    txtinput = zeros(len,7263,'single');
    for k=1:len
        txtinput(k,content(k))=1;
    end
    %transfer it to different location
    win = 57-len;
    input = zeros(56,7263,win,'single');
    for kk = 1:win
        input(kk:kk+len-1,:,kk) = txtinput;
    end
    input = reshape(input,1,56,7263,[]);
    f = getFeature2(net,input,[],'data2','fc5_2bn');
    f = sum(f,4);
    size4 = size(f,4);
    f = reshape(f,[],size4)';
    f = norm_zzd(f);
    ff = cat(1,ff,f);
end
save('./resnet_cuhk_txt.mat','ff','-v7.3');

evaluate;

author: layumi

repo: Image-Text-Embedding


text:
% Base class for all estimators in KAFBOX.
%
% This file is part of the Kernel Adaptive Filtering Toolbox for Matlab.
% https://github.com/steven2358/kafbox/

classdef base_estimator < matlab.mixin.Copyable
    methods
        % get parameter names for the estimator
        function names = get_param_names(obj)
            names = fieldnames(obj);
        end

        % get parameters
        function params = get_params(obj)
            params = struct;
            for fn = fieldnames(obj)'
                params.(fn{1}) = obj.(fn{1});
            end
        end

        % set parameters
        function set_params(obj,params)
            if (nargin > 0) % copy valid parameters
                for fn = fieldnames(params)'
                    if ismember(fn{1},fieldnames(obj))
                        values = params.(fn{1});
                        obj.(fn{1}) = values;
                    else
                        warning('Unknown parameter: %s.',fn{1});
                    end
                end
            end
        end
    end
end

repo: kafbox


text:
%--------------------------------------------------------
% FIR filter design via local moving window LS fitting  -
% A magic smooth and derivative formula generator.      -
% By Dr Yangquan Chen		08-07-1999                   -
% Email=<yqchen@ieee.org>; URL=http://www.crosswinds.net/~yqchen/
% -------------------------------------------------------
% Purpose: general FIR design via local LS fitting.
% total taps = nL+nR+1
% Format: function [c]=sgfilter(nL,nR,M,id)
%          -nL      -nL+1           -1                nR
% FIR=c(1)z   +c(2)z     +...+c(nL)z  +...+c(nL+nR+1)z
%
%       nR
%     ------
%     \                   j
%      >       c(nL+1+j) z
%     /
%     ------
%      j=-nL
% M: the order of LS fit at the moving window of [-nL, ... , nR]
% id: index for the derivative order
%		0: smooth filter,
%		1: 1st order differentiator,
% 		2: 2nd order differentiator, ...
% NOTE: M>=id, set M=(2~4)*(id+1) for reliably results.
%		  to do LS fit, M<nL+nR+1.
%---------------------------------------------------------------
function [c]=sgfilter(nL,nR,M,id)
% Savitzky-Golay smoothing filter.
if (id>M)
   disp('Error in id! (id<M)');return;
end
if (M>(nL+nR))
   disp('Error in M! (M<=nL+nR)');return;
end

A=zeros(nL+nR+1,M+1);
for i=-nL:nR;
   for j=0:M;
      A(i+nL+1,j+1)=i^j;
   end
end
h=zeros(M+1,1);
%h(1)=1;
h(id+1)=1;
b=inv(A'*A)*h;
c=zeros(nL+nR+1,1);
for n=-nL:nR
   nm=n.^[0:M];
   c(n+nL+1)=nm*b;
end
% coefficient for smoothing
return

author: Sable

repo: mcbench-benchmarks


text:
classdef PTKDicomMetadata < handle
    % PTKDicomMetadata. Returns tags for creating Dicom files
    %
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2012.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    properties (Constant)
        DicomVersion = '0.1'
        DicomName = 'TD Pulmonary Toolkit'
        DicomManufacturer = 'www tomdoel com'
        DicomStudyDescription = 'TD Pulmonary Toolkit exported images'
    end
end


text:
function [us,optimState] = searchHedge(searchFcns,u,gpstruct,LB,UB,optimState,options)
%SEARCHWCM Weighted covariance matrix search step (inspired by CMA-ES).

if nargin < 2 || isempty(u) || isempty(gpstruct)
    try
        us = feval(searchFcns{optimState.hedge.chosen}{:});
    catch
        us = 'hedge';
    end
    return;
end

% SUMRULE not specified, all arguments shifted by one
if nargin < 7
    options = optimState;
    optimState = UB;
    UB = LB;
    LB = gpstruct;
    gpstruct = u;
    u = sumrule;
    sumrule = [];
end

hedge = optimState.hedge;

nh = numel(searchFcns);

% Initialize hedge struct
if isempty(hedge)
    hedge.g = zeros(1,nh);
    hedge.g(1) = 10;
    for i = 1:nh; hedge.str{i} = feval(searchFcns{i}{:}); end
    hedge.n = numel(hedge.g);
    hedge.count = 0;
    hedge.gamma = options.HedgeGamma;
    hedge.beta = options.HedgeBeta;
    hedge.decay = options.HedgeDecay;
    % hedge.string = {'mpi','mei'};
end

hedge.count = hedge.count + 1;
% gammaHedge = min(1/nHedge, sqrt(log(nHedge)/(nHedge*hedge.count)));

hedge.p = exp(hedge.beta*(hedge.g - max(hedge.g)))./sum(exp(hedge.beta*(hedge.g - max(hedge.g))));
hedge.p = hedge.p*(1-hedge.n*hedge.gamma) + hedge.gamma;
% hedge.p
hedge.chosen = find(rand() < cumsum(hedge.p),1);
if isempty(hedge.chosen)
    hedge.chosen = randi(nh);
    if options.gpWarnings
        warning('bads:searchHedgeFail', ['Cannot determine best search function in SEARCHHEDGE (P=' numarray2str(hedge.p) '). Attempting to continue.']);
    end
end
if hedge.gamma == 0
    hedge.phat = ones(size(hedge.g));
else
    hedge.phat = Inf(size(hedge.p));
    hedge.phat(hedge.chosen) = hedge.p(hedge.chosen);
end

us = feval(searchFcns{hedge.chosen}{:}, ...
    u, ...
    gpstruct, ...
    LB, ...
    UB, ...
    optimState, ...
    options);

optimState.hedge = hedge;

end

author: acerbilab

repo: bads

save_path: github-repos/MATLAB/acerbilab-bads