text:
function pagerankdemo (steps)
% PAGERANKDEMO draw a 6-node web and compute its pagerank
%
% PAGERANKDEMO draws the 6-node "tiny web" in Section 2.11 of "Numerical
% Computing with MATLAB", by Cleve Moler, SIAM, 2004.  It then simulates the
% computation of Google's PageRank algorithm, by randomly selecting links to
% traverse.  If a link is traversed, the edge and the target node are displayed
% in red.  If the "random surfer" jumps to an arbitrary page, the target node
% is displayed in blue.  The number of hits at each node, and the page rank
% (in %) are displayed % on each node.  Note that after a large number of
% steps, the PageRanks (in percentages) converge to the values given in Section
% 2.11 of Moler (alpha: .321, sigma: .2007, beta: .1705, delta: .1368,
% gamma: .1066, rho: .0643).  See http://www.mathworks.com/moler for more
% details (the pagerank M-file, in particular).
%
% Note that this method is NOT how the PageRank is actually computed.  Instead
% the eigenvalue problem A*x=x is solved for x, where A is the Markov
% transition matrix, A = p*G*D + e*z', where G is the binary matrix used here.
% The method here is a simplistic random-hopping demonstration of the Markov
% process, to motivate the A*x=x formulation of the problem.  In this example,
% A does control how the transitions are made, but the matrix A is not formed
% explicitly.
%
% This demo only operates on a single graph.  It is meant as a simple demo
% only, suitable for in-class use.  To compute the PageRanks for an arbitrary
% graph, use pagerank.m, or the power method (repeat x=A*x until convergence,
% where A is the Markov transition matrix of the web).
%
% Example:
%   pagerankdemo
%   pagerankdemo (1000)   % run 1000 steps with no user input, then quit
%
% See also pagerank
%
% I suggest single-stepping a dozen times or so to see the link traversal in
% process, and then type "1000".  Hit control-C to quit.
%
% Copyright 2007, Tim Davis, University of Florida

% Initial graph
Graph = graphinit ;
rand ('state', 0) ;
n = size (Graph.G, 1) ;

help pagerankdemo

% initialize the page counts
hits = zeros (1,n) ;
oldwhere = 1 ;
where = 1 ;
hits (where) = 1 ;
set (Graph.node (where), 'FaceColor', [0 0 1]) ;

p = 0.85 ;		% probability a link will be followed
c = sum (Graph.G) ;	% outgoing degree

links = cell (1,n) ;
for k = 1:n
    links {k} = find (Graph.G (:,k)) ;
end

follow_link = 0 ;
if (nargin < 1)
    input ('hit enter to start at node alpha: ') ;
end

% write the stats to the figure
set (Graph.nodelabel (where), 'string', ...
	sprintf ('%s %d (%3.1f%%)', Graph.nodes {where}, hits (where), ...
	100 * hits (where) / sum (hits))) ;

if (nargin < 1)
    input ('hit enter to take one step: ') ;
    steps = 1 ;
end

% repeat
while (1)

    % clear the old color and old arrow
    set (Graph.node (where), 'FaceColor', [0 1 0]) ;
    if (follow_link)
	set (Graph.arrows (where,oldwhere), 'LineWidth', 2) ;
	set (Graph.arrows (where,oldwhere), 'Color', [0 0 0]) ;
    end

    % determine where to go to next
    oldwhere = where ;
    if (c (where) == 0 || rand > p)
	% no outgoing links, or ignore the links
	follow_link = 0 ;
	where = floor (n * rand + 1) ;
	set (Graph.node (where), 'FaceColor', [0 0 1]) ;
    else
	% move along the link
	follow_link = 1 ;
	where = links{where}(floor (c (where) * rand + 1)) ;
	set (Graph.node (where), 'FaceColor', [1 0 0]) ;
	set (Graph.arrows (where,oldwhere), 'LineWidth', 5) ;
	set (Graph.arrows (where,oldwhere), 'Color', [1 0 0]) ;
    end

    % increment the hit count
    hits (where) = hits (where) + 1 ;

    % write the stats to the figure
    for k = 1:n
	set (Graph.nodelabel (k), 'string', ...
	sprintf ('%s %d (%3.1f%%)', Graph.nodes {k}, hits (k), ...
	100 * hits (k) / sum (hits))) ;
    end

    drawnow

    % go the next step
    steps = steps - 1 ;
    if (steps <= 0)
        if (nargin > 0)
            break ;
        end
        steps = input ...
	    ('number of steps to make (default 1, control-C to quit): ') ;
	if (steps == 0)
	    break ;
	end
	if (isempty (steps))
	    steps = 1 ;
	end
    end

end

%-------------------------------------------------------------------------------

function Graph = graphinit
% GRAPHINIT create the tiny-web example in Moler, section 2.11, and draw it.
% Example
%   G = graphinit ;

figure (1)
clf

nodes = { 'alpha', 'beta', 'gamma', 'delta', 'rho', 'sigma' } ;

xy = [
0 4
1 3
1 2
2 4
2 0
0 0
] ;

x = xy (:,1) ;
y = xy (:,2) ;

% scale x and y to be in the range 0.1 to 0.9
x = 0.8 * x / 2 + .1 ;
y = 0.8 * y / 4 + .1 ;
xy = [x y] ;

xy_delta = [
 .08 .04 0
-.03 -.02 -1
 .04  0   0
-.05  .04 -1
-.03  0  -1
 .03  0   0
] ;

xd = xy_delta (:,1) ;
yd = xy_delta (:,2) ;
tjust = xy_delta (:,3) ;

G = [
0 0 0 1 0 1
1 0 0 0 0 0
0 1 0 0 0 0
0 1 1 0 0 0
0 0 1 0 0 0
1 0 1 0 0 0 ] ;

clf

n = size (G,1) ;

axes ('Position', [0 0 1 1], 'Visible', 'off') ;

node = zeros (n,1) ;
nodelabel = zeros (n,1) ;
for k = 1:n
    node (k) = annotation ('ellipse', [x(k)-.025 y(k)-.025 .05 .05]) ;
    set (node (k), 'LineWidth', 2) ;
    set (node (k), 'FaceColor', [0 1 0]) ;
    nodelabel (k) = text (x (k) + xd (k), y (k) + yd (k), nodes {k}, ...
	'Units', 'normalized', 'FontSize', 16) ;
    if (tjust (k) < 0)
	set (nodelabel (k), 'HorizontalAlignment', 'right') ;
    end
end

axis off

% Yes, I realize that this is overkill; arrows should be sparse.
% This example is not meant for large graphs.
arrows = zeros (n,n) ;

[i j] = find (G) ;
for k = 1:length (i)
    % get the center of the two nodes
    figx = [x(j(k)) x(i(k))] ;
    figy = [y(j(k)) y(i(k))] ;
%   [figx figy] = dsxy2figxy (gca, axx, axy);
    % shorten the arrows by s units at each end
    s = 0.03 ;
    len = sqrt (diff (figx)^2 + diff (figy)^2) ;
    fy (1) = diff (figy) * (s/len) + figy(1) ;
    fy (2) = diff (figy) * (1-s/len) + figy(1) ;
    fx (1) = diff (figx) * (s/len) + figx(1) ;
    fx (2) = diff (figx) * (1-s/len) + figx(1) ;
    arrows (i(k),j(k)) = annotation ('arrow', fx, fy) ;
    set (arrows (i(k),j(k)), 'LineWidth', 2) ;
    set (arrows (i(k),j(k)), 'HeadLength', 20) ;
    set (arrows (i(k),j(k)), 'HeadWidth', 20) ;
end

Graph.G = G ;
Graph.nodes = nodes ;
Graph.node = node ;
Graph.xy = xy ;
Graph.xy_delta = xy_delta ;
Graph.nodelabel = nodelabel ;
Graph.arrows = arrows ;

author: fieldtrip


text:
function [y]=tt_mvk4(A, x, tol, varargin)

% step_dpow = 0.1; % stepsize for d-power in truncations
% min_dpow = 1; % Minimal d-power for truncation
%
% bot_conv = 0.1; % bottom convergence factor - if better, we can decrease dpow and drank
% top_conv = 0.99; % top convergence factor - if worse, we have to increase dpow and drank

nswp=10;

als_tol_low = 5;
als_iters = 4;

rmax=1000;

verb=1;
kickrank = 10;
 kicktype = 'mr';
% kicktype = 'rand';

% pcatype = 'svd';
pcatype = 'uchol';

y=[];
block_order = [];

for i=1:2:length(varargin)-1
    switch lower(varargin{i})
        case 'nswp'
            nswp=varargin{i+1};
        case 'rmax'
            rmax=varargin{i+1};
        case 'kicktype'
            kicktype=lower(varargin{i+1});
        case 'pcatype'
            pcatype=lower(varargin{i+1});
        case 'y0'
            y=varargin{i+1};
        case 'verb'
            verb=varargin{i+1};
        case 'kickrank'
            kickrank=varargin{i+1};
        case 'step_dpow'
            step_dpow=varargin{i+1};
        case 'min_dpow'
            min_dpow=varargin{i+1};
        case 'bot_conv'
            bot_conv=varargin{i+1};
        case 'top_conv'
            top_conv=varargin{i+1};
        case 'block_order'
            block_order=varargin{i+1};
        case 'als_tol_low'
            als_tol_low=varargin{i+1};
        case 'als_iters'
            als_iters=varargin{i+1};

        otherwise
            error('Unrecognized option: %s\n',varargin{i});
    end
end

rx = x.r;
d = x.d;
tailranks = [false,false];
if (rx(1)~=1)
    e1 = tt_tensor;
    e1.d=1;
    e1.n = [rx(1)];
    e1.r = [1; rx(1)];
    e1.ps = [1; rx(1)^2+1];
    e1.core = reshape(eye(rx(1)), rx(1)^2, 1);
    x = kron(e1, x);

    A = kron(tt_matrix(eye(rx(1))), A);
    tailranks(1)=true;
end;
if (rx(d+1)~=1)
    e1 = tt_tensor;
    e1.d=1;
    e1.n = [rx(d+1)];
    e1.r = [rx(d+1); 1];
    e1.ps = [1; rx(d+1)^2+1];
    e1.core = reshape(eye(rx(d+1)), rx(d+1)^2, 1);
    x = kron(x, e1);

    A = kron(A, tt_matrix(eye(rx(d+1))));
    tailranks(2)=true;
end;

d = x.d;
rx = x.r;
ra = A.r;
n = A.n;
m = A.m;
if (isempty(y))
%     y = tt_rand(n, A.d, kickrank);
    y = tt_ones(n, A.d);
end;

if (isempty(block_order))
    block_order = [+(d), -(d)];
end;

ry = y.r;

cry = core2cell(y);
crA = core2cell(A);
crx = core2cell(x);

phia = cell(d+1,1); phia{1}=1; phia{d+1}=1;

