text:
% surrogdistrib - Build surrogate distribution
%
% surrog = surrogdistrib(data, varargin);
%
% Inputs:
%  data   - [cell] data arrays for which to compute a surrogate
%           distribution.
%
% Optional inputs:
%  'method'  - ['bootstrap'|'perm'] use either 'bootstrap' or 'permutation'
%              method. Bootstrap performs draws with replacement and
%              permutation performs draws without replacement. Default
%              is 'perm'.
%  'pairing' - ['on'|'off'] pair the data arrays.
%  'naccu'   - [integer] number of surrogate. Default is 1.
%  'precomp' - cell array containing precomputed value for speeding up
%              mulitple calls
%
% Output:
%  surrog  - surrogate distribution
%  precomp - cell array containing precomputed value for speeding up
%            mulitple calls
%
% Author: Arnaud Delorme, SCCN/INC/UCSD, La Jolla, 2005-

% Copyright (C) Arnaud Delorme
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

function [res precomp ] = surrogdistrib(data, varargin)

if nargin < 1
    help surrogdistrib,
    return;
end;

if ~strcmpi(varargin{1}, 'precomp')
    opt = finputcheck(varargin, { 'naccu'      'integer'   [1 Inf]             1;
                                  'method'     'string'    { 'perm','permutation','bootstrap' }  'perm';
                                  'pairing'    'string'    { 'on','off' }      'on';
                                  'precomp'    'cell'      {} {} }, 'surrogdistrib');
    if isstr(opt), error(opt); end;
    if strcmpi(opt.method, 'permutation'), opt.method = 'perm'; end;
    if strcmpi(opt.method, 'bootstrap'), bootflag = 1;
    else                                 bootflag = 0;
    end;
    if strcmpi(opt.pairing, 'on')
         pairflag = 1;
    else pairflag = 0;
    end;
else
    opt.precomp = varargin{2};
end;

% concatenate data
% ----------------
if isempty(opt.precomp)
    [ datavals datalen datadims ] = concatdata( data );
    precomp = { datavals datalen datadims bootflag pairflag opt.naccu};
else
    precomp   = opt.precomp;
    datavals  = precomp{1};
    datalen   = precomp{2};
    datadims  = precomp{3};
    bootflag  = precomp{4};
    pairflag  = precomp{5};
    opt.naccu = precomp{6};
end;

% compute surrogate distribution
% ------------------------------
if opt.naccu > 1
    res = supersurrogate( datavals, datalen, datadims, bootflag, pairflag, opt.naccu);
else
    res = surrogate( datavals, datalen, datadims, bootflag, pairflag);
end;

function res = supersurrogate(dat, lens, dims, bootstrapflag, pairedflag, naccu); % for increased speed only shuffle half the indices

    % recompute indices in set and target cell indices
    % ------------------------------------------------
    ncond = length(lens)-1;
    nsubj = lens(2);
    if bootstrapflag
        if pairedflag
             indswap  = mod( repmat([1:lens(end)],[naccu 1]) + ceil(rand(naccu,lens(end))*length(lens))*lens(2)-1, lens(end) )+1;
        else indswap  = ceil(rand(naccu,lens(end))*lens(end));
        end;
    else
        if pairedflag
            [tmp idx] = sort(rand(naccu,nsubj,ncond),3);
            indswap   = ((idx)-1)*nsubj + repmat( repmat([1:nsubj], [naccu 1 1]),[1 1 ncond]);
            indswap   = reshape(indswap, [naccu lens(end)]);
        else
            [tmp indswap] = sort(rand(naccu, lens(end)),2);
        end;
    end;

    for i = 1:length(lens)-1
        if myndims(dat) == 1
             res{i} = reshape(dat(indswap(:,lens(i)+1:lens(i+1))), naccu, lens(i+1)-lens(i));
        else res{i} = reshape(dat(:,indswap(:,lens(i)+1:lens(i+1))), size(dat,1), naccu, lens(i+1)-lens(i));
        end;
    end;
    res = reshape(res, dims);

function res = surrogate(dataconcat, lens, dims, bootstrapflag, pairedflag); % for increased speed only shuffle half the indices

    % recompute indices in set and target cell indices
    % ------------------------------------------------
    if bootstrapflag
        if pairedflag
             indswap  = mod( [1:lens(end)]+ ceil(rand(1,lens(end))*length(lens))*lens(2)-1, lens(end) )+1;
        else indswap  = ceil(rand(1,lens(end))*lens(end));
        end;
    else
        if pairedflag
            indswap  = [1:lens(end)];
            indswap  = reshape(indswap, [lens(2) length(lens)-1]);
            for i = 1:size(indswap,1) % shuffle each row
                [tmp idx] = sort(rand(1,size(indswap,2)));
                indswap(i,:) = indswap(i,idx);
            end;
            indswap  = reshape(indswap, [1 lens(2)*(length(lens)-1)]);
        else
            oriindices = [1:lens(end)]; % just shuffle indices
            [tmp idx] = sort(rand(1,length(oriindices)));
            indswap   = oriindices(idx);
        end;
    end;

    res = {};
    for i = 1:length(lens)-1
        if myndims(dataconcat) == 1
             res{i} = dataconcat(indswap(lens(i)+1:lens(i+1)));
        else res{i} = dataconcat(:,indswap(lens(i)+1:lens(i+1)));
        end;
    end;
    res = reshape(res, dims);

function val = myndims(a)
    if ndims(a) > 2
        val = ndims(a);
    else
        if size(a,1) == 1,
            val = 2;
        elseif size(a,2) == 1,
            val = 1;
        else
            val = 2;
        end;
    end;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function Rx=correlmx(x,p,Rxtype);
% Rx=correlmx(x,p,Rxtype) correlation matrix of a signal
%
%        Rx : correlation matrix (p+1) x (p+1)
%         x : analyzed signal
%         p : last autocorrelation lag
%    Rxtype : computation algorithm (default : 'fbhermitian')
%             possible values : 'hermitian', 'fbhermitian', 'burg' or 'fbburg'
%
% example :
%
% N=100; sig=real(fmconst(N,0.1))+0.4*randn(N,1);
% Rx=correlmx(sig,2,'burg'); [v,d] = eig(Rx), acos(-0.5*v(2,1)/v(1,1))/(2*pi)
% Rx=correlmx(sig,2,'hermitian'); [v,d] = eig(Rx), acos(-0.5*v(2,1)/v(1,1))/(2*pi)

% F. Auger, july 1998.

if (nargin<2),
 error('At least two parameters required');
elseif (nargin==2),
 Rxtype='fbhermitian';
end;

[L,xcol]=size(x);
if xcol>1,
 error('x must be a column vector');
elseif p>L,
 error('L must be greater than p');
elseif p<1,
 error('p must be greater than 0');
end;

Rxtype=upper(Rxtype);
if strcmp(Rxtype,'HERMITIAN')|strcmp(Rxtype,'FBHERMITIAN'),
 vector=x(p+1-(0:p)); Rx=conj(vector) * vector.';
 for t=p+2:L,
  vector=x(t-(0:p)); Rx=Rx+conj(vector) * vector.';
 end;
 Rx=Rx/(L-p);

elseif strcmp(Rxtype,'BURG')|strcmp(Rxtype,'FBBURG'),
 R0=sum(abs(x).^2)/L; % variance
 Rpos=zeros(1,p); Rneg=zeros(1,p);
 for n=1:p,
  Rpos(n)=sum(x(n+1:L).*conj(x(1:L-n)))/(L-n);
  Rneg(n)=sum(x(1:L-n).*conj(x(n+1:L)))/(L-n);
 end;
 Rx=toeplitz([R0 Rpos],[R0 Rneg]);
else error(['unknown algorithm name' Rxtype]);
end;

if strcmp(Rxtype(1:2),'FB'),
 Rx=0.5*(Rx+Rx');
end;
%  This program is free software; you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published by
%  the Free Software Foundation; either version 2 of the License, or
%  (at your option) any later version.
%
%  This program is distributed in the hope that it will be useful,
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
%
%  You should have received a copy of the GNU General Public License
%  along with this program; if not, write to the Free Software
%  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA

author: HeLiangHIT


text:
function [img, h] = rfDensityPlot(pol, ecc, sigma, varargin);
%
% [img, plotHandle] = rfDensityPlot(pol, ecc, sigma, [options]);
%
% Plot the density of estimated population receptive fields
% in the visual field.
%
% ras, 10/2007.
if nargin < 2, error('Not enough input arguments.');    end
if notDefined('sigma'), sigma = ones(size(pol));        end

%% Params
xRng = -15:.1:15;  % sampling rate along X axis
yRng = -15:.1:15;  % sampling rate along Y axis
plotFlag = 1;    % flag to plot the image
polRadians = 0;  % flag: if 1, polar angle is in radians CCW from 3-o-clock
                 % (as per mathematical measures of angle); if 0,
                 % polar angle is degrees CW from 12-o-clock (as is
                 % convenient for describing retinotopy)

%% Parse options
for i = 1:2:length(varargin)
    eval( sprintf('%s = %s', varargin{i}, num2str(varargin{i+1})) );
end

%% Remove NaNs and Infs
ok = find( ~isnan(pol) & ~isinf(pol) & ~isnan(ecc) & ~isinf(ecc) );
pol = pol(ok);
ecc = ecc(ok);
sigma = sigma(:,ok);

%% Make the image
% Get sampling grid of visual field
[X Y] = meshgrid(xRng, yRng);

% initialize a blank image to match the sampling grid
img = zeros(size(X));

% convert polar angle and eccentricity into Cartesian coords
if polRadians==1
    [x0 y0] = pol2cart(pol, ecc);
else
    [x0 y0] = pol2cart( deg2rad(90-pol), ecc );
end

hwait = mrvWaitbar(0, 'Generating RF Density Plot...');

for v = 1:length(x0)
    % compute an RF for this voxel
    RF = rfGaussian2d(X(:), Y(:), sigma(v), sigma(v), 0, x0(v), y0(v));
    RF = reshape(RF, size(img));

    % rescale RF so that area = 1
    RF = RF ./ sum( RF(:) );

    % positive Y should map to up instead of down (MATLAB convention):
    RF = flipud(RF);

    % add to density map
    img = img + RF;

	mrvWaitbar(v/length(x0), hwait);
end

close(hwait)

%% Show the image
if plotFlag==1
    mu = mean(img(:));
    sig = std(img(:));
    clim = [mu-sig mu+sig];
    h = imagesc(image, clim);
    colormap gray;
else
    h = [];
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function test16
%TEST16 test cholmod2 on a large matrix
% Example:
%   test16
% See also cholmod_test

% Copyright 2006-2007, Timothy A. Davis, University of Florida

fprintf ('=================================================================\n');
fprintf ('test16: test cholmod2 on a large matrix\n') ;

rand ('state',1) ;
randn ('state',1) ;

Prob = UFget (936)							    %#ok
A = Prob.A ;
% tic
% [L,s,p] = lchol (A) ;
% toc
% norm (L,1)

n = size (A,1) ;
b = rand (n,1) ;
tic
x = cholmod2(A,b) ;
t = toc ;
fprintf ('time %g\n', t) ;
err = norm (A*x-b) ;

if (err > 1e-5)
    error ('!') ;
end

fprintf ('test16 passed\n') ;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [net, stats] = our_cnn_train(net, imdb, getBatch, varargin)
%CNN_TRAIN  An example implementation of SGD for training CNNs
%    CNN_TRAIN() is an example learner implementing stochastic
%    gradient descent with momentum to train a CNN. It can be used
%    with different datasets and tasks by providing a suitable
%    getBatch function.
%
%    The function automatically restarts after each training epoch by
%    checkpointing.
%
%    The function supports training on CPU or on one or more GPUs
%    (specify the list of GPU IDs in the `gpus` option).

% Copyright (C) 2014-16 Andrea Vedaldi.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).
addpath(fullfile(vl_rootnn, 'examples'));

opts.expDir = fullfile('data','exp') ;
opts.continue = true ;
opts.batchSize = 256 ;
opts.numSubBatches = 1 ;
opts.train = [] ;
opts.val = [] ;
opts.gpus = [] ;
opts.epochSize = inf;
opts.prefetch = false ;
opts.numEpochs = 300 ;
opts.learningRate = 0.001 ;
opts.weightDecay = 0.0005 ;

