text:
%% Read the images.

% J = load('/home/mayank/Work/FlyBowl/newfeatures_walk_2.jab','-mat');
%
% count = zeros(1,2);
% for ndxx = 1:numel(J.x.labels.t0s)
%   for bout = 1:numel(J.x.labels.t0s{ndxx})
%     if J.x.labels.t1s{ndxx}(bout)>1000, continue;end
%     fly = J.x.labels.flies(ndxx);
%     frames = J.x.labels.t0s{ndxx}(bout):J.x.labels.t1s{ndxx}(bout)-1;
%     bothim = {};
%     for tt = 1:2
%       if tt ==1
%         expdir = '/home/mayank/Work/FlySpaceTime/walkMovies/SS03500_test';
%         ftrname = 'DSs';
%       else
%         expdir = '/home/mayank/Work/FlySpaceTime/walkMovies/SS03500';
%         ftrname = 'hs_sups';
%       end

%%
clear tt;
expdir = '/home/mayank/Work/FlySpaceTime/walkMovies/SS03500_test';
ftrname = 'DSs';
% fly = 2;
% frames = 51:200;
fly = 4;
frames = 3342:3353;
% frames = 3532:3543;
% selftrs = [6 7 1; 2 4 1; 6 3 1];
selftrs = [3 7 1; 6 5 1; 3 3 1];
savefigure = false;

%% params.
params = getParams;
npatches = params.npatches;
psize = params.psize;
nbins = params.nbins;
patchsz = params.patchsz;
scale = params.scale;

optflowwinsig = params.optflowwinsig ;
optflowsig = params.optflowsig ;
optreliability = params.optreliability ;

maxflow = 5;
fly_thres = 90;
flow_thres = sqrt(2)/4;

[~,expname] = fileparts(expdir);
moviefilestr = 'movie.ufmf';
moviefile = fullfile(expdir,moviefilestr);
trxfilestr = 'trx.mat';
trxfile = fullfile(expdir,trxfilestr);

[readframe,nframes] = get_readframe_fcn(moviefile);
td = load(trxfile);
tracks = td.trx;
dx = []; dy = [];
theta = [];
rdx = [];rdy = [];
im1 = []; im2 = []; im3 = [];
for i = 1:numel(frames),
  curf = readframe(frames(i));
  nextf = readframe(frames(i)+1);

  trackndx = frames(i) - tracks(fly).firstframe + 1;
  locy = round(tracks(fly).y(trackndx));
  locx = round(tracks(fly).x(trackndx));
  curpatch = extractPatch(curf,...
    locy,locx,tracks(fly).theta(trackndx),patchsz);
  locy = round(tracks(fly).y(trackndx+1));
  locx = round(tracks(fly).x(trackndx+1));
  curpatch2 = extractPatch(nextf,...
    locy,locx,tracks(fly).theta(trackndx+1),patchsz);

  curpatch3 = extractPatch(nextf,...
    locy,locx,tracks(fly).theta(trackndx),patchsz);

  im1(:,:,i) = curpatch;
  im2(:,:,i) = curpatch2;
  im3(:,:,i) = curpatch3;
  dx(i) = round(tracks(fly).x(trackndx+1))-round(tracks(fly).x(trackndx));
  dy(i) = round(tracks(fly).y(trackndx+1))-round(tracks(fly).y(trackndx));
  rdx(i) = tracks(fly).x(trackndx+1)-tracks(fly).x(trackndx)-dx(i);
  rdy(i) = tracks(fly).y(trackndx+1)-tracks(fly).y(trackndx)-dy(i);
  dtheta(i) = tracks(fly).theta(trackndx+1)-tracks(fly).theta(trackndx);
  theta(i) = tracks(fly).theta(trackndx);
end

[nr,nc,~] = size(im1);
%
% figure out the histogram bins
tmptheta = (0:179)*pi/180;
res = nan(nbins,numel(tmptheta));
m = single(ones(10,10));
o = single(zeros(10,10));
for i = 1:numel(tmptheta),
%   fprintf('i = %d, theta = %f\n',i,tmptheta(i));
  o(:) = single(tmptheta(i));
  rescurr = gradientHist(m,o,1,nbins,1);
  res(:,i) = rescurr(1,1,:);
end

bincenters = nan(1,nbins);
for i = 1:nbins,
  bincenters(i) = tmptheta(argmax(res(i,:)));
end

% this seems to be what the centers correspond to
bincenters = linspace(0,pi,nbins+1);
bincenters = bincenters(1:nbins);

% mayank divides by 2
bincenters2 = bincenters*2;
dt = mean(diff(bincenters2));
binedges2 = [bincenters2(1)-dt/2,(bincenters2(1:end-1)+bincenters2(2:end))/2,bincenters2(end)+dt/2];

%% Read the features;

ftrVal = [];
for ndx = 1:size(selftrs,1)
  curname = sprintf('%s_%02d_%02d_%d.mat',ftrname,...
    selftrs(ndx,1),selftrs(ndx,2),selftrs(ndx,3));
  Q = load(fullfile(expdir,'perframe',curname));
  startndx = frames(1)-tracks(fly).firstframe+1;
  endndx  = frames(end)-tracks(fly).firstframe+1;
  ftrVal(ndx,:) = Q.data{fly}(startndx:endndx);
end
maxval = 2;
ftrVal(ftrVal>maxval) = maxval;
%% Visualize

minImg = min(im1,[],3);

nftrs = size(selftrs,1);
padsz = 50; margin = 4;
tsz = (padsz-2*margin)/nftrs;
padImg = padgrab(minImg,255,1,nc+padsz,1,nr,1,size(minImg,3));

if exist('tt','var')
  ff = figure(104+tt);
else
  ff = figure;
end
imshow(uint8(imresize(padImg,scale)));hold on;

xidx = linspace(margin,nc-margin,numel(frames));
cc = hsv(nftrs);
for ndx = 1:nftrs
  ybase = nc + (padsz-margin)/nftrs*ndx;

  plot(xidx*scale,(ybase-ftrVal(ndx,:)*tsz/maxval)*scale,'Color',cc(ndx,:));
  px = [0 psize psize 0 0] + (selftrs(ndx,2)-1)*psize;
  py = [0 0 psize psize 0] + (selftrs(ndx,1)-1)*psize;
  patch(px*scale,py*scale,cc(ndx,:),'FaceColor',cc(ndx,:),'FaceAlpha',0.075,'LineStyle','none');
end

for ndx = 1:npatches-1
  plot([0.5 nr*scale+0.5],[ndx*psize ndx*psize]*scale,'y');
  plot([ndx*psize ndx*psize]*scale,[0.5 nr*scale+0.5],'y');

end

% fd = getframe(ff);
% bothim{tt} = fd.cdata;
%%
%     end
%
%     if savefigure,
%       outimgname = sprintf('%s_%s_Fly%d_From%dTo%d_%s.png',...
%         J.x.labels.names{ndxx}{bout},...
%         expname,fly,frames(1),frames(end),...
%         datestr(now,'yyyymmdd'));
%       outname = fullfile('..','figures',outimgname);
%       imwrite(cat(2,bothim{1},bothim{2}),outname);
%     end
%
%   end
% end

author: kristinbranson

repo: JAABA


text:
function A = GregK316(n);
%GREGK316     Example from Gregory/Karney
%
% Gregory/Karney test case 3.16, n>=2
%
% (  2  -1               )
% ( -1   2  -1           )
% (     -1   2  -1       )
% (        ........      )
% (         -1   2   1   )
% (             -1   2   )
%
%   A = GregK316(n);
%

% written   3/01/95     S.M. Rump
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
%

  if n<2
    error('dimension too small')
  end
  e = ones(n,1);
  A = spdiags( [-e 2*e -e], -1:1, n, n);

repo: OpenMAS


text:
function preconfun = preconhessiansolve(problem, options)
% Preconditioner based on the inverse Hessian, by solving linear systems.
%
% function preconfun = preconhessiansolve(problem)
% function preconfun = preconhessiansolve(problem, options)
%
% Input:
%
% A Manopt problem structure (already containing the manifold and enough
% information to compute the Hessian of the cost) and an options structure
% (optional, currently ignored). Notice that if the Hessian is not positive
% definite, then its inverse is not positive definite either and this
% preconditioner is not suitable.
%
% If the Hessian cannot be computed on 'problem', a warning is issued. An
% approximation of the Hessian will be used instead, and the present
% preconditioner will attempt to invert that (although it may not be a
% linear operator). If no approximate Hessian is provided either, a generic
% approximation is used. Behavior is unspecified.
%
% Output:
%
% Returns a function handle, encapsulating a generic preconditioner of the
% Hessian based on solving linear systems of the form:
%   Hessian(x)[preconfun(x, xdot)] = xdot,
% where x is the point on the manifold, xdot is the input to the
% preconditioner (a tangent vector) and preconfun(x, xdot) is returned
% (also a tangent vector). The solve may be approximate.
%
% The returned preconfun has this calling pattern:
%
%   function precxdot = preconfun(x, xdot)
%   function precxdot = preconfun(x, xdot, storedb)
%   function precxdot = preconfun(x, xdot, storedb, key)
%
% x is a point on the manifold problem.M, xdot is a tangent vector to that
% manifold at x, storedb is a StoreDB object, and key is the StoreDB key to
% point x.
%
% Usage:
%
% Typically, the user will set problem.M and other fields to define the
% cost, the gradient and the Hessian (typically, problem.cost, problem.grad
% and problem.hess, or problem.egrad and problem.ehess). Then, to use this
% generic purpose Hessian preconditioner:
%
%   problem.precon = preconhessiansolve(problem, options);
%
% Passing that problem structure to the conjugategradients solver
% (which uses preconditioning) configured in steepest descent mode results
% in a type of Riemannian Newton method.
%
% See also: conjugategradients

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, April 9, 2015.
% Contributors:
% Change log:

    % Check availability of the Hessian, or at least of an approximation.
    if ~canGetHessian(problem) && ~canGetApproxHessian(problem)
        % Note: we do not give a warning if an approximate Hessian is
        % explicitly given in the problem description, as in that case the
        % user seems to be aware of the issue.
        warning('manopt:getHessian:approx', ...
               ['No Hessian provided. Using an FD approximation instead.\n' ...
                'To disable this warning: warning(''off'', ''manopt:getHessian:approx'')']);
        problem.approxhess = approxhessianFD(problem);
    end

    % Set local defaults here, and merge with user options, if any.
    localdefaults = struct();
    if ~exist('options', 'var') || isempty(options)
        options = struct();
    end
    options = mergeOptions(localdefaults, options);

    % Build and return the function handle here. This extra construct via
    % funhandle makes it possible to make storedb and key optional.
    preconfun = @funhandle;
    function precxdot = funhandle(x, xdot, storedb, key)
        % Allow omission of the key, and even of storedb.
        if ~exist('key', 'var')
            if ~exist('storedb', 'var')
                storedb = StoreDB();
            end
            key = storedb.getNewKey();
        end
        precxdot = hessiansolvehelper(options, problem, x, xdot, ...
                                      storedb, key);
    end

end

function precxdot = hessiansolvehelper(options, problem, x, xdot, storedb, key)
% This function does the actual work.

    % Exclude the case where xdot is zero
    norm_xdot = problem.M.norm(x, xdot);
    if norm_xdot < eps
        precxdot = problem.M.zerovec(x);
        return;
    end

    % Get a shorthand for the Hessian of the cost on M at x.
    hessian = @(u) getHessian(problem, x, u, storedb, key);

    % Setup an optimization problem on the tangent space to problem.M at x.
    M = problem.M;
    tgtspace = tangentspacefactory(M, x);
    prblm.M = tgtspace;
    prblm.cost = @cost;
    prblm.grad = @grad;
    prblm.hess = @(u, udot) 2*hessian(hessian(udot))/norm_xdot;

    function [f, store] = cost(u, store)
        if ~isfield(store, 'residue')
            Hu = hessian(u);
            store.residue = M.lincomb(x, 1, Hu, -1, xdot);
        end
        f = M.norm(x, store.residue).^2 / norm_xdot;
    end
    function [g, store] = grad(u, store)
        if ~isfield(store, 'residue')
            Hu = hessian(u);
            store.residue = M.lincomb(x, 1, Hu, -1, xdot);
        end
        g = 2 * hessian(store.residue) / norm_xdot;
    end

    % checkgradient(prblm); pause;
    % checkhessian(prblm); pause;

    localdefaults.solver = @trustregions;
    localdefaults.verbosity = 0;
    % Merge local defaults with user options, if any.
    if ~exist('options', 'var') || isempty(options)
        options = struct();
    end
    options = mergeOptions(localdefaults, options);

    % Solve the linear system by solving the optimization problem.
    precxdot = manoptsolve(prblm, M.zerovec(), options);

end

author: NicolasBoumal

repo: manopt


text:
%evaluate on the val2 set from imagenet2013
%use the default parameters as used in the demo:
%% load pre-trained edge detection model and set opts (see edgesDemo.m)
addpath(genpath('../'));

model=load('models/forest/modelBsds'); model=model.model;
model.opts.multiscale=0; model.opts.sharpen=2; model.opts.nThreads=4;

%% set up opts for edgeBoxes (see edgeBoxes.m)
opts = edgeBoxes;
opts.alpha = .65;     % step size of sliding window search
opts.beta  = .75;     % nms threshold for object proposals
opts.minScore = .01;  % min score of boxes to detect
opts.maxBoxes = 1e4;  % max number of boxes to detect

%config specifies the image location location to save output etc.
config=createConfig();
imageLoc=config.path.imageLoc;
saveLoc=config.path.outputLoc;
ext=config.opts.imageExt;

imageExt=config.opts.imageExt;
images=dir([imageLoc '*' imageExt]);