% Orthogonalization
for i=d:-1:2
    cr = cry{i};
    cr = reshape(cr, ry(i), n(i)*ry(i+1));
    [cr, rv]=qr(cr.', 0);
    cr2 = cry{i-1};
    cr2 = reshape(cr2, ry(i-1)*n(i-1), ry(i));
    cr2 = cr2*(rv.');
    ry(i) = size(cr, 2);
    cr = reshape(cr.', ry(i), n(i), ry(i+1));
%     psy(i) = psy(i+1)-ry(i)*n(i)*ry(i+1);
%     cry(psy(i):(psy(i+1)-1))=cr(:);
%     psyo(i) = psyo(i-1)+ry(i-1)*n(i-1)*ry(i);
%     cry(psyo(i-1):(psyo(i)-1)) = cr2(:);
    cry{i-1} = reshape(cr2, ry(i-1), n(i-1), ry(i));
    cry{i} = cr;

    phia{i} = compute_next_Phi(phia{i+1}, cr, crA{i}, crx{i}, 'rl');
%     phia{i} = compute_next_Phi(phia{i+1}, cr, reshape(cra(psa(i):(psa(i+1)-1)), ra(i), n(i), n(i), ra(i+1)), reshape(crx(psx(i):(psx(i+1)-1)), rx(i), n(i), rx(i+1)), 'rl');

end;
% cry = cry(psy(1):(psy(d+1)-1));

last_sweep = false;
swp = 1;
i = 1;

% dy_old = ones(d,1);
dy = zeros(d,1);
max_dy = 0;
% For extra-rank addition
% dpows = ones(d,1)*min_dpow;
% dranks = zeros(d,1);

cur_order = block_order;
order_index = 1;
dir = sign(cur_order(order_index));

% DMRG sweeps
while (swp<=nswp)

    % Extract elements - matrix
    Phi1 = phia{i}; Phi2 = phia{i+1};
%     A1 = reshape(cra(psa(i):(psa(i+1)-1)), ra(i), n(i), n(i), ra(i+1));
%     x1 = reshape(crx(psx(i):(psx(i+1)-1)), rx(i), n(i), rx(i+1));
    x1 = crx{i};
    A1 = crA{i};

    newy = reshape(Phi1, ry(i)*ra(i), rx(i));
    newy = newy*reshape(x1, rx(i), m(i)*rx(i+1));
    newy = reshape(newy, ry(i), ra(i)*m(i), rx(i+1));
    newy = permute(newy, [2, 1, 3]);
    newy = reshape(newy, ra(i)*m(i), ry(i)*rx(i+1));
    A1 = permute(A1, [2, 4, 1, 3]);
    A1 = reshape(A1, n(i)*ra(i+1), ra(i)*m(i));
    newy = A1*newy;
    newy = reshape(newy, n(i), ra(i+1), ry(i), rx(i+1));
    newy = permute(newy, [3, 1, 2, 4]);
    newy = reshape(newy, ry(i)*n(i), ra(i+1)*rx(i+1));
    if (dir>0)&&(strcmp(kicktype, 'mr'))
        newy_save = newy;
    end;
    newy = newy*(reshape(Phi2, ry(i+1), ra(i+1)*rx(i+1)).');
    newy = reshape(newy, ry(i)*n(i)*ry(i+1), 1);

    oldy = reshape(cry{i}, ry(i)*n(i)*ry(i+1), 1);

    dy(i) = norm(newy-oldy)/norm(newy);
    max_dy = max(max_dy, dy(i));
%
%     % The new core does not converge - increase rank
%     if (dy(i)/dy_old(i)>top_conv)&&(dy(i)>tol)
%         dranks(i)=dranks(i)+1;
%         dpows(i)=dpows(i)+step_dpow;
%     end;
%     % The new core converges well - try to decrease rank
%     if (dy(i)/dy_old(i)<bot_conv)||(dy(i)<tol)
%         dranks(i)=max(dranks(i)-1, 0);
%         dpows(i)=max(dpows(i)-step_dpow, min_dpow);
%     end;
%
%     if (last_sweep)
%         dpows(i)=0.5;
%         dranks(i)=0;
%     end;

    % Truncation
    if (dir>0) % left-to-right
        newy = reshape(newy, ry(i)*n(i), ry(i+1));
    else
        newy = reshape(newy, ry(i), n(i)*ry(i+1));
    end;

    if (kickrank>=0)
    [u,s,v]=svd(newy, 'econ');
    s = diag(s);
    r = my_chop2(s, tol/sqrt(d)*norm(s));
%     % Artificial rank increasing
%     r = r+dranks(i);
    r = min(r, numel(s));
    r = min(r, rmax);
    else
        if (dir>0)
            [u,v]=qr(newy, 0);
            v=v';
            r = size(u,2);
            s = ones(r,1);
        else
            [v,u]=qr(newy.', 0);
            v=conj(v);
            u=u.';
            r = size(u,2);
            s = ones(r,1);
        end;
    end;

    if (verb>1)
	fprintf('=mvk4=   block %d{%d}, dy: %3.3e, r: %d\n', i, dir, dy(i), r);
%     elseif (verb>0)
%         fprintf('                                                                          \r')
%         fprintf('=mvk4=   block %d{%d}, dy: %3.3e, r: %d\r', i, dir, dy(i), r);
    end;

    if (dir>0)&&(i<d) % left-to-right, kickrank, etc
        u = u(:,1:r);
        v = v(:,1:r)*diag(s(1:r));

        % kick
        radd = 0; rv = 1;
        if (~last_sweep)&&(kickrank>0)
%             newy = reshape(Phi1, ry(i)*ra(i), rx(i));
%             newy = newy*reshape(crx{i}, rx(i), n(i)*rx(i+1));
%             newy = reshape(newy, ry(i), ra(i)*n(i), rx(i+1));
%             newy = permute(newy, [2, 1, 3]);
%             newy = reshape(newy, ra(i)*n(i), ry(i)*rx(i+1));
%             A1 = permute(crA{i}, [2, 4, 1, 3]);
%             A1 = reshape(A1, n(i)*ra(i+1), ra(i)*n(i));
%             newy = A1*newy;
%             newy = reshape(newy, n(i), ra(i+1), ry(i), rx(i+1));
%             newy = permute(newy, [3, 1, 2, 4]);
%             newy = reshape(newy, ry(i)*n(i), ra(i+1)*rx(i+1));
            if (strcmp(kicktype, 'mr'))
                leftresid = [reshape(u*v', ry(i)*n(i), ry(i+1)), -newy_save];
                if (strcmp(pcatype, 'svd'))
                    [uk,sk,vk]=svd(leftresid, 'econ');
                    uk = uk(:,1:min(kickrank, size(uk,2)));
                else
                    uk = uchol(leftresid.', kickrank*2);
                    uk = uk(:,size(uk,2):-1:max(size(uk,2)-kickrank+1,1));
                end;

%                 if (i+1<d)
%                     newy = reshape(phia{i+3}, ry(i+3)*ra(i+3), rx(i+3));
%                     newy = newy*reshape(crx{i+2}, rx(i+2)*n(i+2), rx(i+3)).';
%                     newy = reshape(newy, ry(i+3), ra(i+3), rx(i+2), n(i+2));
%                     newy = permute(newy, [2, 4, 1, 3]);
%                     newy = reshape(newy, ra(i+3)*n(i+2), ry(i+3)*rx(i+2));
%                     A1 = permute(crA{i+2}, [2, 1, 4, 3]);
%                     A1 = reshape(A1, n(i+2)*ra(i+2), ra(i+3)*n(i+2));
%                     newy = A1*newy;
%                     newy = reshape(newy, n(i+2), ra(i+2), ry(i+3), rx(i+2));
%                     newy = permute(newy, [1, 3, 2, 4]);
%                     newy = reshape(newy, n(i+2)*ry(i+3), ra(i+2)*rx(i+2));
%
%                     oldy = reshape(cry{i+2}, ry(i+2), n(i+2)*ry(i+3)).';
%                     rightresid = [oldy, -newy];
%                     vk = uchol(rightresid.', kickrank*2);
%                     vk = vk(:,size(vk,2):-1:max(size(vk,2)-kickrank+1,1));
%
%                     [oldy,rv]=qr([oldy,vk], 0);
%                     newy = reshape(cry{i+1}, ry(i+1)*n(i+1), ry(i+2));
%                     newy = [newy, zeros(ry(i+1)*n(i+1), size(vk,2))];
%                     newy = newy*(rv.');
%                     ry(i+2) = size(oldy,2);
%                     cry{i+1} = reshape(newy, ry(i+1), n(i+1), ry(i+2));
%                     cry{i+2} = permute(reshape(oldy, n(i+2), ry(i+3), ry(i+2)), [3, 1, 2]);
%                     phia{i+2} = compute_next_Phi(phia{i+3}, cry{i+2}, crA{i+2}, crx{i+2}, 'rl');
%                 end;
            else
                uk = rand(ry(i)*n(i), kickrank);
            end;
            [u,rv]=qr([u,uk], 0);
            radd = size(uk,2);
        end;
%         radd = size(u, 2)+size(uk,2)-r;
        v = [v, zeros(ry(i+1), radd)];
        v = rv*(v');
        r = size(u,2);
%         r = r+radd;
        cr2 = cry{i+1};
        cr2 = reshape(cr2, ry(i+1), n(i+1)*ry(i+2));
        v = v*cr2; % size r+radd, n2, r3

%         y{i} = reshape(u, ry(i), n(i), r);
%         y{i+1} = reshape(v, r, n(i+1), ry(i+2));
        ry(i+1) = r;

%         % Add residual component
%         rescomp = als_rank_one(x, 'A', A, 'f', y);
%         % We have to update ALL phis
%         Rl = [1,1];
%         for j=1:i
%             addbas = rescomp{i};
%             addbas = addbas/norm(addbas);
%             newbas = zeros(ry(i)+1, n(i), ry(i+1)+1);
%             newbas(1:ry(i), :, 1:ry(i+1))=y{i};
%             newbas(ry(i)+1, :, ry(i+1)+1)=addbas;
%             newbas = reshape(newbas, ry(i)+1, n(i)*(ry(i+1)+1));
%             newbas = Rl*newbas;
%             ry(i) = size(Rl, 1);
%             newbas = reshape(newbas, ry(i)*n(i), ry(i+1)+1);
%             % We have to orthogonalize it. approximately, if j<i
%             proj = newbas(:,ry(i+1)+1) - newbas(:,1:ry(i+1))*(newbas(:,1:ry(i+1))'*newbas(:,ry(i+1)+1));
%             if (i<j)&&(norm(proj)<tol/(d^dpows(i)))
%                 % Remove the new vector
%                 newbas2 = newbas(:, 1:ry(i+1));
%                 Rl = [eye(ry(i+1)); zeros(1,ry(i+1))];
%             else
%                 % Reort the new vector to the old one
%                 newbas2 = reort(newbas(:,1:ry(i+1)), newbas(:, ry(i+1)+1));
%                 Rl = newbas2'*newbas;
%             end;
%             radd = size(newbas2, 2)-ry(i+1);
%             y{i} = reshape(newbas2, ry(i), n(i), ry(i+1)+radd);
%             if (radd>0) % update phi
%                 addbas = reshape(newbas2(:, ry(i+1)+1), ry(i), n(i), 1);
%                 newphia = zeros(ry(i+1)+1, ra(i+1), rx(i+1));
%                 newphia(1:ry(i+1),:,:)=phia{i+1};
%                 newphia(ry(i+1)+1,:,:) = compute_next_Phi(phia{i}, addbas, A{i}, x{i}, 'lr');
%                 phia{i+1} = newphia;
%             end;
%         end;
%
%         Rr = [1;1];
%         for j=d:-1:i+2
%             addbas = rescomp{i};
%             addbas = addbas/norm(addbas);
%             newbas = zeros(ry(i)+1, n(i), ry(i+1)+1);
%             newbas(1:ry(i), :, 1:ry(i+1))=y{i};
%             newbas(ry(i)+1, :, ry(i+1)+1)=addbas;
%             newbas = reshape(newbas, (ry(i)+1)*n(i), (ry(i+1)+1));
%             newbas = newbas*Rr;
%             ry(i+1) = size(Rr, 2);
%             newbas = reshape(newbas, ry(i)+1, n(i)*ry(i+1));
%             newbas = newbas.';
%             % We have to orthogonalize it. approximately, if j<i
%             proj = newbas(:, ry(i)+1) - newbas(:,1:ry(i))*(newbas(:,1:ry(i))'*newbas(:,ry(i)+1));
%             if (i<j)&&(norm(proj)<tol/(d^dpows(i)))
%                 % Remove the new vector
%                 newbas2 = newbas(:, 1:ry(i));
%                 Rr = [eye(ry(i)), zeros(ry(i), 1)];
%             else
%                 % Reort the new vector to the old one
%                 newbas2 = reort(newbas(:,1:ry(i)), newbas(:, ry(i)+1));
%                 Rr = newbas2'*newbas;
%                 Rr = Rr.';
%             end;
%             radd = size(newbas2, 2)-ry(i);
%             y{i} = reshape(newbas2.', ry(i)+radd, n(i), ry(i+1));
%             if (radd>0) % update phi
%                 addbas = reshape(newbas2(:, ry(i)+1).', 1, n(i), ry(i+1));
%                 newphia = zeros(ry(i)+1, ra(i), rx(i));
%                 newphia(1:ry(i),:,:)=phia{i};
%                 newphia(ry(i)+1,:,:) = compute_next_Phi(phia{i+1}, addbas, A{i}, x{i}, 'lr');
%                 phia{i} = newphia;
%             end;
%         end;

        u = reshape(u, ry(i), n(i), r);
        v = reshape(v, r, n(i+1), ry(i+2));

        % Recompute phi. Left ones, so permute them appropriately
        phia{i+1} = compute_next_Phi(phia{i}, u, crA{i}, crx{i}, 'lr');

        % Stuff back

        cry{i} = u;
        cry{i+1} = v;
    elseif (dir<0)&&(i>1) % right-to-left
        u = u(:,1:r)*diag(s(1:r));
        v = conj(v(:,1:r));
        % kick
        radd = 0; rv = 1;
%          if (~last_sweep)&&(strcmp(kicktype, 'rand'))&&(kickrank>0)
%  %             newy = reshape(Phi2, ry(i+1)*ra(i+1), rx(i+1));
%  %             newy = newy*reshape(crx{i}, rx(i)*n(i), rx(i+1)).';
%  %             newy = reshape(newy, ry(i+1), ra(i+1), rx(i), n(i));
%  %             newy = permute(newy, [2, 4, 1, 3]);
%  %             newy = reshape(newy, ra(i+1)*n(i), ry(i+1)*rx(i));
%  %             A1 = permute(crA{i}, [1, 2, 4, 3]);
%  %             A1 = reshape(A1, ra(i)*n(i), ra(i+1)*n(i));
%  %             newy = A1*newy;
%  %             newy = reshape(newy, ra(i), n(i), ry(i+1), rx(i));
%  %             newy = permute(newy, [1, 4, 2, 3]);
%  %             newy = reshape(newy, ra(i)*rx(i), n(i)*ry(i+1));
%  %             resid = [reshape(u*v', ry(i), n(i)*ry(i+1)); -newy];
%  %             uk = uchol(resid, kickrank+1);
%  %             uk = uk(:,end:-1:max(end-kickrank+1,1));
%              uk = randn(n(i)*ry(i+1), kickrank);
%              [v,rv]=qr([v,uk], 0);
%              radd = size(uk,2);
%  %             v = reort(v, uk);
%          end;

%  %         radd = size(v, 2)+size(uk,2)-r;
        u = [u, zeros(ry(i), radd)];
        u = u*(rv.');
        cr2 = cry{i-1};
        cr2 = reshape(cr2, ry(i-1)*n(i-1), ry(i));
        u = cr2*u;

%         r = r+radd;
        r = size(v,2);

        u = reshape(u, ry(i-1), n(i-1), r);
        v = reshape(v.', r, n(i), ry(i+1));

        % Recompute phi. Here are right phis
        phia{i} = compute_next_Phi(phia{i+1}, v, crA{i}, crx{i}, 'rl');

        % Stuff back
        ry(i) = r;
        cry{i-1} = u;
        cry{i} = v;
    elseif ((dir>0)&&(i==d))||((dir<0)&&(i==1))
        % Just stuff back the last core
        newy = u(:,1:r)*diag(s(1:r))*(v(:,1:r)');
        newy = reshape(newy, ry(i), n(i), ry(i+1));
        cry{i} = newy;
    end;

    i = i+dir;

    % Reversing, residue check, etc
    cur_order(order_index) = cur_order(order_index) - dir;
    % New direction
    if (cur_order(order_index)==0)
        order_index = order_index+1;

        if (verb>0)
%  	    fprintf('                                                   \r')
            fprintf('=mvk4= sweep %d{%d}, max_dy: %3.3e, erank: %g\n', swp, order_index-1, max_dy, sqrt(ry(1:d)'*(n.*ry(2:d+1))/sum(n)));
        end;

        if (last_sweep)
            break;
        end;

        if (kickrank<0)
            kickrank=kickrank-1;
        end;

        if (max_dy<tol*als_tol_low)&&(kickrank>=0)
            kickrank=-1;
        end;

        if (max_dy<tol)&&(kickrank<=-als_iters)
            last_sweep=true;
            kickrank=0;
        end;

        if (order_index>numel(cur_order)) % New global sweep
            cur_order = block_order;
            order_index = 1;
            %residue
            if (last_sweep)
                cur_order = d-1;
            end;

            max_dy = 0;

%             dy_old = dy;
            swp = swp+1;
        end;

        dir = sign(cur_order(order_index));
        i = i+dir;
    end;
end

if (tailranks(1))
    cry{2} = reshape(cry{1}, n(1), ry(2))*reshape(cry{2}, ry(2), n(2)*ry(3));
    cry{2} = reshape(cry{2}, n(1), n(2), ry(3));
end;
if (tailranks(2))
    cry{d-1} = reshape(cry{d-1}, ry(d-1)*n(d-1), ry(d))*reshape(cry{d}, ry(d), n(d));
    cry{d-1} = reshape(cry{d-1}, ry(d-1), n(d-1), n(d));
    cry = cry(1:d-1);
end;
if (tailranks(1))
    cry = cry(2:end);
end;
y = cell2core(y, cry);

end

function [Phi] = compute_next_Phi(Phi_prev, x, A, y, direction)
% Performs the recurrent Phi (or Psi) matrix computation
% Phi = Phi_prev * (x'Ay).
% If direction is 'lr', computes Psi
% if direction is 'rl', computes Phi
% A can be empty, then only x'y is computed.

if (strcmp(direction, 'rl'))
  % Revert ranks to perform the right-to-left recursion
  x = permute(x, [3, 2, 1]);
  y = permute(y, [3, 2, 1]);
  if (~isempty(A))
    A = permute(A, [4, 2, 3, 1]);
  end
end

rx1 = size(x,1); n = size(x,2); rx2 = size(x,3);
ry1 = size(y,1); m = size(y,2); ry2 = size(y,3);
if (~isempty(A))
  ra1 = size(A,1); ra2 = size(A,4);
else
  ra1 = 1; ra2 = 1;
end

Phi = reshape(Phi_prev, [rx1*ra1, ry1]);
y = reshape(y, [ry1, m*ry2]);
Phi = Phi*y;	% complexity §\mcommentfont$\mathcal{O}(n  r_x r_A r_y^2)$§
Phi = reshape(Phi, [rx1, ra1, m, ry2]);
Phi = permute(Phi, [2, 3, 1, 4]);
if (~isempty(A))
  Phi = reshape(Phi, [ra1*m, rx1*ry2]);
  A = permute(A, [4, 2, 1, 3]);
  A = reshape(A, [ra2*n, ra1*m]);
  Phi = A*Phi;	% complexity §\mcommentfont$\mathcal{O}(n^2  r_x r_A^2 r_y)$§
  Phi = reshape(Phi, [ra2, n, rx1, ry2]);
end
Phi = permute(Phi, [3, 2, 1, 4]);
Phi = reshape(Phi, [rx1*n, ra2*ry2]);
x = reshape(x, [rx1*n, rx2]);
Phi = (x')*Phi;	% complexity §\mcommentfont$\mathcal{O}(n  r_x^2 r_A r_y)$§
if (~isempty(A))
  Phi = reshape(Phi, [rx2, ra2, ry2]);
end
end

author: oseledets

repo: TT-Toolbox


text:
% MCMLT - makes matrix of MCMC runs lower triangular
% Copyright (c) 1998, Harvard University. Full copyright in the file Copyright
%
%   [ Alt ] = MCMCLT(A)
%
% A = a chain of matricies, typically covariances
%
% Alt = A with all elements above the diagonal changed to NaNs
%
% Good for removing redundant trace plots and values,
% especially for samples of covariance matricies.
%
% MCMCTRACE will not plot chains that begin with NaNs.
%
% See also: LTVEC, VECLT

function [Alt] = mcmclt(A)

dd = size(A) ;
ll = length(dd) ;
d1 = dd(1) ;
d2 = dd(2) ;

if (ll==2),
  Alt = A ;
  for i1 = 1:d1,
  for i2 = 1:d2,
    if i2>i1,
      Alt(i1,i2)=NaN ;
    end
  end
  end
else
  Alt = A ;
  for i1 = 1:d1,
  for i2 = 1:d2,
    if i2>i1,
      Alt(i1,i2,:)=NaN ;
    end
  end
  end
end

author: Sable

repo: mcbench-benchmarks


text:
function [rVect, vVect] = getStateAtTime_alg(time, sma, ecc, inc, raan, arg, mean, epoch, gmu)
%getStateAtTime Summary of this function goes here
%   Detailed explanation goes here
    numTimes = length(time);

    oneArray = (zeros(1, numTimes)+1);

    sma = sma * oneArray;
    ecc = ecc * oneArray;
    inc = AngleZero2Pi(deg2rad(inc)) * oneArray;
    raan = AngleZero2Pi(deg2rad(raan)) * oneArray;
    arg = AngleZero2Pi(deg2rad(arg)) * oneArray;
    M0 = deg2rad(mean) * oneArray;

    n = computeMeanMotion(sma, gmu);
    deltaT = time - epoch;
    M = (M0(:) + n(:).*deltaT(:))';
    tru = computeTrueAnomFromMean(M, ecc);

    if(length(tru) > 1)
        [rVect,vVect] = vect_getStatefromKepler(sma, ecc, inc, raan, arg, tru, gmu, false);
    else
        [rVect, vVect] = getStatefromKepler_Alg(sma, ecc, inc, raan, arg, tru, gmu);
    end

    rVect(isnan(rVect)) = 0;
    vVect(isnan(vVect)) = 0;
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function vn = normalizeVector3d(v)
%NORMALIZEVECTOR3D Normalize a 3D vector to have norm equal to 1.
%
%   V2 = normalizeVector3d(V);
%   Returns the normalization of vector V, such that ||V|| = 1. Vector V is
%   given as a row vector.
%
%   If V is a N-by-3 array, normalization is performed for each row of the
%   input array.
%
%   See also:
%   vectors3d, vectorNorm3d
%

%   ---------
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 29/11/2004.
%

% HISTORY
% 2005-11-30 correct a bug
% 2009-06-19 rename as normalizeVector3d
% 2010-11-16 use bsxfun (Thanks to Sven Holcombe)

vn   = bsxfun(@rdivide, v, sqrt(sum(v.^2, 2)));


text:
function [ projdata ] = fiff_read_proj(fid,node)

%
% [ projdata ] = fiff_read_proj(fid,node)
%
% Read the SSP data under a given directory node
%

%
%   Author : Matti Hamalainen, MGH Martinos Center
%   License : BSD 3-clause
%
%
%   Revision 1.7  2006/05/16 00:39:32  msh
%   Fixed error in mne_ex_read_raw: All projection items were not activated.
%   List the initial states of projection items when they are loaded.
%
%   Revision 1.6  2006/04/23 15:29:40  msh
%   Added MGH to the copyright
%
%   Revision 1.5  2006/04/21 15:06:03  msh
%   Do not list projection item status because it might be confusing.
%
%   Revision 1.4  2006/04/21 14:43:59  msh
%   Report projection operators.
%   Some more checks in raw data reading.
%
%   Revision 1.3  2006/04/18 20:44:46  msh
%   Added reading of forward solution.
%   Use length instead of size when appropriate
%
%   Revision 1.2  2006/04/12 10:29:02  msh
%   Made evoked data writing compatible with the structures returned in reading.
%
%   Revision 1.1  2006/04/10 23:26:54  msh
%   Added fiff reading routines
%
%

global FIFF;
if isempty(FIFF)
    FIFF = fiff_define_constants();
end

me='MNE:fiff_read_proj';

if nargin ~= 2
    error(me,'Incorrect number of arguments');
end

projdata = struct('kind',{},'active',{},'desc',{},'data',{});
%
%   Locate the projection data
%
nodes = fiff_dir_tree_find(node,FIFF.FIFFB_PROJ);
if length(nodes) == 0
    return;
end
tag = find_tag(nodes(1),FIFF.FIFF_NCHAN);
if ~isempty(tag)
    global_nchan = tag.data;
end
items = fiff_dir_tree_find(nodes(1),FIFF.FIFFB_PROJ_ITEM);
for i = 1:length(items)
    %
    %   Find all desired tags in one item
    %
    item = items(i);
    tag = find_tag(item,FIFF.FIFF_NCHAN);
    if ~isempty(tag)
        nchan = tag.data;
    else
        nchan = global_nchan;
    end
    tag = find_tag(item,FIFF.FIFF_DESCRIPTION);
    if ~isempty(tag)
        desc = tag.data;
    else
        tag = find_tag(item,FIFF.FIFF_NAME);
        if ~isempty(tag)
            desc = tag.data;
        else
            error(me,'Projection item description missing');
        end
    end
    tag = find_tag(item,FIFF.FIFF_PROJ_ITEM_CH_NAME_LIST);
    if ~isempty(tag)
        namelist = tag.data;
    else
        error(me,'Projection item channel list missing');
    end
    tag = find_tag(item,FIFF.FIFF_PROJ_ITEM_KIND);
    if ~isempty(tag)
        kind = tag.data;
    else
        error(me,'Projection item kind missing');
    end
    tag = find_tag(item,FIFF.FIFF_PROJ_ITEM_NVEC);
    if ~isempty(tag)
        nvec = tag.data;
    else
        error(me,'Number of projection vectors not specified');
    end
    tag = find_tag(item,FIFF.FIFF_PROJ_ITEM_CH_NAME_LIST);
    if ~isempty(tag)
        names = fiff_split_name_list(tag.data);
    else
        error(me,'Projection item channel list missing');
    end
    tag = find_tag(item,FIFF.FIFF_PROJ_ITEM_VECTORS);
    if ~isempty(tag)
        data = tag.data;
    else
        error(me,'Projection item data missing');
    end
    tag = find_tag(item,FIFF.FIFF_MNE_PROJ_ITEM_ACTIVE);
    if ~isempty(tag)
        active = tag.data;
    else
        active = false;
    end
    if size(data,2) ~= length(names)
        error(me,'Number of channel names does not match the size of data matrix');
    end
    one.kind           = kind;
    one.active         = active;
    one.desc           = desc;
    %
    %   Use exactly the same fields in data as in a named matrix
    %
    one.data.nrow      = nvec;
    one.data.ncol      = nchan;
    one.data.row_names = [];
    one.data.col_names = names;
    one.data.data      = data;
    %
    projdata(i) = one;
end

if length(projdata) > 0
    fprintf(1,'\tRead a total of %d projection items:\n', ...
        length(projdata));
    for k = 1:length(projdata)
        fprintf(1,'\t\t%s (%d x %d)',projdata(k).desc, ...
            projdata(k).data.nrow,projdata(k).data.ncol);
        if projdata(k).active
            fprintf(1,' active\n');
        else
            fprintf(1,' idle\n');
        end
    end
end

return;

    function [tag] = find_tag(node,findkind)

        for p = 1:node.nent
            if node.dir(p).kind == findkind
                tag = fiff_read_tag(fid,node.dir(p).pos);
                return;
            end
        end
        tag = [];
    end

end

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function out = isempty(v)
% overloads isempty

out = isempty(v.x);

author: mtex-toolbox

repo: mtex

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


text:
function pass = test_coeffs
% TAD, 10 Jan 2014

%% Building blocks
dom = [ - 2  - 0.5 1 2];
I = operatorBlock.eye(dom);
D = operatorBlock.diff(dom);
Z = operatorBlock.zeros(dom);
x = chebfun('x', dom);
u = sin(x.^2);
U = operatorBlock.mult(u);

%% Coefficient realizations
% Quasimatrix [ 1 ]
Ic = toCoeff(I);
err(1) = norm(Ic{1} - x.^0);
% Quasimatrix [ 1 0 ]
Dc = toCoeff(D);
err(2) = norm(Dc{1} - x.^0);
err(3) = norm(Dc{2} - 0);
% Quasimatrix [ x.^2 ]
Uc = toCoeff(U);
err(4) = norm(Uc{1} - u);

%%
A = D*(I+U*D);
Ac = toCoeff(A);
err(5) = norm(Ac{1} - u);
err(6) = norm(Ac{2} - (diff(u)+1));
err(7) = norm(Ac{3} - 0);

pass = err < 1e-15;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function m = curveMoment(curve, p, q)
%CURVEMOMENT  Compute inertia moment of a 2D curve
%   M = curveMoment(CURVE, P, Q)
%
%   Example
%   curveMoment
%
%   See also
%   polygons2d, curveCMoment, curveCSMoment
%
%   Reference
%   Based on ideas and references in:
%   "Affine curve moment invariants for shape recognition"
%   Dongmin Zhao and Jie Chen
%   Pattern Recognition, 1997, vol. 30, pp. 865-901
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2009-03-25,    using Matlab 7.7.0.471 (R2008b)
% Copyright 2009 INRA - Cepia Software Platform.

% coordinate of vertices
px  = curve(:,1);
py  = curve(:,2);

% compute centroids of line segments
cx  = (px(1:end-1)+px(2:end))/2;
cy  = (py(1:end-1)+py(2:end))/2;

% compute length of each line segment
dl  = hypot(px(2:end)-px(1:end-1), py(2:end)-py(1:end-1));

% compute moment
m = zeros(size(p));
for i=1:length(p(:))
    m(i) = sum(cx(:).^p(i) .* cy(:).^q(i) .* dl(:));
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function c = tapas_softmax_config
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Contains the configuration for the softmax observation model for multinomial responses
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2013 Christoph Mathys, TNU, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

% Config structure
c = struct;

% Is the decision based on predictions or posteriors? Comment as appropriate.
c.predorpost = 1; % Predictions
%c.predorpost = 2; % Posteriors

% Model name
c.model = 'softmax';

% Sufficient statistics of Gaussian parameter priors

% Beta
c.logbemu = log(1);
c.logbesa = 4^2;

% Gather prior settings in vectors
c.priormus = [
    c.logbemu,...
         ];

c.priorsas = [
    c.logbesa,...
         ];

% Model filehandle
c.obs_fun = @tapas_softmax;

% Handle to function that transforms observation parameters to their native space
% from the space they are estimated in
c.transp_obs_fun = @tapas_softmax_transp;

return;

author: translationalneuromodeling

repo: tapas


text:
function wt = wfbtput(d,k,w,wt,forceStr)
%WFBTPUT  Put node to the filterbank tree
%   Usage:  wt = wfbtput(d,k,w,wt);
%           wt = wfbtput(d,k,w,wt,'force');
%
%   Input parameters:
%           d   : Level in the tree (0 - root).
%           k   : Index (array of indexes) of the node at level *d* (starting at 0).
%           w   : Node, basic wavelet filterbank.
%           wt  : Wavelet filterbank tree structure (as returned from
%                 |wfbtinit|).
%
%   Output parameters:
%           wt : Modified filterbank structure.
%
%   `wfbtput(d,k,w,wt)` puts the basic filterbank *w* to the filter
%   tree structure *wt* at level *d* and index(es) *k*. The output is a
%   modified tree structure. *d* and *k* have to specify unconnected output
%   of the leaf node. Error is issued if *d* and *k* points to already
%   existing node. For possible formats of parameter *w* see help of |fwt|.
%   Parameter *wt* has to be a structure returned by |wfbtinit|.
%
%   `wfbtput(d,k,w,wt,'force')` does the same but replaces node at *d* and *k*
%   if it already exists. If the node to be replaced has any children,
%   the number of outputs of the replacing node have to be equal to number of
%   outputs of the node beeing replaced.
%
%   Examples:
%   ---------
%
%   This example shows magnitude frequency responses of a tree build from
%   the root:::
%
%      % Initialize empty struct
%      wt = wfbtinit();
%      % Put root node to the empty struct
%      wt1 = wfbtput(0,0,'db8',wt);
%      % Connect a different nodes to both outputs of the root
%      wt2 = wfbtput(1,[0,1],'db10',wt1);
%      % Connect another nodes just to high-pass outputs of nodes just added
%      wt3 = wfbtput(2,[1,3],'db10',wt2);
%      % Add another node at level 3
%      wt4 = wfbtput(3,1,'db16',wt3);
%
%      % Create identical filterbanks
%      [g1,a1] = wfbt2filterbank(wt1,'freq');
%      [g2,a2] = wfbt2filterbank(wt2,'freq');
%      [g3,a3] = wfbt2filterbank(wt3,'freq');
%      [g4,a4] = wfbt2filterbank(wt4,'freq');
%
%      % Plot frequency responses of the growing tree. Linear scale
%      % (both axis) is used and positive frequencies only are shown.
%      subplot(4,1,1);
%      filterbankfreqz(g1,a1,1024,'plot','linabs','posfreq');
%      subplot(4,1,2);
%      filterbankfreqz(g2,a2,1024,'plot','linabs','posfreq');
%      subplot(4,1,3);
%      filterbankfreqz(g3,a3,1024,'plot','linabs','posfreq');
%      subplot(4,1,4);
%      filterbankfreqz(g4,a4,1024,'plot','linabs','posfreq');
%

% AUTHOR: Zdenek Prusa

if nargin<4
   error('%s: Too few input parameters.',upper(mfilename));
end

%if isfield(wt,'dualnodes')
%    error('%s: Cannot modify the dual-tree struct.',upper(mfilename));
%end

do_force = 0;
if nargin==5
    if ~ischar(forceStr)
        error('%s: Fifth parameter should be a string.',upper(mfilename));
    end
    if strcmpi(forceStr,'force')
        do_force = 1;
    end
end

% This was replaced. Calling ltfatargheler was too slow.
%definput.flags.force = {'noforce','force'};
%[flags,kv]=ltfatarghelper({},definput,varargin);

node = fwtinit(w);

oldnodecount = numel(wt.nodes);
nodeschanged = [];

[nodeNoArray,nodeChildIdxArray] = depthIndex2NodeNo(d,k,wt);

for ii=1:numel(nodeNoArray)
 nodeNo = nodeNoArray(ii);
 nodeChildIdx = nodeChildIdxArray(ii);
if(nodeNo==0)
    % adding root
    if(~isempty(find(wt.parents==0,1)))
        if(do_force)
           rootId = find(wt.parents==0,1);
           % if root has children, check if the new root has the same
           % number of them
           if(~isempty(find(wt.children{rootId}~=0,1)))
              if(length(w.g)~=length(wt.nodes{rootId}.g))
                 error('%s: The replacing root have to have %d filters.',mfilename,length(wt.nodes{rootId}.g));
              end
           end
        else
            error('%s: Root already defined. Use FORCE option to replace.',mfilename);
        end
        wt.nodes{rootId} = node;
        nodeschanged(end+1) = rootId;

        if isfield(wt,'dualnodes')
            wt.dualnodes{rootId} = node;
        end
        continue;
    end
    wt.nodes{end+1} = node;
    wt.parents(end+1) = nodeNo;
    wt.children{end+1} = [];

    if isfield(wt,'dualnodes')
        wt.dualnodes{end+1} = node;
    end
    continue;
end

childrenIdx = find(wt.children{nodeNo}~=0);
found = find(childrenIdx==nodeChildIdx,1);
if(~isempty(found))
   if(do_force)
     %check if childrenIdx has any children
     tmpnode = wt.children{nodeNo}(found);
     if(~isempty(find(wt.children{tmpnode}~=0, 1)))
         if length(node.g)~=length(wt.nodes{tmpnode}.g)
            error('%s: The replacing node must have %d filters.',mfilename,length(wt.nodes{tmpnode}.g));
         end
     end
     wt.nodes{tmpnode} = node;
     nodeschanged(end+1) = tmpnode;
     if isfield(wt,'dualnodes')
         wt.dualnodes{tmpnode} = node;
     end
     % Since we are replacing a node, all links are already correct
     continue;
   else
       error('%s: Such node (depth=%d, idx=%d) already exists. Use FORCE option to replace.',mfilename,d,k);
   end
end

wt.nodes{end+1} = node;
wt.parents(end+1) = nodeNo;
wt.children{end+1} = [];
wt.children{nodeNo}(nodeChildIdx) = numel(wt.parents);

if isfield(wt,'dualnodes')
    wt.dualnodes{end+1} = node;
end

end

% We have to correctly shuffle filters in the just added (or modified) filters
% if the tree was already defined as frequency ordered.
if wt.freqOrder
      wt = nat2freqOrder(wt,[nodeschanged,oldnodecount+1:numel(wt.nodes)]);
end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function slsharedisp_decindent(nsteps)
%SLSHAREDISP_DECINDENT Decreases the indent of the displayer
%
% $ Syntax $
%   - slsharedisp_decindent()
%   - slsharedisp_decindent(nsteps)
%
% $ Description $
%   - slsharedisp_decindent() decreases the indent by one step.
%
%   - slsharedisp_decindent(nsteps) decreases the indent by specified
%     number of steps.
%
% $ History $
%   - Created by Dahua Lin, on Aug 29, 2006
%

global GLOBAL_SHARE_DISPLAYER;

if isempty(GLOBAL_SHARE_DISPLAYER)
    error('sltoolbox:gdisperr', ...
        'The global displayer is not open');
end

if nargin == 0
    nsteps = 1;
end

GLOBAL_SHARE_DISPLAYER(end).indent = ...
    GLOBAL_SHARE_DISPLAYER(end).indent - nsteps;

author: lmthang


text:
% nufft_table_test.m

%
% 3D test
%
if 1
	if ~isvar('s.t'), disp 'precompute structure'
%		ktype = 'minmax:kb';
		ktype = 'kaiser';
		Jd = [3 4 3];
		Nd = [10 8 9];
		Ld = [2^8 2^7 2^7];
		Kd = 2 * Nd;
%		n_shift = zeros(size(Nd)); disp 'easy: 0 shift'
		n_shift = [3 7 2]; % stress it

		gam = 2*pi ./ Kd;
		k1 = linspace(-2*Kd(1), 2*Kd(1), 21)';
		k2 = linspace(-2*Kd(2), 2*Kd(2), 18)';
		k3 = linspace(-2*Kd(3), 2*Kd(3), 23)';
		[kk1 kk2 kk3] = ndgrid(k1, k2, k3);
		om = [gam(1)*kk1(:) gam(2)*kk2(:) gam(3)*kk3(:)];

		tic
		s.p = nufft_init(om, Nd, Jd, Kd, n_shift, ktype);
		printf('pre time init %g', toc)

		tic
		s.t = nufft_init(om, Nd, Jd, Kd, n_shift, 'table', Ld, ktype);
		printf('tab time init %g', toc)
	end

	% 3D forward direction
	if 1, disp '3D forward'
%		x = [1:Nd(1)]'*triang(Nd(2))';
		rand('state', 0);
		x = rand(Nd);
		x = zeros(Nd);
		x(4, 2, 3) = 1;

		tic
		Y.d = dtft(x, om, n_shift); % exact
		printf('dtft3 time %g', toc)

		tic
		Y.p = nufft(x, s.p);
		printf('pre time nufft %g', toc)
		printf('pre max%%diff = %g', max_percent_diff(Y.d, Y.p))

		tic
		Y.t = nufft(x, s.t);
		printf('tab time nufft %g', toc)
		printf('tab max%%diff = %g', max_percent_diff(Y.d, Y.t))
		printf('tab vs pre max%%diff = %g', max_percent_diff(Y.p, Y.t))

		if 1
			clf
			im(231, k1, k2, abs(Y.p - Y.d), 'pre'), cbar
			im(232, k1, k2, abs(Y.t - Y.d), 'tab'), cbar
			im(233, k1, k2, abs(Y.t - Y.p), 'tab-pre'), cbar
		end
	end

	% 3D adjoint
	if 1, disp '3D adjoint'
		X = [1:s.t.M]' + 1i * ones(s.t.M,1);

		tic
		c.d = dtft_adj(X, om, Nd, n_shift); % exact
		printf('dtft time adj %g', toc)
		tic
		c.p = nufft_adj(X, s.p);
		printf('pre time nufft adj %g', toc)
		tic
		c.t = nufft_adj(X, s.t);
		printf('tab time nufft adj %g', toc)

		printf('pre max%%diff = %g', max_percent_diff(c.d, c.p))
		printf('tab max%%diff = %g', max_percent_diff(c.d, c.t))
		printf('tab vs pre max%%diff = %g', max_percent_diff(c.p, c.t))

		if 1
			im(234, abs(c.p - c.d), 'pre adj err'), cbar
			im(235, abs(c.t - c.d), 'tab adj err'), cbar
			im(236, abs(c.t - c.p), 'tab-pre adj'), cbar
		end
	end
return
end

%
% 2D test
%
if 1
	if ~isvar('s.t'), disp 'precompute structure'
%		ktype = 'minmax:kb';
		ktype = 'kaiser';
		Jd = [4 5];
		Nd = [20 16];
		Ld = [2^12 2^13];
		Kd = 2 * Nd;
%		n_shift = zeros(size(Nd));
		n_shift = [3 7]; % stress it

		gam = 2*pi ./ Kd;
		k1 = linspace(-2*Kd(1), 2*Kd(1), 51)';
		k2 = linspace(-2*Kd(2), 2*Kd(2), 81)';
		[kk1 kk2] = ndgrid(k1, k2);
		om = [gam(1)*kk1(:) gam(2)*kk2(:)];

		tic
		s.p = nufft_init(om, Nd, Jd, Kd, n_shift, ktype);
		printf('pre time init %g', toc)

		tic
		s.t = nufft_init(om, Nd, Jd, Kd, n_shift, 'table', Ld, ktype);
		printf('tab time init %g', toc)
	end

	% forward direction
	if 1, disp 'test forward'
		x = [1:Nd(1)]'*triang(Nd(2))';

		tic
		Y.d = dtft2(x, om, n_shift); % exact
		printf('dtft2 time %g', toc)

		tic
		Y.p = nufft(x, s.p);
		printf('pre time nufft %g', toc)
		printf('pre max%%diff = %g', max_percent_diff(Y.d, Y.p))

		tic
		Y.t = nufft(x, s.t);
		printf('tab time nufft %g', toc)
		printf('tab max%%diff = %g', max_percent_diff(Y.d, Y.t))
		printf('tab vs pre max%%diff = %g', max_percent_diff(Y.p, Y.t))

		clf
		im(231, k1, k2, abs(Y.p - Y.d), 'pre'), cbar
		im(232, k1, k2, abs(Y.t - Y.d), 'tab'), cbar
		im(233, k1, k2, abs(Y.t - Y.p), 'tab-pre'), cbar
	end

	% 2D adjoint
	if 1, disp '2D adjoint'
		X = [1:s.t.M]' + 1i * ones(s.t.M,1);

		tic
		c.d = dtft2_adj(X, om, Nd(1), Nd(2), n_shift); % exact
		printf('dtft2 time adj %g', toc)
		tic
		c.p = nufft_adj(X, s.p);
		printf('pre time nufft adj %g', toc)
		tic
		c.t = nufft_adj(X, s.t);
		printf('tab time nufft adj %g', toc)

		printf('pre max%%diff = %g', max_percent_diff(c.d, c.p))
		printf('tab max%%diff = %g', max_percent_diff(c.d, c.t))
		printf('tab vs pre max%%diff = %g', max_percent_diff(c.p, c.t))

		im(234, abs(c.p - c.d), 'pre adj err'), cbar
		im(235, abs(c.t - c.d), 'tab adj err'), cbar
		im(236, abs(c.t - c.p), 'tab-pre adj'), cbar
	end
return
end

%
% 1D test
%
if ~isvar('s.t')
%	ktype = 'minmax:kb';
	ktype = 'kaiser';
	Jd = 4;
	Nd = 32;
	Ld = 2^14; % table over-sampling
	Kd = 2 * Nd;
	n_shift = zeros(size(Nd));
	% n_shift = Nd/2;

	gam = 2*pi ./ Kd;
	kv = 2*Kd;
	kv = linspace(-kv, kv, 1001)';
	om = gam * kv;

	tic
	s.p = nufft_init(om, Nd, Jd, Kd, n_shift, ktype);
	printf('pre time init %g', toc)

	tic
	s.t = nufft_init(om, Nd, Jd, Kd, n_shift, 'table', Ld, ktype);
	printf('tab time init %g', toc)
end

if 0
	ramp = [1:Nd]';

	Y.dr = dtft(ramp, om, n_shift);
	Y.de = dtft(eye(Nd), om, n_shift);
	printf('dtft max%%diff = %g', max_percent_diff(Y.dr,Y.de*ramp))

	Y.pr = nufft(ramp, s.p);
	Y.pe = nufft(eye(Nd), s.p);
	printf('pre max%%diff = %g', max_percent_diff(Y.pr,Y.pe*ramp))

	printf('max%%diff = %g', max_percent_diff(Y.dr,Y.pr))
	printf('max%%diff = %g', max_percent_diff(Y.de,Y.pe))

	e = Y.pr - Y.dr;
	clf, plot(kv, abs(e), 'c-')
%	clf, plot(kv, sum(abs(e),2), 'c-')
return
end

% forward direction
if 1
	x = [1:Nd]';
	% x = unitv(Nd, ii);

	Y.d = dtft(x, om, n_shift); % exact
	tic
	Y.p = nufft(x, s.p);
	printf('pre time nufft %g', toc)
	printf('pre max%%diff = %g', max_percent_diff(Y.d,Y.p))

	tic
	Y.t = nufft(x, s.t);
	printf('tab time nufft %g', toc)
	printf('tab max%%diff = %g', max_percent_diff(Y.d, Y.t))

	clf
	subplot(211)
	plot(kv, abs(Y.p - Y.d), 'c-'), title 'pre'
	subplot(212)
	plot(kv, abs(Y.t - Y.d), 'y-'), title 'tab'
prompt
end

% adjoint
if 1
	X = [1:s.t.M]' + 1i * ones(size(om));

	c.d = dtft2_adj(X, [om 0*om], Nd(1), 1, [n_shift 0]); % exact
	tic
	c.p = nufft_adj(X, s.p);
	printf('pre time nufft adj %g', toc)
	tic
	c.t = nufft_adj(X, s.t);
	printf('tab time nufft adj %g', toc)

	printf('pre nufft adj vs dtft adj max%%diff = %g', max_percent_diff(c.d,c.p))
	printf('tab nufft adj vs dtft adj max%%diff = %g', max_percent_diff(c.d,c.t))

	clf
	nn = 0:Nd(1)-1;
	subplot(211)
	plot(nn, abs(c.p - c.d) / mean(abs(c.d)), 'c-'), title 'pre'
	subplot(212)
	plot(nn, abs(c.t - c.d) / mean(abs(c.d)), 'y-'), title 'tab'
end

author: thomaskuestner


text:
function [data] = avgref(data, sel)

% AVGREF computes the average reference in each column
%   [data] = avgref(data)
%
% or it computes the re-referenced data relative to the
% average over the selected channels
%   [data] = avgref(data, sel)

% Copyright (C) 1998-2002, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

% determine the dimension of the data
if length(size(data))==3
  % multiple epochs
  dim=3;
else
  % single epoch with multiple channels
  dim=2;
end

if nargin==1
  % default is to use all channels for average referencing
  if dim==3
    sel=1:size(data,2);
  else
    sel=1:size(data,1);
  end
end

if dim==3
  % the data contains multiple epochs
  for epoch=1:size(data,1)
    reference = mean(data(epoch,sel,:), 2);
    data(epoch,:,:) = data(epoch,:,:) - repmat(reference, [1 size(data,2) 1]);
  end
else
  % the data contains a single epoch
  reference = mean(data(sel,:), 1);
  data = data - repmat(reference, size(data,1), 1);
end


text:
function [Model, Info] = linear_sparse_space_sw(X,Y,Model,parm)
%  Estimate linear weight matrix for input-output mapping
%     Automatic Relevance Prior for each input dimension
%     in original input space (not for delay embedding space)
%     is imposed to get sparse weight matrix
%  Trial data with variable sample number is supported
%     by setting sample number of each trial
%  Notice !!
%     Delay embedding is done in this module,
%     then input vector should be original input data without embedding
%
%     Covariance is calculated in original input space,
%     then increse of embedding dimension make no problem in computation
%
%   [Model, Info] = linear_sparse_space_sw(X,Y,Model,parm)
% --- Input
%  Y  : Output data ( N x T x Ntrial )
%  X  : Input data  ( M x (T + (Dtau-1)*Tau)  x Ntrial)
%  N  =  # of output
%  M  =  # of input (original input space dimension)
%  T  =  # of time sample
%
%  Estimate the following model
%    Y(t) = W * [X(:, t + (Dtau-1)*Tau ); ...; X(:, t)]
%
%  Model : Structure for estimated model
%  if Model is empty, initialization is done before training
%  if Model is previous training result, re-training us done
%
%  parm  : Structure for learning parameter
%  parm.Npre_train :  # of VB-update in initial training
%  parm.Ntrain :  # of training
%  parm.Nskip  :  skip # for print
%  parm.a_min  :  Min value for pruning small variance component
%  parm.Prune  :  = 1 : Prune small variance & irrelevant input dimension
%
%  parm.Tau       = Lag time
%  parm.Dtau      = Number of embedding dimension
%    Total input dimension in embedding space is (M * parm.Dtau)
%  parm.Trial : Valid_sample number for each trial [Ntrial x 1]
%    if parm.Trial = [15 10] and T = 15,
%       last 5 samples in 2nd trial is not used for estimation
% --- Output
%  Model : Structure for estimated model
%  Model.SY  :  Noise variance         ( 1 x 1 )
%  Model.W   :  Weight matrix          ( N x M*D ) , D = parm.Dtau
%  Model.A   :  Prior weight variance  ( 1 x M*D )
%  Model.ix_act : Active index for W after pruning
%
%  Info  : Structure for learning process history
%  Info.FE  = LP + H : Free energy
%  Info.LP  = Log likelihood
%  Info.H   = - Model entropy
%
% 2009-10-18 Made by M. Sato

MINVAL = 1.0e-15;
MinCond = 1e-8;

fprintf('Linear sparse space (variable trial length) start\n')

% Dimension
[N ,T  ,Ntrial]  = size(Y); % N = # of output
[M ,Tx ,Ntrialx] = size(X); % M = # of input without embedding

% Number of samples for each trial
if isfield(parm,'Trial'),
	Trial = parm.Trial;
	if length(Trial) ~= Ntrial, error('Trial sample seting is wrong');end
	if max(Trial) > T, error('Trial sample seting is wrong');end
else
	Trial = repmat(T, [Ntrial 1]);
end
% Set valid sample flag for each trial
trial_sw = zeros(N,T,Ntrial);
for n = 1:Ntrial
	trial_sw(:,1:Trial(n),n) = 1;
end

Tall = sum(Trial);

if Ntrial~=Ntrialx, error('# of trial is different in X and Y'); end;

% # of total training iteration
Ntrain = parm.Ntrain;

Nskip  = 100;   % skip steps for display info
a_min  = 1e-10; % Minimum value for weight pruning
Fdiff  = 1e-10; % Threshold for convergence
Ncheck = 100;   % Minimum number of training iteration
Fstep  = 5;     % Free energy convergence check step
Prune  = 1;     % Prune mode

if isfield(parm,'Nskip'), Nskip  = parm.Nskip; end;
if isfield(parm,'Fdiff'), Fdiff   = parm.Fdiff; end;
if isfield(parm,'a_min'), a_min   = parm.a_min ; end;
if isfield(parm,'Prune'), Prune = parm.Prune; end;
if isfield(parm,'Ncheck'), Ncheck = parm.Ncheck; end
if isfield(parm,'Fstep'),Fstep  = parm.Fstep ; end

% # of embedding dimension
if isfield(parm,'Dtau')
	D    = parm.Dtau;
	tau  = parm.Tau;
else
	D    = 1;
	tau  = 1;
end

if Ntrain < 1, Info = []; return; end;

if isfield(parm,'Npre_train')
	Npre_train = parm.Npre_train;
else
%	Npre_train = Ntrain;
	if Tall >= 2*M*D
		Npre_train = 0;
	elseif Tall >= M*D
		Npre_train = fix(Ntrain/2);
	else
		Npre_train = Ntrain;
	end
end

if Npre_train > Ntrain, Npre_train = Ntrain; end;

fprintf('--- Output Dimension  = %d\n',N)
fprintf('--- Input  Dimension  = %d\n',M)
fprintf('--- Embedding  Dimension  = %d\n',D)
fprintf('--- Number of trials  = %d\n',Ntrial)
fprintf('--- Number of training sample = %d\n',Tall)
fprintf('--- Total update iteration    = %d (%d)\n',Ntrain,Npre_train)

% Original input dimension
Xdim  = M;
M_ALL = Xdim*D;

%
% --- Initialization
%  W , A : Initial variable to use 1st update
%

% Input/Output variance
sx = mean((X(:) - mean(X(:))).^2);
sy = mean((Y(:) - mean(Y(:))).^2);

SY0 = mean(sy);
A0  = 1./mean(sx);

if isfield(parm,'Ta0') && parm.Ta0 > 0,
	Ta0 = parm.Ta0;
	a0  = parm.a0 * A0;
else
	Ta0 = 0;
	a0  = 1;
end

% 'linear_sparse_space'
%    ARD term = alpha * W^2  : A = W^2
if isfield(Model,'ix_act')
	SY  = mean(Model.SY);  % 1 x 1

	% Recover old estimate of A in full embedding space
	ix_act = Model.ix_act;
	Aold = zeros(1,M_ALL); % 1 x M_ALL
	Wold = zeros(N,M_ALL);

	Aold(ix_act) = sum(Model.A,1);
	Wold(:,ix_act) = Model.W;

	if M_ALL ~= Xdim*D,
		fprintf('M_ALL=%d,Xdim=%d,D=%d\n',M_ALL ,Xdim ,D)
		error('M_ALL ~= Xdim*D')
	end
	% Summation over time delay component
	A  = mean(reshape(Aold,[Xdim,D]),2);	% Xdim x 1

	% Active index in input space without embedding
	IX_act = find(A > 0);
	A  = A(IX_act);
	M  = length(IX_act);

	id  = repmat( (0:(D-1))* Xdim ,[M 1]) + repmat(IX_act, [1 D]);
	W   = Wold(:,id(:)) ;  % N x (M*D)

	X  = X(IX_act,:,:); % M x T
else
	A = repmat(A0, [Xdim, 1]);
	W = zeros(N,M_ALL);
	SY  = SY0;
	% Save original input index for pruning
	IX_act = (1:Xdim)';
end

rx  = 1/(D + 1);
SX  = SY*rx;

fprintf('a_min = %g\n', a_min)
fprintf('SY0 = %g\n', SY0)
fprintf('SY  = %g\n', SY)

% Active index in input space without embedding
ix_act  = (1:M)';

% Delay embedding index for W
Wid = [(0:(D-1))'* M + 1 , (1:D)'* M];

% Input covariance (Spatial dimension)
Trial = Trial+(D-1)*tau;

XX = zeros(M,M);
for n=1:Ntrial
	XX = XX + X(:,1:Trial(n), n) * X(:,1:Trial(n), n)';
end
XX  = XX / sum(Trial);

SW = inv( Tall.* XX + diag(1./A) );

% Working variable
dY  = zeros(N,T,Ntrial);% N x T x Ntrial
dYX = zeros(N,M*D);     % N x M*D
G_A = zeros(M,1);       % M x 1
WW  = sum(reshape(sum(W.^2,1),[M,D]),2);	% M x 1

% Free energy histry
FE  = zeros(Ntrain,1);
LP  = zeros(Ntrain,1);
H   = zeros(Ntrain,1);
Err = zeros(Ntrain,1);
Mhist = zeros(Ntrain,1);

% ARD hyper param. history
if isfield(parm,'Debug') && ~isempty(parm.Debug) && parm.Debug > 0
	Debug = 1;
	A_tmp = zeros(M_ALL, ceil(Ntrain/Nskip));
else
	Debug = 0;
end
% recover all component
A_all  = zeros(M,1);

k_save  = 0;

%%%%%% Learning Loop %%%%%%
for k=1:Ntrain
	% Ainv = alpha , A = 1/alpha
	Ainv    = 1./A;

	%    dY  = Y - W * X
	%    dYX = (dY * X')/T;
	dY  = error_delay_time_sw(X, Y, W, tau, Trial); % N x T x Ntrial
	dY  = dY .* trial_sw;                           % mask invalid samples
    dYX = error_corr_delay_sw(X, dY, D, tau, Trial);% N x M*D

    dYY = sum(dY(:).^2)/(N*Tall);
    dYX = dYX/(Tall);

    % Noise variance update
    WWA = sum(WW .*Ainv);
    SX  = (rx) .* dYY + WWA/(N*Tall) ;
	%  SX  = rx .* dYY + sum(SW(:) .* XX(:)) ./rx ;

    % Prevent zero variance
    SX  = max( SX, MINVAL);
	SY  = SX/rx;

    % Log variance
    SWA     = diag(SW) .* Ainv;
    log_sw  = (N*D)*(log_det(SW) + sum(log(Ainv)) - sum(SWA) + M);
    log_sy  = N * sum( log(SY) );
    log_a   = Ta0*sum(log(Ainv) - a0.*Ainv + 1);

    % Free energy
    LP(k)  = - (0.5*Tall) * log_sy ;
    H(k)   = 0.5*( log_sw + log_a );
%    H(k)   = 0.5*( log_sw + log_a - WWA );
    FE(k)  = LP(k) + H(k);
    Err(k) = sum(dYY)/(SY0);

    % Weight variance
    % SW = inv( (T./SX) .* XX + diag(Ainv) );   ( M x M )
    SW  = XX + diag( Ainv./Tall );

   	if rcond(SW) > MinCond,
		SW  = inv( SW );
	else
		SW  = pinv( SW );
	end

    % Weight update
    % W   = (T/SX)*( W * XX + rx .* dYX ) * SW; ( N x MD )
	for j=1:D
    	W(:,Wid(j,1):Wid(j,2)) = ( W(:,Wid(j,1):Wid(j,2)) * XX ...
    	                       + rx * dYX(:,Wid(j,1):Wid(j,2)) ) * SW;
    end

	%  ARD for each input (average over output & delay)
	WW  = sum(reshape(sum(W.^2,1),[M,D]),2);		% M x 1

    % Hyper parameter for weight variance (ARD)
    SW  = SW /Tall; % = inv(Tall*XX + diag(Ainv))
%  	G_A = 1 - diag(SW).*Ainv ;
%  	    = diag( ((Tall*XX + Ainv) - Ainv) * SW )
%  	    = diag( (Tall*XX) * SW )
  	G_A = Tall*sum(XX.*SW,2);
	G_A = max((G_A), MINVAL);

	if k <= Npre_train,
		% VB update rule (Stable)
		% (N*D)*A  = WW./SX + (N*D) * diag(SW)
		% (N*D)*A*(1 - diag(SW)./A)  = WW./SX
		% Modefied VB update rule (Stable)
		A  = (WW./SX + (N*D) * diag(SW) + 2*Ta0*a0)/( N*D + 2*Ta0 );
%	    A  = sqrt(A.*(WW./SX)./(G_A*N*D));
	else
	    % Accelerated update rule
	    A  = sqrt(A.*(WW./SX)./(G_A*N*D));
%	    A  = ((WW./SX) + 2*Ta0*a0)./(G_A*N*D + 2*Ta0);
	end

    % Prune small variance
    if Prune > 0
	    % Find active input dimension
	    ix_act_old = ix_act;

	    % Recover all component
	    switch	Prune
	    case	1
		    A_all(ix_act) = WW/max(WW);    % Prune by Weight
	    case	2
		    A_all(ix_act) = A /max(A);	   % Prune by Alpha
	    case	3
		    A_all(ix_act) = A * (1/SX);    % Prune by Alpha
	    end

	    % Find active input dimension (absolute index)
	    ix_act = find( A_all > a_min ); % effective indices
	    Mnew   = length(ix_act);  		% # of effective input

	    if Mnew < M,
		    % convert to relative index
		    jx_act = trans_index(ix_act,ix_act_old,Xdim);

		    A   = A(jx_act) ;  % 1 x M
		    SW  = SW(jx_act,jx_act);  % M x M

			% Effective delay embedding index for W
			id  = repmat( (0:(D-1))* M ,[Mnew 1]) + repmat(jx_act, [1 D]);
		    W   = W(:,id(:)) ;  % N x MD

			% New delay embedding index for W
		    M   = Mnew;
			Wid = [(0:(D-1))'* M + 1 , (1:D)'* M];

			WW  = WW(jx_act);
		    X   = X(jx_act,:,:);  	 % M x T x Ntrial
		    XX  = XX(jx_act,jx_act); % M x M
		end
    else
		A = max(A,MINVAL);
    end
    % END of if Prune == 1
	Mhist(k) = M;

    if mod(k, Nskip)==0
        % Save history
		if Debug == 1
        	k_save = k_save + 1;
        	A_tmp(:,k_save) = A_all(:);
		end

        fprintf('Iter=%4d, M =%4d, err=%g, F=%g, H=%g, SY=%g, Wmin=%g\n',...
        		        k, M, Err(k), FE(k), -2*H(k)/log(Tall), ...
        		        SY, min(WW)/max(WW));
    end

    % Convergence check
	if k > Ncheck,
		Fdif = (FE(k) - FE(k-Fstep))/(abs(FE(k))+eps);
	else
		Fdif = Fdiff + 1;
	end

	if (Fdiff > abs(Fdif)),
		fprintf('Converged : Free energy change = %g\n',Fdif)
		break;
	end;

end

% convert to relative index
ix_act = IX_act(ix_act);

% Delay embedding index for W in full input space
id  = repmat( (0:(D-1))* Xdim ,[M 1]) + repmat(ix_act, [1 D]);
% Recover A in full embedding input space
A   = repmat(A(:), [1 D]);

% Active index
Model.ix_act = id(:);
Model.M_all  = M_ALL ;

% Save output variable
Model.A    = A(:)' ; % 1 x M*D
Model.W    = W ;
Model.SY   = SY;
Model.SW   = SW; % = inv(Tall*XX + diag(Ainv))

Model.method = 'linear_sparse_space';
Model.mode   = 'cov';
Model.sparse = 'sparse';

% Save history
Info.FE  = FE(1:k);
Info.LP  = LP(1:k);
Info.H   = H(1:k) ;
Info.Err = Err(1:k);
Info.M   = Mhist(1:k);

if exist('A_tmp','var')
	Info.A   = A_tmp(:,1:k_save) ;
end

%%% ---- Index transformation from old active_index to current active_index
function	jx = trans_index(ix,ix_old,M)

N = length(ix_old);
Itrans = zeros(M,1);
Itrans(ix_old) = 1:N;

jx = Itrans(ix);

author: KamitaniLab

repo: GenericObjectDecoding


text:
function a = r82vec_sort_quick_a ( n, a )

%*****************************************************************************80
%
%% R82VEC_SORT_QUICK_A ascending sorts an R82VEC using quick sort.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 February 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of entries in the array.
%
%    Input, real A(2,N), the array to be sorted.
%
%    Output, real A(2,N), the sorted array.
%
  level_max = 25;
  dim_num = 2;

  if ( n < 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R82VEC_SORT_QUICK_A - Fatal error!\n' );
    fprintf ( 1, '  N < 1.\n' );
    error ( 'R82VEC_SORT_QUICK_A - Fatal error!' );
  elseif ( n == 1 )
    return
  end

  level = 1;
  rsave(level) = n + 1;
  base = 1;
  n_segment = n;

  while ( 1 )
%
%  Partition the segment.
%
    [ a(1:2,base:base+n_segment-1), l_segment, r_segment ] = ...
      r82vec_part_quick_a ( n_segment, a(1:2,base:base+n_segment-1) );
%
%  If the left segment has more than one element, we need to partition it.
%
    if ( 1 < l_segment )

      if ( level_max < level )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'R82VEC_SORT_QUICK_A - Fatal error!\n' );
        fprintf ( 1, '  Exceeding recursion maximum of %d\n', level_max );
        error ( 'R82VEC_SORT_QUICK_A - Fatal error!' );
      end

      level = level + 1;
      n_segment = l_segment;
      rsave(level) = r_segment + base - 1;
%
%  The left segment and the middle segment are sorted.
%  Must the right segment be partitioned?
%
    elseif ( r_segment < n_segment )

      n_segment = n_segment + 1 - r_segment;
      base = base + r_segment - 1;
%
%  Otherwise, we back up a level if there is an earlier one.
%
    else

      while ( 1 )

        if ( level <= 1 )
          return
        end

        base = rsave(level);
        n_segment = rsave(level-1) - rsave(level);
        level = level - 1;

        if ( 0 < n_segment )
          break
        end

      end

    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function spdata = APPgetSpData(image, doogFilters, textonFilters, imsegs)
% spdata = APPgetSpData(image, doogFilters, textonFilters, imsegs)
%
% Gets the features and data corresponding to the superpixels given by
% imsegs.
% features:
%   1-12: mean absolute filter response of diff of oriented Gaussian filters
%   13: mean of 1-12
%   14: argmax of 1-12
%   15: max-median of 1-12
%   possibly texton features (15+1:15+nTextons+2)
%   16-18: mean rgb values
%   19-21: hsv conversion from mean rgb values
%   22-23: mean x and y locations
%
% Input:
%   image: rgb image to be analyzed
%   doogFilters: filters for texture (could be empty)
%   textonFilters: filters for texture (could be empty)
%   imsegs: superpixel structure
% Output:
%   spdata: structure for data corresponding to each superpixel
%
% Copyright(C) Derek Hoiem, Carnegie Mellon University, 2005
% Permission granted to non-commercial enterprises for
% modification/redistribution under GNU GPL.
% Current Version: 1.0  09/30/2005

nAngles = size(doogFilters, 3);

% compute texture responses
grayim = rgb2gray(image);

%disp(num2str(size(image, 1)*size(image, 2)*nAngles /1024/1024*8))
orientationImages = zeros(size(image, 1), size(image, 2), nAngles);
% get and record the filter responses for each angle
for (i=1:nAngles)
    orientationImages(:, :, i) = abs(conv2(grayim,doogFilters(:, :, i),'same')-grayim);
end
%spdata.orientationImages = orientationImages;

textonImages = compute_texton_response(grayim, textonFilters);
nTextons = numel(textonFilters);

[height, width, nb] = size(image);

% features: angle response peaks, mean angle response, angle with
% largest response, dominance of largest filter, mean abs texton
% repsonses, largest texton response, dominance of largest response,
% mean of rgb hsv x y
nfeatures = nAngles + 3 + nTextons + 6 + 2;
if nTextons > 0
    nfeatures = nfeatures + 2;
end

% for each segmentation
for i = 1:length(imsegs)

    nseg = imsegs(i).nseg;
    spdata(i).npixels = imsegs(i).npixels;
    spdata(i).adjmat = imsegs(i).adjmat;

    features = zeros(nseg, nfeatures);

    % oriented filters
    maar = zeros(nseg, nAngles);
    for a = 1:nAngles
        maar(:, a) = APPgetSpMeans(imsegs(i), orientationImages(:, :, a));
    end
    features(:, 1:nAngles) = maar;
    features(:, nAngles+1) = mean(maar, 2);
    [maxval, features(:, nAngles+2)] = max(maar, [], 2);
    features(:, nAngles+3) = maxval - median(maar, 2);
    spdata(i).orientation = maar;
    spdata(i).edginess = mean(orientationImages, 3);
    clear orientationImages;

    % textons
    cf = nAngles+3;
    for t = 1:nTextons
        features(:, cf+t) = APPgetSpMeans(imsegs(i), textonImages(:, :, t));
    end
    cf = cf + nTextons;
    if nTextons>0
        [maxval, features(:, cf+1)] = max(features(:, cf+1:cf+nTextons), [], 2);
        features(:, cf+2) = maxval - median(features(:, cf+1:cf+nTextons), 2);
        clear textonImages;
        cf = cf + 2;
    end

    % color
    rgb = zeros(nseg, 3);
    for b = 1:3
        rgb(:, b) = APPgetSpMeans(imsegs(i), image(:, :, b));
    end
    features(:, cf + (1:3)) = rgb;
    features(:, cf + (4:6)) = rgb2hsv(rgb);

    % location
    cf = cf + 6;
    yim = 1-repmat(((0:height-1)/(height-1))', 1, width);
    xim = repmat(((0:width-1)/(width-1)), height, 1);
    features(:, cf+1) = APPgetSpMeans(imsegs(i), yim);
    features(:, cf+2) = APPgetSpMeans(imsegs(i), xim);
    spdata(i).features = features;
end

author: zouchuhang

repo: LayoutNet


text:
function f=ref_igdgt(coef,g,a,M,c_t,c_f,c_w)
%REF_GDGT  Reference generalized DGT
%   Usage:  f=ref_igdgt(c,g,a,M,c_t,c_f,c_w);
%
%   Linear algebra version of the algorithm. Create big matrix
%   containing all the basis functions and multiply with the transpose.

N=size(coef,1)/M;
L=N*a;

F=zeros(L,M*N);

l=(0:L-1).';

if length(g)<L
  g=fir2long(g,L);
end;

for n=0:N-1
  for m=0:M-1
    F(:,M*n+m+1)=exp(2*pi*i*(m+c_f)*(l+c_t)/M).*circshift(g,n*a+c_w);
  end;
end;

f=F*coef;


text:
function visualizemodel(model, components, layers)
% Visualize a mixture of star models.
%   visualizemodel(model)
%
% Arguments
%   model       Model to visualize
%   components  Which components to draw

% AUTORIGHTS
% -------------------------------------------------------
% Copyright (C) 2011-2012 Ross Girshick
% Copyright (C) 2008, 2009, 2010 Pedro Felzenszwalb, Ross Girshick
%
% This file is part of the voc-releaseX code
% (http://people.cs.uchicago.edu/~rbg/latent/)
% and is available under the terms of an MIT-like license
% provided in COPYING. Please retain this notice and
% COPYING if you use this file (or a portion of it) in
% your project.
% -------------------------------------------------------

clf;
if nargin < 2
  components = 1:length(model.rules{model.start});
end

if nargin < 3
  layers = 1;
end

k = 1;
for i = components
  for layer = layers
    visualizecomponent(model, i, length(layers)*length(components), k, layer);
    k = k+1;
  end
end

function visualizecomponent(model, c, nc, k, layer)

rhs = model.rules{model.start}(c).rhs;
root = -1;
parts = [];
defs = {};
anchors = {};
% assume the root filter is first on the rhs of the start rules
if model.symbols(rhs(1)).type == 'T'
  % handle case where there's no deformation model for the root
  root = model.symbols(rhs(1)).filter;
else
  % handle case where there is a deformation model for the root
  root = model.symbols(model.rules{rhs(1)}(layer).rhs).filter;
end
for i = 2:length(rhs)
  defs{end+1} = model_get_block(model, model.rules{rhs(i)}(layer).def);
  anchors{end+1} = model.rules{model.start}(c).anchor{i};
  fi = model.symbols(model.rules{rhs(i)}(layer).rhs).filter;
  parts = [parts fi];
end
% make picture of root filter
pad = 2;
bs = 20;
w = foldHOG(model_get_block(model, model.filters(root)));
scale = max(w(:));
im = HOGpicture(w, bs);
im = imresize(im, 2);
im = padarray(im, [pad pad], 0);
im = uint8(im * (255/scale));

% draw root
numparts = length(parts);
if numparts > 0
  subplot(nc,3,1+3*(k-1));
else
  subplot(nc,1,k);
end
imagesc(im)
colormap gray;
axis equal;
axis off;

% draw parts and deformation model
if numparts > 0
  def_im = zeros(size(im));
  def_scale = 500;
  for i = 1:numparts
    % part filter
    w = model_get_block(model, model.filters(parts(i)));
    p = HOGpicture(foldHOG(w), bs);
    p = padarray(p, [pad pad], 0);
    p = uint8(p * (255/scale));
    % border
    p(:,1:2*pad) = 128;
    p(:,end-2*pad+1:end) = 128;
    p(1:2*pad,:) = 128;
    p(end-2*pad+1:end,:) = 128;
    % paste into root
    x1 = (anchors{i}(1))*bs+1;
    y1 = (anchors{i}(2))*bs+1;
    x2 = x1 + size(p, 2)-1;
    y2 = y1 + size(p, 1)-1;
    im(y1:y2, x1:x2) = p;

    % deformation model
    probex = size(p,2)/2;
    probey = size(p,1)/2;
    for y = 2*pad+1:size(p,1)-2*pad
      for x = 2*pad+1:size(p,2)-2*pad
        px = ((probex-x)/bs);
        py = ((probey-y)/bs);
        v = [px^2; px; py^2; py];
        p(y, x) = defs{i} * v * def_scale;
      end
    end
    def_im(y1:y2, x1:x2) = p;
  end

  % plot parts
  subplot(nc,3,2+3*(k-1));
  imagesc(im);
  colormap gray;
  axis equal;
  axis off;

  % plot deformation model
  subplot(nc,3,3+3*(k-1));
  imagesc(def_im);
  colormap gray;
  axis equal;
  axis off;
end

set(gcf, 'Color', 'white')

author: rbgirshick

repo: voc-dpm

save_path: github-repos/MATLAB/rbgirshick-voc-dpm


text:
function output = F_inner_product_normalized(input_layers)

input1 = input_layers{1}.a;
input2 = input_layers{2}.a;

output = sum(input1 .* input2);

output = output / size(input1,1);

end

author: singaxiong

repo: SignalGraph


text:
function [ n_data, n, fn ] = i4_factorial_values ( n_data )

%*****************************************************************************80
%
%% I4_FACTORIAL_VALUES returns values of the factorial function.
%
%  Discussion:
%
%    0! = 1
%    I! = Product ( 1 <= J <= I ) I
%
%    In Mathematica, the function can be evaluated by:
%
%      n!
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz and Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Wolfram Media / Cambridge University Press, 1999.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, integer N, the argument of the function.
%
%    Output, integer FN, the value of the function.
%
  n_max = 13;

  fn_vec = [ ...
            1, ...
            1, ...
            2, ...
            6, ...
           24, ...
          120, ...
          720, ...
         5040, ...
        40320, ...
       362880, ...
      3628800, ...
     39916800, ...
    479001600 ];

  n_vec = [ ...
     0,  1,  2,  3, ...
     4,  5,  6,  7, ...
     8,  9, 10, 11, ...
    12 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    n = 0;
    fn = 0;
  else
    n = n_vec(n_data);
    fn = fn_vec(n_data);
  end

  return
end


text:
% Excersize 1 - Calculate the equilibrium surface temperature of an icy
% body (e.g. a comet) at a given distance from the Sun.

% This sample program is a script that balances the equation of heat
% transfer on the surface of a ball made of ice, to find an equilibrium
% temperature. This script demonstrates how physical units can be used in
% calculations.

%% Set up the workspace with needed variables
si=setUnits;
% We begin by calling the |setUnits| function in order to get the interface
% structure that holds all predefined units.

LSun=3.8e26*si.joule/si.second;
% Total brightness of the Sun. This is the energy radiated by the Sun every
% second. (Tip: si.s is shorthand for si.second.)

dSun=2.5*1.496e11*si.meter;
% Distance from the Sun to the body whose temperature we must calculate. I
% have chosen a distance of 2.5 astronomical units (Earth-Sun distances)
% that puts our "comet" betwen Mars and Jupiter. (Tip: si.astronomical_unit
% is already defined, and si.AU is even easier. Also, si.m is shorthand for
% si.meter.)

KB=1.3806e-23*si.joule/si.K;
% Boltzmann's constant. (Tip: si.boltzmann is defined.)

stefan=5.67e-8*si.joule/si.m^2/si.s/si.K^4;
% The Stefan-Boltzmann constant. (Tip: Guess what, si.stefan_boltzmann.)

Hice=2.8345e6*si.joule/si.kg;
% Latent heat for sublimation of ice.

A=3.56e12*si.newton/si.m^2; B=6141.667*si.kelvin;
% These parameters are used in an empirical formula for the saturation
% vapor pressure. (Tip: si.N and si.K are shorter.)

avogad=6.022e23*si.mole^(-1);
% Avogadro's number. (Tip: Try si.avogadro.)

mH2O=18*si.g/si.mole;
% The molecular weight of water.

%% Now, Lion-hunt for the temperature of equilibrium
% When we write down the equation:
% radiant energy from the Sun = (radiant + sublimation)energy from the body
% we essentially get an implicit function of temperature vs. distance from
% the Sun. This function however cannot be inverted so we need to hunt for
% the temperature that will balance the equation.

tol=1e-3;
% We'll use this as the relative accuracy.

Tlo=1*si.K; Thi=400*si.K;
% The desired temperature is definitely bracketed by these values.

while (Thi-Tlo)/Thi>tol
    T=(Tlo+Thi)/2;
    fT=(stefan*T^4)+(Hice*sqrt(mH2O/(2*avogad*pi*KB*T))*A*exp(-(B/T)))-...
        (0.25*(LSun/(4*pi*dSun^2)));
    if double(fT)>0
        Thi=T;
    else
        Tlo=T;
    end
end

%% That's it. The desired temperature is
fprintf('The equilibrium temperature is\n')
T
% Try to introduce an easy-to-make-hard-to-detect error in one of the above
% statements, and run the script again. Instead of getting a wrong answer
% you will probably get an informative error message that will help you
% trace the "bug".


text:
function [sP, isNew] = newSphericalPlot(v,varargin)
% split plot in upper and lower hemisphere
%
% 1: axis given -> no sphericalRegion stored -> compute sphericalRegion -> finish
% 2: axis is hold and has sphericalRegion -> use multiplot
% 3: new multiplot

% case 1: predefined axis
% -----------------------
if check_option(varargin,'parent')

  ax = get_option(varargin,'parent');

  % axis is already a spherical plot
  if isappdata(ax(end),'sphericalPlot') && ishold(ax(end))

    for i = 1:length(ax)
      sP(i) = getappdata(ax(i),'sphericalPlot');
    end
    isNew = false;

  else % set up a new spherical axes if required

    % extract spherical region
    % TODO: it might happen that the spherical region needs two axes
    sR = getPlotRegion(v,varargin{:});

    % extract projection
    proj = getProjection(sR,varargin{:});

    % create a new spherical plot
    sP = sphericalPlot(ax,proj(1),varargin{:});
    isNew = true;

  end
  return;
end

% create a new mtexFigure or get a reference to it
[mtexFig,isNew] = newMtexFigure(varargin{:});

%~check_option(varargin,'add2all') ||

if isNew || ~isappdata(mtexFig.currentAxes,'sphericalPlot')

  % get spherical region
  sR = getPlotRegion(v,varargin{:});

  % extract projection(s)
  % this might return two projections for upper and lower hemisphere
  proj = getProjection(sR,varargin{:});
  holdState = getHoldState(mtexFig.gca);

  for i = 1:numel(proj)

    % create a new axis
    if i>1, mtexFig.nextAxis; end
    hold(mtexFig.gca,holdState);

    % display upper/lower if needed
    if numel(proj)>1
      if ~proj(i).sR.isUpper
        tr = {'TR','lower'};
      elseif ~proj(i).sR.isLower
        tr = {'TR','upper'};
      else
        tr = {};
      end
    else
      tr = {};
    end

    % create a new spherical plot
    sP(i) = sphericalPlot(mtexFig.gca,proj(i),tr{:},varargin{:});         %#ok<AGROW>

  end
  mtexFig.drawNow(varargin{:});
  isNew = true;

elseif check_option(varargin,'add2all') % add to or overide existing axes

  for i = 1:numel(mtexFig.children)

    sP(i) = getappdata(mtexFig.children(i),'sphericalPlot'); %#ok<AGROW>

  end

else

  sP = getappdata(mtexFig.currentAxes,'sphericalPlot');

end

end

% ---------------------------------------------------------
function sR = getPlotRegion(varargin)
% returns spherical region to be plotted

% default values from the vectors to plot
if isa(varargin{1},'vector3d')
  sR = varargin{1}.region(varargin{2:end});
else
  sR = sphericalRegion;
end
sR = getClass(varargin,'sphericalRegion',sR);

% check for simple options
if check_option(varargin,'complete')
  sR = sphericalRegion;
end
if check_option(varargin,'upper')
  sR = sR.restrict2Upper;
elseif check_option(varargin,'lower')
  sR = sR.restrict2Lower;
end

% extract antipodal
sR.antipodal = check_option(varargin,'antipodal') || ...
  (isa(varargin{1},'vector3d') && varargin{1}.antipodal);

% for antipodal symmetry reduce to halfsphere
if sR.antipodal && sR.isUpper && sR.isLower &&...
    ~check_option(varargin,'complete')
  sR = sR.restrict2Upper;
end

end
% ---------------------------------------------------------
function proj = getProjection(sR,varargin)

proj = get_option(varargin,'projection','earea');

if ~isa(proj,'sphericalProjection')

  switch proj
    case 'plain', proj = plainProjection(sR);

    case {'stereo','eangle'}, proj = eangleProjection(sR); % equal angle

    case 'edist', proj = edistProjection(sR); % equal distance

    case {'earea','schmidt'}, proj = eareaProjection(sR); % equal area

    case 'orthographic',  proj = orthographicProjection(sR);

    case 'square',  proj = squareProjection(sR);

    case 'gnonomic', proj = gnonomicProjection(sR);

    otherwise

      error('%s\n%s','Unknown projection specified! Valid projections are:',...
        'plain, stereo, eangle, edist, earea, schmidt, orthographic','square')
  end
end

if ~isa(proj,'plainProjection') && sR.isUpper && sR.isLower
  proj = [proj,proj];
  proj(1).sR = proj(1).sR.restrict2Upper;
  proj(2).sR = proj(2).sR.restrict2Lower;
end

end


text:
classdef ImdbDetectionFullSupervision < ImdbMatbox
    properties(SetAccess = protected, GetAccess = public)
        negOverlapRange = [0.1 0.5];
        boxesPerIm = 64;
        boxRegress = true;
        instanceWeighting = false;
        maxImageSize = 600;
    end
    methods
        function obj = ImdbDetectionFullSupervision(imageDir, imExt, matboxDir, filenames, datasetIdx, meanIm)
            obj@ImdbMatbox(imageDir, imExt, matboxDir, filenames, datasetIdx, meanIm);
        end

        function [batchData, numElements] = getBatch(obj, batchInds, net, ~)
            if length(batchInds) > 1
                error('Error: Only supports subbatches of 1!');
            end

            if nargin == 2
                gpuMode = false;
            else
                gpuMode = strcmp(net.device, 'gpu');
            end

            % Load image. Make correct size. Subtract average im.
            [image, oriImSize] = obj.LoadImage(batchInds, gpuMode);

            % Sample boxes
            gStruct = obj.LoadGStruct(batchInds);

            % Flip the image and boxes at training time
            % Note: flipLR alternates between true and false in ImdbMatbox.initEpoch()
            if obj.flipLR && strcmp(obj.datasetMode, 'train')
                currImT = fliplr(image);
                currBoxesT = gStruct.boxes;
                currBoxesT(:,3) = oriImSize(2) - gStruct.boxes(:,1) + 1;
                currBoxesT(:,1) = oriImSize(2) - gStruct.boxes(:,3) + 1;
                gStruct.boxes = currBoxesT;
                image = currImT;
            end

            if ismember(obj.datasetMode, {'train', 'val'})
                [boxes, labels, ~, overlapScores, regressionFactors] = obj.SamplePosAndNegFromGstruct(gStruct, obj.boxesPerIm);

                % Assign elements to cell array for use in training the network
                numElements = obj.boxesPerIm;

                numBatchFields = 2 * (4 + obj.boxRegress + obj.instanceWeighting);
                batchData = cell(numBatchFields, 1);
                idx = 1;
                batchData{idx} = 'input';       idx = idx + 1;
                batchData{idx} = image;         idx = idx + 1;
                batchData{idx} = 'label';       idx = idx + 1;
                batchData{idx} = labels';       idx = idx + 1;
                batchData{idx} = 'boxes';       idx = idx + 1;
                batchData{idx} = boxes';        idx = idx + 1;
                batchData{idx} = 'oriImSize';   idx = idx + 1;
                batchData{idx} = oriImSize;     idx = idx + 1;

                if obj.boxRegress
                    batchData{idx} = 'regressionTargets';   idx = idx + 1;
                    batchData{idx} = regressionFactors';    idx = idx + 1;
                end
                if obj.instanceWeighting
                    instanceWeights = overlapScores;
                    instanceWeights(labels == 1) = 1;
                    instanceWeights = reshape(instanceWeights, [1 1 1 length(instanceWeights)]); % VL-Feat way :-S
                    batchData{idx} = 'instanceWeights';     idx = idx + 1;
                    batchData{idx} = instanceWeights;       %idx = idx + 1;
                end

            else
                % Test set. Get all boxes
                numElements = size(gStruct.boxes,1);
                batchData{6} = oriImSize;
                batchData{5} = 'oriImSize';
                batchData{4} = gStruct.boxes';
                batchData{3} = 'boxes';
                batchData{2} = image;
                batchData{1} = 'input';
            end

        end

        function [image, oriImSize] = LoadImage(obj, batchIdx, gpuMode)
            % image = LoadImage(obj, batchIdx)
            % Loads an image from disk, resizes it, and subtracts the mean image
            imageT = single(imread([obj.imageDir obj.data.(obj.datasetMode){batchIdx} obj.imExt]));
            oriImSize = double(size(imageT));
            if numel(obj.meanIm) == 3
                for colourI = 1:3
                    imageT(:,:,colourI) = imageT(:,:,colourI) - obj.meanIm(colourI);
                end
            else
                imageT = imageT - imresize(obj.meanIm, [oriImSize(1) oriImSize(2)]); % Subtract mean im
            end

            % Resize image to be at most x pixels in each dimension
            % If you notice a "CUDA_ERROR_ILLEGAL_ADDRESS" error, use a
            % smaller maximum image size.
            resizeFactor = obj.maxImageSize / max(oriImSize(1:2));
            if gpuMode
                image = gpuArray(imageT);
                image = imresize(image, resizeFactor);
            else
                image = imresize(imageT, resizeFactor, 'bilinear', 'antialiasing', false);
            end
        end

        function [boxes, labels, keys, overlapScores, regressionTargets] = SamplePosAndNegFromGstruct(obj, gStruct, numSamples)
            % Get positive, negative, and true GT keys
            [maxOverlap, classOverlap] = max(gStruct.overlap, [], 2);

            posKeys = find(maxOverlap >= 0.5 & gStruct.class == 0);
            negKeys = find(maxOverlap < obj.negOverlapRange(2) & maxOverlap >= obj.negOverlapRange(1) & gStruct.class == 0);
            gtKeys = find(gStruct.class > 0);

            % Get correct number of positive and negative samples
            numExtraPos = numSamples * obj.posFraction - length(gtKeys);
            numExtraPos = min(numExtraPos, length(posKeys));
            if numExtraPos > 0
                posKeys = posKeys(randperm(length(posKeys), numExtraPos));
            else
               numExtraPos = 0;
            end
            numNeg = numSamples - numExtraPos - length(gtKeys);
            numNeg = min(numNeg, length(negKeys));
            negKeys = negKeys(randperm(length(negKeys), numNeg));

            % Concatenate for final keys and labs
            keys = cat(1, gtKeys, posKeys, negKeys);
            labels = cat(1, gStruct.class(gtKeys), classOverlap(posKeys), zeros(numNeg, 1));
            labels = single(labels + 1); % Add 1 for background class
            boxes = gStruct.boxes(keys,:);

            overlapScores = cat(1, ones(length(gtKeys),1), maxOverlap(posKeys), maxOverlap(negKeys));

            % Calculate regression targets.
            % Jasper: I simplify Girshick by implementing regression through four
            % scalars which scale the box with respect to its center.
            if nargout == 5
                % Create NaN array: nans represent numbers which will not be active
                % in regression
                regressionTargets = nan([size(boxes,1) 4 * obj.numClasses], 'like', boxes);

                % Get scaling factors for all positive boxes
                gtBoxes = gStruct.boxes(gtKeys,:);
                for bI = 1:length(gtKeys)+length(posKeys)
                    % Get current box and corresponding GT box
                    currPosBox = boxes(bI,:);
                    [~, gtI] = BoxBestOverlapFastRcnn(gtBoxes, currPosBox);
                    currGtBox = gtBoxes(gtI,:);

                    % Get range of regression target based on the label of the gt box
                    targetRangeBegin = 4 * (labels(bI)-1)+1;
                    targetRange = targetRangeBegin:(targetRangeBegin+3);

                    % Set regression targets
                    regressionTargets(bI, targetRange) = BoxRegressionTargetGirshick(currGtBox, currPosBox);
                end
            end
        end

        function SetBoxRegress(obj, doRegress)
            obj.boxRegress = doRegress;
        end

        function SetInstanceWeighting(obj, doInstanceWeighting)
            obj.instanceWeighting = doInstanceWeighting;
        end
    end % End methods
end % End classdef

author: nightrome

repo: matconvnet-calvin


text:
function h = glm_plotTcs(model,voxNames,conds,active,control);
%
% h = glm_plotTcs(model,voxNames,conds,active,control);
%
% Make a figure showing the deconvolved responses
% for each 'voxel'/roi. If voxNames is supplied,
% will title each subplot w/ those names. If active and
% control are supplied, will perform a contrast
% and put up the stats for the peak data point.
%
% if conds is supplied, will only plot the
% specified conditions.
%
% ras 02/05
if ieNotDefined('voxNames')
    voxNames = [];
end

if ieNotDefined('active')
    active = [];
end

if ieNotDefined('control')
    control = [];
end

if ieNotDefined('conds')
    conds = 1:size(model.h_bar,2);
end

h = figure('Color','w',...
           'Units','Normalized',...
           'Position',[0 .2 .4 .6]);

Y = model.h_bar;
E = model.h_bar_sem;

peak = find(Y(:,2,1)==max(Y(:,2,1)))
if ~isempty(active) & ~isempty(control)
    sig = er_contrast(model,active,control,'weights',peak);
else
    sig = [];
end

nvoxels = size(Y,3);

nrows = ceil(sqrt(nvoxels));
ncols = ceil(nvoxels/nrows);

for i = 1:nvoxels
    subplot(nrows,ncols,i);
    errorbar(Y(:,conds,i),E(:,conds,i));
    if ~isempty(voxNames)
        title(voxNames{i});
    end
    if ~isempty(sig)
        AX = axis;
        msg = sprintf('p = %1.3f',10^(-1*sig(i)));
        text(AX(1)+0.63*(AX(2)-AX(1)), AX(3)+0.8*(AX(4)-AX(3)), msg);
    end
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function r = times(a,b)
%TIMES        Gradient multiplication  a .* b
%

% written  10/16/98     S.M. Rump
% modified 11/30/98     S.M. Rump  array/scalar operations
% modified 04/04/04     S.M. Rump  set round to nearest for safety
%                                    improved performance
%                                    complete redesign
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 12/05/05     S.M. Rump  improved performance (thanks to J. Kubitz)
% modified 09/28/08     S.M. Rump  check for rounding to nearest improved
% modified 08/26/12     S.M. Rump  global variables removed
%

  e = 1e-30;
  if 1+e==1-e                           % fast check for rounding to nearest
    rndold = 0;
  else
    rndold = getround;
    setround(0)
  end

  N = getappdata(0,'INTLAB_GRADIENT_NUMVAR');

  if ~isa(a,'gradient')
    m = prod(size(a));
    if m==1                     % non-gradient scalar .* gradient
      r.x = a * b.x;
      r.dx = a * b.dx;
    else                        % non-gradient array .* gradient
      n = prod(size(b.x));
      if n==1                   % non-gradient array .* gradient scalar
        r.x = a * b.x;
        ax = a(:);
        if issparse(b.dx)
          ax = sparse(ax);
        end
        r.dx = ax * b.dx;
      else                      % non-gradient array .* gradient array
        if ~isequal(size(a),size(b))
          error('gradient .* : dimensions not compatible')
        end
        r.x = a .* b.x;
        if issparse(b.dx)
          [ib,jb,sb] = find(b.dx);
          r.dx = sparse(ib,jb,reshape(a(ib),size(sb)).*sb,n,N);
        else
          r.dx = b.dx .* repmat(a(:),1,N);
        end
      end
    end
  elseif ~isa(b,'gradient')      % gradient times non-gradient
    r = b .* a;
    if rndold
      setround(rndold)
    end
    return
  else                          % both factors gradient
    m = prod(size(a.x));
    n = prod(size(b.x));
    sparse_ = issparse(a.dx) | issparse(b.dx);
    if m==1                     % scalar gradient .* gradient
      if n==1                   % scalar gradient .* scalar gradient
        r.x = a.x * b.x;
        r.dx = a.dx * b.x + a.x * b.dx;
      else                      % scalar gradient .* array gradient
        r.x = a.x * b.x;
        if sparse_
          r.dx = sparse(b.x(:)) * a.dx + b.dx * a.x;
        else
          r.dx = b.x(:) * a.dx + b.dx * a.x;
        end
      end
    else                        % array gradient .* gradient
      if n==1                   % array gradient .* scalar gradient
        r.x = a.x * b.x;
        if sparse_
          r.dx = sparse(a.x(:)) * b.dx + a.dx * b.x;
        else
          r.dx = a.x(:) * b.dx + a.dx * b.x;
        end
      else                      % array gradient .* array gradient
        if ~isequal(size(a),size(b))
          error('dimensions not compatible for gradient .*');
        end
        r.x = a.x .* b.x;
        if issparse(a.dx) | issparse(b.dx)
          r.dx = sparse(1:n,1:n,b.x(:))*a.dx + sparse(1:n,1:n,a.x(:))*b.dx;
        else
          r.dx = repmat(b.x(:),1,N) .* a.dx + b.dx .* repmat(a.x(:),1,N);
        end
      end
    end
  end

  r = class(r,'gradient');

  if rndold
    setround(rndold)
  end

repo: OpenMAS


text:
% DATA = MOHSST5_LOADDATA()

% Last modified 2011-01-11
% Copyright (c) Jaakko Luttinen (jaakko.luttinen@tkk.fi)

function data = mohsst5_loaddata()

data = load('mohsst5_data.mat');

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
classdef PSO < ALGORITHM
% <single> <real/integer> <large/none> <constrained/none>
% Particle swarm optimization
% W --- 0.4 --- Inertia weight

%------------------------------- Reference --------------------------------
% R. Eberhart and J. Kennedy, A new optimizer using particle swarm theory,
% Proceedings of the International Symposium on Micro Machine and Human
% Science, 1995, 39-43.
%------------------------------- 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".
%--------------------------------------------------------------------------

    methods
        function main(Algorithm,Problem)
            %% Parameter setting
            W = Algorithm.ParameterSet(0.4);

            %% Generate random population
            Population = Problem.Initialization();
            Pbest      = Population;
            [~,best]   = min(FitnessSingle(Pbest));
            Gbest      = Pbest(best);

            %% Optimization
            while Algorithm.NotTerminated(Population)
                Population     = OperatorPSO(Problem,Population,Pbest,Gbest,W);
                replace        = FitnessSingle(Pbest) > FitnessSingle(Population);
                Pbest(replace) = Population(replace);
                [~,best]       = min(FitnessSingle(Pbest));
                Gbest          = Pbest(best);
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function isgn = i4vec2_compare ( n, a1, a2, i, j )

%*****************************************************************************80
%
%% I4VEC2_COMPARE compares two I4VEC2's.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 February 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of data items.
%
%    Input, integer A1(N), A2(N), contain the two components of each item.
%
%    Input, integer I, J, the items to be compared.
%
%    Output, integer ISGN, the results of the comparison:
%    -1, item I < item J,
%     0, item I = item J,
%    +1, item J < item I.
%
  isgn = 0;

       if ( a1(i) < a1(j) )

    isgn = -1;

  elseif ( a1(i) == a1(j) )

         if ( a2(i) < a2(j) )
      isgn = -1;
    elseif ( a2(i) < a2(j) )
      isgn = 0;
    elseif ( a2(j) < a2(i) )
      isgn = +1;
    end

  elseif ( a1(j) < a1(i) )

    isgn = +1;

  end

  return
end


text:
function asa006_test02 ( )

%*****************************************************************************80
%
%% TEST02 demonstrates the use of CHOLESKY.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 February 2008
%
%  Author:
%
%    John Burkardt
%
  n_max = 15;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST02:\n' );
  fprintf ( 1, '  CHOLESKY computes the Cholesky factorization\n' );
  fprintf ( 1, '  of a positive definite symmetric matrix.\n' );
  fprintf ( 1, '  A compressed storage format is used.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Here we look at the Hilbert matrix\n' );
  fprintf ( 1, '  A(I,J) = 1 / ( I + J - 1 )\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  We expect errors to grow quickly with N.\n' );

  for n = 1 : n_max

    nn = ( n * ( n + 1 ) ) / 2;
%
%  Set A to the Hilbert matrix.
%
    k = 0;
    for i = 1 : n
      for j = 1 : i
        k = k + 1;
        a(k) = 1.0 / ( i + j - 1 );
      end
    end

    [ u, nullty, ifault ] = cholesky ( a, n, nn );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Matrix order N = %d\n', n );
    fprintf ( 1, '  Maxtrix nullity NULLTY = %d\n', nullty );

    k = 0;
    for j = 1 : n
      for i = 1 : j
        k = k + 1;
        ufull(i,j) = u(k);
      end
      for i = j + 1 : n
        ufull(i,j) = 0.0;
      end
    end
%
%  Compute U' * U and compare to A.
%
    k = 0;
    diff = 0.0;
    for i = 1 : n
      for j = 1 : i
        k = k + 1;
        utu = 0.0;
        for l = 1 : n
          utu = utu + ufull(l,i) * ufull(l,j);
        end
        diff = diff + ( a(k) - utu ) * ( a(k) - utu );
      end
    end

    diff = sqrt ( diff );

    fprintf ( 1,  '  RMS ( A - U''*U ) = %f\n', diff );
  end

  return
end