opts.solver = [] ;  % Empty array means use the default SGD solver
[opts, varargin] = vl_argparse(opts, varargin) ;
if ~isempty(opts.solver)
  assert(isa(opts.solver, 'function_handle') && nargout(opts.solver) == 2,...
    'Invalid solver; expected a function handle with two outputs.') ;
  % Call without input arguments, to get default options
  opts.solverOpts = opts.solver() ;
end

opts.momentum = 0.9 ;
opts.saveSolverState = true ;
opts.nesterovUpdate = false ;
opts.randomSeed = 0 ;
opts.memoryMapFile = fullfile(tempdir, 'matconvnet.bin') ;
opts.profile = false ;
opts.parameterServer.method = 'mmap' ;
opts.parameterServer.prefix = 'mcn' ;

opts.conserveMemory = false ;
opts.backPropDepth = +inf ;
opts.sync = false ;
opts.cudnn = true ;

opts.errorFunction ='binary';

if(strcmp(net.layers{end}.type,'ourloss_softmaxlog'))
    opts.errorFunction ='multiclass';
end

opts.errorLabels = {} ;
opts.plotDiagnostics = false ;
opts.plotStatistics = true;
opts.postEpochFn = [] ;  % postEpochFn(net,params,state) called after each epoch; can return a new learning rate, 0 to stop, [] for no change
opts = vl_argparse(opts, varargin) ;

if ~exist(opts.expDir, 'dir'), mkdir(opts.expDir) ; end
if isempty(opts.train), opts.train = find(imdb.images.set==1) ; end
if isempty(opts.val), opts.val = find(imdb.images.set==2) ; end
if isscalar(opts.train) && isnumeric(opts.train) && isnan(opts.train)
  opts.train = [] ;
end
if isscalar(opts.val) && isnumeric(opts.val) && isnan(opts.val)
  opts.val = [] ;
end

% -------------------------------------------------------------------------
%                                                            Initialization
% -------------------------------------------------------------------------

net = vl_simplenn_tidy(net); % fill in some eventually missing values
net.layers{end-1}.precious = 1; % do not remove predictions, used for error
vl_simplenn_display(net, 'batchSize', opts.batchSize) ;

evaluateMode = isempty(opts.train) ;
if ~evaluateMode
  for i=1:numel(net.layers)
    J = numel(net.layers{i}.weights) ;
    if ~isfield(net.layers{i}, 'learningRate')
      net.layers{i}.learningRate = ones(1, J) ;
    end
    if ~isfield(net.layers{i}, 'weightDecay')
      net.layers{i}.weightDecay = ones(1, J) ;
    end
  end
end

% setup error calculation function
hasError = true ;
if isstr(opts.errorFunction)
  switch opts.errorFunction
    case 'none'
      opts.errorFunction = @error_none ;
      hasError = false ;
    case 'multiclass'
      opts.errorFunction = @error_multiclass ;
      if isempty(opts.errorLabels), opts.errorLabels = {'top1err'} ; end
    case 'binary'
      opts.errorFunction = @error_binary ;
      if isempty(opts.errorLabels), opts.errorLabels = {'binerr'} ; end
    otherwise
      error('Unknown error function ''%s''.', opts.errorFunction) ;
  end
end

state.getBatch = getBatch ;
stats = [] ;

% -------------------------------------------------------------------------
%                                                        Train and validate
% -------------------------------------------------------------------------

modelPath = @(ep) fullfile(opts.expDir, sprintf('net-epoch-%d.mat', ep));
modelFigPath = fullfile(opts.expDir, 'net-train.pdf') ;

start = opts.continue * findLastCheckpoint(opts.expDir) ;
if start >= 1
  fprintf('%s: resuming by loading epoch %d\n', mfilename, start) ;
  [net, state, stats] = loadState(modelPath(start)) ;
else
  state = [] ;
end