for i=1:length(images)

	imageName=images(i).name;
	im=imread([imageLoc imageName]);
	if(size(im, 3) == 1)
		im=repmat(im,[1,1,3]);
	end
	bbs=edgeBoxes(im,model,opts);
	if(isfield((config.opts),'numProposals'))
		numProposals=config.opts.numProposals;
	        if(size(bbs,1)>=numProposals)
        	        bbs=bbs(1:numProposals);
        	else
                	fprintf('Only %d proposals were generated for image: %s\n',size(bbs,1),imageName);
        	end
	end
 	%edges boxes produces baoxes as "[x y w, h]"
	%we convert to [x y x+w y+h]==[xmin ymin xmax ymax]
        boxes=bbs(:,1:4);
	boxes=[boxes(:,1) boxes(:,2) boxes(:,1)+ boxes(:,3) boxes(:,2)+boxes(:,4)];
	proposals.boxes= boxes;
	proposalFileName=strrep(imageName,imageExt,'.mat');

	save([saveLoc proposalFileName], 'proposals');
end

author: Cloud-CV

repo: object-proposals


text:
function a = oto ( m, n )

%*****************************************************************************80
%
%% OTO returns the OTO matrix.
%
%  Example:
%
%    N = 5
%
%    2  1  .  .  .
%    1  2  1  .  .
%    .  1  2  1  .
%    .  .  1  2  1
%    .  .  .  1  2
%
%  Properties:
%
%    A is banded, with bandwidth 3.
%
%    A is tridiagonal.
%
%    Because A is tridiagonal, it has property A (bipartite).
%
%    A is integral: int ( A ) = A.
%
%    A is Toeplitz: constant along diagonals.
%
%    A is symmetric: A' = A.
%
%    Because A is symmetric, it is normal.
%
%    Because A is normal, it is diagonalizable.
%
%    A is persymmetric: A(I,J) = A(N+1-J,N+1-I).
%
%    A is weakly diagonally dominant, but not strictly diagonally dominant.
%
%    The family of matrices is nested as a function of N.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns of A.
%
%    Output, real A(M,N), the matrix.
%
  a = zeros ( m, n );

  for i = 1 : m
    for j = 1 : n

      if ( j == i - 1 )
        a(i,j) = 1.0;
      elseif ( j == i )
        a(i,j) = 2.0;
      elseif ( j == i + 1 )
        a(i,j) = 1.0;
      else
        a(i,j) = 0.0;
      end

    end
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ point_coord, edge_point, face_order, face_point ] = ...
  icos_shape ( point_num, edge_num, face_num, face_order_max )

%*****************************************************************************80
%
%% ICOS_SHAPE describes an icosahedron.
%
%  Discussion:
%
%    The input data required for this routine can be retrieved from
%    ICOS_SIZE.
%
%    The vertices lie on the unit sphere.
%
%    The dual of an icosahedron is a dodecahedron.
%
%    The data has been rearranged from a previous assignment.
%    The STRIPACK program refuses to triangulate data if the first
%    three nodes are "collinear" on the sphere.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 July 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer POINT_NUM, the number of points (12).
%
%    Input, integer EDGE_NUM, the number of edges (30).
%
%    Input, integer FACE_NUM, the number of faces (20).
%
%    Input, integer FACE_ORDER_MAX, the maximum number of vertices per face (3).
%
%    Output, real POINT_COORD(3,POINT_NUM); the points.
%
%    Output, integer EDGE_POINT(2,EDGE_NUM), the points that make up each
%    edge, listed in ascending order of their indexes.
%
%    Output, integer FACE_ORDER(FACE_NUM), the number of vertices per face.
%
%    Output, integer FACE_POINT(FACE_ORDER_MAX,FACE_NUM); FACE_POINT(I,J)
%    is the index of the I-th point in the J-th face.  The points are listed
%    in the counter-clockwise direction defined by the outward normal at the
%    face.  The nodes of each face are ordered so that the lowest index
%    occurs first.  The faces are then sorted by nodes.
%
  dim_num = 3;
%
%  Set point coordinates.
%
  phi = 0.5 * ( sqrt ( 5.0 ) + 1.0 );
  b = 1.0 / sqrt ( 1.0 + phi * phi );
  a = phi * b;
  z = 0.0;
%
%  Set the points.
%
  point_coord(1:dim_num,1:point_num) = [ ...
      a,  b,  z; ...
      a, -b,  z; ...
      b,  z,  a; ...
      b,  z, -a; ...
      z,  a,  b; ...
      z,  a, -b; ...
      z, -a,  b; ...
      z, -a, -b; ...
     -b,  z,  a; ...
     -b,  z, -a; ...
     -a,  b,  z; ...
     -a, -b,  z ]';
%
%  Set the edges.
%
  edge_point(1:2,1:edge_num) = [ ...
     1,  2; ...
     1,  3; ...
     1,  4; ...
     1,  5; ...
     1,  6; ...
     2,  3; ...
     2,  4; ...
     2,  7; ...
     2,  8; ...
     3,  5; ...
     3,  7; ...
     3,  9; ...
     4,  6; ...
     4,  8; ...
     4, 10; ...
     5,  6; ...
     5,  9; ...
     5, 11; ...
     6, 10; ...
     6, 11; ...
     7,  8; ...
     7,  9; ...
     7, 12; ...
     8, 10; ...
     8, 12; ...
     9, 11; ...
     9, 12; ...
    10, 11; ...
    10, 12; ...
    11, 12 ]';
%
%  Set the face orders.
%
  face_order(1:face_num) = [ ...
    3, 3, 3, 3, 3, 3, 3, 3, 3, 3, ...
    3, 3, 3, 3, 3, 3, 3, 3, 3, 3 ]';
%
%  Set the faces.
%
  face_point(1:face_order_max,1:face_num) = [ ...
     1,  2,  4; ...
     1,  3,  2; ...
     1,  4,  6; ...
     1,  5,  3; ...
     1,  6,  5; ...
     2,  3,  7; ...
     2,  7,  8; ...
     2,  8,  4; ...
     3,  5,  9; ...
     3,  9,  7; ...
     4,  8, 10; ...
     4, 10,  6; ...
     5,  6, 11; ...
     5, 11,  9; ...
     6, 10, 11; ...
     7,  9, 12; ...
     7, 12,  8; ...
     8, 12, 10; ...
     9, 11, 12; ...
    10, 12, 11 ]';

  return
end


text:
args = {};
if strcmp(jobName, 'edges_to_ucms')
  % Uses the precomptued edge cues to do edge detection at multiple scales
  % This edges at multiple scales are used to generate a multiscale UCM.
  p = get_paths();
  modelfile = fullfile_ext(p.contours_model_dir, 'forest', 'modelNyuRgbd-3', 'mat');
  K = 40;
  for i = 1:K,
    args{i} = {sprintf('all_%d_%d_end', i, K), modelfile, p.ucm_dir, p.contours_cues_dir, cropCamera(getCameraParam('color'))};
  end
  jobParam = struct('numThreads', 1, 'codeDir', pwd(), 'preamble', '', 'matlabpoolN', 0, 'globalVars', {{}}, 'fHandle', @wrapper_rgbd_to_ucm, 'numOutputs', 0);
  resourceParam = struct('mem', 8, 'hh', 20, 'numJobs', K, 'ppn', 1, 'nodes', 1, 'logDir', '/work4/sgupta/pbsBatchDir/', 'queue', 'psi', 'notif', false, 'username', 'sgupta', 'headNode', 'zen');
  [jobId jobDir] = jobParallel(jobName, resourceParam, jobParam, args);
  while ~collectJob(jobDir), pause(120); end
end

if strcmp(jobName, 'benchmark_multi_ucm')
  %% Testing code
  p = get_paths();
  c = benchmarkPaths();
  jobParam = struct('numThreads', 1, 'codeDir', pwd(), 'preamble', '', 'matlabpoolN', 1, 'globalVars', {{}}, 'fHandle', @empty, 'numOutputs', 1);
  resourceParam = struct('mem', 2, 'hh', 1, 'numJobs', 40, 'ppn', 1, 'nodes', 1, 'logDir', '/work4/sgupta/pbsBatchDir/', 'queue', 'psi', 'notif', false, 'username', 'sgupta', 'headNode', 'zen');
  pDistrEval={'type', 'psi', 'pLaunch', struct('jobParam', jobParam, 'resourceParam', resourceParam)};
  edgesEvalDir('resDir', fullfile(p.ucm_dir, 'multi-png', filesep()), ...
    'gtDir', fullfile(c.benchmarkGtDir, filesep()), ...
    'pDistr', pDistrEval, 'cleanup', 0, 'thrs', 99, 'maxDist', 0.011, ...
    'ids', getImageSet('test'));
end

if strcmp(jobName, 'pareto')
  p = get_paths();
  params = nyud_params('root_cache_dir', p.cache_dir, 'feature_id', 'depth', 'depth_features', true, 'camera_matrix', cropCamera(getCameraParam('color')));
  num_regions = 20000;
  pareto_choose_point(params, num_regions);
end

if strcmp(jobName, 'cache-mcg')
  % Cache the base features for MCG
  jobParam = struct('numThreads', 1, 'codeDir', pwd(), 'preamble', '', 'matlabpoolN', 0, 'globalVars', {{}}, 'fHandle', @cache_mcg_features, 'numOutputs', 0);
  resourceParam = struct('mem', 4, 'hh', 20, 'numJobs', 80, 'ppn', 1, 'nodes', 1, 'logDir', '/work4/sgupta/pbsBatchDir/', 'queue', 'psi', 'notif', false, 'username', 'sgupta', 'headNode', 'zen');

  params = nyud_params('root_cache_dir', p.cache_dir, 'feature_id', 'depth', 'depth_features', true, 'camera_matrix', cropCamera(getCameraParam('color')));

  imNames = getImageSet('trainval');
  for i = 1:length(imNames), args{i} = {params, [], imNames{i}, params.cacheDir}; end

  [jobId jobDir] = jobParallel(jobName, resourceParam, jobParam, args);
  while ~collectJob(jobDir), pause(120); end
end

if strcmp(jobName, 'rank_training'),
  p = get_paths();
  params = nyud_params('root_cache_dir', p.cache_dir, 'feature_id', 'depth', 'depth_features', true, 'camera_matrix', cropCamera(getCameraParam('color')));
  rank_training(params, params.cacheDir);
end

if strcmp(jobName, 'region-detect')
  p = get_paths();
  params = nyud_params('root_cache_dir', p.cache_dir, 'feature_id', 'depth', 'depth_features', true, 'camera_matrix', cropCamera(getCameraParam('color')));
  imlist = getImageSet('all');
  rf = loadvar(params.files.trained_classifier,'rf');
  n_cands = loadvar(params.files.pareto_point,'n_cands');

  parfor i = 1:length(imlist),
    % Load the UCMs here
    imname = imlist{i};
    ucms = {};
    for ii = 1:length(params.hier_dirs),
      tmp = load(fullfile_ext(params.hier_dirs{ii}, imname, 'mat'));
      ucms{ii} = tmp.ucm2;
    end
    D = getImage(imname, 'depth');
    RD = getImage(imname, 'rawdepth');
    mcg_cache_obj = cache_mcg_features(params, ucms, imname, params.cacheDir);
    candidates = compute_mcg_cands(params, rf, n_cands, mcg_cache_obj, D, RD);
    % Save regions here!
    out_file = fullfile_ext(params.files.output_dir, imname, 'mat');
    parsave(out_file, 'scores', candidates.scores, 'bboxes', candidates.bboxes, 'superpixels', candidates.superpixels, 'labels', candidates.labels);
  end
  stats = eval_labels(params.files.output_dir, 'nyud40Obj', 'test');
  save(sprintf('%s-eval.mat', params.files.output_dir), 'stats', 'params');
end

if strcmp(jobName, 'show_region_results')
  p = get_paths();
  params = nyud_params('root_cache_dir', p.cache_dir, 'feature_id', 'depth', 'depth_features', true, 'camera_matrix', cropCamera(getCameraParam('color')));
  method{1} = params.mcg_id; method{2} = 'old-depth'; method{3} = 'old-nodepth';
  % params = nyud_params('feature_id', 'no_depth', 'depth_features', false, 'camera_matrix', cropCamera(getCameraParam('color')));
  % method{2} = params.mcg_id;
  legends = {
    ['Our (MCG (RGBD edges,' char(10) 'RGBD feats.)) [RGBD]'], ...
    ['old - depth'], ...
    ['old - nodepth'], ...
  };
  show_region_results(fullfile(p.output_dir, 'regions'), method, legends);
end

author: s-gupta

repo: rcnn-depth

save_path: github-repos/MATLAB/s-gupta-rcnn-depth


text:
function varargout = createCube()
%CREATECUBE Create a 3D mesh representing the unit cube
%
%   [V, E, F] = createCube
%   Create a unit cube, as a polyhedra representation.
%   c has the form [V E F], where V is a 8-by-3 array with vertices
%   coordinates, E is a 12-by-2 array containing indices of neighbour
%   vertices, and F is a 6-by-4 array containing vertices array of each
%   face.
%
%   [V, F] = createCube;
%   Returns only the vertices and the face vertex indices.
%
%   MESH = createCube;
%   Returns the data as a mesh structure, with fields 'vertices', 'edges'
%   and 'faces'.
%
%   Example
%   [n, e, f] = createCube;
%   drawMesh(n, f);
%
%   See also
%   meshes3d, drawMesh
%   createOctahedron, createTetrahedron, createDodecahedron
%   createIcosahedron, createCubeOctahedron
%

% ---------
% author : David Legland
% e-mail: david.legland@inra.fr
% INRA - TPV URPOI - BIA IMASTE
% created the 10/02/2005.