for epoch=start+1:opts.numEpochs

  % Set the random seed based on the epoch and opts.randomSeed.
  % This is important for reproducibility, including when training
  % is restarted from a checkpoint.

  rng(epoch + opts.randomSeed) ;
  prepareGPUs(opts, epoch == start+1) ;

  % Train for one epoch.
  params = opts ;
  params.epoch = epoch ;
  params.learningRate = opts.learningRate(min(epoch, numel(opts.learningRate))) ;
  params.train = opts.train(randperm(numel(opts.train))) ; % shuffle
  params.train = params.train(1:min(opts.epochSize, numel(opts.train)));
  params.val = opts.val(randperm(numel(opts.val))) ;
  params.imdb = imdb ;
  params.getBatch = getBatch ;

  if numel(params.gpus) <= 1
    [net, state] = processEpoch(net, state, params, 'train') ;
    [net, state] = processEpoch(net, state, params, 'val') ;
    if ~evaluateMode
      saveState(modelPath(epoch), net, state) ;
    end
    lastStats = state.stats ;
  else
    spmd
      [net, state] = processEpoch(net, state, params, 'train') ;
      [net, state] = processEpoch(net, state, params, 'val') ;
      if labindex == 1 && ~evaluateMode
        saveState(modelPath(epoch), net, state) ;
      end
      lastStats = state.stats ;
    end
    lastStats = accumulateStats(lastStats) ;
  end

  stats.train(epoch) = lastStats.train ;
  stats.val(epoch) = lastStats.val ;
  clear lastStats ;
  if ~evaluateMode
    saveStats(modelPath(epoch), stats) ;
  end

  if params.plotStatistics
    switchFigure(1) ; clf ;
    plots = setdiff(...
      cat(2,...
      fieldnames(stats.train)', ...
      fieldnames(stats.val)'), {'num', 'time'}) ;
    for p = plots
      p = char(p) ;
      values = zeros(0, epoch) ;
      leg = {} ;
      for f = {'train', 'val'}
        f = char(f) ;
        if isfield(stats.(f), p)
          tmp = [stats.(f).(p)] ;
          values(end+1,:) = tmp(1,:)' ;
          leg{end+1} = f ;
        end
      end
      subplot(1,numel(plots),find(strcmp(p,plots))) ;
      plot(1:epoch, values','o-') ;
      xlabel('epoch') ;
      title(p) ;
      legend(leg{:}) ;
      grid on ;
    end
    drawnow ;
    print(1, modelFigPath, '-dpdf') ;
  end

  if ~isempty(opts.postEpochFn)
    if nargout(opts.postEpochFn) == 0
      opts.postEpochFn(net, params, state) ;
    else
      lr = opts.postEpochFn(net, params, state) ;
      if ~isempty(lr), opts.learningRate = lr; end
      if opts.learningRate == 0, break; end
    end
  end
end

% With multiple GPUs, return one copy
if isa(net, 'Composite'), net = net{1} ; end

% -------------------------------------------------------------------------
function err = error_multiclass(params, labels, res)
% -------------------------------------------------------------------------
predictions = gather(res(end-1).x) ;
[~,predictions] = sort(predictions, 3, 'descend') ;
[~,labels]=max(labels,[],3);
error = ~bsxfun(@eq, predictions, labels) ;
err=sum(sum(sum(error(:,:,1,:),1),2),4) ;

% -------------------------------------------------------------------------
function err = error_binary(params, labels, res)
% -------------------------------------------------------------------------
predictions = gather(res(end-1).x) ;
error = bsxfun(@times, predictions, labels) < 0 ;
err = sum(error(:));
err = err/size(labels,3);

% -------------------------------------------------------------------------
function err = error_none(params, labels, res)
% -------------------------------------------------------------------------
err = zeros(0,1) ;

% -------------------------------------------------------------------------
function [net, state] = processEpoch(net, state, params, mode)
% -------------------------------------------------------------------------
% Note that net is not strictly needed as an output argument as net
% is a handle class. However, this fixes some aliasing issue in the
% spmd caller.

% initialize with momentum 0
if isempty(state) || isempty(state.solverState)
  for i = 1:numel(net.layers)
    state.solverState{i} = cell(1, numel(net.layers{i}.weights)) ;
    state.solverState{i}(:) = {0} ;
  end
end

% move CNN  to GPU as needed
numGpus = numel(params.gpus) ;
if numGpus >= 1
  net =our_vl_simplenn_move(net, 'gpu') ;
  for i = 1:numel(state.solverState)
    for j = 1:numel(state.solverState{i})
      s = state.solverState{i}{j} ;
      if isnumeric(s)
        state.solverState{i}{j} = gpuArray(s) ;
      elseif isstruct(s)
        state.solverState{i}{j} = structfun(@gpuArray, s, 'UniformOutput', false) ;
      end
    end
  end
end
if numGpus > 1
  parserv = ParameterServer(params.parameterServer) ;
  vl_simplenn_start_parserv(net, parserv) ;
else
  parserv = [] ;
end

% profile
if params.profile
  if numGpus <= 1
    profile clear ;
    profile on ;
  else
    mpiprofile reset ;
    mpiprofile on ;
  end
end

subset = params.(mode) ;
num = 0 ;
stats.num = 0 ; % return something even if subset = []
stats.time = 0 ;
adjustTime = 0 ;
res = [] ;
error = [] ;

if(size(params.imdb.images.labels,1)>1)
    density=mean(params.imdb.images.labels>0,2);
else
    density=0;
end

start = tic ;
for t=1:params.batchSize:numel(subset)
  fprintf('%s: epoch %02d: %3d/%3d:', mode, params.epoch, ...
          fix((t-1)/params.batchSize)+1, ceil(numel(subset)/params.batchSize)) ;
  batchSize = min(params.batchSize, numel(subset) - t + 1) ;
  for s=1:params.numSubBatches
    % get this image batch and prefetch the next
    batchStart = t + (labindex-1) + (s-1) * numlabs ;
    batchEnd = min(t+params.batchSize-1, numel(subset)) ;
    batch = subset(batchStart : params.numSubBatches * numlabs : batchEnd) ;
    num = num + numel(batch) ;
    if numel(batch) == 0, continue ; end

    [im, labels] = params.getBatch(params.imdb, batch) ;

    if params.prefetch
      if s == params.numSubBatches
        batchStart = t + (labindex-1) + params.batchSize ;
        batchEnd = min(t+2*params.batchSize-1, numel(subset)) ;
      else
        batchStart = batchStart + numlabs ;
      end
      nextBatch = subset(batchStart : params.numSubBatches * numlabs : batchEnd) ;
      params.getBatch(params.imdb, nextBatch) ;
    end

    if numGpus >= 1
      im = gpuArray(im) ;
    end

    if strcmp(mode, 'train')
      dzdy = 1 ;
      evalMode = 'normal' ;
    else
      dzdy = [] ;
      evalMode = 'test' ;
    end
    net.layers{end}.class = labels ;
    net.layers{end}.iter=params.epoch;
    net.layers{end}.density=density;
    res=our_vl_simplenn(net, im, dzdy, res, ...
                      'accumulate', s ~= 1, ...
                      'mode', evalMode, ...
                      'conserveMemory', params.conserveMemory, ...
                      'backPropDepth', params.backPropDepth, ...
                      'sync', params.sync, ...
                      'cudnn', params.cudnn, ...
                      'parameterServer', parserv, ...
                      'holdOn', s < params.numSubBatches) ;

    % accumulate errors
    error = sum([error, [...
      sum(double(gather(res(end).x))) ;
      reshape(params.errorFunction(params, labels, res),[],1) ; ]],2) ;
  end

  % accumulate gradient
  if strcmp(mode, 'train')
    if ~isempty(parserv), parserv.sync() ; end
    [net, res, state] = accumulateGradients(net, res, state, params, batchSize, parserv) ;

    labelNum=size(net.layers{end}.class,3);
    for lay=1:numel(net.layers) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        if(strcmp(net.layers{lay}.type,'conv_mask'))
            if(isempty(net.layers{lay}.sliceMag))
                net.layers{lay}.sliceMag=zeros(size(res(lay).x,3),labelNum);
            end
            for lab=1:labelNum
                tmp=gather(max(max(res(lay).x(:,:,:,net.layers{end}.class(:,:,lab,:)==1),[],1),[],2));
                if(~isempty(tmp))
                    %meantmp=sum(tmp,4)./max(sum(tmp>0,4),1);
                    meantmp=mean(tmp,4);

                    if(sum(net.layers{lay}.sliceMag(:,lab))==0)
                        net.layers{lay}.sliceMag(:,lab)=max(meantmp(:),0.1);
                    else
                        tmptmp=0.9;
                        meantmp(meantmp==0)=net.layers{lay}.sliceMag(meantmp==0);
                        net.layers{lay}.sliceMag(:,lab)=net.layers{lay}.sliceMag(:,lab).*tmptmp+meantmp(:).*(1-tmptmp);
                    end
                end
            end
        end
    end
  end

  % get statistics
  time = toc(start) + adjustTime ;
  batchTime = time - stats.time ;
  stats = extractStats(net, params, error / num) ;
  stats.num = num ;
  stats.time = time ;
  currentSpeed = batchSize / batchTime ;
  averageSpeed = (t + batchSize - 1) / time ;
  if t == 3*params.batchSize + 1
    % compensate for the first three iterations, which are outliers
    adjustTime = 4*batchTime - time ;
    stats.time = time + adjustTime ;
  end

  fprintf(' %.1f (%.1f) Hz', averageSpeed, currentSpeed) ;
  for f = setdiff(fieldnames(stats)', {'num', 'time'})
    f = char(f) ;
    fprintf(' %s: %.3f', f, stats.(f)) ;
  end
  fprintf('\n') ;

  % collect diagnostic statistics
  if strcmp(mode, 'train') && params.plotDiagnostics
    switchFigure(2) ; clf ;
    diagn = [res.stats] ;
    diagnvar = horzcat(diagn.variation) ;
    diagnpow = horzcat(diagn.power) ;
    subplot(2,2,1) ; barh(diagnvar) ;
    set(gca,'TickLabelInterpreter', 'none', ...
      'YTick', 1:numel(diagnvar), ...
      'YTickLabel',horzcat(diagn.label), ...
      'YDir', 'reverse', ...
      'XScale', 'log', ...
      'XLim', [1e-5 1], ...
      'XTick', 10.^(-5:1)) ;
    grid on ; title('Variation');
    subplot(2,2,2) ; barh(sqrt(diagnpow)) ;
    set(gca,'TickLabelInterpreter', 'none', ...
      'YTick', 1:numel(diagnpow), ...
      'YTickLabel',{diagn.powerLabel}, ...
      'YDir', 'reverse', ...
      'XScale', 'log', ...
      'XLim', [1e-5 1e5], ...
      'XTick', 10.^(-5:5)) ;
    grid on ; title('Power');
    subplot(2,2,3); plot(squeeze(res(end-1).x)) ;
    drawnow ;
  end
end

% Save back to state.
state.stats.(mode) = stats ;
if params.profile
  if numGpus <= 1
    state.prof.(mode) = profile('info') ;
    profile off ;
  else
    state.prof.(mode) = mpiprofile('info');
    mpiprofile off ;
  end
end
if ~params.saveSolverState
  state.solverState = [] ;
else
  for i = 1:numel(state.solverState)
    for j = 1:numel(state.solverState{i})
      s = state.solverState{i}{j} ;
      if isnumeric(s)
        state.solverState{i}{j} = gather(s) ;
      elseif isstruct(s)
        state.solverState{i}{j} = structfun(@gather, s, 'UniformOutput', false) ;
      end
    end
  end
end

net =our_vl_simplenn_move(net, 'cpu') ;

% -------------------------------------------------------------------------
function [net, res, state] = accumulateGradients(net, res, state, params, batchSize, parserv)
% -------------------------------------------------------------------------
numGpus = numel(params.gpus) ;
otherGpus = setdiff(1:numGpus, labindex) ;

for l=numel(net.layers):-1:1
  for j=numel(res(l).dzdw):-1:1

    if ~isempty(parserv)
      tag = sprintf('l%d_%d',l,j) ;
      parDer = parserv.pull(tag) ;
    else
      parDer = res(l).dzdw{j}  ;
    end

    if j == 3 && strcmp(net.layers{l}.type, 'bnorm')
      % special case for learning bnorm moments
      thisLR = net.layers{l}.learningRate(j) ;
      net.layers{l}.weights{j} = vl_taccum(...
        1 - thisLR, ...
        net.layers{l}.weights{j}, ...
        thisLR / batchSize, ...
        parDer) ;
    else
      % Standard gradient training.
      thisDecay = params.weightDecay * net.layers{l}.weightDecay(j) ;
      thisLR = params.learningRate * net.layers{l}.learningRate(j) ;

      if thisLR>0 || thisDecay>0
        % Normalize gradient and incorporate weight decay.
        parDer = vl_taccum(1/batchSize, parDer, ...
                           thisDecay, net.layers{l}.weights{j}) ;

        if isempty(params.solver)
          % Default solver is the optimised SGD.
          % Update momentum.
          state.solverState{l}{j} = vl_taccum(...
            params.momentum, state.solverState{l}{j}, ...
            -1, parDer) ;

          % Nesterov update (aka one step ahead).
          if params.nesterovUpdate
            delta = params.momentum * state.solverState{l}{j} - parDer ;
          else
            delta = state.solverState{l}{j} ;
          end

          % Update parameters.
          net.layers{l}.weights{j} = vl_taccum(...
            1, net.layers{l}.weights{j}, ...
            thisLR, delta) ;

        else
          % call solver function to update weights
          [net.layers{l}.weights{j}, state.solverState{l}{j}] = ...
            params.solver(net.layers{l}.weights{j}, state.solverState{l}{j}, ...
            parDer, params.solverOpts, thisLR) ;
        end
      end
    end

    % if requested, collect some useful stats for debugging
    if params.plotDiagnostics
      variation = [] ;
      label = '' ;
      switch net.layers{l}.type
        case {'conv','convt','conv_mask'}
          if isnumeric(state.solverState{l}{j})
            variation = thisLR * mean(abs(state.solverState{l}{j}(:))) ;
          end
          power = mean(res(l+1).x(:).^2) ;
          if j == 1 % fiters
            base = mean(net.layers{l}.weights{j}(:).^2) ;
            label = 'filters' ;
          else % biases
            base = sqrt(power) ;%mean(abs(res(l+1).x(:))) ;
            label = 'biases' ;
          end
          variation = variation / base ;
          label = sprintf('%s_%s', net.layers{l}.name, label) ;
      end
      res(l).stats.variation(j) = variation ;
      res(l).stats.power = power ;
      res(l).stats.powerLabel = net.layers{l}.name ;
      res(l).stats.label{j} = label ;
    end
  end
end

% -------------------------------------------------------------------------
function stats = accumulateStats(stats_)
% -------------------------------------------------------------------------

for s = {'train', 'val'}
  s = char(s) ;
  total = 0 ;

  % initialize stats stucture with same fields and same order as
  % stats_{1}
  stats__ = stats_{1} ;
  names = fieldnames(stats__.(s))' ;
  values = zeros(1, numel(names)) ;
  fields = cat(1, names, num2cell(values)) ;
  stats.(s) = struct(fields{:}) ;

  for g = 1:numel(stats_)
    stats__ = stats_{g} ;
    num__ = stats__.(s).num ;
    total = total + num__ ;

    for f = setdiff(fieldnames(stats__.(s))', 'num')
      f = char(f) ;
      stats.(s).(f) = stats.(s).(f) + stats__.(s).(f) * num__ ;

      if g == numel(stats_)
        stats.(s).(f) = stats.(s).(f) / total ;
      end
    end
  end
  stats.(s).num = total ;
end

% -------------------------------------------------------------------------
function stats = extractStats(net, params, errors)
% -------------------------------------------------------------------------
stats.lossForClassification = errors(1) ;
for i = 1:numel(params.errorLabels)
  stats.(params.errorLabels{i}) = errors(i+1) ;
end

% -------------------------------------------------------------------------
function saveState(fileName, net, state)
%if((fileName(end-4)=='5')||(fileName(end-4)=='0'))
% -------------------------------------------------------------------------
save(fileName, 'net', 'state') ;
%end

% -------------------------------------------------------------------------
function saveStats(fileName, stats)
%if((fileName(end-4)=='5')||(fileName(end-4)=='0'))
% -------------------------------------------------------------------------
if exist(fileName)
  save(fileName, 'stats', '-append') ;
else
  save(fileName, 'stats') ;
end
%end

% -------------------------------------------------------------------------
function [net, state, stats] = loadState(fileName)
% -------------------------------------------------------------------------
load(fileName, 'net', 'state', 'stats') ;
net = vl_simplenn_tidy(net) ;
if isempty(whos('stats'))
  error('Epoch ''%s'' was only partially saved. Delete this file and try again.', ...
        fileName) ;
end

% -------------------------------------------------------------------------
function epoch = findLastCheckpoint(modelDir)
% -------------------------------------------------------------------------
list = dir(fullfile(modelDir, 'net-epoch-*.mat')) ;
tokens = regexp({list.name}, 'net-epoch-([\d]+).mat', 'tokens') ;
epoch = cellfun(@(x) sscanf(x{1}{1}, '%d'), tokens) ;
epoch = max([epoch 0]) ;

% -------------------------------------------------------------------------
function switchFigure(n)
% -------------------------------------------------------------------------
if get(0,'CurrentFigure') ~= n
  try
    set(0,'CurrentFigure',n) ;
  catch
    figure(n) ;
  end
end

% -------------------------------------------------------------------------
function clearMex()
% -------------------------------------------------------------------------
%clear vl_tmove vl_imreadjpeg ;
disp('Clearing mex files') ;
clear mex ;
clear vl_tmove vl_imreadjpeg ;

% -------------------------------------------------------------------------
function prepareGPUs(params, cold)
% -------------------------------------------------------------------------
numGpus = numel(params.gpus) ;
if numGpus > 1
  % check parallel pool integrity as it could have timed out
  pool = gcp('nocreate') ;
  if ~isempty(pool) && pool.NumWorkers ~= numGpus
    delete(pool) ;
  end
  pool = gcp('nocreate') ;
  if isempty(pool)
    parpool('local', numGpus) ;
    cold = true ;
  end
end
if numGpus >= 1 && cold
  fprintf('%s: resetting GPU\n', mfilename) ;
  clearMex() ;
  if numGpus == 1
    disp(gpuDevice(params.gpus)) ;
  else
    spmd
      clearMex() ;
      disp(gpuDevice(params.gpus(labindex))) ;
    end
  end
end

repo: interpretableCNN


text:
function [x,mn,mx]=melbankm(p,n,fs,fl,fh,w)
if nargin < 6
    w='tz';
    if nargin < 5
        fh=0.5;
        if nargin < 4
            fl=0;
        end
    end
end
f0=700/fs;
fn2=floor(n/2);
lr=log((f0+fh)/(f0+fl))/(p+1);
bl=n*((f0+fl)*exp([0 1 p p+1]*lr)-f0);
b2=ceil(bl(2));
b3=floor(bl(3));
if any(w=='y')
    pf=log((f0+(b2:b3)/n)/(f0+fl))/lr;
    fp=floor(pf);
    r=[ones(1,b2) fp fp+1 p*ones(1,fn2-b3)];
    c=[1:b3+1 b2+1:fn2+1];
    v=2*[0.5 ones(1,b2-1) 1-pf+fp pf-fp ones(1,fn2-b3-1) 0.5];
    mn=1;
    mx=fn2+1;
else
    b1=floor(bl(1))+1;
    b4=min(fn2,ceil(bl(4)))-1;
    pf=log((f0+(b1:b4)/n)/(f0+fl))/lr;
    fp=floor(pf);
    pm=pf-fp;
    k2=b2-b1+1;
    k3=b3-b1+1;
    k4=b4-b1+1;
    r=[fp(k2:k4) 1+fp(1:k3)];
    c=[k2:k4 1:k3];
    v=2*[1-pm(k2:k4) pm(1:k3)];
    mn=b1+1;
    mx=b4+1;
end
if any(w=='n')
    v=1-cos(v*pi/2);
elseif any(w=='m')
    v=1-0.92/1.08*cos(v*pi/2);
end
if nargout > 1
    x=sparse(r,c,v);
else
    x=sparse(r,c+mn-1,v,p,1+fn2);
end

author: decouples


text:
function i = labindex()
i = 1 ;

author: willard-yuan

repo: cnn-for-image-retrieval


text:
function [success] = writeHDF5Data(inputData,fileSavePath,varargin)
	% Saves input data to a HDF5 file, tries to preserve datatype.
	% Biafra Ahanonu
	% started: 2013.11.01
	%
	% inputs
		% inputData: matrix, [x y frames] preferred.
		% fileSavePath: str, path where HDF5 file should be saved.
	% outputs
		% success = 1 if successful save, 0 if error.
	% options
		% datasetname = HDF5 hierarchy where data should be stored

	[success] = ciapkg.hdf5.writeHDF5Data(inputData,fileSavePath,'passArgs', varargin);
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
addpath /home/fowlkes/gpb_owt/grouping/lib/
outFile = '/home/fowlkes/gpb_owt/grouping/data/101087_gPb.mat_pbs.mat';
simFile = '/home/fowlkes/gpb_owt/grouping/data/101087_gPb.mat_pbs.mat_sim.tmp';
imgFile = '/home/fowlkes/gpb_owt/grouping/data/101087_gPb.mat_pbs.mat_img.jpg';
latFile = '/home/fowlkes/gpb_owt/grouping/data/101087_gPb.mat_pbs.mat_lat.tmp';
load /home/fowlkes/gpb_owt/grouping/data/mPb.mat  %load in mPb data.

%
% old code
%
l{1} = zeros(size(mPb, 1) + 1, size(mPb, 2));
l{1}(2:end, :) = mPb;
l{2} = zeros(size(mPb, 1), size(mPb, 2) + 1);
l{2}(:, 2:end) = mPb;
latFile = strcat(outFile, '_lat.tmp');
write_array(latFile, l);

system(sprintf('/home/fowlkes/gpb_owt/grouping/lib/segment -image %s -smatrixfile %s -readlattice true -latticefile %s -dthresh 5', imgFile, simFile, latFile));
Wold = read_smatrix(simFile);

%
% new code
%

% mex buildW.cpp -Iutil smatrix.cc ic.cc affinity.cc util/libutil.a

[x,y,z] = buildW(l{1},l{2});
Wnew = sparse(x,y,z);

sum(sum( (Wold-Wnew).^2 ))

author: s-gupta

repo: rcnn-depth

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


text:
function [parent_pop] = osy(parent_pop)
%   This procedure implements Osyczka and Kundu's function.
%   The canonical zdt1 function is defined as below --
%   f_1 = -25.0 * ((x_1 - 2)^2 - (x_2 - 2)^2 - (x_3 - 1)^2
%                           - (x_4 - 4)^2 - (x_5 - 1)^2)
%   f_2 = sum_{i=1}^n x_i^2.0
%   s.t.
%       c_1(x) = x_1 + x_2 - 2 >= 0
%       c_2(x) = 6 - x_1 - x_2 >= 0
%       c_3(x) = 2 - x_2 + x_1 >= 0
%       c_4(x) = 2 - x_1 + 3x_2 >= 0
%       c_5(x) = 4 - (x_3 - 3)^2 - x_4 >= 0
%       c_6(x) = (x_5 - 3)^2 + x_6 - 4 >= 0
%   where,
%       0 <= x_1 <= 10.0
%       0 <= x_2 <= 10.0
%       1 <= x_3 <= 5.0
%       0 <= x_4 <= 6.0
%       1 <= x_5 <= 5.0
%       0 <= x_6 <= 10.0

global nreal ;
global ncon ;
global nobj ;

cindex = nreal + nobj + 1 : nreal + nobj + ncon ;

% disp(parent_pop(1,:))

x = parent_pop(:,1:nreal);
f1 = -(25.0 .* ((x(:,1) - 2.0) .^ 2.0) + ((x(:,2) - 2.0) .^ 2.0) ...
             + ((x(:,3) - 1.0) .^ 2.0) + ((x(:,4) - 4.0) .^ 2.0) ...
             + ((x(:,5) - 1.0) .^ 2.0));
f2 = sum((x .^ 2.0), 2) ;

c = parent_pop(:,cindex) ;
% disp(c(1,:))
c(:,1) = ((x(:,1) + x(:,2)) ./ 2.0) - 1.0 ;
c(:,2) = 1.0 - ((x(:,1) + x(:,2)) ./ 6.0);
c(:,3) = 1.0 - (x(:,2) ./ 2.0) + (x(:,1) ./ 2.0);
c(:,4) = 1.0 - (x(:,1) ./ 2.0) + (3.0 .* (x(:,2) ./ 2.0));
c(:,5) = 1.0 - (((x(:,3) - 3.0) .^ 2.0) ./ 4.0) - (x(:,4) ./ 4.0);
c(:,6) = (((x(:,5) - 3.0) .^ 2.0) ./ 4.0) + (x(:,6) ./ 4.0) - 1.0;
parent_pop(:, (nreal+1)) = f1;
parent_pop(:, (nreal+2)) = f2;
parent_pop(:,cindex) = c ;
end

author: chudur-budur


text:
function [Q_gppca, Q_vbpca] = mohsst5_experiment_periodic(data, D, maxiter, ...
                                                  densityW, densityX, ...
                                                  updatepseudow, updatepseudox, ...
                                                  vbinit, covfuncset)

seed = ceil(sum(clock))
randn('state', seed)
rand('state', seed)

if nargin < 1
  maxiter = 10;
end
if nargin < 6
  updatepseudow = true;
end
if nargin < 7
  updatepseudox = true;
end
if nargin < 8
  vbinit = false;
end
if nargin < 9
  covfuncset = 3;
end

[M,N] = size(data.observations);

% Divide data into test and train sets
testsize = 20; % testset size in percents
Itest = load(sprintf('/share/bayes/data/jaakko/mohsst5/ind%dtest.mat',testsize));
Itrain = load(sprintf('/share/bayes/data/jaakko/mohsst5/ind%dtrain.mat',testsize));
Ytest = data.observations;
Ytest(Itrain.Itrain) = nan;
Yobs = data.observations;
Yobs(Itest.Itest) = nan;

datestring = datestr(now, 'yyyymmdd');

% Weight observations!
weights = sqrt(cosd(data.coordinates(2,:)));
Yobsw = bsxfun(@times, weights(:), Yobs);
Ytestw = bsxfun(@times, weights(:), Ytest);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

init = [];
switch vbinit
 case 0
  vbinitstring = '';
 case 1
  Qvb = pca_full(Yobsw,D,'maxiters',10,'rotate2pca',true, 'algorithm','vb');
  init.W = zeros(M,D);
  init.X = randn(D,N);
  init.W(:,1:10) = Qvb.A(:,1:10);
  init.X(1:10,:) = Qvb.S(1:10,:);
  vbinitstring = '_vbinit';
 case 2
  Qvb = pca_full(Yobsw,D,'maxiters',10,'rotate2pca',true, 'algorithm','vb');

  init.W = zeros(M,D);
  init.X = randn(D,N);

  samplerate = 1/30; % one sample in month

  % 3 very slow components
  lowpass = 1 ./ (10*365);
  [Af,Sf,S] = dss(Qvb.A, Qvb.S, Qvb.Sv, 2 * lowpass/samplerate, 3);
  init.W(:,1:3) = Af;
  init.X(1:3,:) = S;

  % 5 slow components
  lowpass = 1 ./ [10*365, 2*365];
  [Af,Sf,S] = dss(Qvb.A, Qvb.S, Qvb.Sv, 2 * lowpass/samplerate, 5);
  init.W(:,4:8) = Af;
  init.X(4:8,:) = S;

  % 2 periodical components
  bandpass = 1 ./ [2*365, 0.5*365];
  [Af,Sf,S] = dss(Qvb.A, Qvb.S, Qvb.Sv, 2 * bandpass/samplerate, 2);
  init.W(:,9:10) = Af;
  init.X(9:10,:) = S;

% $$$   % DEBUGGING STUFF
% $$$   Q_gppca = Af;%init.W;
% $$$   Q_vbpca.Xf = Sf;
% $$$   Q_vbpca.X = S;
% $$$   return
% $$$
  vbinitstring = '_vbdssinit';

end

% Covariance functions
% W distances in kilometers
for d=1:D
  switch 1 % select covariance functions
   case 1
    gpcovW = {@gpcovScale, @(logtheta,x1,x2) gpcov(logtheta,x1,x2,@ ...
                                                   sqdistEarth)};
    logthetaW = [log(2); log(1e3)];
    Wcov = 'se';
  end
  initthetaW{d} = logthetaW;
  covfunW{d} = gpcovW;
end

% X distances in days
for d=1:D
  switch covfuncset % select covariance functions
   case 1
    % 80 pediodical
    if d==1
      disp('Using 80 periodical');
    end
    gpcovX = {@gpcovProduct, @gpcov, @gpcovPeriodic};
    logthetaX = log([1e5, 0.5, 365*(D-d+3)]);
    Xcov = 'per';
   case 2
    % 80 piecewise polynomial
    if d==1
      disp('Using 80 piecewise polynomial');
    end
    gpcovX = @gpcovPP;
    logthetaX = log(30*(5+40*rand));
    Xcov = 'cs';
   case 3
    % 5 trend SE + 5 slow SE + 5 periodical + rest faster piecewise
    % polynomial
    if d==1
      disp(['Using 5 trends + 5 slow + 5 periodical + rest faster ' ...
            'piecewise polynomial']);
    end
    if d <= 5
      gpcovX = @gpcov;
      logthetaX = log(365*(10+10*rand)); % 10-20 years
      densityX(d) = 0.05;
    elseif d <= 10
      gpcovX = @gpcov;
      logthetaX = log(365*(2+8*rand)); % 2-10 years
      densityX(d) = 0.2;
    elseif d <= 15
      gpcovX = {@gpcovProduct, @gpcov, @gpcovPeriodic};
      logthetaX = log([200*365, 0.5, 365*(0.5+1.5*rand)]); % 0.5-2 years
      densityX(d) = 0.2;
    else
      gpcovX = @gpcovPP;
      logthetaX = log(30*(6+6*rand)); % 0.5-1 years
      densityX(d) = 1;
    end
    Xcov = 'se-per-cs';
   case 4
    % 3 slow piecewise polynomial + 5 middle + 2 periodical + rest faster
    % piecewise polynomial
    if d==1
      disp(['Using 3 trends + 5 slow + 2 periodical + rest faster ' ...
            'piecewise polynomial']);
    end
    if d <= 3
      gpcovX = @gpcov;
      logthetaX = log(365*(10+10*rand)); % 10-20 years
      densityX(d) = 0.05;
    elseif d <= 8
      gpcovX = @gpcov;
      logthetaX = log(365*(2+8*rand)); % 2-10 years
      densityX(d) = 0.2;
    elseif d <= 10
      gpcovX = {@gpcovProduct, @gpcov, @gpcovPeriodic};
      logthetaX = log([80*365, 0.5, 365*(0.5+1.5*rand)]); % 0.5-2 years
      densityX(d) = 0.2;
    else
      gpcovX = @gpcovPP;
      logthetaX = log(30*(6+6*rand)); % 0.5-1 years
      densityX(d) = 1;
    end
    Xcov = 'se-per-cs';
  end
  initthetaX{d} = logthetaX;
  covfunX{d} = gpcovX;
end

filename = sprintf(['/share/bayes/data/jaakko/mohsst5/' ...
                    'mohsst5_weighted%s_testsize=%d_D=%d_Wcov=%s_Wden=' ...
                    '%d_Xcov=%s_Xden=%d_iter=%d_date=%s'], vbinitstring, ...
                   testsize, D, Wcov, ceil(100*min(densityW)), Xcov, ...
                   ceil(100*min(densityX)), maxiter, datestring)

disp('-- Run GP PCA --')
Q_gppca = vbgppcamv(Yobsw, D,...
                    data.coordinates, data.time,...
                    covfunW, initthetaW,...
                    covfunX,initthetaX, ...
                    'maxiter',maxiter, ...
                    'pseudodensityx', densityX, ...
                    'pseudodensityw', densityW, ...
                    'loglikelihood', true, ...
                    'updatehyper', [1 6 21 51 101 201 501], ...
                    'updatepseudox', updatepseudox, ...
                    'updatepseudow', updatepseudow, ...
                    'maxsearchx', 3, ...
                    'maxsearchw', 3, ...
                    'checkgradx', false, ...
                    'checkgradw', false, ...
                    'init', init, ...
                    'autosavetime', 3600, ...
                    'autosavefile', filename);

% Error measures (using weights!)
Yh_gppca = Q_gppca.W * Q_gppca.X;
trainrmse_gppca = rmse(Yobsw, Yh_gppca)
testrmse_gppca = rmse(Ytestw, Yh_gppca)

% Huuuuuge matrices... Don't save them..
Q_gppca.CovXp = [];
Q_gppca.CovWp = [];

% INVERSE WEIGHTS!!!
Q_gppca.W = bsxfun(@rdivide, Q_gppca.W, weights(:));
Q_gppca.varW = bsxfun(@rdivide, Q_gppca.varW, weights(:).^2);
for d=1:D
  lat = Q_gppca.pseudoW{d}(2,:);
  Q_gppca.Wp{d} = Q_gppca.Wp{d} ./ sqrt(abs(cosd(lat(:))));
end

% Save
Q_gppca.filename = filename;
Q_gppca.seed = seed;
save(filename, '-struct', 'Q_gppca');
fprintf('Saved results in %s\n', filename);

% $$$ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% $$$ disp('-- Run PCA --')
% $$$
% $$$ % Alexander's vbpca
% $$$ disp('version by Alex')
% $$$ Q_vbpca = pca_full(Yobs,D,'maxiters',maxiter,'rotate2pca',true,
%'algorithm','vb');
% $$$ Yh_vbpca = Q_vbpca.A * Q_vbpca.S + repmat(Q_vbpca.Mu,1,N);
% $$$ vb_filename = sprintf(['/share/bayes/data/jaakko/mohsst5/' ...
% $$$
%'mohsst5_VBPCA_testsize=%d_D=%d_iter=%d_date=%s'], ...
% $$$                       testsize, D, maxiter, datestring);
% $$$ save(vb_filename, '-struct', 'Q_vbpca');
% $$$
% $$$ % My vbpca
% $$$ %disp('version by Jaakko')
% $$$ %Q_vbpca = vbpcamv(Yobs,D,'maxiters',maxiter);
% $$$ %Yh_vbpca = Q_vbpca.W * Q_vbpca.X + repmat(Q_vbpca.mu,1,N);
% $$$
% $$$ trainrmse_vbpca = rmse(Yobs, Yh_vbpca)
% $$$ testrmse_vbpca = rmse(Ytest, Yh_vbpca)

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function [merged unused] = propval(propvals, defaults, varargin)

% Create a structure combining property-value pairs with default values.
%
% [MERGED UNUSED] = PROPVAL(PROPVALS, DEFAULTS, ...)
%
% Given a cell array or structure of property-value pairs
% (i.e. from VARARGIN or a structure of parameters), PROPVAL will
% merge the user specified values with those specified in the
% DEFAULTS structure and return the result in the structure
% MERGED.  Any user specified values that were not listed in
% DEFAULTS are output as property-value arguments in the cell array
% UNUSED.  STRICT is disabled in this mode.
%
% ALTERNATIVE USAGE:
%
% [ ARGS ] = PROPVAL(PROPVALS, DEFAULTS, ...)
%
% In this case, propval will assume that no user specified
% properties are meant to be "picked up" and STRICT mode will be enforced.
%
% ARGUMENTS:
%
% PROPVALS - Either a cell array of property-value pairs
%   (i.e. {'Property', Value, ...}) or a structure of equivalent form
%   (i.e. struct.Property = Value), to be merged with the values in
%   DEFAULTS.
%
% DEFAULTS - A structure where field names correspond to the
%   default value for any properties in PROPVALS.
%
% OPTIONAL ARGUMENTS:
%
% STRICT (default = true) - Use strict guidelines when processing
%   the property value pairs.  This will warn the user if an empty
%   DEFAULTS structure is passed in or if there are properties in
%   PROPVALS for which no default was provided.
%
% EXAMPLES:
%
% Simple function with two optional numerical parameters:
%
% function [result] = myfunc(data, varargin)
%
%   defaults.X = 5;
%   defaults.Y = 10;
%
%   args = propvals(varargin, defaults)
%
%   data = data * Y / X;
%
% >> myfunc(data)
%    This will run myfunc with X=5, Y=10 on the variable 'data'.
%
% >> myfunc(data, 'X', 0)
%    This will run myfunc with X=0, Y=10 (thus giving a
%    divide-by-zero error)
%
% >> myfunc(data, 'foo', 'bar') will run myfunc with X=5, Y=10, and
%    PROPVAL will give a warning that 'foo' has no default value,
%    since STRICT is true by default.
%

% License:
%=====================================================================
%
% This is part of the Princeton MVPA toolbox, released under
% the GPL. See http://www.csbmb.princeton.edu/mvpa for more
% information.
%
% The Princeton MVPA toolbox is available free and
% unsupported to those who might find it useful. We do not
% take any responsibility whatsoever for any problems that
% you have related to the use of the MVPA toolbox.
%
% ======================================================================

% Backwards compatibility
pvdef.ignore_missing_default = false;
pvdef.ignore_empty_defaults = false;

% check for the number of outputs
if nargout == 2
  pvdef.strict = false;
else
  pvdef.strict = true;
end

pvargs = pvdef;

% Recursively process the propval optional arguments (possible
% because we only recurse if optional parameters are given)
if ~isempty(varargin)
  pvargs = propval(varargin, pvdef);
end

% NOTE: Backwards compatibility with previous version of propval
if pvargs.ignore_missing_default | pvargs.ignore_empty_defaults
  pvargs.strict = false;
end

% check for a single cell argument; assume propvals is that argument
if iscell(propvals) && numel(propvals) == 1
  propvals = propvals{1};
end

% check for valid inputs
if ~iscell(propvals) & ~isstruct(propvals)
  error('Property-value pairs must be a cell array or a structure.');
end

if ~isstruct(defaults) & ~isempty(defaults)
  error('Defaults struct must be a structure.');
end

% check for empty defaults structure
if isempty(defaults)
  if pvargs.strict & ~pvargs.ignore_missing_default
   error('Empty defaults structure passed to propval.');
  end
  defaults = struct();
end

defaultnames = fieldnames(defaults);
defaultvalues = struct2cell(defaults);

% prepare the defaults structure, but also prepare casechecking
% structure with all case stripped
defaults = struct();
casecheck = struct();

for i = 1:numel(defaultnames)
  defaults.(defaultnames{i}) = defaultvalues{i};
  casecheck.(lower(defaultnames{i})) = defaultvalues{i};
end

% merged starts with the default values
merged = defaults;
unused = {};
used = struct();

properties = [];
values = [];

% To extract property value pairs, we use different methods
% depending on how they were passed in
if isstruct(propvals)
  properties = fieldnames(propvals);
  values = struct2cell(propvals);
else
  properties = { propvals{1:2:end} };
  values = { propvals{2:2:end} };
end

if numel(properties) ~= numel(values)
  error(sprintf('Found %g properties but only %g values.', numel(properties), ...
                numel(values)));
end

% merge new properties with defaults
for i = 1:numel(properties)

  if ~ischar(properties{i})
    error(sprintf('Property %g is not a string.', i));
  end

  % convert property names to lower case
  properties{i} = properties{i};

  % check for multiple usage
  if isfield(used, properties{i})
    error(sprintf('Property %s is defined more than once.\n', ...
                  properties{i}));
  end

  % Check for case errors
  if isfield(casecheck, lower(properties{i})) & ...
    ~isfield(merged, properties{i})
    error(['Property ''%s'' is equal to a default property except ' ...
           'for case.'], properties{i});
  end

  % Merge with defaults
  if isfield(merged, properties{i})
    merged.(properties{i}) = values{i};
  else
    % add to unused property value pairs
    unused{end+1} = properties{i};
    unused{end+1} = values{i};

    % add to defaults, just in case, if the user isn't picking up "unused"
    if (nargout == 1 & ~pvargs.strict)
      merged.(properties{i}) = values{i};
    end

    if pvargs.strict
      error('Property ''%s'' has no default value.', properties{i});
    end

  end

  % mark as used
  used.(properties{i}) = true;
end

author: jiajunwu


text:
function w = composite_weights ( order )

%*****************************************************************************80
%
%% COMPOSITE_WEIGHTS computes weights for a composite midpoint rule.
%
%  Discussion:
%
%     Right now, our composite rule is assumed to be a composition of
%     1 point rules, applied to the interval [-1,1], with equal spacing.
%
%     Except for the rule of order 1, the endpoints are included
%     as abscissas.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 August 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer ORDER, the order of the rule.
%
%    Output, real W(ORDER), the weights of the rule.
%
  a = -1.0;
  b = +1.0;

  if ( order == 1 )
    w(1) = b - a;
  else
    w(1)         = ( b - a ) / ( order - 1 ) / 2.0;
    w(2:order-1) = ( b - a ) / ( order - 1 );
    w(order)     = ( b - a ) / ( order - 1 ) / 2.0;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
classdef RMMEDA_F4 < PROBLEM
% <multi> <real> <large/none>
% Benchmark MOP for testing RM-MEDA

%------------------------------- Reference --------------------------------
% Q. Zhang, A. Zhou, and Y. Jin, RM-MEDA: A regularity model-based
% multiobjective estimation of distribution algorithm, IEEE Transactions on
% Evolutionary Computation, 2008, 12(1): 41-63.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        %% Default settings of the problem
        function Setting(obj)
            obj.M = 3;
            if isempty(obj.D); obj.D = 30; end
            obj.lower    = zeros(1,obj.D);
            obj.upper    = ones(1,obj.D);
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,X)
            g = sum((X(:,3:end)-repmat(X(:,1),1,size(X,2)-2)).^2,2);
            PopObj(:,1) = cos(pi/2*X(:,1)).*cos(pi/2*X(:,2)).*(1+g);
            PopObj(:,2) = cos(pi/2*X(:,1)).*sin(pi/2*X(:,2)).*(1+g);
            PopObj(:,3) = sin(pi/2*X(:,1)).*(1+g);
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(obj,N)
            R = UniformPoint(N,3);
            R = R./repmat(sqrt(sum(R.^2,2)),1,3);
        end
        %% Generate the image of Pareto front
        function R = GetPF(obj)
            a = linspace(0,pi/2,10)';
            R = {sin(a)*cos(a'),sin(a)*sin(a'),cos(a)*ones(size(a'))};
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function T = spm_type(x, arg)
% Translate data type specifiers between SPM & MATLAB representations
% FORMAT T = spm_type(x, arg)
% x    - specifier
% T    - type
% arg  - optional string argument, can be:
%         - 'maxval'  - return maximum allowed value.
%         - 'minval'  - return minimum allowed value.
%         - 'nanrep'  - return 1 if there is a NaN representation.
%         - 'bits'    - return the number of bits per voxel.
%         - 'intt'    - return 1 if values rounded to nearest integer.
%         - 'conv'    - return conversion function handle.
%__________________________________________________________________________
%
% Format specifiers are based on NIFTI-1.
% If the input is a number then the corresponding MATLAB string is
% returned by default.
% If the input is a string then the appropriate TYPE is returned.
% However, if the optional arg argument is supplied then other
% information will be returned instead.
%
% With no arguments, a list of data types is returned.
%__________________________________________________________________________
% Copyright (C) 1996-2015 Wellcome Trust Centre for Neuroimaging

% John Ashburner & Andrew Holmes
% $Id: spm_type.m 6656 2015-12-24 16:49:52Z guillaume $

prec   = {'uint8','int16','int32','float32','float64','int8','uint16','uint32'};
conv   = {@uint8,@int16,@int32,@single,@double,@int8,@uint16,@uint32};
types  = [    2      4      8   16   64   256    512    768];
maxval = [2^8-1 2^15-1 2^31-1  Inf  Inf 2^7-1 2^16-1 2^32-1];
minval = [    0  -2^15  -2^31 -Inf -Inf  -2^7      0      0];
nanrep = [    0      0      0    1    1     0      0      0];
bits   = [    8     16     32   32   64     8     16     32];
intt   = [    1      1      1    0    0     1      1      1];

if ~nargin
    T = types;
    return;
end

if ischar(x)
    sel = find(strcmpi(prec,deblank(x)));
else
    if numel(x) == 1
        sel = find(types == x);
    else
        sel = find(ismember(types,x));
    end
end
if nargin == 1
    if ischar(x)
        if isempty(sel), T = NaN;
        else T = types(sel); end
    else
        if isempty(sel), T = 'unknown';
        else T = char(prec(sel)); end
    end
elseif isempty(sel)
    T = NaN;
else
    switch lower(arg)
        case 'maxval', T = maxval(sel);
        case 'minval', T = minval(sel);
        case 'nanrep', T = nanrep(sel);
        case 'bits',   T = bits(sel);
        case 'intt',   T = intt(sel);
        case 'conv',   T = conv(sel);
        otherwise,     T = NaN;
    end
end

author: spm


text:
classdef SubcellsMesher_1D < SubcellsMesher_Interior_Abstract

    methods (Access = public)

        function obj = SubcellsMesher_1D(cParams)
           obj.init(cParams);
        end

    end

    methods (Access = protected)

        function computeAllPossibleSubcellsInCell(obj)
            obj.subcells_connec = [1 3; 2 3];
        end

    end
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
%
% Overview: Test ECGBeatFitter program for mean ECG beat extraction and
% parameter optimization. The script does the following:
% 1. Reads the ECG data in the sample_ecg_datam file in the 'Physionet-Cardiovascular-Signal-Toolbox/Tools/ECG_Analysis_Tools/MV/testdata/' subfolder.
% 2. Estimates the average ECG beat for the first fifteen seconds.
% 3. The ECGBeatFitter algorithm is used to estimate the parameters of
% Gaussian functions that accurately estimate the shape of the average ECG
% beat. The number of Gaussians is specified by selecting points along the
% average beat in the ECGBeatFitter GUI.
% 4. The output of the script are the parameters for the Gaussians stored in the following variables:
%
%   ai: contains the amplitudes of the Gaussians.
%   bi: contains the standard deviations of the Gaussians.
%   tetai: contains the phase of each Gaussian function.
%
% 5. These parameters may be estimated for each of the x, y and z
% components of a VCG and substituted into the generate_resp_modulated_ecg
% function in the Physionet-Cardiovascular-Signal-Toolbox to generate an artificial
% VCG with a morphology similar to VCG used for estimating the parameters.
% The Dower transform may be applied to generate 12 lead ECG from the VCG.
%
%
% ORIGINAL SOURCE AND AUTHORS:
% Open Source ECG Toolbox, version 1.0, November 2006
% Released under the GNU General Public License
% Copyright (C) 2006  Reza Sameni
% Sharif University of Technology, Tehran, Iran -- LIS-INPG, Grenoble, France
% reza.sameni@gmail.com
% editted by Ismail Sadiq on 11/29/2020.

% This program is free software; you can redistribute it and/or modify it
% under the terms of the GNU General Public License as published by the
% Free Software Foundation; either version 2 of the License, or (at your
% option) any later version.
% This program is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
% Public License for more details. You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.

clc
clear
close all;

% Load test ecg data
addpath(genpath('../../../../../PhysioNet-Cardiovascular-Signal-Toolbox-master/')); % add all dependencies from cardiovascular signal toolbox
sig = load('./testdata/mvm/sample_ecg_datam');
siginfo = readheader('./testdata/mvm/sample_ecg_datam.hea');
fs = siginfo.freq;
ecg = (sig.val - siginfo.adczero)./siginfo.gain; data=ecg(1:15000);% Adjust signal according to gain and dc offset
clear siginfo sig ecg;

t = (0:length(data)-1)/fs;

f = 1;                                          % approximate R-peak frequency

bsline = LPFilter(data,.7/fs);                  % baseline wander removal (may be replaced by other approaches)
%bsline = BaseLineKF(data,.5/fs);                % baseline wander removal (may be replaced by other approaches)

data1 = data-bsline;

%//////////////////////////////////////////////////////////////////////////
% Making the data noisy
SNR = 20;
SignalPower = mean(data1.^2);
NoisePower = SignalPower / 10^(SNR/10);
x = data1 + sqrt(NoisePower)*randn(size(data1));
%//////////////////////////////////////////////////////////////////////////

peaks = PeakDetection(x,f/fs);                  % peak detection

[phase, phasepos] = PhaseCalculation(peaks);     % phase calculation

teta = 0;                                       % desired phase shift
pphase = PhaseShifting(phase,teta);             % phase shifting

bins = round(fs/3);                                     % number of phase bins
[ECGmean,ECGsd,meanphase] = MeanECGExtraction(x,pphase,bins,1); % mean ECG extraction

OptimalParams = ECGBeatFitter(ECGmean,ECGsd,meanphase);               % ECG beat fitter GUI

% display the optimal parameters
L = length(OptimalParams)/3;
ai = OptimalParams(1:L)
bi = OptimalParams(L+1:2*L)
tetai = OptimalParams(2*L+1:3*L)

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox


text:
clear
cleanFile = 'DataSet/TIMIT_10_TEST-White+5db.wav';
enhancedFile_1 = 'Enhanced/out_TIMIT_10_TEST_IRMDNN_0308_10.wav';
[snr_mean, segsnr_mean]= comp_snr(cleanFile, enhancedFile_1);

author: vipchengrui

repo: traditional-speech-enhancement


text:
function npartitions = npart_enum ( n, npart )

%*****************************************************************************80
%
%% NPART_ENUM enumerates the number of partitions of N with NPART parts.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 January 2011
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Donald Kreher, Douglas Simpson,
%    Combinatorial Algorithms,
%    CRC Press, 1998,
%    ISBN: 0-8493-3988-X,
%    LC: QA164.K73.
%
%  Parameters:
%
%    Input, integer N, the integer to be partitioned.
%    Normally N must be positive, but for this routine any
%    N is allowed.
%
%    Input, integer NPART, the number of parts of the partition.
%    Normally, 1 <= NPART <= N is required,
%    but for this routine any value of NPART is allowed.
%
%    Output, integer NPARTITIONS is the number of partitions of N
%    with NPART parts.
%
  if ( n <= 0 )

    npartitions = 0;

  elseif ( npart <= 0 || n < npart )

    npartitions = 0;

  else

    offset = 1;

    p = npart_table ( n, npart );

    npartitions = p(n+offset,npart+offset);

  end

  return
end


text:
% created: Zoya Bylinskii, Aug 2014
% Code for computing degrees of visual angle

% This code has been written to make it easier to report visual angle in
% papers, and to ensure a standardized calculation is used.

% This code can take variable input (e.g. horizontal, vertical, or diagonal
% screen measurements in any units), and will calculate any number of
% degrees of visual angle in pixels, as well as degrees of visual angle
% spanned by the computer screen or an image.

% The code is provided on: http://saliency.mit.edu/

% Resources used: http://osdoc.cogsci.nl/miscellaneous/visual-angle/,
%                 http://en.wikipedia.org/wiki/Visual_angle
% Thank you to Lavanya Sharan and Melissa Vo for helpful discussion.

function computeDegVisAngle_gui()

% ------ INPUT ARGUMENTS: ------
% necessary:
% 'distance'         distance to screen with units (e.g. 56 cm)

% provide at least one of:
% 'height'           height of screen with units
% 'width'            width of screen with units
% 'diagonal'         diagonal of screen with units

% provide at least one of:
% 'hres'             horizontal resolution of screen (provide this if also provided 'height')
% 'vres'             vertical resolution of screen (provide this if also provided 'width')
% provide both if only provided 'diagonal'

% optional:
%'outputUnit'        output units to measure screen, distance in (default: cm)
%'numDegs'           number of degres of visual angle (in pixels) to calculate (default numDegs=1)
%'imageWidth'        width of image in pixels
%'imageHeight'       height of image in pixels

% NOTE: if you provide the screen diagonal and resolution, then the screen
% height and width will be estimated under the assumption that the aspect
% ratio of the screen dimensions matches that of the set resolution.
% If these estimates look wrong, either:
% (1) you measured the monitor dimensions, instead of the screen dimensions
% (2) the screen resolution aspect ratio does not match the screen
%     dimensions
% (3) explicitely enter the screen height and width instead of the diagonal

% Common mistake: remember that the number of pixels per unit size of screen
% or image is a linear function, but number of pixels per degree of visual
% angle is NOT.

% Remember: when writing a paper, include enough information for visual
% angle to be calculated - providing screen resolution without screen
% dimensions or vice versa is insufficient! Also, remember to measure the
% screen itself, rather than the whole monitor.

% units required for distance, height, width, and diagonal values
ui=inputdlg({'distance (with units)','height (with units)',...
    'width (with units)','diagonal (with units)',...
    'horizontal resolution','vertical resolution',...
    'output units (optional)','number of degrees (optional)',...
    'image width in pixels (optional)','image height in pixels (optional)'},...
    'Visual Angle Calculator',1,...
    {'','','','','','','cm','1','',''});

dist = ui{1};
height = ui{2};
width = ui{3};
diagonal = ui{4};
hres = str2num(ui{5});
vres = str2num(ui{6});
outputUnit = ui{7};
numDegs = str2num(ui{8});
imageWidth = str2num(ui{9});
imageHeight = str2num(ui{10});

computeDegVisAngle_helper(dist,height,width,diagonal,hres,vres,outputUnit,numDegs,imageWidth,imageHeight);

author: cvzoya

repo: saliency

save_path: github-repos/MATLAB/cvzoya-saliency


text:
function str = mean2str(m,s,n,varargin)

%mean2str - Convert mean or median and SEM (or confidence interval) to string.
%
%  USAGE
%
%    s = mean2str(m,s,n,<options>)
%
%    m              mean or median
%    s              standard error (see sem or semedian) or confidence interval
%    n              optional number of observations
%    <options>      optional list of property-value pairs (see table below)
%
%    =========================================================================
%     Properties    Values
%    -------------------------------------------------------------------------
%     'precision'   number of digits (default = 3)
%     'split'       split the output in two cells, e.g. for small figures
%                   (default = 'off')
%    =========================================================================

% Copyright (C) 2013 by Michaël Zugaro
%
% 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.

% Default values
precision = 3;
split = 'off';

% Check parameters
if nargin < 2,
  error('Incorrect number of parameters (type ''help mean2str'' for details).');
end
if ~isdscalar(m),
  error('Incorrect mean or median (type ''help mean2str'' for details).');
end
if ~isdscalar(s) && ~isdvector(s,'<=','#2'),
  error('Incorrect SEM or confidence interval (type ''help mean2str'' for details).');
end

% Optional number of observations
if nargin < 3,
	n = [];
elseif ischar(n),
	varargin = {n,varargin{:}};
	n = [];
elseif ~isiscalar(n,'>0'),
  error('Incorrect number of observations (type ''help mean2str'' for details).');
end

% Parse parameter list
for i = 1:2:length(varargin),
  if ~ischar(varargin{i}),
    error(['Parameter ' num2str(i+2) ' is not a property (type ''help mean2str'' for details).']);
  end
  switch(lower(varargin{i})),
    case 'precision',
      precision = varargin{i+1};
      if ~isiscalar(precision,'>0'),
        error('Incorrect value for property ''precision'' (type ''help mean2str'' for details).');
      end
    case 'split',
      split = lower(varargin{i+1});
      if ~isstring_FMAT(split,'on','off'),
        error('Incorrect value for property ''split'' (type ''help mean2str'' for details).');
      end
    otherwise,
      error(['Unknown property ''' num2str(varargin{i}) ''' (type ''help mean2str'' for details).']);
  end
end

format = ['%.' int2str(precision) 'f'];
if isdscalar(s),
	str = [sprintf(format,m) ' +- '  sprintf(format,s)];
else
	str = [sprintf(format,m) ' [' sprintf(format,s(1)) ',' sprintf(format,s(2)) ']'];
end
if ~isempty(n),
	if strcmp(split,'on'),
		str = {str,['(N=' int2str(n) ')']};
	else
		str = [str ' (N=' int2str(n) ')'];
	end
end

author: buzsakilab


text:
function perform_frame_selection(params)
  %% Load images
  clc;
  sequence_name = params.sequences_name(params.current_sequence).name;
  disp(['Loading images from sequence: ' sequence_name]);
  V = load_images(params); % show_4dvideo(V);
  disp(size(V));

  %%% Frame difference
  disp('Performing frame difference');
  O = zeros(size(V,1),size(V,2),size(V,4));
  vDiffSum = zeros(size(V,4),1);
  myfilter = fspecial('gaussian',[3 3], 0.5);
  I0 = rgb2gray(V(:,:,:,1));
  switch(sequence_name)
    case 'Toscana'
      beta = 1e-2;
    otherwise
      beta = 1e-3;
  end
  for i = 2 : size(V,4)
    if(mod(i,10) == 0) fprintf('%d ',i); end

    I = rgb2gray(V(:,:,:,i));
    I = imfilter(I, myfilter, 'replicate');
    %I = medfilt2(I,[5,5]);
    %I = imgaussfilt(I,2);

    switch(sequence_name)
      case {'Candela_m1.10','HallAndMonitor'}
        if(i == 2) disp('Using first image as reference...'); end
        I_old = I0;
      otherwise
        if(i == 2) disp('Using the previous frame as reference...'); end
        I_old = rgb2gray(V(:,:,:,i-1));
    end
    I_old = imfilter(I_old, myfilter, 'replicate');

    I_mag = sqrt(power(I-I_old,2));
    I_res = I_mag;

    I_thr = 0.5*I_res.^2 > beta;
    O(:,:,i) = I_thr;
    vDiffSum(i) = sum(I_thr(:));
  end
  disp(i);
  vDiffSum = vDiffSum./max(vDiffSum(:));
  % show_3dvideo(O);
  clear i I I0 I_old I_mag I_res I_thr myfilter beta;

  %%% Average
%   switch(sequence_name)
%     case 'HallAndMonitor'
%       avgDiffSum = medfilt1(vDiffSum,5,'zeropad');
%     otherwise
       avgDiffSum = medfilt1(vDiffSum,3,'zeropad');
%   end
  if(params.debug)
    clf,plot(vDiffSum),hold on,plot(avgDiffSum),hold off;
    pause(1);
  end

  %% Gradient
  disp('Calculating difference gradient');

  %S = vDiffSum;
  S = avgDiffSum;

  Y = gradient(S);
  %Y = diff(S);
  Y = (Y-min(Y))/(max(Y)-min(Y)); % normalize
  %Y = medfilt1(Y,3,'truncate');
  %Y = medfilt1(Y,3,'zeropad');

  switch(sequence_name)
    case {'HallAndMonitor','HighwayII'}
      Tt = 0.075;
    case {'CAVIAR1','Candela_m1.10','HighwayI','IBMtest2','PeopleAndFoliage'}
      Tt = 0.1;
    case {'Foliage'}
      Tt = 0.2;
    case {'HumanBody2'}
      Tt = 0.05;
    case {'Board','CAVIAR2','CaVignal','Snellen','Toscana'}
      Tt = 0.125;
  end

  T1 = (Y > (mean(Y) + Tt));
  T2 = (Y < (mean(Y) - Tt));
  T = double(or(T1,T2));

  if(params.debug)
    clf;
    plot(S);
    hold on;
    plot(Y);
    plot(T.*.05,'s');
    hold off;
    title('Frame Selection');
    xlabel('Frames'),ylabel({'Difference between';'consecutive frames'});
    legend('vector d normalized ','derivative of vector d',...
      'selected frames',...
      'Location','southoutside',... %northoutside southoutside
      'Orientation','horizontal');
    pause(1);
  end
  disp(['Number of selected frames: ' num2str(sum(T))]);

  %% Frame selection
  disp('Performing frame selection');
  switch(sequence_name)
    case {'Toscana'}
      V2 = V;
      O2 = O;
    otherwise
      V2 = V(:,:,:,logical(T));
      O2 = O(:,:,logical(T));
  end
  % size(Vnew), size(Onew)

  if(params.debug)
    clf,show_4dvideo(V2);
    %clf,show_3dvideo(O2);slice3(O2);
  end

  %% Save
  clear T Tt S T1 T2 Y avgDiffSum;
  disp(['Saving at: ' fullfile(params.sequences_mat,[sequence_name '.mat'])]);
  save(fullfile(params.sequences_mat,[sequence_name '.mat']));
end

author: andrewssobral


text:
function s=AAFT(x,c)
%Syntax: s=AAFT(x,c)
%___________________
%
% Makes c Amplitude Adjusted Fourier Transformed (AAFT) surrogates of a time
% series x.
%
% s is the AAFT time series.
% x is the original time series.
% c is the number of surrogates.
%
%
% References:
%
% Theiler J, Galdrikian B, Longtin A, Eubank S, Farmer D J (1992): Using
% Surrogate Data to Detect Nonlinearity in Time Series. In Nonlinear Modeling
% and Forecasting, eds. Casdagli M & Eubank S. 163-188. Addison-Wesley
%
% Theiler J, Eubank S,Galdrikian B, Longtin A,  Farmer D J (1992): Testing
% for nonlinearity in time series: the method of surrogate data. Physica D
% 58: 77-94
%
%
% Alexandros Leontitsis
% Department of Education
% University of Ioannina
% 45110 - Dourouti
% Ioannina
% Greece
%
% University e-mail: me00743@cc.uoi.gr
% Lifetime e-mail: leoaleq@yahoo.com
% Homepage: http://www.geocities.com/CapeCanaveral/Lab/1421
%
% 12 Apr 2001.

if nargin<1 | isempty(x)==1
   error('You should provide a time series.');
else
   % x must be a vector
   if min(size(x))>1
      error('Invalid time series.');
   end
   x=x(:);
end

if nargin<2 | isempty(c)==1
   c=1;
else
   % c must be scalar
   if sum(size(c))>2
      error('c must be scalar.');
   end
   % c must be greater or equal than 1
   if c<1
      error('c must be greater or equal than 1.');
   end
end

for i=1:c
    % Initialize
    y=x;
    % Make n normal random devaiates
    normal=sort(randn(size(y)));
    % Sort y and extract the ranks
    [y,T]=sort(y);
    [T,r]=sort(T);
    % Assign the ranks of y to the normal deviates and apply the phase
    %  randomization
    normal=phaseran(normal(r));
    % Extract the ranks of the phase randomized normal deviates
    [normal,T]=sort(normal);
    [T,r]=sort(T);
    % Assign the ranks of the phase randomized normal deviates to y and
    %  obtain the AAFT surrogates
    s(:,i)=y(r);
end

author: Sable

repo: mcbench-benchmarks


text:
function [ alphas, betas, scaling, finalLikelihood] = CCRF_training_bfgs( num_seqs, thresholdX, thresholdFun, x, y, yUnnormed, alphas, betas, lambda_a, lambda_b, similarityFNs, Precalc_Bs, Precalc_Bs_flat, Precalc_yBys, varargin)
%GRADIENTDESCENTCCRF Performs CCRF gradient descen given the initial state
%and gradient descent parameters
%   Detailed explanation goes here

    % if these are not provided calculate them, TODO this might be

    % It is possible to predefine the component B^(k) required
    % to compute B term and partial derivatives, also can predefine yB^(k)y,
    % as they also do not change through the iterations
    if(sum(strcmp(varargin,'PrecalcBs')) && sum(strcmp(varargin,'PrecalcBsFlat'))...
             && sum(strcmp(varargin,'Precalc_yBy')))

        ind = find(strcmp(varargin,'PrecalcBs')) + 1;
        Precalc_Bs = varargin{ind};

        ind = find(strcmp(varargin,'PrecalcBsFlat')) + 1;
        Precalc_Bs_flat = varargin{ind};

        ind = find(strcmp(varargin,'Precalc_yBys')) + 1;
        Precalc_yBys = varargin{ind};
    else
        % if these are not provided calculate them
        [ ~, Precalc_Bs, Precalc_Bs_flat, Precalc_yBys ] = CalculateSimilarities( num_seqs, x, similarityFNs, y);
    end

    params = [alphas; betas];

    objectiveFun = @(params)objectiveFunction(params, numel(alphas), lambda_a, lambda_b, Precalc_Bs, x, y, Precalc_yBys, Precalc_Bs_flat);

    options = optimset('Algorithm','interior-point','GradObj','on', 'TolX', thresholdX, 'TolFun', thresholdFun, 'Hessian', 'bfgs', 'display','off', 'useParallel', 'Always');

    if(sum(strcmp(varargin,'max_iter')))
        options.MaxIter = varargin{find(strcmp(varargin,'max_iter')) + 1};
    end

    params = fmincon(objectiveFun, params, [], [],[],[], zeros(numel(params),1), Inf(numel(params), 1), [], options);
    alphas = params(1:numel(alphas));
    betas = params(numel(alphas)+1:end);

    finalLikelihood = LogLikelihoodCCRF(y, x, alphas, betas, lambda_a, lambda_b, Precalc_Bs_flat);
%     fprintf('Final log likelihood at iteration; logL %f, learning rate\n', finalLikelihood);

    % establish the scaling
    scaling = getScaling2(alphas, betas, x, yUnnormed, Precalc_Bs);

end

function [loss, gradient] = objectiveFunction(params, numAlpha, lambda_a, lambda_b, PrecalcBs, x, y, Precalc_yBys, PrecalcBsFlat)

    alphas = params(1:numAlpha);
    betas = params(numAlpha+1:end);
    [gradient, SigmaInvs, CholDecomps, Sigmas] = gradientCCRFFull(params, lambda_a, lambda_b, PrecalcBs, x, y, Precalc_yBys, PrecalcBsFlat);
    % as bfgs does gradient descent rather than ascent, negate the results
    gradient = -gradient;
    loss = -LogLikelihoodCCRF(y, x, alphas, betas, lambda_a, lambda_b, PrecalcBsFlat, SigmaInvs, CholDecomps, Sigmas);
end

author: TadasBaltrusaitis

repo: OpenFace


text:
function [E,J]=SynthMeasWatsonSHStickTortIsoVIsoDot_B0(x, protocol, fibredir)
% Substrate: Impermeable sticks (cylinders with zero radius) in a homogeneous
% background.
% Orientation distribution: Watson's distribution with SH approximation
% Signal approximation: Not applicable
% Notes: This version estimates the hindered diffusivity from the free diffusivity
% and packing density using Szafer et al's tortuosity model for randomly
% packed cylinders.
% This version includes an isotropic diffusion compartment with its own
% diffusivity.
% This version includes a stationary water compartment.
% Includes a free parameter for the measurement at b=0.
%
% [E,J]=SynthMeasWatsonSHStickTortIsoV_B0(x, protocol, fibredir)
% returns the measurements E according to the model and the Jacobian J of the
% measurements with respect to the parameters.  The Jacobian does not
% include derivates with respect to the fibre direction.
%
% x is the list of model parameters in SI units:
% x(1) is the volume fraction of the intracellular space.
% x(2) is the free diffusivity of the material inside and outside the cylinders.
% x(3) is the concentration parameter of the Watson's distribution.
% x(4) is the volume fraction of the isotropic compartment.
% x(5) is the diffusivity of the isotropic compartment.
% x(6) is the volume fraction of the isotropic restriction.
% x(7) is the measurement at b=0.
%
% protocol is the object containing the acquisition protocol.
%
% fibredir is a unit vector along the symmetry axis of the Watson's
% distribution.  It must be in Cartesian coordinates [x y z]' with size [3 1].
%
% author: Gary Hui Zhang (gary.zhang@ucl.ac.uk)
%

xcyl=[x(1) x(2) 0 x(3) x(4) x(5) x(6) x(7)];

if(nargout == 1)
    E=SynthMeasWatsonSHCylSingleRadTortIsoVIsoDot_GPD_B0(xcyl, protocol, fibredir, 0);
else
    [E,Jcyl]=SynthMeasWatsonSHCylSingleRadTortIsoVIsoDot_GPD_B0(xcyl, protocol, fibredir, 0);
end

if(nargout>1)
    J(:,1) = Jcyl(:,1);
    J(:,2) = Jcyl(:,2);
    J(:,3) = Jcyl(:,4);
    J(:,4) = Jcyl(:,5);
    J(:,5) = Jcyl(:,6);
    J(:,6) = Jcyl(:,7);
    J(:,7) = Jcyl(:,8);
end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
classdef DehomogenisationPrinter < handle

    properties (Access = private)
        mesh
        dataRes
        alpha
        density
    end

    properties (Access = private)
        fileName
        folderPath
    end

    methods (Access = public)

        function obj = DehomogenisationPrinter(cParams)
            obj.init(cParams);
        end

        function print(obj)
            obj.wrapResMeshData();
            obj.projectAlphaToNodes();
            obj.projectDensityToNodes();
            obj.printDehomogResFile();
            obj.printDehomogMshFile();
        end

    end

    methods (Access = private)

        function init(obj,cParams)
            obj.fileName   = cParams.fileName;
            obj.folderPath = cParams.folderPath;
        end

        function wrapResMeshData(obj)
            s.fileName   = obj.fileName;
            s.folderPath = obj.folderPath;
            w = WrapperMshResFiles(s);
            w.compute();
            obj.dataRes = w.dataRes;
            obj.mesh    = w.mesh;
        end

        function projectAlphaToNodes(obj)
            alphaN = obj.dataRes.AlphaGauss;
            alpha1 = obj.projectInNodesScalarPieceWiseFunction(alphaN(:,1));
            alpha2 = obj.projectInNodesScalarPieceWiseFunction(alphaN(:,2));
            normA = sqrt(alpha1.^2 + alpha2.^2);
            obj.alpha = [alpha1,alpha2]./normA;
        end

        function projectDensityToNodes(obj)
            alphaN = obj.dataRes.DensityGauss;
            dens   = obj.projectInNodesScalarPieceWiseFunction(alphaN(:,1));
            obj.density = dens;
        end

        function alpha = projectInNodesScalarPieceWiseFunction(obj,fValues)
            s.mesh   = obj.mesh;
            s.fValues = fValues;
            f = P0Function(s);
            alpha = f.projectToLinearNodalFunction();
        end

        function printDehomogResFile(obj)
            fD = [obj.fileName,'DehomogRes.txt'];
            fOutName = fullfile(obj.folderPath,fD);
            s.fileName = fOutName;
            s.values(:,1) = obj.dataRes.DesignVar1;
            s.values(:,2) = obj.dataRes.DesignVar2;
            s.values(:,3) = obj.alpha(:,1);
            s.values(:,4) = obj.alpha(:,2);
            s.values(:,5) = obj.dataRes.SuperEllipseExponent;
            s.values(:,6) = obj.density;
            p = DehomogenisationOutputPrinter(s);
            p.print();
        end

        function printDehomogMshFile(obj)
            fD = [obj.fileName,'DehomogMesh.txt'];
            fOutName = fullfile(obj.folderPath,fD);
            s.fileName = fOutName;
            s.mesh     = obj.mesh;
            p = DehomogenisationMeshPrinter(s);
            p.print();
        end

    end

end


text:
function [scan3dd] = elec_load_scan_3ddasc(filename)

% elec_load_scan_3ddasc - read ascii export of Neuroscan 3DD file
%
% Load an ascii electrode file and return
% the electrode labels and coordinates.  This
% function can read and return Cartesian (x,y,z)
% and/or spherical (theta,phi,r) coordinates.
%
% [scan3dd] = elec_load_scan_3ddasc(filename)
%
% where:
%
% file = '<path><filename>' with row format '%s %d %f %f %f'
%
% The file format is that of NeuroScan 3Dspace ascii
% export files.  Each row of the file comprises an electrode label,
% an electrode type code (see below), and the x,y,z coordinates (cm).
% Each field is separated by spaces.  For example,
%
% Fz    69  Xcm Ycm Zcm
%
% Example result:
%
% scan3dd =
%
%      label: {1x128 cell}
%          x: [128x1 double]
%          y: [128x1 double]
%          z: [128x1 double]
%        hsp: [1617x3 double]
%     nasion: [-0.0333 9.0341 0]
%        lpa: [-6.9128 0 0]
%        rpa: [6.9128 0 0]
%     origin: [0 0 0]
%        ref: [0 9.9341 0]
%
% All coordinates are in centimeters; +X is right, +Y is anterior
% and +Z is superior.
%
% Notes:
%
% i) Type is defined as follows (from Neuroscan 3Dspace ascii export):
%
%           Electrode               Type
%           ---------               ----
%           Nasion                   110 (or 78)
%           Left                     108 (or 76)
%           Right                    114 (or 82)
%           electrodes                69
%           Centroid (origin)         99 (or 67; eg <0,0,0>)
%           Ref                      120 (or 88)
%           Scalp Points (hsp)        32
%
%

% $Revision: 1.1 $ $Date: 2009-04-28 22:13:55 $

% Licence:  GNU GPL, no express or implied warranties
% History:  08/1999, Darren.Weber_at_radiology.ucsf.edu
%           08/2003, Darren.Weber_at_radiology.ucsf.edu
%                    complete rewrite, replacing elec_load.m
%                    new version of 3Dspace exports different type numbers ;-)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ver = '$Revision: 1.1 $';
fprintf('\nELEC_LOAD_SCAN_3DDASC [v %s]\n',ver(11:15));

tic;

[path,name,ext] = fileparts(filename);
file = fullfile(path,[name ext]);

fprintf('...loading electrodes from:\n\t%s\n', file);

scan3dd = read_3dd(file);

fprintf('...loaded %d electrodes\n', size(scan3dd.x,1));

t = toc; fprintf('...done (%6.2f sec).\n\n',t);

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function scan3dd = read_3dd(file),

scan3dd.label = [];
scan3dd.x = [];
scan3dd.y = [];
scan3dd.z = [];
scan3dd.hsp = [];

fid = fopen(file);
if fid < 0,
  msg = sprintf('cannot open file: %s\n',file);
  error(msg);
end

% 3DD ascii files contain position information for
% fiducial, sensor, and head shape fields.

% First get the fiducial points by reading the whole file
% (clumsy, but exhaustive search for field indicators).

% Fiducial points are required for MRI registration. They are
% the nasion, left and right preauricular points, eg:

%      Nasion	78	-0.033340	9.034104	0.000000
%      Left	76	-6.912818	-0.000000	0.000000
%      Right	82	6.912818	0.000000	-0.000000

fprintf('...searching for fiducials...');

n = 0;
while n < 5,
  tmp = fgetl(fid);
  if tmp < 0, break; end

  tmp = lower(tmp);

  if strfind(tmp,'nasion'),
    tmp = sscanf(tmp,'%s %d %f %f %f');
    scan3dd.nasion = [tmp(end-2) tmp(end-1) tmp(end)];
    n = n + 1;
    continue;
  end
  if strfind(tmp,'left'),
    tmp = sscanf(tmp,'%s %d %f %f %f');
    scan3dd.lpa = [tmp(end-2) tmp(end-1) tmp(end)];
    n = n + 1;
    continue;
  end
  if strfind(tmp,'right'),
    tmp = sscanf(tmp,'%s %d %f %f %f');
    scan3dd.rpa = [tmp(end-2) tmp(end-1) tmp(end)];
    n = n + 1;
    continue;
  end
  if strfind(tmp,'centroid'),
    tmp = sscanf(tmp,'%s %d %f %f %f');
    scan3dd.origin = [tmp(end-2) tmp(end-1) tmp(end)];
    n = n + 1;
    continue;
  end
  if strfind(tmp,'ref'),
    tmp = sscanf(tmp,'%s %d %f %f %f');
    scan3dd.ref = [tmp(end-2) tmp(end-1) tmp(end)];
    n = n + 1;
    continue;
  end
end

frewind(fid);
fprintf('done\n');

fprintf('...searching for electrodes...');

ok = 1;
while ok,
  tmp = fgetl(fid);
  if tmp < 0, break; end

  if strfind(lower(tmp),'nasion'),   continue; end
  if strfind(lower(tmp),'left'),     continue; end
  if strfind(lower(tmp),'right'),    continue; end
  if strfind(lower(tmp),'centroid'), continue; end
  if strfind(lower(tmp),'ref'),      continue; end

  tmp = sscanf(tmp,'%s %d %f %f %f');

  if tmp(end-3) == 69,
    scan3dd.label{end+1} = char(tmp(1:end-4))';
    scan3dd.x(end+1,1) = tmp(end-2);
    scan3dd.y(end+1,1) = tmp(end-1);
    scan3dd.z(end+1,1) = tmp(end);
    continue;
  end
  if strfind(char(tmp(1:end-4))','32'),
    break;
    % found a head shape point
    scan3dd.hsp(end+1,:) = [tmp(end-2),tmp(end-1),tmp(end)];
    continue;
  end
end

frewind(fid);
fprintf('done\n');

fprintf('...searching for head shape points...');

ok = 1;
while ok,
  tmp = fgetl(fid);
  if tmp < 0, break; end

  if strfind(lower(tmp),'nasion'),   continue; end
  if strfind(lower(tmp),'left'),     continue; end
  if strfind(lower(tmp),'right'),    continue; end
  if strfind(lower(tmp),'centroid'), continue; end
  if strfind(lower(tmp),'ref'),      continue; end
  if strfind(lower(tmp),'69'),       continue; end

  tmp = sscanf(tmp,'%d %f %f %f');

  if tmp(end-3) == 32,
    % found a head shape point
    scan3dd.hsp(end+1,:) = [tmp(end-2),tmp(end-1),tmp(end)];
    continue;
  end
end

frewind(fid);
fprintf('done\n');

fclose(fid);

return

author: PatternRecognition

repo: OpenBMI


text:
function results = vl_test_binsearch(varargin)
% VL_TEST_BINSEARCH
vl_test_init ;

function test_inf_bins()
x = [-inf -1 0 1 +inf] ;
vl_assert_equal(vl_binsearch([],          x), [0 0 0 0 0]) ;
vl_assert_equal(vl_binsearch([-inf 0],    x), [1 1 2 2 2]) ;
vl_assert_equal(vl_binsearch([-inf],      x), [1 1 1 1 1]) ;
vl_assert_equal(vl_binsearch([-inf +inf], x), [1 1 1 1 2]) ;

function test_empty()
vl_assert_equal(vl_binsearch([], []), []) ;

function test_bnd()
vl_assert_equal(vl_binsearch([], [1]),    [0]) ;
vl_assert_equal(vl_binsearch([], [-inf]), [0]) ;
vl_assert_equal(vl_binsearch([], [+inf]), [0]) ;

vl_assert_equal(vl_binsearch([1], [.9]),   [0]) ;
vl_assert_equal(vl_binsearch([1], [1]),    [1]) ;
vl_assert_equal(vl_binsearch([1], [-inf]), [0]) ;
vl_assert_equal(vl_binsearch([1], [+inf]), [1]) ;

function test_basic()
vl_assert_equal(vl_binsearch(-10:10, -10:10), 1:21) ;
vl_assert_equal(vl_binsearch(-10:10, -11:10), 0:21) ;
vl_assert_equal(vl_binsearch(-10:10, [-inf, -11:10, +inf]), [0 0:21 21]) ;

function test_frac()
vl_assert_equal(vl_binsearch(1:10, 1:.5:10), floor(1:.5:10))
vl_assert_equal(vl_binsearch(1:10, fliplr(1:.5:10)), ...
                fliplr(floor(1:.5:10))) ;

function test_array()
a = reshape(1:100,10,10) ;
b = reshape(1:.5:100.5, 2, []) ;
c = floor(b) ;
vl_assert_equal(vl_binsearch(a,b), c) ;

author: Cloud-CV

repo: object-proposals


text:
function csv = readCSV(filename)
FID = fopen(filename);
currLine = fgetl(FID);
idx = find(currLine==',');
field{1} = currLine(1:idx(1)-1);
csv.(field{1}) = {};
for i=2:length(idx)
    field{end+1} = regexprep(currLine(idx(i-1)+1:idx(i)-1), ' ', '_');
    csv.(field{end}) = {};
end
field{end+1} = regexprep(currLine(idx(end)+1:end), ' ', '_');
csv.(field{end}) = {};
nField = length(field);

line_cnt = 0;
while 1
    currLine = fgetl(FID);
    line_cnt = line_cnt + 1;
    PrintProgress(line_cnt, 1e5, 5000);

    if currLine == -1
        break;
    end

    idx = find(currLine==',');
    for i=1:nField
        if i==1
            idx1 = 1;
        else
            idx1 = idx(i-1)+1;
        end
        if i==nField
            idx2 = length(currLine);
        else
            idx2 = idx(i)-1;
        end

        csv.(field{i}){end+1} = currLine(idx1:idx2);
    end
end

fclose(FID);

author: singaxiong

repo: SignalGraph


text:
%% Author: epokh
%% Website: www.epokh.org/drupy
%% This software is under GPL

%%Function description:
%%Transformation of the end effector frame following the Euler angles
%%in closed form
%%input: rotation angles along y,z,y axis expressed in degrees

function [EulerMat]=Euler(fi1,theta,fi2)

EulerMat=[cosd(fi1)*cosd(theta)*cosd(fi2)-sind(fi1)*sind(fi2),-cosd(fi1)*cosd(theta)*sind(fi2)-sind(fi1)*cosd(fi2),cosd(fi1)*sind(theta),0;
        sind(fi1)*cosd(theta)*cosd(fi2)+cosd(fi1)*sind(fi2),-sind(fi1)*cosd(theta)*sind(fi2)+cosd(fi1)*cosd(fi2),sind(fi1)*sind(theta),0;
        -sind(theta)*cosd(fi2),sind(theta)*sind(fi2),cosd(theta),0;
        0                      ,0                   ,0          ,1;];

end