%   HISTORY
%   04/01/2007: remove unused variables

x0 = 0; dx= 1;
y0 = 0; dy= 1;
z0 = 0; dz= 1;

nodes = [...
    x0 y0 z0; ...
    x0+dx y0 z0; ...
    x0 y0+dy z0; ...
    x0+dx y0+dy z0; ...
    x0 y0 z0+dz; ...
    x0+dx y0 z0+dz; ...
    x0 y0+dy z0+dz; ...
    x0+dx y0+dy z0+dz];

edges = [1 2;1 3;1 5;2 4;2 6;3 4;3 7;4 8;5 6;5 7;6 8;7 8];

% faces are oriented such that normals point outwards
faces = [1 3 4 2;5 6 8 7;2 4 8 6;1 5 7 3;1 2 6 5;3 7 8 4];

% format output
varargout = formatMeshOutput(nargout, nodes, edges, faces);

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
% TOPOIMAGE - plot concatenated multichannel time/frequency images
%               in a topographic format
%               Uses a channel location file with the same format as TOPOPLOT
%               or else plots data on a rectangular grid of axes.
%               Click on individual images to examine separately.
%
% Usage:
%    >> topoimage(data,'chan_locs',ntimes,limits);
%    >> topoimage(data,[rows cols],ntimes,limits);
%    >> topoimage(data,'chan_locs',ntimes,limits,title,...
%                                 channels,axsize,colors,ydir,rmbase)
%
% Inputs:
%   data       = data consisting of nchans images, each size (rows,ntimes*chans)
%  'chan_locs' = file of channel locations as in >> topoplot example
%                Else [rows cols] matrix of locations. Example: [6 4]
%   ntimes     = columns per image
%  [limits]    = [mintime maxtime minfreq maxfreq mincaxis maxcaxis]
%                Give times in msec {default|0 (|both caxis 0) -> use data limits)
%  'title'     = plot title {0 -> none}
%   channels   = vector of channel numbers to plot & label {0 -> all}
%   axsize     = [x y] axis size {default [.08 .07]}
%  'colors'    = file of color codes, 3 chars per line
%                 ( '.' = space) {0 -> default color order}
%   ydir       = y-axis polarity (pos-up = 1; neg-up = -1) {def -> pos-up}
%   rmbase     = if ~=0, remove the mean value for times<=0 for each freq {def -> no}
%
% Author: Scott Makeig, SCCN/INC/UCSD, La Jolla, 12-10-1999
%
% See also: TOPOPLOT, TIMEF

% Copyright (C) 12-10-99 Scott Makeig, SCCN/INC/UCSD, scott@sccn.ucsd.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

% 1-16-00  debugged help msg and improved presentation -sm
% 3-16-00  added AXCOPY -sm
% 8-07-00  added LOGIMAGESC via 'LOGIT' -sm
% 9-02-00  added RMBASE option below, plus colorbar to key image -sm ???
% 1-25-02  reformated help & license, added link -ad

function topoimage(data,loc_file,times,limits,plottitle,channels,axsize,colors,ydr,rmbas)

% Options:
% LOGIT = 1; % comment out for non-log imaging
% YVAL = 10;          % plot horizontal lines at 10 Hz (comment to omit)
% RMBASE = 0;  % remove <0 mean for each image row

MAXCHANS = 256;
DEFAULT_AXWIDTH  = 0.08;
DEFAULT_AXHEIGHT = 0.07;
DEFAULT_SIGN = 1;   % Default - plot positive-up
LINEWIDTH = 2.0;
FONTSIZE = 14;      % font size to use for labels
CHANFONTSIZE = 10;  % font size to use for channel names
TICKFONTSIZE=10;    % font size to use for axis labels
TITLEFONTSIZE = 16;

PLOT_WIDTH = 0.75;  % width and height of plot array on figure!
PLOT_HEIGHT = 0.81;
ISRECT = 0;         % default

if nargin < 1,
    help topoimage
    return
end

if nargin < 4,
    help topoimage
    error('topoimage(): needs four arguments');
end

if times <0,
    help topoimage
    return
elseif times==1,
    fprintf('topoimage: cannot plot less than 2 times per image.\n');
    return
else
    freqs = 0;
end

axcolor= get(0,'DefaultAxesXcolor'); % find what the default x-axis color is
plotfile = 'topoimage.ps';
ls_plotfile = 'ls -l topoimage.ps';

%
%%%%%%%%%%%%%%%%%%%%%%%%%% Substitute defaults for missing parameters %%%%%
%
SIGN = DEFAULT_SIGN;
if nargin < 10
   rmbas = 0;
end
if nargin < 9
   ydr = 0;
end
if ydr == -1
   SIGN = -1;
end

if nargin < 8
    colors = 0;
end

if nargin < 7,
  axwidth  = DEFAULT_AXWIDTH;
  axheight = DEFAULT_AXHEIGHT;
elseif all(size(axsize) == [1 1]) && axsize(1) == 0
  axwidth  = DEFAULT_AXWIDTH;
  axheight = DEFAULT_AXHEIGHT;
elseif size(axsize) == [1 2]
  axwidth  = axsize(1);
  axheight = axsize(2);
  if axwidth > 1 || axwidth < 0 || axheight > 1 || axwidth < 0
    help topoimage
    return
  end
else
  help topoimage
  return
end

[freqs,framestotal]=size(data);             % data size
chans = framestotal/times;

fprintf('\nPlotting data using axis size [%g,%g]\n',axwidth,axheight);
if nargin < 6
   channels = 0;
end
if channels == 0
   channels = 1:chans;
end
if nargin < 5
    plottitle = 0; %CJH
end
limitset = 0;
if nargin < 4,
    limits = 0;
elseif length(limits)>1
    limitset = 1;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%% Test parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
  icadefs; % read MAXPLOTDATACHANS constant from icadefs.m

  if max(channels) > chans
    fprintf('topoimage(): max channel index > %d channels in data.\n',...
                       chans);
    return
  end
  if min(channels) < 1
    fprintf('topoimage(): min channel index (%g) < 1.\n',...
                       min(channels));
    return
  end
  if length(channels)>MAXPLOTDATACHANS,
    fprintf('topoimage(): not set up to plot more than %d channels.\n',...
                       MAXPLOTDATACHANS);
    return
  end
%
%%%%%%%%%%%%% Extend the size of the plotting area in the window %%%%%%%%%%%%
%
  curfig = gcf;
  h=figure(curfig);
  set(h,'Color',BACKCOLOR); % set the background color
  set(h,'PaperUnits','normalized'); % use percentages to avoid US/A4 difference
  set(h,'PaperPosition',[0.0235308 0.0272775 0.894169 0.909249]); % equivalent
  % orient portrait
  axis('normal');
%
%%%%%%%%%%%%%%%%%%%% Read the channel names %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
  if ischar(channels) == 0,
    % channames = zeros(MAXPLOTDATACHANS,4);
    % for c=1:length(channels),
    %     channames(c,:)= sprintf('%4d',channels(c));
    % end
    channames = num2str(channels(:));                   %%CJH
  else,
    if ~ischar(channels)
       fprintf('topoimage(): channel file name must be a string.\n');
       return
    end
    chid = fopen(channels,'r');
    if chid <3,
        fprintf('topoimage(): cannot open file %s.\n',channels);
        return
    end
    channames = fscanf(chid,'%s',[4 MAXPLOTDATACHANS]);
    channames = channames';
       [r c] = size(channames);
    for i=1:r
        for j=1:c
            if channames(i,j)=='.',
                channames(i,j)=' ';
            end
        end
    end
  end; % setting channames
%
%%%%%%%%%%%%%%%%%%%%%%%%% Plot and label specified channels %%%%%%%%%%%%%%%%%%
%
data = matsel(data,times,0,0,channels);
chans = length(channels);
%
%%%%%%%%%%%%%%%%%%%%%%%%% Read the color names %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
  if colors ~=0,
    if ~ischar(colors)
       fprintf('topoimage(): color file name must be a string.\n');
       return
    end
    cid = fopen(colors,'r');
    % fprintf('cid = %d\n',cid);
    if cid <3,
        fprintf('topoimage: cannot open file %s.\n',colors);
        return
    end
    colors = fscanf(cid,'%s',[3 MAXPLOTDATAEPOCHS]);
    colors = colors';
       [r c] = size(colors);
    for i=1:r
        for j=1:c
            if colors(i,j)=='.',
                colors(i,j)=' ';
            end
        end
    end
  else % use default color order (no yellow!)
     colors =['r  ';'b  ';'g  ';'c  ';'m  ';'r  ';'b  ';'g  ';'c  ';'m  ';'r  ';'b  ';'g  ';'c  ';'m  ';'r  ';'b  ';'g  ';'c  ';'m  ';'r  ';'b  ';'g  ';'c  ';'m  ';'r  ';'b  ';'g  ';'c  ';'m  ';'r  ';'b  ';'g  ';'c  ';'m  '];
     colors = [colors; colors];  % make > 64 available
  end
  for c=1:length(colors)   % make white traces black unless axis color is white
    if colors(c,1)=='w' && axcolor~=[1 1 1]
         colors(c,1)='k';
    end
  end
%
%%%%%%%%%%%%%%%%%%%%%%% Read and adjust limits %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
  if limits==0,      % == 0 or [0 0 0 0]
    xmin=min(times);
    xmax=max(times);
    ymin=min(freqs);
    ymax=max(freqs);
  else
    if length(limits)~=6,
      fprintf( ...
       'topoimage: limits should be 0 or an array [xmin xmax ymin ymax zmin zmax].\n');
      return
    end
    xmin = limits(1);
    xmax = limits(2);
    ymin = limits(3);
    ymax = limits(4);
    zmin = limits(5);
    zmax = limits(6);
  end

  if xmax == 0 && xmin == 0,
    x = [0:1:times-1];
    xmin = min(x);
    xmax = max(x);
  else
    dx = (xmax-xmin)/(times-1);
    x=xmin*ones(1,times)+dx*(0:times-1); % compute x-values
    xmax = xmax*times/times;
  end
  if xmax<=xmin,
      fprintf('topoimage() - xmax must be > xmin.\n')
      return
  end

  if ymax == 0 && ymin == 0,
      y=[1:1:freqs];
      ymax=freqs;
      ymin=1;
  else
    dy = (ymax-ymin)/(freqs-1);
    y=ymin*ones(1,freqs)+dy*(0:freqs-1); % compute y-values
    ymax = max(y);
  end
  if ymax<=ymin,
      fprintf('topoimage() - ymax must be > ymin.\n')
      return
  end

  if zmax == 0 && zmin == 0,
      zmax=max(max(data));
      zmin=min(min(data));
      fprintf('Color axis limits [%g,%g]\n',zmin,zmax);
  end
  if zmax<=zmin,
      fprintf('topoimage() - zmax must be > zmin.\n')
      return
  end

  xlabel = 'Time (ms)';
  ylabel = 'Hz';

%
%%%%%%%%%%%%%%%%%%%%%%%% Set up plotting environment %%%%%%%%%%%%%%%%%%%%%%%%%
%
  h = gcf;
  % set(h,'YLim',[ymin ymax]);       % set default plotting parameters
  % set(h,'XLim',[xmin xmax]);
  % set(h,'FontSize',18);
  % set(h,'DefaultLineLineWidth',1); % for thinner postscript lines
%
%%%%%%%%%%%%%%%%%%%%%%%%%% Print plot info %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
  % clf;   % clear the current figure

  % print plottitle over (left) subplot 1
  if plottitle==0,
    plottitle = '';
  end
  h=gca;title(plottitle,'FontSize',TITLEFONTSIZE); % title plot and
  hold on
  msg = ['\nPlotting %d traces of %d frames with colors: '];

  msg = [msg ' -> \n'];    % print starting info on screen . . .
  fprintf(...
    '\nlimits: [xmin,xmax,ymin,ymax] = [%4.1f %4.1f %4.2f %4.2f]\n',...
                xmin,xmax,ymin,ymax);

  set(h,'YLim',[ymin ymax]);            % set default plotting parameters
  set(h,'XLim',[xmin xmax]);
  set(h,'FontSize',FONTSIZE);            % choose font size

  set(h,'FontSize',FONTSIZE);            % choose font size
  set(h,'YLim',[ymin ymax]);            % set default plotting parameters
  set(h,'XLim',[xmin xmax]);

  axis('off')

%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Read chan_locs %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if size(loc_file,2) == 2 % plot in a rectangular grid
   ISRECT = 1;
   ht = loc_file(1);
   wd = loc_file(2);
   if chans > ht*wd
      fprintf(...
         '\ntopoimage(): (d%) channels to be plotted > grid size [%d %d]\n\n',...
                           chans,ht,wd);
      return
   end
   hht = (ht-1)/2;
   hwd = (wd-1)/2;
   xvals = zeros(ht*wd,1);
   yvals = zeros(ht*wd,1);
   dist  = zeros(ht*wd,1);
   for i=1:wd
    for j=1:ht
      xvals(i+(j-1)*wd) = -hwd+(i-1);
      yvals(i+(j-1)*wd) =  hht-(j-1);
       dist(i+(j-1)*wd) =  sqrt(xvals(j+(i-1)*ht).^2+yvals(j+(i-1)*ht).^2);
    end
   end
   maxdist = max(dist);
   for i=1:wd
    for j=1:ht
      xvals(i+(j-1)*wd) = 0.499*xvals(i+(j-1)*wd)/maxdist;
      yvals(i+(j-1)*wd) = 0.499*yvals(i+(j-1)*wd)/maxdist;
    end
   end
   channames = repmat(' ',ht*wd,4);
   for i=1:ht*wd
     channum = num2str(i);
     channames(i,1:length(channum)) = channum;
   end

else % read chan_locs file
  fid = fopen(loc_file);
  if fid<1,
    fprintf('topoimage(): cannot open eloc_file "%s"\n',loc_file)
    return
  end
  A = fscanf(fid,'%d %f %f %s',[7 MAXCHANS]);
  fclose(fid);
  A = A';

 if length(channels) > size(A,1),
   error('topoimage(): data channels must be <= chan_locs channels')
 end

  channames = setstr(A(channels,4:7));
  idx = find(channames == '.');                       % some labels have dots
  channames(idx) = setstr(abs(' ')*ones(size(idx)));  % replace them with spaces

  Th = pi/180*A(channels,2);                          % convert degrees to rads
  Rd = A(channels,3);
  % ii = find(Rd <= 0.5); % interpolate on-head channels only
  % Th = Th(ii);
  % Rd = Rd(ii);

  [yvals,xvals] = pol2cart(Th,Rd);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
xvals = 0.5+PLOT_WIDTH*xvals;   % controls width of plot array on page!
yvals = 0.5+PLOT_HEIGHT*yvals;  % controls height of plot array on page!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Plot traces %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
  xdiff=xmax-xmin;
  rightmost = max(xvals);
  basetimes = find(x<=0);
P=0;
  Axes = [];
    fprintf('\ntrace %d: ',P+1);
    for I=1:chans,%%%%%%%%%% for each data channel %%%%%%%%%%%%%%%%%%%%%%%%%%
      if P>0
        axes(Axes(I))
        hold on;                      % plot down left side of page first
        axis('off')

      else % P <= 0

      %
      %%%%%%%%%%%%%%%%%%%%%%% Plot data images %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      %
        xcenter = xvals(I);
        ycenter = yvals(I);
        Axes = [Axes axes('Units','Normal','Position', ...
              [xcenter-axwidth/2 ycenter-axheight/2 axwidth axheight])];
        axes(Axes(I))
        imageaxes = gca;
        axislcolor = get(gca,'Xcolor');   %%CJH

        dataimage = matsel(data,times,0,0,I);
        if rmbas~=0                  % rm baseline
          dataimage = dataimage ...
             - repmat(mean(matsel(data,times,basetimes,0,I)')',1,times);
        end

        if exist('LOGIT')
          logimagesc(x,y,dataimage); % <---- plot logfreq image
          if exist('YVAL')
             YVAL = log(YVAL);
          end
        else
          imagesc(x,y,dataimage);    % <---- plot image
        end
        hold on

      curax = axis;
      xtk  = get(gca,'xtick'); % use these for cal axes below
      xtkl = get(gca,'xticklabel');
      ytk  = get(gca,'ytick');
      ytkl = get(gca,'yticklabel');

      set(gca,'tickdir','out');
      set(gca,'ticklength',[0.02 0.05]);
      set(gca,'xticklabel',[]);
      set(gca,'yticklabel',[]);
      set(gca,'ydir','normal');
      caxis([zmin zmax]);

        if exist('YVAL') && YVAL>=curax(3) && YVAL<=curax(4)
            hold on
            hp=plot([xmin xmax],[YVAL YVAL],'r-');%,'color',axislcolor);
                                                       % draw horizontal axis
            set(hp,'Linewidth',1.0)
        end

        if xmin<0 && xmax>0
          hold on
          vl= plot([0 0],[curax(3) curax(4)],'color',axislcolor);  % draw vert axis
          set(vl,'linewidth',2);
        end

       % if xcenter == rightmost
       %    colorbar
       %    rightmost = Inf;
       % end

       % secondx = 200;                              % draw second vert axis
       % axis('off');plot([secondx secondx],[ymin ymax],'color',axislcolor);

       %
       %%%%%%%%%%%%%%%%%%%%%%% Print channel names %%%%%%%%%%%%%%%%%%%%%%%%%%
       %
       NAME_OFFSET = 0.01;
       if channels~=0,                               % print channames
          if ~ISRECT    % print before topographically arrayed image

             % axis('off');
             hold on
             h=text(xmin-NAME_OFFSET*xdiff,(curax(4)+curax(3))*0.5,[channames(I,:)]);
             set(h,'HorizontalAlignment','right');
             set(h,'FontSize',CHANFONTSIZE);                % choose font size

          else % print before rectangularly arrayed image
            if xmin<0
               xmn = 0;
            else
               xmn = xmin;
            end
            % axis('off');
            h=text(xmin-NAME_OFFSET*xdiff,ymax,[channames(I,:)]);
            set(h,'HorizontalAlignment','right');
            set(h,'FontSize',TICKFONTSIZE);                % choose font size
          end
       end; % channels
      end; % P=0

      % if xcenter == rightmost
      %    colorbar
      %    rightmost = Inf;
      % end

      % if xmin<0 & xmax>0
       % axes(imageaxes);
       % hold on; plot([0 0],[curax(3) curax(4)],'k','linewidth',2);
      % end

      drawnow
      fprintf(' %d',I);
    end; % %%%%%%%%%%%%%%% chan I %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  fprintf('\n');

%
%%%%%%%%%%%%%%%%%%%%% Make time and freq cal axis %%%%%%%%%%%%%%%%%%%%%%%%%
%
      ax = axes('Units','Normal','Position', ...
                         [0.80 0.1 axwidth axheight]);
      axes(ax)
      axis('off');
      imagesc(x,y,zeros(size(dataimage))); hold on   % <---- plot green
      caxis([zmin zmax]);
      set(gca,'ydir','normal');
      if xmin <=0
          py=plot([0 0],[curax(3) curax(4)],'color','k'); % draw vert axis at time zero
      else
          py=plot([xmin xmin],[curax(3) curax(4)],'color','k'); % vert axis at xmin
      end
      hold on
      if exist('YVAL') && YVAL>=curax(3)
          px=plot([xmin xmax],[YVAL YVAL],'color',axislcolor);
                                                      % draw horiz axis at YVAL
      else
          px=plot([xmin xmax],[curax(3) curax(4)],'color',axislcolor);
                                                      % draw horiz axis at ymin
      end

      axis(curax);
      set(gca,'xtick',xtk); % use these for cal axes
      set(gca,'xticklabel',xtkl);
      set(gca,'ytick',ytk);
      set(gca,'yticklabel',ytkl);
      set(gca,'ticklength',[0.02 0.05]);
      set(gca,'tickdir','out');
      h = colorbar;
      cbp = get(h,'position');
      set(h,'position',[cbp(1) cbp(2) 2*cbp(3) cbp(4)]);
      caxis([zmin zmax]);

      % secondx = 200;                    % draw second vert axis
      % axis('off');plot([secondx secondx],[curax(3) curax(4)],'color',axislcolor);

      %
      %%%%%%%%%%%%%%%%%%%%% Plot axis values %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      %

    if 0 % DETOUR
          signx = xmin-0.15*xdiff;
          axis('off');h=text(signx,SIGN*curax(3),num2str(curax(3),3));
            set(h,'FontSize',TICKFONTSIZE);         % choose font size
            set(h,'HorizontalAlignment','right','Clipping','off');

          textx = xmin-0.6*xdiff;
          axis('off');h=text(textx,(curax(3)+curax(4))/2,ylabel);        % text Hz
            set(h,'Rotation',90);
            set(h,'FontSize',TICKFONTSIZE);         % choose font size
            set(h,'HorizontalAlignment','center','Clipping','off');

        % axis('off');h=text(signx,SIGN*ymax,['+' num2str(ymax,3)]);  % text +ymax
          axis('off');h=text(signx,SIGN*ymax,[    num2str(ymax,3)]);  % text  ymax
            set(h,'FontSize',TICKFONTSIZE);         % choose font size
            set(h,'HorizontalAlignment','right','Clipping','off');

          ytick = curax(3)-0.3*(curax(4)-curax(3));
          tick = [int2str(xmin)]; h=text(xmin,ytick,tick);            % text xmin
            set(h,'FontSize',TICKFONTSIZE);         % choose font size
            set(h,'HorizontalAlignment','center',...
                      'Clipping','off');  % center text

           h=text(xmin+xdiff/2,ytick-0.5*(curax(4)-curax(3)),xlabel);% text Times
             set(h,'FontSize',TICKFONTSIZE);         % choose font size
             set(h,'HorizontalAlignment','center',...
                       'Clipping','off');  % center text

          tick = [int2str(xmax)]; h=text(xmax,ytick,tick);            % text xmax
            set(h,'FontSize',TICKFONTSIZE);         % choose font size
            set(h,'HorizontalAlignment','center',...
                      'Clipping','off');  % center text
          axis on
          set(ax,'xticklabel','');
          set(ax,'yticklabel','');
          set(ax,'ticklength',[0.02 0.05]);
          set(ax,'tickdir','out');
          caxis([zmin zmax]);
          hc=colorbar;
          cmapsize = size(colormap,1);
          set(hc,'ytick',[1 cmapsize]); %

          minlabel = num2str(zmin,3);
          while (length(minlabel)<4)
             if ~contains(minlabel,'.')
                  minlabel = [minlabel '.'];
             else
                  minlabel = [minlabel '0'];
             end
          end
          maxlabel = num2str(zmax,3);
          if zmin<0 && zmax>0
             maxlabel = ['+' maxlabel];
          end
          while (length(maxlabel)<length(minlabel))
             if ~contains(maxlabel,'.')
                  maxlabel = [maxlabel '.'];
             else
                  maxlabel = [maxlabel '0'];
             end
          end
          while (length(maxlabel)>length(minlabel))
             if ~contains(minlabel,'.')
                  minlabel = [minlabel '.'];
             else
                  minlabel = [minlabel '0'];
             end
          end
          set(hc,'yticklabel',[minlabel;maxlabel]);
          set(hc,'Color',BACKCOLOR);
          set(hc,'Zcolor',BACKCOLOR);
      end % DETOUR

axcopy(gcf); % turn on pop-up axes
%
%%%%%%%%%%%%%%%%%% Make printed figure fill page %%%%%%%%%%%%%%%%%%%%%%%%%%%
%
  % orient tall
  % curfig = gcf;
  % h=figure(curfig);
  % set(h,'PaperPosition',[0.2 0.3 7.6 10]); % stretch out the plot on the page

function [returnval] = contains(strng,chr)
   returnval=0;
   for i=1:length(strng)
     if strng(i)==chr
         returnval=1;
         break
     end
   end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function gX = mlpKernGradX(kern, X, X2)

% MLPKERNGRADX Gradient of MLP kernel with respect to input locations.
% FORMAT
% DESC computes the gradident of the multi-layer perceptron
% kernel with respect to the input positions where both the row
% positions and column positions are provided separately.
% ARG kern : kernel structure for which gradients are being
% computed.
% ARG x1 : row locations against which gradients are being computed.
% ARG x2 : column locations against which gradients are being computed.
% RETURN g : the returned gradients. The gradients are returned in
% a matrix which is numData2 x numInputs x numData1. Where numData1 is
% the number of data points in X1, numData2 is the number of data
% points in X2 and numInputs is the number of input
% dimensions in X.
%
% SEEALSO mlpKernParamInit, kernGradX, mlpKernDiagGradX
%
% COPYRIGHT : Neil D. Lawrence, 2004, 2005, 2006, 2009

% KERN

gX = zeros(size(X2, 1), size(X2, 2), size(X, 1));
for i = 1:size(X, 1);
  gX(:, :, i) = mlpKernGradXpoint(kern, X(i, :), X2);
end

function gX = mlpKernGradXpoint(kern, x, X2)

% MLPKERNGRADXPOINT Gradient with respect to one point of x.

innerProd = X2*x';
numer = innerProd*kern.weightVariance + kern.biasVariance;
vec1 = sum(x.*x, 2)*kern.weightVariance + kern.biasVariance + 1;
vec2 = sum(X2.*X2, 2)*kern.weightVariance + kern.biasVariance + 1;
denom = sqrt(vec2*vec1');
arg = numer./denom;
gX = zeros(size(X2));
twooverpi = 2/pi;
for j = 1:size(X2, 2)
  gX(:, j)=X2(:, j)./denom - vec2.*x(:, j).*numer./denom.^3;
  gX(:, j) = twooverpi*kern.weightVariance*kern.variance*gX(:, j)./sqrt(1-arg.*arg);
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function tet_mesh_test006 ( )

%*****************************************************************************80
%
%% TET_MESH_TEST006 tests TET_MESH_TET_NEIGHBORS.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 August 2009
%
%  Author:
%
%    John Burkardt
%
  tet_order = 4;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TET_MESH_TEST006\n' );
  fprintf ( 1, '  TET_MESH_TET_NEIGHBORS computes the 4 neighboring\n' );
  fprintf ( 1, '  tetrahedrons of each tetrahedron in a tet mesh.\n' );
  fprintf ( 1, '  containing a point.\n' );
%
%  Set up the example tetrahedron mesh.
%
  [ node_num, tet_num ] = tet_mesh_order4_example_size ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  This mesh has tetrahedron order %d\n', tet_order );
  fprintf ( 1, '  The number of tetrahedrons is   %d\n', tet_num );

  [ node_xyz, tet_node ] = tet_mesh_order4_example_set ( node_num, tet_num );
%
%  Print the tets.
%
  i4mat_transpose_print_some ( tet_order, tet_num, tet_node, ...
    1, 1, tet_order, 10, '  First 10 Tets:' );
%
%  The TET_NEIGHBOR array is needed by TET_MESH_DELAUNAY_SEARCH.
%
  tet_neighbor = tet_mesh_neighbor_tets ( tet_order, tet_num, tet_node );

  i4mat_transpose_print_some ( 4, tet_num, tet_neighbor, ...
    1, 1, 4, 10, '  First 10 Tet Neighbors:' );

  return
end


text:
function outcell = parse_char_to_cell(invar,sepval)
% take a row of characters and separate into cells, breaking at either
% spaces or tabs
%
% :Usage:
% :Usage:
% ::
%
%     outcell = parse_char_to_cell(invar,sepval)
%
% :Examples:
% ::
%
%    % copy from Excel as row, then parse:
%    disorder = ['SAD	PTSD	PTSD	PTSD	PTSD	PTSD	SP	SAD	PTSD	PTSD	PTSD	PTSD	PTSD	SAD	SAD	PTSD	PTSD	SP	SP	SP	PTSD	PTSD	PTSD	PTSD	SAD	SAD	SAD	SP	SAD	SAD	SAD	SP	SP	SP	SAD	SP	PTSD	SP	PTSD	PTSD'];
%    disorder = parse_char_to_cell(disorder, 'tab');
%
% ..
%    tor wager
% ..

switch sepval
    case 'space'
whsep = find(invar == ' ');
    case 'tab'
whsep = find(invar == sprintf('\t'));
    otherwise
        error('Enter space or tab')
end

st = [1 whsep+1];
en = [st(2:end) - 1 length(invar)+1];
n = length(st);

for i = 1:n
    outcell{i} = invar(st(i):en(i)-1);
end

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
% [im,stamp] = sbfmfreadframe(f,fid,frame2file,bgcenter)
function [im,stamp] = sbfmfreadframe(f,fid,frame2file,bgcenter)

if fseek(fid,frame2file(f),'bof') == 0,
  [im,stamp] = sbfmf_read_frame(fid,bgcenter);
else
  im = -1;
  stamp = -1;
end


text:
function ds_dsm = cosmo_dissimilarity_matrix_measure(ds, varargin)
% Compute a dissimilarity matrix measure
%
% ds_dsm = cosmo_dissimilarity_matrix_measure(ds[, varargin])
%
% Inputs:
%  dataset            dataset struct with fields .samples (PxQ) and
%                     .sa.targets (Px1) for P samples and Q features.
%                     each target should occur exactly once
%  args               optional struct:
%      .metric        a string with the name of the distance
%                     metric to be used by pdist (default: 'correlation')
%      .center_data   If true, then data is centered before the pair-wise
%                     distances are computed. The default is false; when
%                     used with the 'correlation' metric, it is recommended
%                     to use center_data, 'true'.
%
%   Returns
%
% Output:
%    ds_sa            Struct with fields:
%      .samples       Nx1 flattened lower triangle of a dissimilarity
%                     matrix as returned by [cosmo_]pdist, where
%                     N=P*(P-1)/2 is the number of pairwise distances
%                     between all samples in the dataset.
%      .a.sdim.labels Set to
%      .sa            Struct with field:
%        .targets1    } Nx1 vectors indicating the pairs of indices in the
%        .targets2    } lower part of the square form of the dissimilarity
%                       matrix. if .dsm_pairs(k,:)==[i,j] then .samples(k)
%                     the dissimlarity between the i-th and j-th sample
%                     target.
%
%
% Example:
%     % ds is a dataset struct with ds.sa.targets=(11:16)';
%     ds=struct();
%     ds.samples=[1 2 3; 1 2 3; 1 0 1; 1 1 2; 1 1 2];
%     ds.sa.targets=(11:15)';
%     %
%     % compute dissimilarity with centered data
%     dsm_ds=cosmo_dissimilarity_matrix_measure(ds,'center_data',true);
%     cosmo_disp(dsm_ds);
%     %|| .sa
%     %||   .targets1
%     %||     [ 2
%     %||       3
%     %||       4
%     %||       :
%     %||       4
%     %||       5
%     %||       5 ]@10x1
%     %||   .targets2
%     %||     [ 1
%     %||       1
%     %||       1
%     %||       :
%     %||       3
%     %||       3
%     %||       4 ]@10x1
%     %|| .a
%     %||   .sdim
%     %||     .labels
%     %||       { 'targets1'  'targets2' }
%     %||     .values
%     %||       { [ 11    [ 11
%     %||           12      12
%     %||           13      13
%     %||           14      14
%     %||           15 ]    15 ] }
%     %|| .samples
%     %||   [         0
%     %||             2
%     %||             2
%     %||         :
%     %||      1.11e-16
%     %||      1.11e-16
%     %||     -2.22e-16 ]@10x1
%     %
%     % map results to matrix. values of 0 mean perfect correlation
%     [samples, labels, values]=cosmo_unflatten(dsm_ds,1,...
%                                           'set_missing_to',NaN);
%     cosmo_disp(samples)
%     %|| [ NaN       NaN       NaN       NaN       NaN
%     %||     0       NaN       NaN       NaN       NaN
%     %||     2         2       NaN       NaN       NaN
%     %||     2         2  1.11e-16       NaN       NaN
%     %||     2         2  1.11e-16 -2.22e-16       NaN ]
%     %
%     cosmo_disp(labels)
%     %|| { 'targets1'  'targets2' }
%     %
%     cosmo_disp(values)
%     %|| { [ 11    [ 11
%     %||     12      12
%     %||     13      13
%     %||     14      14
%     %||     15 ]    15 ] }
%
%     % Searchlight using this measure
%     ds=cosmo_synthetic_dataset('ntargets',6,'nchunks',1);
%     % (in this toy example there are only 6 voxels, and the radius
%     %  of the searchlight is 1 voxel. Real-life examples use larger
%     %  datasets and a larger radius)
%     nbrhood=cosmo_spherical_neighborhood(ds,'radius',1,'progress',false);
%     opt=struct();
%     opt.progress=false;          % do not show progress
%     opt.metric='euclidean'; % (instead of default 'correlation')
%     measure=@cosmo_dissimilarity_matrix_measure;
%     sl_ds=cosmo_searchlight(ds, nbrhood, measure, opt);
%     cosmo_disp(sl_ds);
%     %|| .a
%     %||   .fdim
%     %||     .labels
%     %||       { 'i'  'j'  'k' }
%     %||     .values
%     %||       { [ 1         2         3 ]  [ 1         2 ]  [ 1 ] }
%     %||   .vol
%     %||     .mat
%     %||       [ 2         0         0        -3
%     %||         0         2         0        -3
%     %||         0         0         2        -3
%     %||         0         0         0         1 ]
%     %||     .dim
%     %||       [ 3         2         1 ]
%     %||     .xform
%     %||       'scanner_anat'
%     %||   .sdim
%     %||     .labels
%     %||       { 'targets1'  'targets2' }
%     %||     .values
%     %||       { [ 1    [ 1
%     %||           2      2
%     %||           3      3
%     %||           4      4
%     %||           5      5
%     %||           6 ]    6 ] }
%     %|| .fa
%     %||   .nvoxels
%     %||     [ 3         4         3         3         4         3 ]
%     %||   .radius
%     %||     [ 1         1         1         1         1         1 ]
%     %||   .center_ids
%     %||     [ 1         2         3         4         5         6 ]
%     %||   .i
%     %||     [ 1         2         3         1         2         3 ]
%     %||   .j
%     %||     [ 1         1         1         2         2         2 ]
%     %||   .k
%     %||     [ 1         1         1         1         1         1 ]
%     %|| .samples
%     %||   [  3.1      3.68      3.56      1.47      2.96      2.27
%     %||     6.06      6.39      3.29      6.54      4.43       4.1
%     %||     5.85       4.2      2.11      6.47      6.18      3.47
%     %||       :         :        :          :         :         :
%     %||     4.62      3.18     0.829      5.53      5.53      3.14
%     %||      3.7      3.08      1.75      4.71      4.95      3.39
%     %||     1.23      0.83      1.48      1.03      1.75      1.31 ]@15x6
%     %|| .sa
%     %||   .targets1
%     %||     [ 2
%     %||       3
%     %||       4
%     %||       :
%     %||       5
%     %||       6
%     %||       6 ]@15x1
%     %||   .targets2
%     %||     [ 1
%     %||       1
%     %||       1
%     %||       :
%     %||       4
%     %||       4
%     %||       5 ]@15x1
%     %||
%
%     % limitation: cannot have repeated targets
%     ds=cosmo_synthetic_dataset('nchunks',2,'ntargets',3);
%     cosmo_dissimilarity_matrix_measure(ds);
%     %|| error('...')
%
%     % averaging the samples for each unique target resolves the issue of
%     % repeated targets
%     ds=cosmo_synthetic_dataset('nchunks',2,'ntargets',3);
%     ds_avg=cosmo_fx(ds,@(x)mean(x,1),'targets');
%     ds_dsm=cosmo_dissimilarity_matrix_measure(ds_avg);
%     cosmo_disp(ds_dsm);
%     ||.sa
%||  .targets1
%||    [ 2
%||      3
%||      3 ]
%||  .targets2
%||    [ 1
%||      1
%||      2 ]
%||.a
%||  .sdim
%||    .labels
%||      { 'targets1'  'targets2' }
%||    .values
%||      { [ 1    [ 1
%||          2      2
%||          3 ]    3 ] }
%||.samples
%||  [  1.68
%||     1.71
%||    0.711 ]
%
% Notes:
%   - it is recommended to set the 'center_data' to true when using
%     the default 'correlation' metric, as this removes a main effect
%     common to all samples; but note that this option is disabled by
%     default due to historical reasons.
%  -  [cosmo_]pdist defaults to 'euclidean' distance, but correlation
%     distance is preferable for neural dissimilarity matrices, hence it
%     is used as the default here
%  -  Results from this function, when used with the default 'correlation'
%     metric, should *not* be Fisher transformed (using atanh) because
%     the output ranges from 0 to 2 (=one minus Pearson correlation)
%     and the Fisher transform of a value >1 is complex (non-real). This is
%     generally a Bad Thing.
%
% See also: cosmo_pdist, pdist
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #

    % check input
    check_input(ds);
    args=get_args(varargin);

    % make new dataset
    ds_dsm=struct();

    % if center_data, then subtract the mean first
    samples=ds.samples;
    if args.center_data
        samples=bsxfun(@minus,samples,mean(samples,1));
    end

    % compute pair-wise distances between all samples using cosmo_pdist,
    % then store them as samples in ds_dsm
    % >@@>
    dsm = cosmo_pdist(samples, args.metric)';

    % store dsm
    ds_dsm=get_sample_attributes(ds.sa.targets);
    ds_dsm.samples=dsm;
    % <@@<

function check_input(ds)
    if ~(isstruct(ds) && ...
                isfield(ds,'samples') && ...
                isfield(ds,'sa') && ...
                isfield(ds.sa,'targets'))
        error(['require dataset structure with fields '...
                    '.samples and .sa.targets']);
    end

function args=get_args(varargin)
    persistent cached_varargin;
    persistent cached_args;

    if ~isequal(varargin, cached_varargin)
        cached_args=cosmo_structjoin('metric','correlation',...
                            'center_data',false,...
                            varargin);
        cached_varargin=varargin;
    end

    args=cached_args;

function ds_skeleton=get_sample_attributes(targets)
    persistent cached_targets;
    persistent cached_ds_skeleton;

    if ~isequal(targets, cached_targets)
        ntargets=numel(targets);

        % unique targets
        classes=unique(targets);
        nclasses=numel(classes);

        % each should occur exactly once
        if nclasses~=ntargets
            error(['.sa.targets should be permutation of unique targets; '...
                    'to average samples with the same targets, consider '...
                    'ds_mean=cosmo_fx(ds,@(x)mean(x,1),''targets'')'],...
                        nclasses);
        end

        % store single sample attribute: the pairs of sample attribute indices
        % used to compute the dsm.
        [i,j]=find(triu(repmat(1:nclasses,nclasses,1),1)');
        cached_ds_skeleton.sa=struct();
        cached_ds_skeleton.sa.targets1=i;
        cached_ds_skeleton.sa.targets2=j;

        % set sample dimensions
        add_labels={'targets1','targets2'};
        add_values={targets, targets};

        cached_ds_skeleton.a.sdim=struct();
        cached_ds_skeleton.a.sdim.labels=add_labels;
        cached_ds_skeleton.a.sdim.values=add_values;

        cached_targets=targets;
    end

    ds_skeleton=cached_ds_skeleton;

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
function out = combineCells(x,dim,d)
%combines a collection of cells either horizontally or vertically,
%depending on the initial orientation

    x = x(returnCellLengths(x) > 0);

    if isempty(x)

        out = [];

    else

        L = length(x);
        lengths = returnCellLengths(x);
        [~,maxIdx] = max(lengths);
        s = size(x{maxIdx});

        if nargin < 2 || isempty(dim)
            dim = argmax(s);
        end

        if dim == 1

            if nargin  < 3 || isempty(d)
                d = s(2);
            end
            lengths = returnCellLengths(x) ./ d;
            s2 = size(lengths);
            if s2(2) > s2(1)
                lengths = lengths';
            end
            cVals = [0; cumsum(lengths)];

            out = zeros(cVals(end),d);

            for i=1:L
                if ~isempty(x{i})
                    out(cVals(i)+1:cVals(i+1),:) = x{i};
                end
            end

        else

            if nargin  < 3 || isempty(d)
                d = s(1);
            end
            lengths = returnCellLengths(x) ./ d;
            s2 = size(lengths);
            if s2(1) > s2(2)
                lengths = lengths';
            end
            cVals = [0 cumsum(lengths)];

            out = zeros(d,cVals(end));

            for i=1:L;
                if ~isempty(x{i})
                    out(:,cVals(i)+1:cVals(i+1)) = x{i};
                end
            end

        end

    end

author: gordonberman

repo: MotionMapper


text:
function imgOut = ConvertLinearSpace(img, mtx)
%
%       imgOut = ConvertLinearSpace(img, mtx)
%
%
%        Input:
%           -img: image to convert into a new color space
%           -mtx: a 3x3 matrix that defines a linear color transformation
%
%        Output:
%           -imgOut: converted image
%
%     Copyright (C) 2011  Francesco Banterle
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
%

%Is it a three color channels image?
check3Color(img);

%Is it a 3x3 matrix?
[r, c]=size(mtx);
if(r ~= 3 || c ~= 3)
    error('The matrix for color transformation is not 3x3.');
end

imgOut = zeros(size(img));
for i=1:3
    imgOut(:,:,i) = img(:,:,1) * mtx(i,1) + img(:,:,2) * mtx(i,2) + img(:,:,3) * mtx(i,3);
end

end

author: banterle


text:
function mrAnatAverageAcpcAnalyze(fileNameList, outFileBaseName, alignLandmarks, newMmPerVox, weights, bb, clipVals, showFigs);
%
% mrAnatAverageAcpcAnalyze(fileNameList, outFileBaseName, [alignLandmarks=[]], ...
%                          [newMmPerVox=[1 1 1]], [weights=ones(size(fileNameList))], ...
%                          [bb], [clipVals], [showFigs=1])
%
% Reslices the first analyze file to ac-pc space at
% 1x1x1mm resolution and then aligns all the rest of the analyze
% files to that one and averages them all together.
%
% You can specify the ac-pc landmarks as a 3x3 matrix of the form:
%  [ acX, acY, acZ; pcX, pcY, pcZ; midSagX, midSagY, midSagZ ]
% Where:
%  X is the (left -> right) location of a point,
%  Y is the (posterior -> anterior) location of a point,
%  Z is the (inferior -> superior) location of a point,
% ac is the anterior commissure, pc is the posterior commissure, and midSag
% is another poit in the mid-sagittal plane that is somewhat distant from
% the ac-pc line. These 3 points define the rotation/translation into ac-pc
% space. If the ac is properly set in the image header, then you can just
% pass in the pc and midSag coords, specified as the offset from the ac.
%
% If these landmarks aren't provided, then the xform from the image header
% is used. Note that if the image origin isn't close to the AC, then the
% default bounding box ([-90,90; -126,90; -72,108]') won't be appropriate.
%
% weights specifies the weighting factor to be applied to each of the input
% images (fileNameList).
%
% REQUIRES:
%  * Stanford anatomy tools (eg. /usr/local/matlab/toolbox/mri/Anatomy)
%  * spm2 tools (eg. /usr/local/matlab/toolbox/mri/spm2)
%
% HISTORY:
% 2004.11.10 RFD (bob@white.stanford.edu) wrote it, based on averageAnalyze
if ~exist('showFigs','var') | isempty(showFigs), showFigs = 1; end

if (~exist('fileNameList','var') | isempty(fileNameList) | ...
        ~exist('outFileBaseName','var') | isempty(outFileBaseName))
    help(mfilename);
    return;
end
if (~exist('newMmPerVox','var') | isempty(newMmPerVox))
    newMmPerVox = [1 1 1];
end
if ~exist('saveIntermediate','var') | isempty(saveIntermediate)
  saveIntermediate = 0;
end

if (~exist('alignLandmarks','var')) alignLandmarks = []; end

if ~exist('weights','var') | isempty(weights)
  weights = ones(size(fileNameList));
end

if(~exist('bb','var') | isempty(bb))
    % Bounding box, in physical space (ie. mm from the origin, which should be
    % at or near the AC).
    bb = [-90,90; -126,90; -72,108]';
end

if (~exist('clipVals','var')) clipVals = []; end

% from spm_bsplins:
% d(1:3) - degree of B-spline (from 0 to 7) along different dimensions
% d(4:6) - 1/0 to indicate wrapping along the dimensions
% not sure what wrapping is, but '7' is the highest quality (but slowest).
bSplineParams = [7 7 7 0 0 0];

% We explicitly initialize the spm_defaults global here, and ensure that
% the analyze_flip option is turned off. (Our analyze files are never
% right-left reversed!)
spm_defaults;
defaults.analyze.flip = 0;

numImages = length(fileNameList);

% Load the first image (the reference)
[refImg,mmPerVox,refImgHdr] = loadAnalyze(fileNameList{1});

if(isempty(clipVals))
  clipVals = repmat([0.4 0.98],numImages,1);
end

refImg = mrAnatHistogramClip(refImg, clipVals(1,1), clipVals(1,2));
%[refImg, lc, uc] = mrAnatHistogramClipOptimal(refImg, 99);
%fprintf('\nClipped reference image at [%0.1f, %0.1f].\n', lc, uc);

if(isempty(alignLandmarks))
    % *** TO DO: allow user to select the landmarks!
    tal2ref = inv(refImgHdr.mat);

else
    if(size(alignLandmarks,1)==2)
        origin = inv(refImgHdr.mat)*[0 0 0 1]'-0.5;
        %origin(3) = size(refImg,3)-origin(3);
        %alignLandmarks(:,3) = -alignLandmarks(:,3);
        origin = origin(1:3)';
        imY = alignLandmarks(1,:); imY = imY./norm(imY);
        imZ = alignLandmarks(2,:); imZ = imZ./norm(imZ);
    else
        %% flip 3rd axis
        %alignLandmarks(:,3) = size(refImg,3)-alignLandmarks(:,3);
        % The first landmark should be the anterior commissure (AC)- our origin
        origin = alignLandmarks(1,:);
        % Define the current image axes by re-centering on the origin (the AC)
        imY = alignLandmarks(2,:)-origin; imY = imY./norm(imY);
        imZ = alignLandmarks(3,:)-origin; imZ = imZ./norm(imZ);
    end
    % x-axis (left-right) is the normal to [ac, pc, mid-sag] plane
    imX = cross(imZ,imY);
    % Make sure the vectors point right, superior, anterior
    if(imX(1)<0) imX = -imX; end
    if(imY(2)<0) imY = -imY; end
    if(imZ(3)<0) imZ = -imZ; end
    % Project the current image axes to the cannonical AC-PC axes. These
    % are defined as X=[1,0,0], Y=[0,1,0], Z=[0,0,1], with the origin
    % (0,0,0) at the AC. Note that the following are the projections
    x = [0 1 imY(3)]; x = x./norm(x);
    y = [1  0 imX(3)]; y = y./norm(y);
    %z = [0  imX(2) 1]; z = z./norm(z);
    z = [0  -imY(1) 1]; z = z./norm(z);
    % Define the 3 rotations using the projections. We have to set the sign
    % of the rotation, depending on which side of the plane we came from.
    rot(1) = sign(x(3))*acos(dot(x,[0 1 0])); % rot about x-axis (pitch)
    rot(2) = sign(y(3))*acos(dot(y,[1 0 0])); % rot about y-axis (roll)
    rot(3) = sign(z(2))*acos(dot(z,[0 0 1])); % rot about z-axis (yaw)

    scale = mmPerVox;

    % Affine build assume that we need to translate before rotating. But,
    % our rotations have been computed about the origin, so we'll pass a
    % zero translation and set it ourselves (below).
    ref2tal = affineBuild([0 0 0], rot, scale, [0 0 0]);
    tal2ref = inv(ref2tal);

    % Insert the translation.
    tal2ref(1:3,4) = [origin+newMmPerVox/2]';
end

% Resample it to 1x1x1
disp('Resampling reference image to ac-pc space, isotropic voxels...');
[newRefImg,xform] = mrAnatResliceSpm(refImg, tal2ref, bb, newMmPerVox, bSplineParams, showFigs);
newOrigin = inv(xform)*[0 0 0 1]'; newOrigin = newOrigin(1:3)'-newMmPerVox/2;
% Reclip in case the interpolation introduced out-of-range values
newRefImg(newRefImg<0) = 0; newRefImg(newRefImg>1) = 1;
if(showFigs)
    o = round(newOrigin);
    figure; set(gcf,'Name',[fileNameList{1} ' (ref)']);
    subplot(1,3,1); imagesc(flipud(squeeze(newRefImg(:,:,o(3)))')); axis image; colormap gray;
    subplot(1,3,2); imagesc(flipud(squeeze(newRefImg(:,o(2),:))')); axis image; colormap gray;
    subplot(1,3,3); imagesc(flipud(squeeze(newRefImg(o(1),:,:))')); axis image; colormap gray;
    %imagesc(makeMontage(refImg,[20:4:size(newRefImg,3)-18]));axis image;colormap gray;
    %title([fileNameList{1} ' (reference image) aligned.']);
    pause(0.1);
end
Vref.uint8 = uint8(round(newRefImg.*255));
Vref.mat = xform;

outImg = newRefImg.*weights(1);
numSamples = zeros(size(outImg));
nans = isnan(outImg);
numSamples(~nans) = weights(1);
outImg(nans) = 0;
for(ii=2:numImages)
    if(saveIntermediate)
        % We do this first since the final result will be saved after this loop.
        fname = sprintf('%s_%d',outFileBaseName,ii-1);
        disp(['writing intermediate result ',fname,'...']);
        img = outImg;
        nz = numSamples>0;
        img(nz) = img(nz)./numSamples(nz);
        img = img-min(img(:));
        img = int16(img.*(32767/max(img(:))));
        saveAnalyze(img, fname, newMmPerVox, ['AVERAGE:' refImgHdr.descrip], newOrigin);
    end
    fprintf('Aligning %s to reference image...\n',fileNameList{ii});
    [img,mmPerVox,hdr] = loadAnalyze(fileNameList{ii});
    img = mrAnatHistogramClip(img, clipVals(ii,1), clipVals(ii,2));
    %[img, lc, uc] = mrAnatHistogramClipOptimal(img, 99);
    %fprintf('Clipped image %d at [%0.1f, %0.1f].\n', ii, lc, uc);
    Vin.uint8 = uint8(round(img.*255));
    Vin.mat = hdr.mat;
    %transRot = spm_coreg(Vin, Vref);
    %xform = spm_matrix(transRot(:)')*Vin.mat\Vref.mat*inv(Vref.mat);
    transRot = spm_coreg(Vref, Vin);
    xform = inv(Vin.mat)*spm_matrix(transRot(:)');
    fprintf('Resampling %s to reference image...\n',fileNameList{ii});
    [img,xform] = mrAnatResliceSpm(img, xform, bb, newMmPerVox, bSplineParams, showFigs);
    % Reclip in case the interpolation introduced out-of-range values
    img(img<0) = 0; img(img>1) = 1;
    if(showFigs)
        o = round(newOrigin);
    	figure; set(gcf,'Name',[fileNameList{ii}]);
        subplot(1,3,1); imagesc(flipud(squeeze(img(:,:,o(3)))')); axis image; colormap gray;
        subplot(1,3,2); imagesc(flipud(squeeze(img(:,o(2),:))')); axis image; colormap gray;
        subplot(1,3,3); imagesc(flipud(squeeze(img(o(1),:,:))')); axis image; colormap gray;
        %figure; imagesc(makeMontage(img,[20:4:size(img,3)-18])); axis image; colormap gray;
        %title([fileNameList{ii} ' aligned.']);
        pause(0.1);
    end
    nans = isnan(img);
    numSamples(~nans) = numSamples(~nans)+weights(ii);
    img(nans) = 0;
    outImg = outImg+img.*weights(ii);
end
% Rescale based on the number of samples at each voxel
nz = numSamples>0;
outImg(nz) = outImg(nz)./numSamples(nz);

if(showFigs)
    o = round(newOrigin);
    figure; set(gcf,'Name',['Average']);
    subplot(1,3,1); imagesc(flipud(squeeze(outImg(:,:,o(3)))')); axis image; colormap gray;
    subplot(1,3,2); imagesc(flipud(squeeze(outImg(:,o(2),:))')); axis image; colormap gray;
    subplot(1,3,3); imagesc(flipud(squeeze(outImg(o(1),:,:))')); axis image; colormap gray;
    figure; imagesc(makeMontage(outImg,[20:4:size(outImg,3)-18])); axis image; colormap gray;
    title(['Average aligned.']);
    pause(0.1);
end

% rescale to 15 bits (0-32767)
outImg = outImg-min(outImg(:));
outImg = int16(outImg.*(32767/max(outImg(:))));

disp(['writing ',outFileBaseName,'...']);
hdr = saveAnalyze(outImg, outFileBaseName, newMmPerVox, ['AVERAGE:' refImgHdr.descrip], newOrigin);

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function convBconv_last()
    global config mem;
    curr_layer_idx = config.misc.current_layer;
    mem.deltas{curr_layer_idx+1} = reshape(accumarray(mem.convBconv{curr_layer_idx+2}, mem.deltas{curr_layer_idx+1}(:)), ...
                                            mem.orig_activation_size{curr_layer_idx+1}(2), mem.orig_activation_size{curr_layer_idx+1}(1))';
                                           %size(mem.layer_inputs{curr_layer_idx+1}, 1), size(mem.layer_inputs{curr_layer_idx+1}, 2));
end

author: jimmy-ren


text:
function [ C ] = gsp_gwft(G,g,f, param )
%GSP_GWFT Generalized windowed Fourier transform
%   Usage:  C = gsp_gwft(G,g,f, param );
%           C = gsp_gwft(G,g,f);
%
%   Input parameters:
%         G     : Graph
%         g     : Window (graph signal or kernel)
%         f     : Graph signal (column vector)
%         param : Structure of optional parameter
%   Output parameters:
%         C     : Coefficient.
%
%   This function compute the graph windowed Fourier transform of a signal
%   *f* with the window *g*. The function returns a matrix of size N^2*N.
%
%   *param* a Matlab structure containing the following fields:
%
%     * *param.verbose* : 0 no log, 1 print main steps, 2 print all steps.
%       By default, it is 1.
%     * *param.lowmemory* : use less memory. By default, it is 1.
%
%   Reference: shuman2013windowed

% Author : Nathanael Perraudin
% Testing: test_gwft

% Optional input arguments
if nargin<4, param=struct; end
if ~isfield(param, 'verbose'), param.verbose=1 ; end
if ~isfield(param, 'lowmemory'), param.lowmemory=1 ; end

if ~isfield(G,'U')
   error(['GSP_GWFT: You need first to compute the Fourier basis. ',...
       'You can do it with the function gsp_compute_fourier_basis']);
end

if sum(G.U(:,1)<eps(2))
    error(['GSP_GWFT: The current implementation of this function is ',...
       'not working for disconnected graphs'])
end

Nf = size(f,2);

if iscell(g)
    g = gsp_igft(G,g{1}(G.e));
end

if isa(g, 'function_handle')
    g = gsp_igft(G,g(G.e));
end

if ~param.lowmemory
    % Compute the Frame into a big matrix
    Frame=gsp_gwft_frame_matrix(G,g,param);

    C=Frame'*f;

    C=reshape(C,G.N,G.N,Nf);
else
    % Compute the translate of g
    ghat=G.U'*g;
    Ftrans=sqrt(G.N)*G.U*(repmat(ghat,1,G.N).*G.U');
    C=zeros(G.N);
    for jj = 1:Nf
        for ii=1:G.N
             C(:,ii,jj)=(repmat(1./G.U(:,1),1,G.N) .* ...
                      G.U.*repmat(Ftrans(:,ii),1,G.N))'*f(:,jj);
        end
    end

end

repo: gspbox


text:
function [lbm] = g2lbm(g)
% Convert mass from grams to pounds-mass.
% Chad Greene 2012
lbm = g*0.002204622621849 ;

author: Sable

repo: mcbench-benchmarks


text:
function g = power(f, b)
%.^   CHEBTECH power.
%   F.^G returns a CHEBTECH F to the scalar power G, a scalar F to the CHEBTECH
%   power G, or a CHEBTECH F to the CHEBTECH power G. F and or G may be complex.
%
%   This function assumes that the curve traced out by F in the complex plane
%   both (1) does not come too close to zero and (2) does not cross over the
%   branch cut in POWER along the negative real axis.  That is, F should not
%   vanish at any point of [-1, 1], and the imaginary part of F should not
%   vanish at any point of (-1, 1) where the real part of F is negative.  If any
%   of these assumptions are violated, garbage may be returned with no warning.
%
%   H = POWER(F, G) is called for the syntax 'F .^ G'.
%
% See also SQRT.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

% If F is complex and the imaginary part of F vanishes exactly at the domain
% boundary, tiny rounding errors in evaluating the imaginary part at the
% boundary points can cause it to appear to change sign there.  If this happens,
% and if the real part of F is negative there, the endpoint function values will
% cross over the branch cut in POWER, creating a jump discontinuity that should
% not be there.  Enabling extrapolation resolves this issue by avoiding these
% function evaluations:

pref = f.techPref();
if ( ~isreal(f) )
    pref.extrapolate = 1;
end

% Simply call the compose function:
if ( isa(f, 'chebtech') && isa(b, 'chebtech') )
    % Both F and G are CHEBTECHs:
	g = compose(f, @power, b, [], [], pref);
elseif ( isa(f, 'chebtech') )
    % F is CHEBTECH and G is constant:
	g = compose(f, @(f) power(f, b), [], [], pref);
else
    % F is constant and G is CHEBTECH:
	g = compose(b, @(b) power(f, b), [], [], pref);
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [m] = in2m(in)
% Convert length from inches to meters.
% Chad Greene 2012
m = in*.0254;


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Spherical Harmonic Modeling and Analysis Toolkit (SPHARM-MAT) is a 3D
% shape modeling and analysis toolkit.
% It is a software package developed at Shenlab in Center for Neuroimaging,
% Indiana University (SpharmMat@gmail.com, http://www.iupui.edu/~shenlab/)
% It is available to the scientific community as copyright freeware
% under the terms of the GNU General Public Licence.
%
% Copyright 2009, 2010, ShenLab, Center for Neuroimaging, Indiana University
%
% This file is part of SPHARM-MAT.
%
% SPHARM-MAT 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.
%
% SPHARM-MAT 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 SPHARM-MAT. If not, see <http://www.gnu.org/licenses/>.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%
%
%

function write_metaNcoef(filename, confs)

load(filename);

if ~exist('faces', 'var') | ~exist('vertices', 'var') | ~exist('fvec', 'var')
    disp('One or more of faces, vertices, and SPHARM coefficients are missing');
    return;
end

[fvec] = read_coef(coef);

[pa, na, ex]=fileparts(filename);

type = size(faces,2);
filename1 = sprintf('%s/%s_surfSPHARM.meta', confs.OutDirectory, na(1:end-4));
filename2 = sprintf('%s/%s_surfSPHARM.coef', confs.OutDirectory, na(1:end-4));

if exist(filename1,'file') & exist(filename2,'file')
    prompt = {'Enter new filenames:'};
    dlg_title = 'New File Names';
    num_lines = 2;
    def = {filename1; filename2};
    answer = inputdlg(prompt,dlg_title,num_lines,def);
    filename1 = answer{1};
    filename2 = answer{2};
end

% Write a meta file
if type == 3
    disp('This function does not support any objects in the general triangular mesh now');
elseif type == 4
    write_meta_quad(filename1, vertices, faces);
%    write_meta_quad(filename2, sph_verts, faces);
end

% Write a coef file
write_coef(filename2, fvec);

return;

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
%% RATE OF CONVERGENCE OF LINEAR ELEMENT FOR POISSON EQUATION
%
% This example is to show the rate of convergence of linear finite element
% approximation of the Poisson equation on the unit cube with the
% following boundary conditions:
%
% - Non-empty Dirichlet boundary condition.
% - Pure Neumann boundary condition.
% - Robin boundary condition.
%
% The basis, data structure and numerical test is summarized in <a
% href="matlab:ifem Poisson3femrate">Poisson3femrate</a>.
%
% See also Poisson3P2femrate, Poissonfemrate
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

clear variables

%% Setting
[node,elem] = cubemesh([0,1,0,1,0,1],0.5);
mesh = struct('node',node,'elem',elem);
option.L0 = 1;
option.maxIt = 4;
option.elemType = 'P1';
option.printlevel = 1;
option.plotflag = 1;

%% Non-empty Dirichlet boundary condition.
fprintf('Mixed boundary conditions. \n');
pde = sincosdata3;
mesh.bdFlag = setboundary3(node,elem,'Dirichlet','~(x==0)','Neumann','x==0');
% bdFlag = setboundary3(node,elem,'Dirichlet');
femPoisson3(mesh,pde,option);

%% Pure Neumann boundary condition.
fprintf('Pure Neumann boundary condition. \n');
option.plotflag = 0;
pde = sincosdata3Neumann;
mesh.bdFlag = setboundary3(node,elem,'Neumann');
femPoisson3(mesh,pde,option);

%% Robin boundary condition.
fprintf('Robin boundary condition. \n');
pde = sincosRobindata3;
mesh.bdFlag = setboundary3(node,elem,'Robin');
femPoisson3(mesh,pde,option);

%% Conclusion
%
% The optimal rate of convergence of the H1-norm (1st order) and L2-norm
% (2nd order) is observed. The 2nd order convergent rate between two
% discrete functions ||DuI-Duh|| is known as superconvergence.
%
% MGCG converges uniformly in all cases.
%
% Note that for Neumann and Robin boundary condition, the order of the error in the
% maximum norm is not close to 2. Neumann problem is around 1.5 and Robin
% is near 1.3.

repo: ifem


text:
function ActiveImportanceSampling = activeimportancesampling_vbmc(vp,gp,acqfun,acqinfo,options)
%ACTIVEIMPORTANCESAMPLING_VBMC Setup importance sampling acquisition functions

% This function samples from the base importance sampling (IS) density in
% three steps:
% 1) Use importance sampling-resampling (ISR) to sample from the
%    base IS density based on a proposal distribuion which is a mixture of
%    a smoothed variational posterior and on box-uniform mixture centered
%    around the current training points
% 2) Optionally use MCMC initialized with the ISR samples to sample from
%    the base IS density
% 3) Compute IS statistics and book-keeping

% Does the importance sampling step use the variational posterior?
isamplevp_flag = isfield(acqinfo,'importance_sampling_vp') ...
    && acqinfo.importance_sampling_vp;

% Do we simply sample from the variational posterior?
onlyvp_flag = isfield(acqinfo,'variational_importance_sampling') ...
    && acqinfo.variational_importance_sampling;

D = size(gp.X,2);

Ns_gp = numel(gp.post); % # GP hyperparameter samples

% Input space bounds and typical scales (for MCMC only)
widths = std(gp.X,[],1);
MaxBnd = 0.5;
diam = max(gp.X) - min(gp.X);
LB = min(gp.X) - MaxBnd*diam;
UB = max(gp.X) + MaxBnd*diam;

ActiveImportanceSampling.lnw = [];
ActiveImportanceSampling.Xa = [];
ActiveImportanceSampling.fs2a = [];

if onlyvp_flag
    %% Step 0: Simply sample from the variational posterior

    if ischar(options.ActiveImportanceSamplingMCMCSamples)
        K = vp.K;
        nvars = D;
        Na = ceil(eval(options.ActiveImportanceSamplingMCMCSamples));
    elseif isscalar(options.ActiveImportanceSamplingMCMCSamples)
        Na = ceil(options.ActiveImportanceSamplingMCMCSamples);
    else
        Na = 0;
    end

    if ~isfinite(Na) || ~isscalar(Na) || Na <= 0
        error('OPTIONS.ActiveImportanceSamplingMCMCSamples should be (or evaluate to) a positive integer.');
    end

    Xa = vbmc_rnd(vp,Na,0);
    [~,~,fmu,fs2] = gplite_pred(gp,Xa,[],[],1,0);

    if isfield(acqinfo,'mcmc_importance_sampling') && acqinfo.mcmc_importance_sampling

        % Compute fractional effective sample size (ESS)
        fESS = fess_vbmc(vp,fmu,Xa);

        % If fESS is less than thresh ==> major mismatch, do MCMC
        if fESS < options.ActiveImportanceSamplingfESSThresh
            Xa_old = Xa;

            if isamplevp_flag
                logpfun = @(x) log_isbasefun(x,acqfun,gp,vp);
            else
                logpfun = @(x) log_isbasefun(x,acqfun,gp,[]);
            end

            % Get MCMC options
            Nmcmc_samples = Na*options.ActiveImportanceSamplingMCMCThin;
            thin = 1;
            burnin = 0;
            sampleopts = get_mcmcopts([],thin,burnin);
            logPfuns = logpfun;
            W = Na;    % # walkers

            % Perform a single MCMC step for all samples
            Xa = eissample_lite(logPfuns,Xa,Nmcmc_samples,W,widths,LB,UB,sampleopts);
            Xa = Xa(end-Na+1:end,:);
            [~,~,fmu,fs2] = gplite_pred(gp,Xa,[],[],1,0);

            if 0
                hold off;
                scatter(Xa_old(:,1),Xa_old(:,2),'b'); hold on;
                scatter(Xa(:,1),Xa(:,2),'k'); hold on;
                drawnow;
            end
        end
    end

    if isamplevp_flag
        vlnpdf = max(vbmc_pdf(vp,Xa,0,1),log(realmin));
        lny = acqfun('islogf1',vlnpdf,[],[],fmu,fs2);
    else
        lny = acqfun('islogf1',[],[],[],fmu,fs2);
    end

    ActiveImportanceSampling.fs2a = fs2;
    ActiveImportanceSampling.lnw = lny';
    ActiveImportanceSampling.Xa = Xa;

else
    %% Step 1: Importance sampling-resampling

    Nvp_samples = options.ActiveImportanceSamplingVPSamples;
    Nbox_samples = options.ActiveImportanceSamplingBoxSamples;
    w_vp = Nvp_samples/(Nvp_samples + Nbox_samples);

    rect_delta = 2*std(gp.X);

    % Smoothed posterior for importance sampling-resampling
    if Nvp_samples > 0
        scale_vec = [0.05,0.2,1];
        %scale = sqrt(mean(var(gp.X,[],1)));

        vp_is = vp;
        for ii = 1:numel(scale_vec)
            vp_is.K = vp_is.K + vp.K;
            vp_is.w = [vp_is.w, vp.w];
            vp_is.mu = [vp_is.mu, vp.mu];
            vp_is.sigma = [vp_is.sigma, sqrt(vp.sigma.^2 + scale_vec(ii)^2)];
        end
        vp_is.w = vp_is.w/sum(vp_is.w);

        % Sample from smoothed posterior
        Xa_vp = vbmc_rnd(vp_is,Nvp_samples,0,0);
        [lnw,fs2a_vp] = activesample_proposalpdf(Xa_vp,gp,vp_is,w_vp,rect_delta,acqfun,vp,isamplevp_flag);
        ActiveImportanceSampling.lnw = [ActiveImportanceSampling.lnw, lnw'];
        ActiveImportanceSampling.Xa = [ActiveImportanceSampling.Xa; Xa_vp];
        ActiveImportanceSampling.fs2a = [ActiveImportanceSampling.fs2a; fs2a_vp];
    else
        vp_is = [];
    end

    % Box-uniform sampling around training inputs
    if Nbox_samples > 0
        jj = randi(size(gp.X,1),[1,Nbox_samples]);
        Xa_box = gp.X(jj,:) + bsxfun(@times,2*rand(numel(jj),D)-1,rect_delta);
        [lnw,fs2a_box] = activesample_proposalpdf(Xa_box,gp,vp_is,w_vp,rect_delta,acqfun,vp,isamplevp_flag);
        ActiveImportanceSampling.lnw = [ActiveImportanceSampling.lnw, lnw'];
        ActiveImportanceSampling.Xa = [ActiveImportanceSampling.Xa; Xa_box];
        ActiveImportanceSampling.fs2a = [ActiveImportanceSampling.fs2a; fs2a_box];
    end

    ActiveImportanceSampling.lnw(~isfinite(ActiveImportanceSampling.lnw)) = -Inf;
    % optimState.w = exp(optimState.lnw - max(optimState.lnw))';
    % optimState.w = optimState.w / sum(optimState.w);
    % 1./sum(optimState.w.^2)

    %% Step 2 (optional): MCMC sample

    Nmcmc_samples = options.ActiveImportanceSamplingMCMCSamples;

    if Nmcmc_samples > 0

        ActiveImportanceSampling_old = ActiveImportanceSampling;

        ActiveImportanceSampling.lnw = zeros(Ns_gp,Nmcmc_samples);
        ActiveImportanceSampling.Xa = zeros(Nmcmc_samples,D,Ns_gp);
        ActiveImportanceSampling.fs2a = zeros(Nmcmc_samples,Ns_gp);

        gp1 = gp;   % Consider only one GP sample at a time

        for s = 1:Ns_gp

            gp1.post = [];  % Assign current GP sample
            gp1.post = gp.post(s);

            if D == 2 && 0
                % We could use a quasi-random grid for D <= 2, but not implemented
                Xa = rand(Na,D).*(UB - LB) + LB;
                [~,~,fmu,optimState.fs2a] = gplite_pred(gp1,Xa,[],[],1,0);
                optimState.lnw = mean(fmu,2)';
                optimState.Xa = Xa;
            else
                if isamplevp_flag
                    logpfun = @(x) log_isbasefun(x,acqfun,gp1,vp);
                else
                    logpfun = @(x) log_isbasefun(x,acqfun,gp1,[]);
                end

                % Get MCMC options
                thin = options.ActiveImportanceSamplingMCMCThin;
                burnin = ceil(thin*Nmcmc_samples/2);
                sampleopts = get_mcmcopts(Nmcmc_samples,thin,burnin);

                logPfuns = logpfun;
                % sampleopts.TransitionOperators = {'transSliceSampleRD'};
                W = 2*(D+1);    % # walkers

                if 0
                    % Take starting points from high posterior density region
                    hpd_frac = 0.5;
                    N = numel(gp1.y);
                    N_hpd = min(N,max(W,round(hpd_frac*N)));
                    [~,ord] = sort(gp1.y,'descend');
                    X_hpd = gp1.X(ord(1:N_hpd),:);
                    x0 = X_hpd(randperm(N_hpd,min(W,N_hpd)),:);
                    x0 = bsxfun(@min,bsxfun(@max,x0,LB),UB);
                else
                    % Use importance sampling-resampling
                    [~,~,fmu,fs2] = gplite_pred(gp1,ActiveImportanceSampling_old.Xa,[],[],1,0);
                    lnw = ActiveImportanceSampling_old.lnw(s,:) + acqfun('islogf2',[],[],[],fmu,fs2)';
                    w = exp(bsxfun(@minus,lnw,max(lnw,[],2)));
                    x0 = zeros(W,D);
                    for ii = 1:W
                        idx = catrnd(w,1);
                        w(idx) = 0;
                        x0(ii,:) = ActiveImportanceSampling_old.Xa(idx,:);
                    end
                end

                [Xa,logp] = eissample_lite(logPfuns,x0,Nmcmc_samples,W,widths,LB,UB,sampleopts);
                [~,~,fmu,fs2] = gplite_pred(gp1,Xa,[],[],1,0);

                % Fixed log weight for importance sampling (log fixed integrand)
                if isamplevp_flag
                    vlnpdf = max(vbmc_pdf(vp,Xa,0,1),log(realmin));
                    lny = acqfun('islogf1',vlnpdf,[],[],fmu,fs2);
                else
                    lny = acqfun('islogf1',[],[],[],fmu,fs2);
                end
                % lny = lny - warpvars_vbmc(Xa,'logp',vp.trinfo);

                ActiveImportanceSampling.fs2a(:,s) = fs2;
                ActiveImportanceSampling.lnw(s,:) = bsxfun(@minus,lny',logp');
                ActiveImportanceSampling.Xa(:,:,s) = Xa;
            end

        end
    end
end

if 0
    hold off;
    scatter(ActiveImportanceSampling.Xa(:,1),ActiveImportanceSampling.Xa(:,2)); hold on;
    % scatter(x0(:,1),x0(:,2),'ro','MarkerFaceColor','r')
    xlim([-2,2]);
    ylim([-2,2]);
    drawnow;
end

%% Step 3: Pre-compute quantities for importance sampling calculations

% Precompute cross-kernel matrix on importance points
Kax_mat = zeros(size(ActiveImportanceSampling.Xa,1),size(gp.X,1),Ns_gp);
Ctmp_mat = zeros(size(gp.X,1),size(ActiveImportanceSampling.Xa,1),Ns_gp);
for s = 1:Ns_gp
    if size(ActiveImportanceSampling.Xa,3) == 1
        Xa = ActiveImportanceSampling.Xa;
    else
        Xa(:,:) = ActiveImportanceSampling.Xa(:,:,s);
    end
    hyp = gp.post(s).hyp;
    L = gp.post(s).L;
    Lchol = gp.post(s).Lchol;
    sn2_eff = 1/gp.post(s).sW(1)^2;
    if gp.covfun(1) == 1    % Hard-coded SE-ard for speed
        ell = exp(hyp(1:D))';
        sf2 = exp(2*hyp(D+1));
        Kax_tmp = sq_dist(Xa*diag(1./ell),gp.X*diag(1./ell));
        Kax_mat(:,:,s) = sf2 * exp(-Kax_tmp/2);
    else
        error('Other covariance functions not supported yet.');
    end

    if Lchol
        Ctmp_mat(:,:,s) = (L\(L'\Kax_mat(:,:,s)'))/sn2_eff;
    else
        Ctmp_mat(:,:,s) = (L*Kax_mat(:,:,s)');
    end
end
ActiveImportanceSampling.Kax_mat = Kax_mat;
ActiveImportanceSampling.Ctmp_mat = Ctmp_mat;

% Precompute integrated mean basis function on importance points
if isfield(gp,'intmeanfun') && gp.intmeanfun > 0
    plus_idx = gp.intmeanfun_var > 0;
    if size(ActiveImportanceSampling.Xa,3) == 1
        Ha = gplite_intmeanfun(ActiveImportanceSampling.Xa,gp.intmeanfun);
        ActiveImportanceSampling.Ha = Ha(plus_idx,:);
    else
        for s = 1:Ns_gp
            Ha = gplite_intmeanfun(ActiveImportanceSampling.Xa(:,:,s),gp.intmeanfun);
            ActiveImportanceSampling.Ha(:,:,s) = Ha(plus_idx,:);
        end
    end
end

end

%--------------------------------------------------------------------------
function [lnw,fs2] = activesample_proposalpdf(Xa,gp,vp_is,w_vp,rect_delta,acqfun,vp,isamplevp_flag)
%ACTIVESAMPLE_PROPOSALPDF Compute importance weights for proposal pdf

[N,D] = size(gp.X);
Na = size(Xa,1);

[~,~,fmu,fs2] = gplite_pred(gp,Xa,[],[],1,0);

Ntot = 1 + N; % Total number of mixture elements

if w_vp < 1; templpdf = zeros(Na,Ntot); end

% Mixture of variational posteriors
if w_vp > 0
    logflag = true;
    templpdf(:,1) = vbmc_pdf(vp_is,Xa,0,logflag) + log(w_vp);
else
    templpdf(:,1) = -Inf;
end

% Fixed log weight for importance sampling (log fixed integrand)
if isamplevp_flag
    vlnpdf = max(vbmc_pdf(vp,Xa,0,1),log(realmin));
    lny = acqfun('islogf1',vlnpdf,[],[],fmu,fs2);
else
    lny = acqfun('islogf1',[],[],[],fmu,fs2);
end
% lny = lny - warpvars_vbmc(Xa,'logp',vp.trinfo);

% Mixture of box-uniforms
if w_vp < 1
    VV = prod(2*rect_delta);

    for ii = 1:N
        templpdf(:,ii+1) = log(all(abs(bsxfun(@minus,Xa,gp.X(ii,:))) < rect_delta,2) / VV / N * (1 - w_vp));
    end

    mmax = max(templpdf,[],2);
    lpdf = log(sum(exp(templpdf - mmax),2));
    lnw = bsxfun(@minus,lny,lpdf + mmax);
else
    lnw = bsxfun(@minus,lny,templpdf);
end

end

%--------------------------------------------------------------------------
function y = log_isbasefun(x,acqfun,gp,vp)
%LOG_ISBASEFUN Base importance sampling proposal log pdf

[fmu,fs2] = gplite_pred(gp,x);
if isempty(vp)
    y = acqfun('islogf',[],[],[],fmu,fs2);
else
    vlnpdf = max(vbmc_pdf(vp,x,0,1),log(realmin));
    y = acqfun('islogf',vlnpdf,[],[],fmu,fs2);
end

end

%--------------------------------------------------------------------------
%SQ_DIST Compute matrix of all pairwise squared distances between two sets
% of vectors, stored in the columns of the two matrices, a (of size n-by-D)
% and b (of size m-by-D).
function C = sq_dist(a,b)

n = size(a,1);
m = size(b,1);
mu = (m/(n+m))*mean(b,1) + (n/(n+m))*mean(a,1);
a = bsxfun(@minus,a,mu); b = bsxfun(@minus,b,mu);
C = bsxfun(@plus,sum(a.*a,2),bsxfun(@minus,sum(b.*b,2)',2*a*b'));
C = max(C,0);

end

%--------------------------------------------------------------------------
function sampleopts = get_mcmcopts(Ns,thin,burnin)
%GET_MCMCOPTS Get standard MCMC options.

if nargin < 2 || isempty(thin); thin = 1; end
if nargin < 3; burnin = []; end

sampleopts.Thin = thin;
if isempty(burnin)
    sampleopts.Burnin = ceil(sampleopts.Thin*Ns/2);
else
    sampleopts.Burnin = burnin;
end
sampleopts.Display = 'off';
sampleopts.Diagnostics = false;
sampleopts.VarTransform = false;
sampleopts.InversionSample = false;
sampleopts.FitGMM = false;

end

%--------------------------------------------------------------------------
function x = catrnd(p,n)
%CATRND Sample from categorical distribution.

maxel = 1e6;
Nel = n*numel(p);
stride = ceil(maxel/numel(p));

cdf(1,:) = cumsum(p);
u = rand(n,1)*cdf(end);

% Split for memory reasons
if Nel <= maxel
    x = sum(bsxfun(@lt, cdf, u),2) + 1;
else
    x = zeros(n,1);
    idx_min = 1;
    while idx_min <= n
        idx_max = min(idx_min+stride-1,n);
        idx = idx_min:idx_max;
        x(idx) = sum(bsxfun(@lt, cdf, u(idx)),2) + 1;
        idx_min = idx_max+1;
    end
end

end

author: acerbilab

repo: vbmc

save_path: github-repos/MATLAB/acerbilab-vbmc


text:
% unique_bc - unique backward compatible with Matlab versions prior to 2013a

function [C,IA,IB] = unique_bc(A,varargin);

errorFlag = error_bc;

v = version;
indp = find(v == '.');
v = str2num(v(1:indp(2)-1));
if v > 7.19, v = floor(v) + rem(v,1)/10; end;

if nargin > 2
    ind = strmatch('legacy', varargin);
    if ~isempty(ind)
        varargin(ind) = [];
    end;
end;

if v >= 7.14
    [C,IA,IB] = unique(A,varargin{:},'legacy');
    if errorFlag
        [C2,IA2] = unique(A,varargin{:});
        if ~isequal(C, C2) || ~isequal(IA, IA2) || ~isequal(IB, IB2)
            warning('backward compatibility issue with call to unique function');
        end;
    end;
else
    [C,IA,IB] = unique(A,varargin{:});
end;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function filum_test02 ( )

%*****************************************************************************80
%
%% TEST02 tests FILE_CHAR_COUNT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 June 2007
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST02\n' );
  fprintf ( 1, '  FILE_CHAR_COUNT counts the characters in a file.\n' );

  filename = 'story.txt';

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Examining file:\n' );
  fprintf ( 1, '    "%s".', filename );
  fprintf ( 1, '\n' );

  char_num = file_char_count ( filename );
  fprintf ( 1, '  Number of characters in file: %d\n', char_num );

  return
end


text:
function i4mat_l1_inverse_test ( )

%*****************************************************************************80
%
%% I4MAT_L1_INVERSE_TEST tests I4MAT_L1_INVERSE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 April 2009
%
%  Author:
%
%    John Burkardt
%

  n = 6;
%
%  Each row of this definition is a COLUMN of the matrix.
%
  a = [ ...
     1,  2,  0,  5,  0, 75; ...
     0,  1,  0,  0,  0,  0; ...
     0,  0,  1,  3,  0,  0; ...
     0,  0,  0,  1,  0,  6; ...
     0,  0,  0,  0,  1,  4; ...
     0,  0,  0,  0,  0,  1 ]';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'I4MAT_L1_INVERSE_TEST\n' );
  fprintf ( 1, '  I4MAT_L1_INVERSE inverts a unit lower triangular matrix.\n' );

  i4mat_print ( n, n, a, '  The original matrix:' );

  b = i4mat_l1_inverse ( n, a );

  i4mat_print ( n, n, b, '  The inverse matrix:' );

  c(1:n,1:n) = a(1:n,1:n) * b(1:n,1:n);

  i4mat_print ( n, n, c, '  The product:' );

  return
end