text:
function report(dir1, suffix1, showcls, do_auc_ap)
% Print scores for all classes.

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

conf = voc_config();

if nargin < 3
  showcls = true;
end

if nargin < 4
  do_auc_ap = false;
end

count1 = 0;
for i=1:length(conf.pascal.VOCopts.classes)
  cls = conf.pascal.VOCopts.classes{i};
  try
    load([dir1 cls suffix1]);
    ap1 = ap;
    if do_auc_ap
      ap2 = xVOCap(recall, prec);
      score2(i) = ap2;
    end
    if showcls
      if do_auc_ap
        fprintf('%12s %.3f %.3f\n', cls, ap1, ap2);
      else
        fprintf('%12s %.3f\n', cls, ap1);
      end
    else
      if do_auc_ap
        fprintf('%.3f %.3f\n', ap1);
      else
        fprintf('%.3f\n', ap1);
      end
    end
    score1(i) = ap1;
  catch
    score1(i) = nan;
    score2(i) = nan;
    if showcls
      fprintf('%12s -\n', cls);
    else
      fprintf('-\n');
    end
  end
end

a1 = nanmean(score1);
if do_auc_ap
  a2 = nanmean(score2);
end
if showcls
  fprintf('%s\n', repmat('-', [1 12]));
  if do_auc_ap
    fprintf('%12s %.3f %.3f\n', 'mAP', a1, a2);
  else
    fprintf('%12s %.3f\n', 'mAP', a1);
  end
else
  if do_auc_ap
    fprintf('%.3f %.3f\n', a1, a2);
  else
    fprintf('%.3f\n', a1);
  end
end

author: rbgirshick

repo: voc-dpm

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


text:
function [M,Vr,Vc] = meanShiftIm( X,sigSpt,sigRng,softFlag,maxIter,minDel )
% Applies the meanShift algorithm to a joint spatial/range image.
%
% See "Mean Shift Analysis and Applications" by Comaniciu & Meer for info.
%
% Assumes X is an MxNxP array, where an X(i,j,:) represents the range data
% at locations (i,j).  This function runs meanShift on each of the MxN data
% points.  It takes advantage of the lattice structure of an image for
% efficiency (it only needs to calculate full distance between two points
% if they are near each other spatially).
%
% In the original formulation of the algorithm, after normalization of the
% data, the search window around each point x has radius 1 (ie
% corresponding to 1 std of the data). That is the search window only
% encloses 2*s+1 pixels, and of those, all which fall within 1 unit from x
% are used to calcluate the new mean.  If softFlag==0 the original
% formulation is used.  If softFlag==1, instead of using a fixed radius,
% each point p is used in the calulation of the mean with points close to x
% given significantly more weight.  Specifically, each point p is given
% weight exp(-dist(x,p)).  So instead of having a fixed cutoff at r, the
% cutoff is 'soft' (same idea as in softmax), and occurs at approximately
% r.  The implementation remains efficient by actually using a hard cutoff
% at points further then 2r spatially from x.
%
% The resulting matrix M is of size MxNx(P+2).  M(i,j,1) represents the
% convergent row location of X(i,j,:) - (which had initial row location i)
% and M(i,j,2) represents the final column location.  M(i,j,p+2) represents
% the convergent value for X(i,j,p).  The optionaly outputs Vr and Vc are
% 2D arrays where Vr(i,j)=M(i,j,1)-i and Vc(i,j)=M(i,j,2)-j.  That is they
% represent the spatial offset between the original location of a point and
% its convergent location.  Display using quiver(Vc,Vr,0).
%
% USAGE
%  [M,Vr,Vc] = meanShiftIm( X,sigSpt,sigRng,[softFlag],[maxIter],[minDel] )
%
% INPUTS
%  X        - MxNxP data array, P may be 1
%  sigSpt   - integer specifying spatial standard deviation
%  sigRng   - value specifying the standard deviation of the range data
%  softFlag - [0]- see above
%  maxIter  - [100] maximum number of iterations per data point
%  minDel   - [.001] minimum amount of spatial change defining convergence
%
% OUTPUTS
%  M        - array of convergent locations [see above]
%  Vr       - spatial motion in row direction
%  Vc       - spatial motion in col direction
%
% EXAMPLE
%  I=double(imread('cameraman.tif'))/255;
%  [M,Vr,Vc] = meanShiftIm( I,5,.2 );
%  figure(1); im(I); figure(2); im( M(:,:,3) );
%  % color image:
%  I=double(imread('hestain.png'))/255;
%  [M,Vr,Vc] = meanShiftIm( I,5,.2 );
%  figure(1); im(I); figure(2); im( M(:,:,3:end) );
%
% See also MEANSHIFT, MEANSHIFTIMEXPLORE
%
% Piotr's Computer Vision Matlab Toolbox      Version 2.0
% Copyright 2014 Piotr Dollar.  [pdollar-at-gmail.com]
% Licensed under the Simplified BSD License [see external/bsd.txt]

[sigSpt,er] = checkNumArgs( sigSpt, 1, 0, 1 ); error(er);
if( nargin<4 || isempty(softFlag)); softFlag = 0; end
if( nargin<5 || isempty(maxIter) ); maxIter = 100; end
if( nargin<6 || isempty(minDel)); minDel = .001; end

[mrows, ncols, p] = size(X); p = p+2;
[gridRs, gridCs] = ndgrid( 1:mrows, 1:ncols );
data = cat( 3, cat( 3, gridRs/sigSpt, gridCs/sigSpt), X/sigRng );

%%% MAIN LOOP
M = data;
ticId = ticStatus('meanShiftIm');  %t0 = clock;  tlast = t0;
if( softFlag ); radius = sigSpt*2; else radius = sigSpt; end
for i=1:mrows; for j=1:ncols; %#ok<ALIGN>
    Mij = data(i,j,:); Mij = Mij(:)';
    itercount = 0; diff = 1;
    while( itercount < maxIter && diff>minDel )

      % get data which is possibly relevant (within spatial range)
      r = round( Mij(1)*sigSpt );  c = round( Mij(2)*sigSpt );
      boundsr = max(1,r-radius):min(mrows,r+radius);
      boundsc = max(1,c-radius):min(ncols,c+radius);
      dataWin = data( boundsr, boundsc, : );
      dataWinF = reshape( dataWin, [], p );

      % get next mean
      MijOld = Mij;
      n = size( dataWinF, 1);
      D = sum( (dataWinF - ones(n,1)*Mij).^2, 2 );
      if( softFlag )
        S = exp( -D ); sumS = sum(S); Srep = S(:,ones(1,p));
        Mij = sum( dataWinF .* Srep, 1 ) / sumS;
      else
        dataWinF = dataWinF( D < 1, : );
        Mij = sum( dataWinF, 1 ) / size( dataWinF,1 );
      end

      % check if Mij changed [only on basis of x,y location]
      diff = sum( (MijOld(1:2)-Mij(1:2)).^2 );
      itercount = itercount+1;

    end
    M(i,j,:) = Mij(:);
    fracdone = ((i-1)*ncols+j) / (mrows*ncols);
    tocStatus( ticId, fracdone );
  end
end
M = cat(3, M(:,:,1:2)*sigSpt, M(:,:,3:end)*sigRng );

%%% Output spatial difference
if( nargout>1 )
  Vr = M(:,:,1)-gridRs;  Vc = M(:,:,2)-gridCs;
end

author: pdollar

repo: toolbox

save_path: github-repos/MATLAB/pdollar-toolbox


text:
function X_rec = recoverData(Z, U, K)
%RECOVERDATA Recovers an approximation of the original data when using the
%projected data
%   X_rec = RECOVERDATA(Z, U, K) recovers an approximation the
%   original data that has been reduced to K dimensions. It returns the
%   approximate reconstruction in X_rec.
%

% You need to return the following variables correctly.
X_rec = zeros(size(Z, 1), size(U, 1));

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the approximation of the data by projecting back
%               onto the original space using the top K eigenvectors in U.
%
%               For the i-th example Z(i,:), the (approximate)
%               recovered data for dimension j is given as follows:
%                    v = Z(i, :)';
%                    recovered_j = v' * U(j, 1:K)';
%
%               Notice that U(j, 1:K) is a row vector.
%

for i = 1:size(Z, 1)
    v = Z(i, :);

    recovered_j = v * U(:, 1:K)';
    X_rec(i, :) = recovered_j;
end

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

end

author: ecmadao

repo: Coding-Guide


text:
% DEMUSPSVARGPLVM2 Demonstrate linear variational GPLVM (Bayesian PCA) on USPS data.

% VARGPLVM

% Fix seeds
randn('seed', 1e5);
rand('seed', 1e5);

dataSetName = 'usps';
experimentNo = 2;
printDiagram = 1;

% load data
[YTrain, lblsTrain, YTest, lblsTest] = lvmLoadData(dataSetName);

% Set up model
options = vargplvmOptions('dtcvar');
options.kern = {'linard2', 'white'};
options.numActive = 20;

options.optimiser = 'scg';
latentDim = 10;
d = size(YTrain, 2);

iters = 1000;
display = 1;

% create a separate vargplvm for each digit
for i=1:10
    %
    Y = YTrain(lblsTrain(:,i)==1,:);

    model = vargplvmCreate(latentDim, d, Y, options);
    %
    model = vargplvmParamInit(model, model.m, model.X);

    model = vargplvmOptimise(model, display, iters);

    varmodel{i} = model;
    %
end

iters = 100;
display = 0;

capName = dataSetName;
capName(1) = upper(capName(1));
modelType = model.type;
modelType(1) = upper(modelType(1));
save(['dem' capName modelType num2str(experimentNo) '.mat'], 'varmodel');

% measure performance on test data
indexPresent = 1:size(YTest,2);
TestError = 0;
for n=1:size(YTest,1)
    %
    % compute the approximate class conditional density for each digit
    for i=1:10
       prob(n,i) = vargplvmProbabilityCompute(varmodel{i}, YTest(n,:), indexPresent)
    end
    [maxP C] = max(prob(n,:));
    if lblsTest(n,C) == 0
        TestError = TestError + 1;
    end
    %
end

%
save(['dem' capName modelType num2str(experimentNo) '.mat'], 'varmodel', 'prob', 'TestError');

author: SheffieldML

repo: vargplvm


text:
% Test file for chebfun constructor (periodic).

function pass = test_constructor_basic_periodic(pref)

if ( nargin == 0 )
    pref = chebfunpref();
end

%% Test construction of some basic functions:

% Some basic test functions:
FF = {@(x) exp(sin(pi*x)), @(x) exp([sin(pi*x), cos(pi*x)]), @(x) exp([sin(pi*x), cos(pi*x), -cos(pi*x).^2])};

for j = 1:numel(FF);
    % Initialise k:
    k = 0;

    % Pick the test function:
    F = FF{j};

    % Test on [-1 1]:
    f = chebfun(F, [-1, 1], pref, 'periodic');
    g = chebfun(F, [-1, 1], pref, 'trig');
    xx = linspace(-1, 1);
    err = norm(feval(f, xx) - F(xx), inf);
    pass(j, k+1) = (err < 10*eps*vscale(f)) && (norm(f-g) < pref.chebfuneps);
    pass(j, k+2) = err < 50*pref.chebfuneps;
    k = k + 2;

    % Test on [-1 1] (no domain passed):
    f = chebfun(F, pref, 'periodic');
    xx = linspace(-1, 1);
    err = norm(feval(f, xx) - F(xx), inf);
    pass(j, k+1) = err < 10*eps*vscale(f);
    pass(j, k+2) = err < 500*pref.chebfuneps;
    k = k + 2;

    % Test on [0 10000]:
    f = chebfun(F, [-100, 100], pref, 'periodic');
    g = chebfun(F, [-100, 100], pref, 'trig');
    xx = linspace(-100, 100);
    err = norm(feval(f, xx) - F(xx), inf);
    pass(j, k+1) = err < 1e3*eps*vscale(f);
    pass(j, k+2) = (err < 100*hscale(f)*pref.chebfuneps) && (norm(f-g) <...
        100*pref.chebfuneps);
    k = k + 2;

end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function w_omega = rotationalPID_velocity(w_R_b,w_R_b_des,w_omega_des,Kp)

    % ROTATIONALPID_VELOCITY implements a controller trajectory tracking inside
    %                        the group SO(3). The angular velocity is assumed
    %                        to be a control input.
    %
    % FORMAT: w_omega = rotationalPID_velocity(w_R_b,w_R_b_des,w_omega_des,Kp)
    %
    % INPUT: - w_R_b       = [3 * 3] rotation matrix
    %        - w_R_b_des   = [3 * 3] desired rotation matrix
    %        - w_omega     = [3 * 1] angular velocity (expressed in the world frame)
    %        - w_omega_des = [3 * 1] desired angular velocity (expressed in the world frame)
    %        - Kp          = [3 * 3] orientation gains
    %
    % OUTPUT: - w_omega    = [3 * 1] input angular velocity
    %
    % Authors: Daniele Pucci, Marie Charbonneau, Gabriele Nava
    %
    %          all authors are with the Italian Istitute of Technology (IIT)
    %          email: name.surname@iit.it
    %
    % Genoa, Dec 2017
    %

    %% --- Initialization ---

    skv      = wbc.skewVee(w_R_b*transpose(w_R_b_des));
    w_omega  = w_omega_des -Kp*skv;
end

author: robotology

repo: whole-body-controllers


text:
function out = proj_halfspace_box(x,a,b,lb,ub)
%PROJ_HALFSPACE_BOX computes the orthogonal projection onto the intersection of a halfspace and a box
%                   {x : <a,x> <= b, lb<=x<=ub}
%
%  Usage:
%  out = PROJ_HALFSPACE_BOX(x,a,b,[lb],[ub])
%  ===========================================
%  Input:
%  x - point to be projected (vector/matrix)
%  a - vector/matrix
%  b - scalar
%  lb - lower bound (vector/matrix/scalar) [default: -inf]
%  ub - upper bound (vector/matrix/scalar) [default: inf]
%  ===========================================
%  Assumptions:
%  The intersection of the halfspace and the box is nonempty
%  ===========================================
%  Output:
%  out - projection vector

% This file is part of the FOM package - a collection of first order methods for solving convex optimization problems
% Copyright (C) 2017 Amir and Nili Beck
%
% 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/>.

%reading the user x and setting defalut values when required.
if (nargin < 3)
    error ('usage: proj_halfspace_box(x,a,b,[lb],[ub]') ;
end

if (nargin < 5)
    %upper bound is not given, setting to defalut value :inf
    ub = inf ;
end
if ((nargin < 4) || (isempty( lb)))
    %lower bound is not given, setting to defalut value :0
    lb = -inf;
end

if  (any(any(lb > ub)))
    error('Set is infeasible') ;
end

if trace(a'*min(max(x,lb),ub)) <= b
    out =  min(max(x,lb),ub) ;
else
    %use proj_hyperplane_box
    out = proj_hyperplane_box(x,a,b,lb,ub) ;
end

author: hiroyuki-kasai

repo: SGDLibrary


text:
function zData = constructMissingSurvey(zDataTrue, Gdraw, PsiAll, dataValues)

%PsiAll=PsiDraw_prop;
% Initialize
Ppsi = dataValues.Ppsi; %load selection matrix for the states (local mean), where survey data is available
PsiZ = PsiAll * Ppsi'; %only use those variables
[T, Mz] = size(zDataTrue);
zData = zDataTrue;

% for each survey variable
for iM = 1:Mz

    % determine time of incoming survey data
    Tz              = find(isfinite(zData(:,iM)),1,'first');
    iNaNs           = find(isnan(zData(:,iM)));
    iNaNs(iNaNs<Tz) = [];

    zDataDraw_i     = PsiZ(:, iM) + sqrt(Gdraw(:,iM)).*randn(T,1);
    % fill missing values (to previous values)
    zData(iNaNs,iM) = zDataDraw_i(iNaNs);
end

end

author: european-central-bank

repo: BEAR-toolbox


text:
function f1 = p03_f1 ( x )

%*****************************************************************************80
%
%% P03_F1 evaluates the first derivative for problem 3.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 February 2002
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the value of the variable.
%
%    Output, real F1, the first derivative of the
%    objective function.
%
  f1 = ( 4.0 * x * x + 4.0 ) * x + 1.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function euler = rotMat2euler(R)
    % from paper: "Adaptive Filter for a Miniature MEMS Based Attitude and
    % Heading Reference System" by Wang et al, IEEE.

    phi = atan2(R(3,2,:), R(3,3,:) );
    theta = -atan(R(3,1,:) ./ sqrt(1-R(3,1,:).^2) );
    psi = atan2(R(2,1,:), R(1,1,:) );

    euler = [phi(1,:)' theta(1,:)' psi(1,:)'];
end

author: xioTechnologies

repo: Gait-Tracking-With-x-IMU


text:
% generates a circle as a test track for tracking control applications
close all;

r = 50;
phi = 0:0.05:2*pi;
x_m = (cos(phi) - 1)*r;
y_m = sin(phi)*r;
s_m = [0, cumsum(sqrt(diff(x_m).^2+diff(y_m).^2))];
kappa_radpm = ones(length(x_m), 1)/r;
ax = zeros(length(x_m), 1);
vx = ones(length(x_m), 1)*20;
banking = zeros(length(x_m), 1);
figure;
plot(x_m, y_m);
axis equal;
grid on;

% write to csv
dlmwrite('Circle.csv', [s_m', x_m', y_m', phi', kappa_radpm, vx, ax, banking],';');

author: TUMFTM


text:
function benchNewAlg(benchdir,pres,alg)
% benchNewAlg(benchdir,pres,alg)
%
% benchdir = root directory of benchmark results
% pres = presentation, 'gray' or 'color'
% alg = algorithm directory containing *.bmp pb files
%
% Things you should do before running this script:
%
% (1) Download and untar BSDS images and human data.
% (2) Download the code (you must have done that if you're
%     reading this!
% (3) Edit Dataset/bsdsRoot.m script to point to your BSDS
%     directory.
% (4) Run 'gmake install' from the code directory.
% (5) Create a benchmark directory e.g. 'bench' and untar the human
%     benchmark data into it.  That should create 'bench/gray'
%     and 'bench/color' sub-directories, each with a 'human'
%     subdirectory.
% (6) Put your algorithm's pb files into the directory under the
%     proper 'gray' or 'color' directory.
% (7) Create the name.txt and about.html files in your algorithm
%     directory.

algdir = fullfile(benchdir,pres,alg);

fprintf(2,'Benchmarking algorithm (takes a while!)...\n');
boundaryBench(algdir,pres);

fprintf(2,'Generating benchmark graphs for this algorithm...\n');
boundaryBenchGraphs(algdir);

fprintf(2,'Generating benchmark graphs to compare algorithms...\n');
boundaryBenchGraphsMulti(benchdir);

fprintf(2,'Generating benchmark web pages...\n');
boundaryBenchHtml(benchdir);

author: Cloud-CV

repo: object-proposals


text:
function OBJ=read_wobj(fullfilename)
% Read the objects from a Wavefront OBJ file
%
% OBJ=read_wobj(filename);
%
% OBJ struct containing:
%
% OBJ.vertices : Vertices coordinates
% OBJ.vertices_texture: Texture coordinates
% OBJ.vertices_normal : Normal vectors
% OBJ.vertices_point  : Vertice data used for points and lines
% OBJ.material : Parameters from external .MTL file, will contain parameters like
%           newmtl, Ka, Kd, Ks, illum, Ns, map_Ka, map_Kd, map_Ks,
%           example of an entry from the material object:
%       OBJ.material(i).type = newmtl
%       OBJ.material(i).data = 'vase_tex'
% OBJ.objects  : Cell object with all objects in the OBJ file,
%           example of a mesh object:
%       OBJ.objects(i).type='f'
%       OBJ.objects(i).data.vertices: [n x 3 double]
%       OBJ.objects(i).data.texture:  [n x 3 double]
%       OBJ.objects(i).data.normal:   [n x 3 double]
%
% Example,
%   OBJ=read_wobj('examples\example10.obj');
%   FV.vertices=OBJ.vertices;
%   FV.faces=OBJ.objects(3).data.vertices;
%   figure, patch(FV,'facecolor',[1 0 0]); camlight
%
% Function is written by D.Kroon University of Twente (June 2010)

verbose=true;

if(exist('fullfilename','var')==0)
    [filename, filefolder] = uigetfile('*.obj', 'Read obj-file');
    fullfilename = [filefolder filename];
end
filefolder = fileparts( fullfilename);
if(verbose),disp(['Reading Object file : ' fullfilename]); end

% Read the DI3D OBJ textfile to a cell array
file_words = file2cellarray( fullfilename);
% Remove empty cells, merge lines split by "\" and convert strings with values to double
[ftype fdata]= fixlines(file_words);

% Vertex data
vertices=[]; nv=0;
vertices_texture=[]; nvt=0;
vertices_point=[]; nvp=0;
vertices_normal=[]; nvn=0;
material=[];

% Surface data
no=0;

% Loop through the Wavefront object file
for iline=1:length(ftype)
    if(mod(iline,10000)==0),
        if(verbose),disp(['Lines processed : ' num2str(iline)]); end
    end

    type=ftype{iline}; data=fdata{iline};

    % Switch on data type line
    switch(type)
        case{'mtllib'}
            if(iscell(data))
                datanew=[];
                for i=1:length(data)
                    datanew=[datanew data{i}];
                    if(i<length(data)), datanew=[datanew ' ']; end
                end
                data=datanew;
            end

            filename_mtl=fullfile(filefolder,data);
            material=readmtl(filename_mtl,verbose);
        case('v') % vertices
            nv=nv+1;
            if(length(data)==3)
                % Reserve block of memory
                if(mod(nv,10000)==1), vertices(nv+1:nv+10001,1:3)=0; end
                % Add to vertices list X Y Z
                vertices(nv,1:3)=data;
            else
                % Reserve block of memory
                if(mod(nv,10000)==1), vertices(nv+1:nv+10001,1:4)=0; end
                % Add to vertices list X Y Z W
                vertices(nv,1:4)=data;
            end
        case('vp')
            % Specifies a point in the parameter space of curve or surface
            nvp=nvp+1;
            if(length(data)==1)
                % Reserve block of memory
                if(mod(nvp,10000)==1), vertices_point(nvp+1:nvp+10001,1)=0; end
                % Add to vertices point list U
                vertices_point(nvp,1)=data;
            elseif(length(data)==2)
                % Reserve block of memory
                if(mod(nvp,10000)==1), vertices_point(nvp+1:nvp+10001,1:2)=0; end
                % Add to vertices point list U V
                vertices_point(nvp,1:2)=data;
            else
                % Reserve block of memory
                if(mod(nvp,10000)==1), vertices_point(nvp+1:nvp+10001,1:3)=0; end
                % Add to vertices point list U V W
                vertices_point(nvp,1:3)=data;
            end
        case('vn')
            % A normal vector
            nvn=nvn+1; if(mod(nvn,10000)==1),  vertices_normal(nvn+1:nvn+10001,1:3)=0; end
            % Add to vertices list I J K
            vertices_normal(nvn,1:3)=data;
        case('vt')
            % Vertices Texture Coordinate in photo
            % U V W
            nvt=nvt+1;
            if(length(data)==1)
                % Reserve block of memory
                if(mod(nvt,10000)==1), vertices_texture(nvt+1:nvt+10001,1)=0; end
                % Add to vertices texture list U
                vertices_texture(nvt,1)=data;
            elseif(length(data)==2)
                % Reserve block of memory
                if(mod(nvt,10000)==1), vertices_texture(nvt+1:nvt+10001,1:2)=0; end
                % Add to vertices texture list U V
                vertices_texture(nvt,1:2)=data;
            else
                % Reserve block of memory
                if(mod(nvt,10000)==1), vertices_texture(nvt+1:nvt+10001,1:3)=0; end
                % Add to vertices texture list U V W
                vertices_texture(nvt,1:3)=data;
            end
        case('l')
            no=no+1; if(mod(no,10000)==1), objects(no+10001).data=0; end
            array_vertices=[];
            array_texture=[];
            for i=1:length(data),
                switch class(data)
                    case 'cell'
                        tvals=str2double(stringsplit(data{i},'/'));
                    case 'string'
                        tvals=str2double(stringsplit(data,'/'));
                    otherwise
                        tvals=data(i);
                end
                val=tvals(1);
                if(val<0), val=val+1+nv; end
                array_vertices(i)=val;
                if(length(tvals)>1),
                    val=tvals(2);
                    if(val<0), val=val+1+nvt; end
                    array_texture(i)=val;
                end
            end
            objects(no).type='l';
            objects(no).data.vertices=array_vertices;
            objects(no).data.texture=array_texture;
        case('f')
            no=no+1; if(mod(no,10000)==1), objects(no+10001).data=0; end
            array_vertices=[];
            array_texture=[];
            array_normal=[];
            for i=1:length(data);
                switch class(data)
                    case 'cell'
                        tvals=str2double(stringsplit(data{i},'/'));
                    case 'string'
                        tvals=str2double(stringsplit(data,'/'));
                    otherwise
                        tvals=data(i);
                end
                val=tvals(1);

                if(val<0), val=val+1+nv; end
                array_vertices(i)=val;
                if(length(tvals)>1),
                    if(isfinite(tvals(2)))
                        val=tvals(2);
                        if(val<0), val=val+1+nvt; end
                        array_texture(i)=val;
                    end
                end
                if(length(tvals)>2),
                    val=tvals(3);
                    if(val<0), val=val+1+nvn; end
                    array_normal(i)=val;
                end
            end

            % A face of more than 3 indices is always split into
            % multiple faces of only 3 indices.
            objects(no).type='f';
            findex=1:min (3,length(array_vertices));

            objects(no).data.vertices=array_vertices(findex);
            if(~isempty(array_texture)),objects(no).data.texture=array_texture(findex); end
            if(~isempty(array_normal)),objects(no).data.normal=array_normal(findex); end
            for i=1:length(array_vertices)-3;
                no=no+1; if(mod(no,10000)==1), objects(no+10001).data=0; end
                findex=[1 2+i 3+i];
                findex(findex>length(array_vertices))=findex(findex>length(array_vertices))-length(array_vertices);
                objects(no).type='f';
                objects(no).data.vertices=array_vertices(findex);
                if(~isempty(array_texture)),objects(no).data.texture=array_texture(findex); end
                if(~isempty(array_normal)),objects(no).data.normal=array_normal(findex); end
            end
        case{'#','$'}
            % Comment
            tline='  %';
            if(iscell(data))
                for i=1:length(data), tline=[tline ' ' data{i}]; end
            else
                tline=[tline data];
            end
            if(verbose), disp(tline); end
        case{''}
        otherwise
            no=no+1;
            if(mod(no,10000)==1), objects(no+10001).data=0; end
            objects(no).type=type;
            objects(no).data=data;
    end
end

% Initialize new object list, which will contain the "collapsed" objects
objects2(no).data=0;

index=0;

i=0;
while (i<no), i=i+1;
    type=objects(i).type;
    % First face found
    if((length(type)==1)&&(type(1)=='f'))
        % Get number of faces
        for j=i:no
            type=objects(j).type;
            if((length(type)~=1)||(type(1)~='f'))
                j=j-1; break;
            end
        end
        numfaces=(j-i)+1;

        index=index+1;
        objects2(index).type='f';
        % Process last face first to allocate memory
        objects2(index).data.vertices(numfaces,:)= objects(i).data.vertices;
        if(isfield(objects(i).data,'texture'))
            objects2(index).data.texture(numfaces,:) = objects(i).data.texture;
        else
            objects2(index).data.texture=[];
        end
        if(isfield(objects(i).data,'normal'))
            objects2(index).data.normal(numfaces,:)  = objects(i).data.normal;
        else
            objects2(index).data.normal=[];
        end
        % All faces to arrays
        for k=1:numfaces
            objects2(index).data.vertices(k,:)= objects(i+k-1).data.vertices;
            if(isfield(objects(i).data,'texture'))
                objects2(index).data.texture(k,:) = objects(i+k-1).data.texture;
            end
            if(isfield(objects(i).data,'normal'))
                objects2(index).data.normal(k,:)  = objects(i+k-1).data.normal;
            end
        end
        i=j;
    else
        index=index+1;
        objects2(index).type=objects(i).type;
        objects2(index).data=objects(i).data;
    end
end

% Add all data to output struct
OBJ.objects=objects2(1:index);
OBJ.material=material;
OBJ.vertices=vertices(1:nv,:);
OBJ.vertices_point=vertices_point(1:nvp,:);
OBJ.vertices_normal=vertices_normal(1:nvn,:);
OBJ.vertices_texture=vertices_texture(1:nvt,:);
if(verbose),disp('Finished Reading Object file'); end

function twords=stringsplit(tline,tchar)
% Get start and end position of all "words" separated by a char
i=find(tline(2:end-1)==tchar)+1; i_start=[1 i+1]; i_end=[i-1 length(tline)];
% Create a cell array of the words
twords=cell(1,length(i_start)); for j=1:length(i_start), twords{j}=tline(i_start(j):i_end(j)); end

function file_words=file2cellarray(filename)
% Open a DI3D OBJ textfile
fid=fopen(filename,'r');
file_text=fread(fid, inf, 'uint8=>char')';
fclose(fid);
file_lines = regexp(file_text, '\n+', 'split');
file_words = regexp(file_lines, '\s+', 'split');

function [ftype fdata]=fixlines(file_words)
ftype=cell(size(file_words));
fdata=cell(size(file_words));

iline=0; jline=0;
while(iline<length(file_words))
    iline=iline+1;
    twords=removeemptycells(file_words{iline});
    if(~isempty(twords))
        % Add next line to current line when line end with '\'
        while(strcmp(twords{end},'\')&&iline<length(file_words))
            iline=iline+1;
            twords(end)=[];
            twords=[twords removeemptycells(file_words{iline})];
        end
        % Values to double

        type=twords{1};
        stringdold=true;
        j=0;
        switch(type)
            case{'#','$'}
                for i=2:length(twords)
                    j=j+1; twords{j}=twords{i};
                end
            otherwise
                for i=2:length(twords)
                    str=twords{i};
                    val=str2double(str);
                    stringd=~isfinite(val);
                    if(stringd)
                        j=j+1; twords{j}=str;
                    else
                        if(stringdold)
                            j=j+1; twords{j}=val;
                        else
                            twords{j}=[twords{j} val];
                        end
                    end
                    stringdold=stringd;
                end
        end
        twords(j+1:end)=[];
        jline=jline+1;
        ftype{jline}=type;
        if(length(twords)==1), twords=twords{1}; end
        fdata{jline}=twords;
    end
end
ftype(jline+1:end)=[];
fdata(jline+1:end)=[];

function b=removeemptycells(a)
j=0; b={};
for i=1:length(a);
    if(~isempty(a{i})),j=j+1; b{j}=a{i}; end;
end

function  objects=readmtl(filename_mtl,verbose)
if(verbose),disp(['Reading Material file : ' filename_mtl]); end
file_words=file2cellarray(filename_mtl);
% Remove empty cells, merge lines split by "\" and convert strings with values to double
[ftype fdata]= fixlines(file_words);

% Surface data
objects.type(length(ftype))=0;
objects.data(length(ftype))=0;
no=0;
% Loop through the Wavefront object file
for iline=1:length(ftype)
  type=ftype{iline}; data=fdata{iline};

    % Switch on data type line
    switch(type)
        case{'#','$'}
            % Comment
            tline='  %';
            if(iscell(data))
                for i=1:length(data), tline=[tline ' ' data{i}]; end
            else
                tline=[tline data];
            end
            if(verbose), disp(tline); end
        case{''}
        otherwise
            no=no+1;
            if(mod(no,10000)==1), objects(no+10001).data=0; end
            objects(no).type=type;
            objects(no).data=data;
    end
end
objects=objects(1:no);
if(verbose),disp('Finished Reading Material file'); end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


function EXP_A = get_RS_matrix_exponential( Q, dt, xi, drifts, vols, psi_J)
N = length(xi);
m_0 = length(drifts);

Qt = dt*Q';
EXP_A = ones(m_0,m_0,N);

drifts = dt*drifts;
vols = (0.5*dt)*vols.^2;

if nargin < 8
    for j = 1:N
        EXP_A(:,:,j) = expm(Qt + diag(drifts*xi(j) - vols*xi(j)^2));  %This incorporates the dt
    end
else
    for j = 1:N
        EXP_A(:,:,j) = expm(Qt + diag(drifts*xi(j) - vols*xi(j)^2 + dt*psi_J(xi(j))));  %This incorporates the dt
    end
end

end


text:
%UDC Uncorrelated normal based quadratic Bayes classifier (BayesNormal_U)
%
%  W = UDC(A)
%  W = A*UDC
%
% INPUT
%  A  input dataset
%
% OUTPUT
%  W   output mapping
%
% DESCRIPTION
% Computation a quadratic classifier between the classes in the
% dataset A assuming normal densities with uncorrelated features.
%
% The use of probabilistic labels is supported.  The classification A*W is
% computed by normal_map.
%
% EXAMPLES
% PREX_DENSITY
%
% SEE ALSO (PRTools Guide)
% MAPPINGS, DATASETS, NMC, NMSC, LDC, QDC, QUADRC, NORMAL_MAP

% Copyright: R.P.W. Duin, duin@ph.tn.tudelft.nl
% Faculty of Applied Physics, Delft University of Technology
% P.O. Box 5046, 2600 GA Delft, The Netherlands

% $Id: udc.m,v 1.6 2007/06/05 12:45:44 duin Exp $

function W = udc(a)

	if nargin == 0
		W = prmapping(mfilename);
		W = setname(W,'Bayes-Normal-U');
		return
	end

	islabtype(a,'crisp','soft');
	isvaldfile(a,2,2); % at least 2 objects per class, 2 classes
	[m,k,c] = getsize(a);

	[U,G] = meancov(a); %computing mean and covariance matrix
	p = getprior(a);
	for j = 1:c
		G(:,:,j) = diag(diag(G(:,:,j)));
	end
	w.mean = +U;
	w.cov = G;
	w.prior = p;
	%W = prmapping('normal_map','trained',w,getlab(U),k,c);
	W = normal_map(w,getlab(U),k,c);
	W = setname(W,'Bayes-Normal-U');
	W = setcost(W,a);

return

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
% Plots of temperature, pressure and density with altitude increments.
sldata;
chart;
A=chart;

A(:,2)=A(:,2)*T1;
A(:,3)=A(:,3)*p1;
A(:,4)=A(:,4)*rho1;

% Absolute variables are separated in different vectors.

h=A(:,1);
T=A(:,2);
p=A(:,3);
rho=A(:,4);

plott;
figure;
plotp;
figure;
plotrho;

author: Sable

repo: mcbench-benchmarks


text:
% NAME:
% test_RANSAC_line_02.m
%
% DESC:
% test to estimate the parameters of a line using real data

close all
clear
% clc

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load LineData

% set RANSAC options
options.epsilon = 1e-6;
options.P_inlier = 0.99;
options.sigma = 1;
options.est_fun = @estimate_line;
options.man_fun = @error_line;
options.mode = 'MSAC';
options.Ps = [];
options.notify_iters = [];
options.min_iters = 100;
options.fix_seed = false;
options.reestimate = true;
options.stabilize = false;

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% RANSAC
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% run RANSAC
[results, options] = RANSAC(X, options);

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Results Visualization
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure;
hold on
ind = results.CS;
plot(X(1, ind), X(2, ind), '.g')
plot(X(1, ~ind), X(2, ~ind), '.r')
xlabel('x')
ylabel('y')
title('RANSAC results for 2D line estimation')
legend('Inliers', 'Outliers')
axis equal tight

author: RANSAC

repo: RANSAC-Toolbox


text:
% Compute derivative k'(x^p,x^q) of a stationary covariance k(d2) (ard or iso)
% w.r.t. to squared distance d2 = (x^p - x^q)'*inv(P)*(x^p - x^q) measure. Here
% P is either diagonal with ARD parameters ell_1^2,...,ell_D^2 where D is the
% dimension of the input space or ell^2 times the unit matrix for isotropic
% covariance.
% The derivatives can only be computed for one of the following eight
% covariance functions: cov{Matern|PP|RQ|SE}{iso|ard}.
%
% Copyright (c) by Hannes Nickisch, 2013-10-28.
%
% See also INFFITC.M, COVFITC.M.

function Kp = cov_deriv_sq_dist(cov,hyp,x,z)
  if nargin<4, z = []; end                                 % make sure, z exists
  xeqz = numel(z)==0; dg = strcmp(z,'diag') && numel(z)>0;      % determine mode

  if iscell(cov), covstr = cov{1}; else covstr = cov; end
  if ~ischar(covstr), covstr = func2str(covstr); end
  if numel([strfind(covstr,'iso'),strfind(covstr,'ard')])==0
    error('Only iso|ard covariances allowed for derivatives w.r.t. xu.')
  elseif numel([strfind(covstr,'covLIN');
                strfind(covstr,'covGabor');
                strfind(covstr,'covPER')])>0
    error('Gabor|LIN|PER covariances not allowed for derivatives w.r.t. xu.')
  end

  [n,D] = size(x);
  if numel(strfind(covstr,'iso')), id = 1:D; else id = 1; end  % *iso covariance
  ell1 = exp(hyp(1));                        % first characteristic length scale

  Kp = feval(cov{:},hyp,x,z,1);              % use derivative w.r.t. log(ell(1))
  % precompute squared distances
  if dg                                                             % vector kxx
    d2 = zeros(n,1);
  else
    if xeqz                                               % symmetric matrix Kxx
      d2 = sq_dist(x(:,id)'/ell1);
    else                                                 % cross covariances Kxz
      d2 = sq_dist(x(:,id)'/ell1,z(:,id)'/ell1);
    end
  end
  Kp = -1/2*Kp./d2; Kp(d2==0) = 0;

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function [alignedSignal] = alignSignal(responseSignal, alignmentSignal,timeSeq,varargin)
	% Aligns values in a response signal (can be multiple response signals) to binary points in an alignment signal (e.g. 1=align to this time-point, 0=don't align).
	% Biafra Ahanonu
	% started: 2013.11.13 [23:47:34]
	% inputs
		% responseSignal = MxN matrix of M signals over N points
		% alignmentSignal = a 1xN vector of 0s and 1s, where 1s will be alignment points
		% timeSeq = 1xN sequence giving time around alignments points to process, e.g. -2:2.
	% options
		% overallAlign = align all response signals to alignmentSignal pts
	% outputs
		% alignedSignal = a matrix of size 1xlength(timeSeq) if sum all signals or Mxlength(timeSeq) if keep the sums for each signal separate

	[alignedSignal] = ciapkg.signal_processing.alignSignal(responseSignal, alignmentSignal,timeSeq,'passArgs', varargin);
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function err = err_fit(params, sample, tes)
% return residual

%size(sample)
%size(signal(params, tes))

err = sample - signal(params, tes);

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function export_protocol(iProtocol, iSubject, OutputFile)
% EXPORT_PROTOCOL: Export a protocol into a zip file.
%
% USAGE:  export_protocol(iProtocol, iSubject, OutputFile)
%         export_protocol(iProtocol, iSubject)        : Ask for the output filename
%         export_protocol(iProtocol)                  : Export all the subjects of protocol, ask for the output filename
%         export_protocol()                           : Export current protocol, ask for the output filename

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2012-2015; Martin Cousineau, 2019

global GlobalData;

%% ===== PARSE INPUTS =====
if (nargin < 3)
    OutputFile = [];
end
if (nargin < 2)
    iSubject = [];
end
if (nargin < 1) || isempty(iProtocol)
    iProtocol = bst_get('iProtocol');
else
    gui_brainstorm('SetCurrentProtocol', iProtocol);
end
% Check protocol indice
if isempty(iProtocol) || (iProtocol == 0)
    bst_error('Invalid protocol indice.', 'Export protocol', 0);
    return
end
% Get protocol description
ProtocolInfo = bst_get('ProtocolInfo');
% Get output filename
if isempty(OutputFile)
    % Get default directories
    LastUsedDirs = bst_get('LastUsedDirs');
	% Default output filename
    if isempty(iSubject)
        OutputFile = bst_fullfile(LastUsedDirs.ExportProtocol, file_standardize([ProtocolInfo.Comment, '.zip']));
    else
        sSubject = bst_get('Subject', iSubject);
        OutputFile = bst_fullfile(LastUsedDirs.ExportProtocol, file_standardize([ProtocolInfo.Comment, '_', sSubject.Name, '.zip']));
    end
    % File selection
    OutputFile = java_getfile('save', 'Export protocol', OutputFile, 'single', 'files', ...
                              {{'.zip'}, 'Zip files (*.zip)', 'ZIP'}, 1);
    if isempty(OutputFile)
        return
    end
    % Save new default export path
    LastUsedDirs.ExportProtocol = bst_fileparts(OutputFile);
    bst_set('LastUsedDirs', LastUsedDirs);
end
% Force saving database
db_save(1);

%% ===== ZIP FILES =====
% Progress bar
bst_progress('start', 'Export protocol', 'Creating zip file...');
% Cd to protocol folder
prevFolder = pwd;
cd(bst_fileparts(ProtocolInfo.SUBJECTS, 1));
% Prefixes to add to the folders
[tmp__, anatFolder] = bst_fileparts(ProtocolInfo.SUBJECTS, 1);
[tmp__, dataFolder] = bst_fileparts(ProtocolInfo.STUDIES, 1);
% Build list of files to zip
if isempty(iSubject)
    % Add the entire subject folder
    ListZip = {anatFolder};
    % List files in studies: add all files
    allFiles = dir(ProtocolInfo.STUDIES);
    for i = 1:length(allFiles)
        if (allFiles(i).name(1) ~= '.')
            ListZip{end+1} = bst_fullfile(dataFolder, allFiles(i).name);
        end
    end
else
    % Get default study for this subject
    sSubject = bst_get('Subject', iSubject, 1);
    sStudy   = bst_get('AnalysisIntraStudy', iSubject);
    % Create a temporary protocol.mat for future database update information
    ProtocolMat = struct();
    ProtocolMat.ProtocolInfo      = GlobalData.DataBase.ProtocolInfo(iProtocol);
    ProtocolMat.ProtocolSubjects  = GlobalData.DataBase.ProtocolSubjects(iProtocol);
    ProtocolMat.ProtocolStudies   = GlobalData.DataBase.ProtocolStudies(iProtocol);
    ProtocolMat.DbVersion         = GlobalData.DataBase.DbVersion;
    ProtocolMat.LastAccessDate    = datestr(now);
    ProtocolMat.LastAccessUserDir = bst_get('UserDir');
    % Remove useless fields
    StudyDir = ProtocolMat.ProtocolInfo.STUDIES;
    ProtocolMat.ProtocolInfo = rmfield(ProtocolMat.ProtocolInfo, 'STUDIES');
    ProtocolMat.ProtocolInfo = rmfield(ProtocolMat.ProtocolInfo, 'SUBJECTS');
    ProtocolMat.ProtocolSubjects.Subject = ProtocolMat.ProtocolSubjects.Subject(iSubject);
    [sStudies, iStudies] = bst_get('StudyWithSubject', ProtocolMat.ProtocolSubjects.Subject.FileName, 'default_study', 'intra_subject');
    ProtocolMat.ProtocolStudies.Study = ProtocolMat.ProtocolStudies.Study(iStudies);
    ProtocolFile = bst_fullfile(dataFolder, 'protocol.mat');
    bst_save(ProtocolFile, ProtocolMat, 'v7');
    % List all files that might be useful for this subject
    ListZip = {bst_fullfile(anatFolder, bst_fileparts(sSubject.FileName)), ...
               bst_fullfile(anatFolder, bst_get('DirDefaultSubject')), ...
               bst_fullfile(dataFolder, bst_fileparts(bst_fileparts(sStudy.FileName))), ...
               bst_fullfile(dataFolder, bst_get('DirDefaultStudy')), ...
               bst_fullfile(dataFolder, bst_get('DirAnalysisInter')), ...
               ProtocolFile};
end
% Zip
zip(OutputFile, ListZip);
% Restore protocol file
if ~isempty(iSubject)
    % Remove temporary protocol file
    file_delete(bst_fullfile(StudyDir, 'protocol.mat'), 1);
    % Save again the original one
    GlobalData.DataBase.isProtocolModified(iProtocol) = 1;
    db_save(1);
end
% Restore initial folder
cd(prevFolder);
% Close progress bar
bst_progress('stop');

author: brainstorm-tools


text:
function C = partitions(M,K)
%PARTITIONS Find all partitions of a set.
% C = PARTITIONS(N), for scalar N, returns a cell array, wherein each cell
% has a various number of other cells representing the partitions of the
% set of numbers {1,2,3,...N}.  The length of C is the Bell number for N.
% C = PARTITIONS(N), for vector N, returns the partitions of the vector
% elements treated as members of a set.
% C = PARTITIONS(N), for cell N, returns the partitions of the cell
% elements treated as members of a set.
% PARTITIONS(N,K) returns only the partitions of which have K members.
%
% K must be less than or equal to N for scalar N, or length(N).
%
%
% Examples:
%
%     C = partitions(4); % Find the partitions of set {1 2 3 4}
%     home  % Makes for nice display for small enough N.  Else use clc.
%     partdisp(C) % File accompanying PARTITIONS
%
%     C = partitions({'peanut','butter','yummy'});
%     home
%     partdisp(C)
%
%     C = partitions([5 7 8 50],3); % partitions of set {5 7 8 50}
%     home
%     partdisp(C,3)
%
%     C = partitions(['a','b','c','d']); % for longer chars, use {}, not []
%     home
%     partdisp(C)
%
% Class support for inputs N,K:
%      float: double, single, char (N only)
%
% See also, nchoosek, perms, npermutek (On the FEX)
%
% Author: Matt Fig,
% Contact: popkenai@yahoo.com
% Date: 5/17/2009

Kflg = false;  % No K passed in.
Cflg = false;  % A set not passed in.

if length(M)>1
    N = length(M); % User passed {'here','there','everywhere'} for example.
    Cflg = true;
else
    if iscell(M)  % User passed {2} for example.
        error('Set arguments must have more than one element.')
    end

    N = M;
end

if nargin>1
    Kflg = true;  % Used in while loop below, K passed in.
    S = stirling(N,K);  % The number of partitions, Stirling number.
    C = cell(S,min(1,K)); % Main Cell.
    cnt = 1;  % Start the while loop counter.
else
    K = 0;  % Since user doesn't want only certain partitions.
    S = ceil(sum((1:2*N).^N./cumprod(1:2*N))/exp(1)); % Bell number.
    C = cell(S,1); % Main Cell.

    if Cflg
        C{1} = {M};
    else
        C{1} = {1:N}; % First one is easy.
    end

    cnt = 2; % Start the while loop counter.
end

if K~=ceil(K) || numel(K)~=1
    error('Second argument must be an integer of length 1.  See help.')
end

if N<0 || K>N || K<0
   error('Arguments must be greater than zero, and K cannot exceed N.')
end

if N==0
    C = {{}}; % Easy case.
    return
end

if K==1
    if Cflg
        C{1} = {M};  % Easy case.
    else
        C{1} = {1:N};  % Easy case.
    end
    return
end

if Cflg
    NV = M;
else
    NV = 1:M; % Vector: base partition, RGF indexes into this guy.
end

NU = 1; % Number of unique indices in current partition.
stp1 = N; % Controls assigning of indices.
RGF = ones(1,N); % Holds the indexes.
BLD = cell(1,N); % Smaller cell array will be used in creating larger.

while cnt<=S
    idx1 = N; % Index into RGF.
    stp2 = RGF(idx1); % Works with stp1.

    while stp1(stp2)==1
        RGF(idx1) = 1;  % Assign value to RGF.
        idx1 = idx1 - 1; % Need to increment idx1 for translation below.
        stp2 = RGF(idx1); % And set this guy for stp1 assign below.
    end

    NU = NU + idx1 - N;  % Get provisional number of unique vals.
    stp1(1) = stp1(1) + N - idx1;

    if stp2==NU % Increment the number of unique elements.
        NU = NU +1;
        stp1(NU) = 0;
    end

    RGF(idx1) = stp2 + 1;  % Increment this position.
    stp1(stp2) = stp1(stp2) - 1;  % Translate indices of these two.
    stp1(stp2+1) = stp1(stp2+1) + 1; % Re-assign stopper.

    if NU==(~Kflg * NU + Kflg * K)
    % We could use:  C{cnt} = accumarray(RGF',NV',[],@(x) {x});   (SLOW!!)
    % or the next lines to C{cnt} = TMP;
        TMP = BLD(1:NU); % Build subcell of correct size.
        TMP{1} = NV(RGF==1); % These first two are always here.... no loop.
        TMP{2} = NV(RGF==2);

        for ii = 3:NU % Build the rest of cell array, if any.
            TMP{ii} = NV(RGF==ii);
        end

        C{cnt} = TMP;  % Assign cell at jj.
        cnt = cnt + 1;
    end
end

function S = stirling(N,K)
% Calculate the Stirling number of the second kind. Subfunc to partitions.

for ii = K:-1:0
    S(ii+1) = (-1)^ii * prod(1:K) / (prod(1:(K-ii)) *...
              (prod(1:ii)))*(K - ii)^N;  %#ok
end

S = 1/prod(1:K) * sum(S);


text:
function results = utl_run_batchjob(opts,d,appname,setnames)
% Internal: the actual processing function of bci_batchtrain.
% Results = ult_run_batchjob(Options,DatasetIndex,ApproachName,DatasetNames)
%
% This function executes one of the individual jobs generated by bci_batchtrain (application of one
% approach to one dataset).
%
% In:
%   Options : bci_batchtrain options struct
%
%   DatasetIndex : index into Options.datasets that determines which dataset to process
%
%   ApproachName : name of a field in Options.approaches that determines which approach to use
%
%   DatasetNames : names of the respective datasets for diagnostics
%
% Out:
%   Results : result data structure
%
% See also:
%   bci_batchtrain
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2011-08-19
dp;

% scaffolding to support custom error handlers (ultimately calls do_run to do the actual work)
results = [];
if strcmp(char(opts.handler),'rethrow')
    results = do_run(opts,d,appname,setnames);
else
    try
        results = do_run(opts,d,appname,setnames);
    catch e
        fprintf('Error processing data set "%s" with approach "%s".',setnames{d},appname);
        if ischar(opts.handler)
            try
                opts.handler = str2func(opts.handler);
            catch e
                error('The given error handler is malformed: %s',hlp_tostring(opts.handler,100));
            end
        end
        opts.handler(e);
    end
end

function results = do_run(opts,d,appname,setnames)
    % run a particular approach on a particular data set
    results = [];

    % determine where results would be saved
    storename = env_translatepath(strrep(strrep(opts.storepatt,'%set',setnames{d}),'%approach',appname));
    if opts.reuse && ~isempty(storename) && exist(storename,'file')

        % attempt to load existing results from disk
        fprintf('Reusing existing result for approach "%s" on set "%s".\n',appname,setnames{d});
        t0=tic; fprintf('Trying to load file %s...',storename);
        if opts.cache_loaded
            persistent result_cache; %#ok<TLEV>
            storetag = ['x' hlp_cryptohash(storename)];
            if isfield(result_cache,storetag)
                res = result_cache.(storetag);
            else
                io_load(storename);
                result_cache.(storetag) = res; %#ok<NODEF>
            end
        else
            io_load(storename);
        end
        fprintf('%.1f seconds (%.1f MB/s).\n',toc(t0),getfield(whos('res'),'bytes')/2^20/toc(t0));

    elseif ~opts.loadonly

        % train a model on the Dataset using bci_train
        [res.loss,res.model,res.stats] = bci_train(opts.trainargs{:}, 'data',opts.datasets{d}, 'approach',opts.approaches.(appname), 'markers',opts.markers);

        % optionally run bci_predict on each of the PredictSets
        if ~isempty(opts.predictsets) && ~isempty(opts.predictsets{d})
            for k=1:length(opts.predictsets{d})
                try
                    [res.pred_predictions{k},res.pred_loss(k),res.pred_stats{k},res.pred_targets{k}] = bci_predict(opts.predictargs{:},'data',opts.predictsets{d}{k},'model',res.model, 'markers',opts.markers);
                catch e
                    [res.pred_predictions{k},res.pred_loss(k),res.pred_stats{k},res.pred_targets{k}] = deal([],NaN,struct(),[]);
                    fprintf('Error computing predictions for set "%s", prediction set #%i with approach "%s".\n',setnames{d},k,appname);
                    opts.handler(e);
                end
            end
            res.pred_loss = res.pred_loss';
        end

        % save results
        if ~isempty(storename)
            t0=tic; fprintf('Saving results as %s...',storename);
            io_save(storename,opts.saveargs{:},'res');
            fprintf('%.1f seconds (%.1f MB/s).\n',toc(t0),getfield(whos('res'),'bytes')/2^20/toc(t0));
        end
    end

    % assign results to result data structure
    if ~isempty(opts.resultpatt)
        try
            statement = [strrep(strrep(opts.resultpatt,'%num',num2str(d)),'%approach',appname) ' res;'];
            eval(statement);
        catch e
            fprintf('Failed to evaluate ResultPattern (statement was "%s") with error: %s.\n',statement,e.message);
        end
    end

    % print diagnostics
    try
        if isfield(res,'pred_loss')
            fprintf('%s@%s: train loss: %.4f, prediction loss: %.4f\n',appname,setnames{d},res.loss,res.pred_loss);
        else
            fprintf('%s@%s: train loss: %.4f\n',appname,setnames{d},res.loss);
        end
    catch e
        fprintf('Failed to display loss estimates with error: %s\n',e.message);
    end

repo: REST


text:
function [w, infos] = hb(problem, in_options)
% Heavy Ball algorithm.
%
% Inputs:
%       problem     function (cost/grad/hess)
%       in_options  options
% Output:
%       w           solution of w
%       infos       information
%
% This file is part of GDLibrary and SGDLibrary.
%
% Created by H.Kasai on Oct. 28, 2020

    % set dimensions and samples
    d = problem.dim;
    n = problem.samples;

    % set local options
    local_options = [];
    local_options.step_alg = 'backtracking';
    local_options.beta = 0.01;

    % merge options
    options = mergeOptions(get_default_options(d), local_options);
    options = mergeOptions(options, in_options);

    % initialise
    iter = 0;
    grad_calc_count = 0;
    w = options.w_init;
    w_old = w;
    prev_step = options.step_init;

    % for stepsize
    if strcmp(options.step_alg, 'fix') || strcmp(options.step_alg, 'no_change')
        if isprop(problem, 'L')
            if problem.L > 0
                if isprop(problem, 'mu')
                    if problem.mu > 0
                        % This casse is L-smooth and mu-strongly convex.
                        cn = problem.L/problem.mu;
                        options.step_init = 4/(sqrt(problem.L)+sqrt(problem.mu))^2;
                        options.beta = ( (sqrt(cn)-1)/(sqrt(cn)+1) )^2;
                    else
                        % This casse is L-smooth
                        options.step_init = 1/problem.L;
                    end
                else
                    % This casse is L-smooth
                    options.step_init = 1/problem.L;
                end
            else
                options.step_alg = 'backtracking';
            end
        end
    end

    % initialize by BB step-size
    if strcmp(options.step_init_alg, 'bb_init')
        options.step_init = bb_init(problem, w);
    end

    % store first infos
    clear infos;
    [infos, f_val, optgap, grad, gnorm] = store_infos(problem, w, options, [], iter, grad_calc_count, 0);
    grad_old = grad;

    % display info
    if options.verbose
        fprintf('HB: Iter = %03d, cost = %.24e, gnorm = %.4e, optgap = %.4e\n', iter, f_val, gnorm, optgap);
    end

    % set start time
    start_time = tic();

    % main loop
    while (optgap > options.tol_optgap) && (gnorm > options.tol_gnorm) && (iter < options.max_epoch)

        options.iter = iter;
        [step, ~] = options.linesearchfun(options.step_alg, problem, w, w_old, grad, grad_old, prev_step, options);

        prev_step = step;
        w_w_old_diff = w - w_old;
        w_old = w;

        % update w
        w = w - step * grad + options.beta * w_w_old_diff;

        % proximal operator
        if ismethod(problem, 'prox')
            w = problem.prox(w, step);
        end

        % store gradient
        grad_old = grad;

        % measure elapsed time
        elapsed_time = toc(start_time);

        % count gradient evaluations
        grad_calc_count = grad_calc_count + n;

        % update iter
        iter = iter + 1;

        % store infos
        [infos, f_val, optgap, grad, gnorm] = store_infos(problem, w, options, infos, iter, grad_calc_count, elapsed_time);

        % display infos
        if options.verbose
            fprintf('HB: Iter = %03d, cost = %.24e, gnorm = %.4e, optgap = %.4e\n', iter, f_val, gnorm, optgap);
        end
    end

    if gnorm < options.tol_gnorm
        fprintf('Gradient norm tolerance reached: tol_gnorm = %g\n', options.tol_gnorm);
    elseif optgap < options.tol_optgap
        fprintf('Optimality gap tolerance reached: tol_optgap = %g\n', options.tol_optgap);
    elseif iter == options.max_epoch
        fprintf('Max iter reached: max_epoch = %g\n', options.max_epoch);
    end

end


text:
function ypp = fppcube ( x )

%*****************************************************************************80
%
%% FPPCUBE sets the value of the second derivative of the cubic function.
%
%  Discussion:
%
%    Y(X) = ( ( X + 2 ) * X + 3 ) * X + 4
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 February 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the argument.
%
%    Output, real YPP, the value of the second derivative of the cubic function.
%
  ypp = 6.0E+00 * x + 4.0E+00;

  return
end


text:
function result = ball_unit_14_3d ( func )

%*****************************************************************************80
%
%% BALL_UNIT_14_3D approximates an integral inside the unit ball in 3D.
%
%  Integration region:
%
%    Points (X,Y,Z) such that:
%
%      X**2 + Y**2 + Z**2 <= 1.
%
%  Discussion:
%
%    A 288 point 14-th degree formula is used, Stroud number S3:14-1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 November 2004
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Reference:
%
%    Arthur H Stroud,
%    Approximate Calculation of Multiple Integrals,
%    Prentice Hall, 1971.
%
%  Parameters:
%
%    Input, external FUNC, the name of the user supplied
%    function which evaluates F(X,Y,Z), of the form
%      function value = func ( x, y, z )
%
%    Output, real RESULT, the approximate integral of the function.
%
  r = [ 0.968160240E+00, 0.836031107E+00, 0.613371433E+00, 0.324253423E+00 ];
  weight = [ 0.076181268E+00, 0.126263673E+00, 0.098048133E+00, 0.032840260E+00 ];
  xtab = [ ...
    -0.151108275E+00, 0.315838353E+00, 0.346307112E+00, ...
    -0.101808787E+00, -0.409228403E+00 ];
  ytab = [ ...
    0.155240600E+00, 0.257049387E+00, 0.666277790E+00, ...
    0.817386065E+00, 0.501547712E+00 ];
  ztab = [ ...
    0.976251323E+00, 0.913330032E+00, 0.660412970E+00, ...
    0.567022920E+00, 0.762221757E+00 ];

  quad = 0.0;

  for m = 1 : 4

    w1 = 125.0 * weight(m) / 3360.0;
    x = 0.525731112 * r(m);
    y = 0.850650808 * r(m);
    z = 0.0;

    for j = 1 : 2
      x = -x;
      for k = 1 : 2
        y = -y;
        for l = 1 : 3
          [ x, y, z ] = r8_swap3 ( x, y, z );
          quad = quad + w1 * feval ( func, x, y, z );
        end
      end
    end

    w2 = 143.0 * weight(m) / 3360.0;

    for n = 1 : 5

      x = xtab(n) * r(m);
      y = ytab(n) * r(m);
      z = ztab(n) * r(m);

      for i = 1 : 3

        temp = x;
        x = z;
        z = -y;
        y = -temp;

        for j = 1 : 3

          [ x, y, z ] = r8_swap3 ( x, y, z );

          quad = quad + w2 * feval ( func, x, y, z );

        end

        y = -y;
        z = -z;
        quad = quad + w2 * feval ( func, x, y, z );

      end

    end

  end

  volume = ball_unit_volume_3d  ( );
  result = quad * volume;

  return
end


text:
function FD = FermiDiract(z,nu)
%% Fermi-Diract function
% Implementation of the Fermi-Diract Statistical function
% in latex words:
%
% $$F_\nu(z)=\frac{1}{\Gamma(\nu)} \int_{0}^{\infty}\frac{x^{\nu-1}}{z^{-1}e^x+1}
%   \approx \sum_{l=1}^{\infty}(-1)^{l-1}\frac{z^l}{l^\nu}$$
%
% As we can notice this function is of the order $\nu$

l   = 1:50; % up to l = 50 to ensure a fair accuaracy
fd  = (-1).^(l-1) .* (z.^l) ./ (l.^nu);
FD  = sum(fd);
return

repo: Aero-matlab


text:
function [ know, x ] = p19_sol ( n )

%*****************************************************************************80
%
%% P19_SOL returns the solution for problem 19.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 March 2000
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the problem.  This value
%    is only needed for those problems with variable N.
%
%    Output, integer KNOW.
%    If KNOW is 0, then the solution is not known.
%    If KNOW is positive, then the solution is known, and is returned in X.
%
%    Output, real X(N), the solution, if known.
%
  know = 1;

  x = [ 1.0, 1.0 ]';

  return
end


text:
% T in K
% es in Pa
function es = calculate_saturation_vapor_pressure_liquid (T, method)
    if (nargin < 2) || isempty(method),  method = 'Murphy&Koop2005';  end
    switch method
    case {'Bolton1980', 'approx'}
        t = T - 273.15;
        es2 = 0.6112 .* exp(17.67 .* t ./ (t + 243.5));
        es = es2 * 1000;
        % The AMS Glossary gives [1] the simplified formula above,
        % accordingly to Bolton (1980, eq. 10).
        %
        % AMS Glossary. <http://amsglossary.allenpress.com/glossary/search?id=clausius-clapeyron-equation1>
        % Bolton, D., 1980: The computation of equivalent potential temperature. Mon. Wea. Rev., 108, 1046-1053. <http://dx.doi.org/10.1175/1520-0493(1980)108<1046:TCOEPT>2.0.CO;2>
    case 'Murphy&Koop2005'
        if ~all(123 < T & T < 332)
            warning('calculate_saturation_vapor_pressure_liquid:outRange', ...
                'Temperature out of range [123-332] K.');
        end
        temp = 54.842763 - 6763.22 ./ T - 4.210 .* log(T) + 0.000367 .* T ...
          + tanh( 0.0415 * (T - 218.8) ) ...
          .* (53.878 - 1331.22 ./ T - 9.44523 .* log(T) + 0.014025 .* T);
        es = exp(temp);
        % D. M. MURPHY and T. KOOP
        % Review of the vapour pressures of ice and supercooled water for
        % atmospheric applications
        % Q. J. R. Meteorol. Soc. (2005), 131, pp. 1539-1565
        % doi: 10.1256/qj.04.94
        % <http://dx.doi.org/10.1256/qj.04.94>
        %
        % "Widely used expressions for water vapour (Goff and Gratch 1946; Hyland and Wexler 1983) are being applied outside the ranges of data used by the original authors for their fits. This work may be the first time that data on the molar heat capacity of supercooled water [i.e., at temperatures below its freezing temperature] have been used to constrain its vapour pressure.
    otherwise
        error('calculate_saturation_vapor_pressure_liquid:methodUnknown', ...
            'Method "%s" unknown.', method);
    end
end

%!test
%! % Table C1. VALUES RECOMMENDED FOR CHECKING COMPUTER CODES
%! temp = [...
%! % Temperature (K), Liquid vapor pressure (Pa)
%! 150     1.562e-5    0.001e-5
%! 180     0.011239    0.000001
%! 210     1.2335      0.0001
%! 240     37.667      00.001
%! 273.15  611.213     000.001
%! 273.16  611.657     000.001
%! 300     3536.8      0000.1
%! ];
%! T = temp(:,1);
%! es = temp(:,2);
%! tol = temp(:,3)/10*5;  % rounding error.
%! es2 = calculate_saturation_vapor_pressure_liquid (T);
%! %[es, es2, es2-es]  % DEBUG
%! myassert(es2, es, -tol);

%!test
%! T = (123:10:332)';
%! es  = calculate_saturation_vapor_pressure_liquid (T, 'Murphy&Koop2005');
%! es2 = calculate_saturation_vapor_pressure_liquid (T, 'Bolton1980');
%! %[T-273.15, es, es2, es2-es, 100*(es2-es)./es]  % DEBUG
%! % worse agreement (in percentage) for supercooled temperatures, as expected.
%! temp = 100*(es2-es)./es;
%! myassert(abs(temp(1)) > abs(temp(end)))


text:
function cmdString = tapas_physio_create_spm_toolbox_link(pathPhysIO)
% Creates a symbolik link of PhysIO/code folder to subfolder SPM/toolbox/PhysIO
% to make toolbox visible to SPM Batch editor
%
%   tapas_physio_create_spm_toolbox_link()
%
% IN
%
% OUT
%   cmdString   string of executed command
%
% EXAMPLE
%   tapas_physio_create_spm_toolbox_link
%
%   See also

% Author: Lars Kasper
% Created: 2018-02-17
% Copyright (C) 2018 TNU, Institute for Biomedical Engineering,
%                    University of Zurich and ETH Zurich.
%
% This file is part of the TAPAS PhysIO Toolbox, which is released under the terms of the GNU General Public
% License (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.
%
cmdString = '';

if nargin < 1
    % 2x fileparts, because parent folder
    pathPhysIO = fileparts(fileparts(mfilename('fullpath')));
end

pathSpm = fileparts(which('spm'));

if isempty(pathSpm)
    warning('SPM folder not found. Could not create symbolic link to PhysIO Toolbox');
else
    pathLinkPhysIOSPM = fullfile(pathSpm, 'toolbox', 'PhysIO');
    if ~exist(pathLinkPhysIOSPM, 'dir')
        % Create Link or hard-copy folder, OS dependent
        if ispc
            % unfortunately, system link does not work for SPM, has to be
            % hard copy
            fprintf('Copying %s to %s, because symlink not sufficient on Windows...\n', pathPhysIO, pathLinkPhysIOSPM);
            cmdString = sprintf('copyfile(''%s'', ''%s'');',pathPhysIO, pathLinkPhysIOSPM);
            eval(cmdString);
            % linking indeed the other way around than in Linux/Mac
            % cmdString = sprintf('mklink /D %s %s', pathLinkPhysIOSPM, pathPhysIO);
        else %unix/Mac
            cmdString = sprintf('ln -s ''%s'' ''%s''', pathPhysIO, pathLinkPhysIOSPM);
            system(cmdString);
        end
    else
        warning('Destination spm/toolbox folder %s already found. Will not overwrite or re-link', pathLinkPhysIOSPM);
    end
end

author: translationalneuromodeling

repo: tapas


text:
function model = rmSliceSet(model,tmp,slice)
% rmSliceSet - put slice info into model struct
% we convert back to double precision here.
%
% model = rmSliceSet(model,tmp,slice);
%
% 2008/01 SOD: extracted from rmGridFit.

% loop over models
for n=1:numel(model),
    % variables may have slightly different names
    ftmp   = {'x0','y0','s','x02','y02','s2','s_major','s_minor','s_theta','rss','rss2','rsspos','rssneg','rawrss','rawrss2', 'exponent'};

    % now get values from model and put in new slice values
    for fn = 1:numel(ftmp),
        % check whether data exists and has data in tmp structure
        if isfield(tmp{n},ftmp{fn}) && ~isempty(tmp{n}.(ftmp{fn}))
            % get data from model structure
            val = rmGet(model{n},ftmp{fn});
            % if no data in model structure make it (zeros)
            if isempty(val)
                val = zeros(size(rmGet(model{n},'x0')));
            end
            % now put data in model structure format
            switch length(size(val))
                case 3 % presumably means INPLANE view, where slice is the third dimension
                    data = reshape(double(tmp{n}.(ftmp{fn})), size(val,1), size(val,2));
                    val(:,:,slice) = data;
                otherwise
                    if size(val,2) == size(tmp{n}.(ftmp{fn}),2)
                        val(slice,:)  = double(tmp{n}.(ftmp{fn}));
                    elseif slice == 1 && size(val,2) ~= size(tmp{n}.(ftmp{fn}),2)
                        val = double(tmp{n}.(ftmp{fn}));
                    end
            end
            % save data in model
            model{n}      = rmSet(model{n},ftmp{fn},val);
        end
    end

    % other params
    if isfield(tmp{n},'df')
        model{n}       = rmSet(model{n},'dfglm',double(tmp{n}.df));
    end
    if isfield(tmp{n},'desc')
        model{n}       = rmSet(model{n},'desc',tmp{n}.desc);
    end

    % distribute beta values
    val = rmGet(model{n},'b');
    switch length(size(val))
        % switch on the number of dimensions, since inplane data has diff
        % dimesnsionality than other views
        case 4 % 4 dimensions means Inplane model:
            %  3 dimensions of coords, and one dimension of beta values.
            %  Hence the dims are x, y, slice, betas.
            nbetas = size(val,4);
            %tmp{n}.b = zeros(nbetas,nvoxelsPerSlice,'single');
            for fn = 1:nbetas,
                data = reshape(double(tmp{n}.b(fn,:)), size(val,1), size(val,2));
                val(:,:,slice,fn) = data;
            end;

        otherwise
            val = val(:,:,1:size(tmp{n}.b,1));
            for fn = 1:size(val,3),
                val(slice,:,fn) = double(tmp{n}.b(fn,:));
            end;
    end
    model{n} = rmSet(model{n},'b',val);
end;
return;
%-----------------------------------


text:
function optimal_set = GMM_spike_sorter(SpikeMat,sdd,REM,INPCA)
% this function sorts detected spikes based on GMM clustering method
% SpikeMat is matrix of detected spikes, sdd is settings, REM contains
% spikes after statistical filtering. INPCA is a logical value which
% determines whether considering noise spikes in computing PCA or not.

% default values for REM and INPCA
if nargin < 3
    REM = [];
    INPCA = true;
end

if nargin < 4
    INPCA = true;
end

% removing REM if it is given in wrong format
if ~islogical(REM) || length(REM) ~= size(SpikeMat,1)
    REM = [];
end

seed = sdd.sort.random_seed;

g_max = sdd.sort.g_max;
g_min = sdd.sort.g_min;

error = sdd.sort.error; % the termination tolerance for the loglikelihood function value.

% removing outliers using given REM
if ~INPCA && ~isempty(REM)
    SpikeMat(REM,:) = [];
end

if isempty(SpikeMat)
    optimal_set = [];
    return
end

[~,SpikeMat] = pca(SpikeMat,'NumComponents',sdd.sort.n_pca);

if ~isempty(REM) && INPCA
    SpikeMat(REM,:) = [];
end

max_iter = sdd.sort.max_iter;

% simple clustering method to determine centers
rng(seed)
options = statset('MaxIter',max_iter, 'TolFun', error);
myfunc = @(X,K)(cluster(fitgmdist(X, K, 'Replicates',5, 'Options', options), X));
eva = evalclusters(SpikeMat,myfunc, 'CalinskiHarabasz', 'KList',[g_min:g_max]);
BestModel = fitgmdist(SpikeMat,eva.OptimalK,'Options',options,'Replicates', 5);
if isempty(REM)
    optimal_set.cluster_index = cluster(BestModel, SpikeMat);
else
    optimal_set.cluster_index(~REM) = cluster(BestModel, SpikeMat);
    optimal_set.cluster_index(REM) = 255; % removed
end

author: ramintoosi

repo: ROSS

save_path: github-repos/MATLAB/ramintoosi-ROSS


text:
function [meanAmp,meanPh,seZ,meanStd] = vectorMean(view,scanNum,ROIcoords)
%
% [meanAmp,meanPh,seZ,meanStd] = vectorMean(view,scanNum,ROIcoords)
%
% Calculates mean amplitude and phase for pixels that are in
% ROIcoords.
% The standard error (std(z)/sqrt(length(z))) can also
% returned.  This is the average distance of the complex values
% z = amp*exp(-i(ph)) from the mean of z.
%
% Can >also< return the mean noise std in the ROI:
% Computed from the coherence and the amplitude.
%
% scanNum: scan number (integer)
% ROIcoords: 3xN array of (y,x,z) coords (e.g., corresponding to
%   the selected ROI).
%
% djh 4/23/98
% bw  2/17/99 Added seZ computation.
%
% Get co and ph (vectors) for the desired scan, within the
% current ROI.
%
subAmp = getCurDataROI(view,'amp',scanNum,ROIcoords);
subPh = getCurDataROI(view,'ph',scanNum,ROIcoords);
subCo = getCurDataROI(view,'co',scanNum,ROIcoords);

% Remove NaNs from subCo and subAmp that may be there if ROI
% includes volume voxels where there is no data.
NaNs = find(isnan(subPh));
if ~isempty(NaNs)
  disp('ROI includes voxels that have no data.  These voxels are being ignored.');
  notNaNs = find(~isnan(subPh));
  subPh = subPh(notNaNs);
  subAmp = subAmp(notNaNs);
  subCo= subCo(notNaNs);

end

% Compute the mean co right here...
meanCo=mean(subCo);

% convert to complex numbers
z = subAmp .* exp(sqrt(-1)*subPh);
if isempty(z)
   disp('Warning: no activity seen for current ROI');
   seZ = 0;
   meanAmp = 0;
   meanPh = 0;
else
	meanZ = mean(z);

    if nargout > 2
        seZ   = std(z)/sqrt(length(z));
    end
    meanAmp = abs(meanZ);
	meanPh = angle(meanZ);
end
% Compute the meanStd right here...
meanStd=meanAmp.*sqrt((1/mean(subCo).^2)-1);

return;

% Debug

coords = INPLANE{1}.ROIs(INPLANE{1}.selectedROI).coords;
scan = 1;
[meanAmp,meanPh] = vectorMean(INPLANE{1},scan,coords)


text:
%                               SYMBADA = getsymbada(At,Ajc,DAt,psdblkstart)
% GETSYMBADA
%   Ajc points to start of PSD-nonzeros per column
%   DAt.q has the nz-structure of ddotA.
%
% ******************** INTERNAL FUNCTION OF SEDUMI ********************
%
% See also sedumi, partitA, getada1, getada2.

function SYMBADA = getsymbada(At,Ablkjc,DAt,psdblkstart)
%
% This file is part of SeDuMi 1.1 by Imre Polik and Oleksandr Romanko
% Copyright (C) 2005 McMaster University, Hamilton, CANADA  (since 1.1)
%
% Copyright (C) 2001 Jos F. Sturm (up to 1.05R5)
%   Dept. Econometrics & O.R., Tilburg University, the Netherlands.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% Affiliation SeDuMi 1.03 and 1.04Beta (2000):
%   Dept. Quantitative Economics, Maastricht University, the Netherlands.
%
% Affiliations up to SeDuMi 1.02 (AUG1998):
%   CRL, McMaster University, Canada.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% 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
%

Alpq = spones(extractA(At,Ablkjc,0,3,1,psdblkstart(1)));
Ablks = findblks(At,Ablkjc,3,[],psdblkstart);
if spars(Ablks)==1 || spars(Alpq)==1 || (~isempty(DAt.q) && spars(DAt.q)==1)
    SYMBADA=sparse(ones(size(At,2),size(At,2)));
else
    SYMBADA = Alpq' * Alpq + Ablks'*Ablks + DAt.q'*DAt.q;
end


text:
function value = r8_log_10 ( x )

%*****************************************************************************80
%
%% R8_LOG_10 returns the logarithm base 10 of |X|.
%
%  Discussion:
%
%    value = Log10 ( |X| )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 June 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the number whose base 2 logarithm is desired.
%    X should not be 0.
%
%    Output, real VALUE, the logarithm base 10 of the absolute
%    value of X.  It should be true that |X| = 10**R8_LOG_10.
%
  if ( x == 0.0 )
    value = -inf;
  else
    value = log10 ( abs ( x ) );
  end

  return
end


text:
function plotQuadCon(Q,l,rl,ru,data)
%PLOTQUADCON Plot Quadratic Constraints on the current figure
%   plotQuadCon(Q,l,r)

%   Copyright (C) 2013 Jonathan Currie (IPL)

xl = data.xl; yl = data.yl;

hold on;

%Colour
dkg = [0.2 0.2 0.2];

%Determine number of quad constraints
if(iscell(Q))
    no = length(Q);
else
    no = 1;
end
%>= 2014b requires more points...
if(~verLessThan('matlab','8.4'))
    nmul = 2;
else
    nmul = 1;
end

%Generate Constraint Surface Points
[x1,x2] = meshgrid(linspace(xl(1),xl(2),data.npts*nmul),linspace(yl(1),yl(2),data.npts*nmul));

%For each quadratic constraint, plot
for i = 1:no
    %Get Constraint Variables
    if(iscell(Q))
        nQ = Q{i}; nl = l(:,i); nrl = rl(i); nru = ru(i);
    else
        nQ = Q; nl = l; nrl = rl; nru = ru;
    end
    %Form Constraint Function
    if(data.ndec==1)
        con = @(x) x(1)*nQ*x(1) + nl*x(1);
    else
        con = @(x) x.'*nQ*x + nl.'*x;
    end
    %Plot Each Quad Con as General Nonlinear Constraint
    plotNLCon(con,[],nrl,nru,x1,x2,dkg,data);
end
hold off;

%OLD CODE
% %Plot Quadratic Inequality Constraints (Inefficient.. need to solve the quadratic)
% [x1,x2] = meshgrid(linspace(xl(1),xl(2),npts),linspace(yl(1),yl(2),npts));
% nox = size(x1);
% noy = size(x2);
% obj = zeros(nox(1),noy(2));
% if(iscell(Q))
%     no = length(Q);
% else
%     no = 1;
% end
% for i = 1:no
%     %get vars
%     if(iscell(Q))
%         nQ = Q{i}; nl = l(:,i); nrl = rl(i); nru = ru(i);
%     else
%         nQ = Q; nl = l; nrl = rl; nru = ru;
%     end
%     eq = false; sense = 'L';
%     % check for double constraint or equality
%     if(~isinf(nrl) && ~isinf(nru))
%         if(nrl == nru)
%             nr = nru;
%             eq = true;
%         else
%             nr = [nrl nru];
%             sense = 'GL';
%         end
%     elseif(~isinf(nrl))
%         nr = nrl;
%     else
%         nr = nru;
%     end
%     % create surface
%     for n = 1:nox(1)
%         for m = 1:noy(2)
%             x = [x1(n,m) x2(n,m)]';
%             obj(n,m) = x.'*nQ*x + nl.'*x;
%         end
%     end
%     if(eq)
%         contour(x1,x2,obj,'color',[0 0 1],'levellist',nr);
%     else
%         for j = 1:length(nr) %each row constraint (max 2)
%             c = contour(x1,x2,obj,'color',dkg,'levellist',nr(j));
%             %Plot Hatch
%             if(~isempty(c))
%                 %See if we have multiple contours (non-convex or sd)
%                 len = size(c,2)-1;
%                 if(c(2,1) ~= len)
%                     %Build contour array
%                     cstrt = 2; cend = []; n = 2; ind = 1;
%                     while(ind <= len)
%                         ind = ind + c(2,ind) + 1;
%                         cend(n-1) = ind-1; %#ok<AGROW>
%                         cstrt(n) = ind+1;  %#ok<AGROW>
%                         n = n + 1;
%                     end
%                 else
%                     cstrt = 2;
%                     cend = len;
%                 end
%                 %Plot each contour hatch
%                 for n = 1:length(cend)
%                     %Get contour vectors
%                     vecx = diff(c(1,cstrt(n):cend(n)));
%                     vecy = diff(c(2,cstrt(n):cend(n)));
%                     if(isempty(vecx) || isempty(vecy))
%                         continue;
%                     end
%                     %Rotate hatch lines based on infeasible region
%                     xt = [c(1,cstrt(n))+vecy(1) c(2,cstrt(n))-vecx(1)]'; %check rotated -90
%                     fval = xt.'*nQ*xt + nl.'*xt ;
%                     if(fval <= nr(j) || sense(j) == 'G') %rotate 90
%                         hvecx = -vecy;
%                         hvecy = vecx;
%                     else %rotate -90
%                         hvecx = vecy;
%                         hvecy = -vecx;
%                     end
%                     %Normalize
%                     av = mean(sqrt(hvecx.^2 + hvecy.^2));
%                     dirs = atan2(hvecy,hvecx);
%                     hvecx = av*cos(dirs);
%                     hvecy = av*sin(dirs);
%                     %Shift origin
%                     hvecx = c(1,cstrt(n):cend(n)-1) + hvecx;
%                     hvecy = c(2,cstrt(n):cend(n)-1) + hvecy;
%                     %Plot
%                     line([c(1,cstrt(n):cend(n)-1)' hvecx']',[c(2,cstrt(n):cend(n)-1)' hvecy']','Color','k')
%                 end
%             else
%                 optiwarn('opti:plot','Cannot plot inequality constraint as contour data is empty!');
%             end
%         end
%     end
% end
%
% hold off;

% CONVEX QC CODE
% %Get contour vectors
% vecx = diff(c(1,2:end));
% vecy = diff(c(2,2:end));
% %Rotate hatch lines based on infeasible region
% xt = [c(1,2)+vecy(1) c(2,2)-vecx(1)]'; %check rotated -90
% fval = xt.'*nQ*xt + nl.'*xt ;
% if(fval <= nr) %rotate 90
%     hvecx = -vecy;
%     hvecy = vecx;
% else %rotate -90
%     hvecx = vecy;
%     hvecy = -vecx;
% end
% %Normalize
% av = mean(sqrt(hvecx.^2 + hvecy.^2));
% dirs = atan2(hvecy,hvecx);
% hvecx = av*cos(dirs);
% hvecy = av*sin(dirs);
% %Shift origin
% hvecx = c(1,2:end-1) + hvecx;
% hvecy = c(2,2:end-1) + hvecy;
% %Plot
% line([c(1,2:end-1)' hvecx']',[c(2,2:end-1)' hvecy']','Color','k')

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [ x, seed ] = normal_truncated_ab_sample ( mu, s, a, b, seed )

%*****************************************************************************80
%
%% NORMAL_TRUNCATED_AB_SAMPLE samples the truncated Normal PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 August 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real MU, S, the mean and standard deviation of the
%    parent Normal distribution.
%
%    Input, real A, B, the lower and upper truncation limits.
%
%    Input/output, integer SEED, a seed for the random number
%    generator.
%
%    Output, real X, a sample of the PDF.
%
  alpha = ( a - mu ) / s;
  beta = ( b - mu ) / s;

  alpha_cdf = normal_01_cdf ( alpha );
  beta_cdf = normal_01_cdf ( beta );

  [ u, seed ] = r8_uniform_01 ( seed );
  xi_cdf = alpha_cdf + u * ( beta_cdf - alpha_cdf );
  xi = normal_01_cdf_inv ( xi_cdf );

  x = mu + s * xi;

  return
end


text: %MAKEFIGS_ANALYTIC has been moved to JLAB_MAKEFIGS.

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
base64 = org.apache.commons.codec.binary.Base64;
tsv_file = 'MsCelebV1-Faces-Cropped.Samples.tsv';

root_folder = 'data';
mkdir(root_folder);

tsv_fid = fopen(tsv_file, 'r');
total_lines = 3481187;
image_id = 1;
while ~feof(tsv_fid)
    line = fgetl(tsv_fid);%MID,EntityNameString,ImageURL,FaceID,FaceRectangle_Base64Encoded,FaceData_Base64Encoded
    C = strsplit(line,'\t');
    folder = fullfile(root_folder,C{1});
    if exist(folder,'dir')==0
        mkdir(folder);
    end;
    filename = fullfile(folder,[C{2} '-' C{5} '.jpg']);
    raw  = base64decode(C{7});
%     img_fid = fopen(filename, 'wb');
%     fwrite(img_fid, raw, 'uint8');
%     fclose(img_fid);
    jImg = javax.imageio.ImageIO.read(java.io.ByteArrayInputStream(raw));
    h = jImg.getHeight;
    w = jImg.getWidth;
    p = typecast(jImg.getData.getDataStorage, 'uint8');
    img = permute(reshape(p, [3 w h]), [3 2 1]);
    img = img(:,:,[3 2 1]);
    imshow(img);
    disp([num2str(image_id) '/' num2str(total_lines) ' ' filename]);
    image_id = image_id + 1;
end;
fclose(tsv_fid);

author: happynear

repo: FaceVerification


text:
function [P_O] = getColourMap(patch, bg_hist, fg_hist, n_bins, grayscale_sequence)
%% GETCOLOURMAP computes pixel-wise probabilities (PwP) given PATCH and models BG_HIST and FG_HIST
    % check whether the patch has 3 channels
    [h, w, d] = size(patch);
    % figure out which bin each pixel falls into
    bin_width = 256/n_bins;
    % convert image to d channels array
    patch_array = reshape(double(patch), w*h, d);
    % to which bin each pixel (for all d channels) belongs to
    bin_indices = floor(patch_array/bin_width) + 1;
    % Get pixel-wise posteriors (PwP)
    P_bg = getP(bg_hist, h, w, bin_indices, grayscale_sequence);
    P_fg = getP(fg_hist, h, w, bin_indices, grayscale_sequence);

    % Object-likelihood map
    P_O = P_fg ./ (P_fg + P_bg);
end

author: bertinetto

repo: staple

save_path: github-repos/MATLAB/bertinetto-staple


text:
%this is a simple example, demonstrating delay estimation. Estimates the
%delay of yn1 with respect to yn0 in the presence of noise.

N_p=256;%number of points;

d=5*rand;%delay

y0=randn(1,N_p);%reference signal

y1=fft_circshift(y0,d);%perform circular shift

yn0=0.05*randn(1,N_p);%noise
yn1=0.05*randn(1,N_p);%noise

y0n=y0+yn0;%add noise to signals
y1n=y1+yn1;

alpha=-1;%tuning parameter for modified gaussian method
%(only used in the case of negative xc if alpha is negative, otherwise
%always use it. See delayest_3point.m)

%estimate delays
d_hat_parabola=delayest_3point(y1n,y0n,'parabola','xc');
d_hat_Gaussian=real(delayest_3point(y1n,y0n,'Gaussian','xc'));
d_hat_modGaussian=delayest_3point(y1n,y0n,'modGaussian','xc',alpha);
d_hat_cosine=real(delayest_3point(y1n,y0n,'cosine','xc'));
d_hat_phase=delayest_fft(y1n,y0n);
d_hat_iterative=delayest_iterative(y1n,y0n);

fprintf('Delay=%f samples\n',d);
fprintf('Estimates (samples):\nparabola %f\nGaussian %f\nmodified Gaussian %f\ncosine %f\nphase %f\niterative %f\n',...
    d_hat_parabola,d_hat_Gaussian,d_hat_modGaussian,d_hat_cosine,d_hat_phase,d_hat_iterative);

figure(1)
plot([y0n;y1n]','.-')


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

g = zeros(size(z));

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

g = sigmoid(z) .* (1-sigmoid(z));

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

end

author: AvaisP


text:
clc;clear all;close all;
%***********************************************%
% This code runs on the Market-1501 dataset.    %
% Please modify the path to your own folder.    %
% We use the mAP and hit-1 rate as evaluation   %
%***********************************************%
% if you find this code useful in your research, please kindly cite our
% paper as,
% Liang Zheng, Liyue Sheng, Lu Tian, Shengjin Wang, Jingdong Wang, and Qi Tian,
% Scalable Person Re-identification: A Benchmark, ICCV, 2015.

% Please download Market-1501 dataset and unzip it in the "dataset" folder.
addpath(genpath('LOMO_XQDA/'));
addpath(genpath('utils/'));
addpath(genpath('KISSME/'));
run('KISSME/toolbox/init.m');

%% network name
netname = 'ResNet_50'; % network: CaffeNet  or ResNet_50

%% test info
galFea = importdata(['feat/IDE_' netname '_test.mat']);
galFea = double(galFea);
probFea = importdata(['feat/IDE_' netname '_query.mat']);
probFea = double(probFea);
label_gallery = importdata('data/testID.mat');
label_query = importdata('data/queryID.mat');
cam_gallery =   importdata('data/testCam.mat');
cam_query =  importdata('data/queryCam.mat');

%% normalize
sum_val = sqrt(sum(galFea.^2));
for n = 1:size(galFea, 1)
    galFea(n, :) = galFea(n, :)./sum_val;
end

sum_val = sqrt(sum(probFea.^2));
for n = 1:size(probFea, 1)
    probFea(n, :) = probFea(n, :)./sum_val;
end

% Instead of pdist2. 50 times faster than pdist2
my_pdist2 = @(A, B) sqrt( bsxfun(@plus, sum(A.^2, 2), sum(B.^2, 2)') - 2*(A*B'));
%% Euclidean

dist_eu = my_pdist2(galFea', probFea');
[CMC_eu, map_eu, ~, ~] = evaluation(dist_eu, label_gallery, label_query, cam_gallery, cam_query);

fprintf(['The IDE (' netname ') + Euclidean performance:\n']);
fprintf(' Rank1,  mAP\n');
fprintf('%5.2f%%, %5.2f%%\n\n', CMC_eu(1) * 100, map_eu(1)*100);

repo: Learning-via-Translation


text:
function a = setvalues(a, v)

%tstoolbox/@achse/setvalues
%   Syntax:
%     * a = setvalues(a, v)
%
% Copyright 1997-2001 DPI Goettingen, License http://www.physik3.gwdg.de/tstool/gpl.txt

a.values = v;


text:
function [prior,binSize] = part_location_prior(expidx)

p = rcnn_exp_params(expidx);
exp_dir = [p.expDir '/' p.shortName];

conf = rcnn_config('sub_dir', '/cachedir/test', 'exp_dir', exp_dir);
saveTo = [conf.cache_dir '/gtHist/'];
if (~exist(saveTo,'dir'))
    mkdir(saveTo);
end

% sc = p.refHeight/200;
sc = 2;
binSize = p.prior_bin_size*sc;
delta = 300*sc;
nBins = round(2*delta/binSize);

try
    assert(false);
%     load(fname, 'prior');
catch
%     load(p.testGT);
    load(p.trainGT);
%     if (exist('single_person_annolist','var'))
%         annolist = single_person_annolist;
%     end

%     fnameDist = [fileparts(p.evalTest) '/distAll'];
%     load(fnameDist, 'keypointsAll');
%     assert(length(keypointsAll) == length(annolist));

    [~,parts] = util_get_parts24();
    nParts = length(p.pidxs);
    prior = zeros(2*delta,2*delta,nParts);
    histLoc = zeros(nBins,nBins,nParts);
    nEx = zeros(nParts,1);
    X0 = size(prior,1)/2;
    Y0 = size(prior,2)/2;

    for imgidx = 1:length(annolist)
        fprintf('.');
        rect = annolist(imgidx).annorect;
        for ridx = 1:length(rect)
            objpos = rect(ridx).objpos;
            points = rect(ridx).annopoints.point;
            for pidx = 1:length(parts)
                pidxLin = find(p.pidxs == parts(pidx).id);
                if (~isempty(pidxLin))
                    assert(parts(pidx).pos(1) == parts(pidx).pos(2));
                    pp = util_get_annopoint_by_id(points,parts(pidx).pos(1));
%                     pp = keypointsAll(imgidx).det(parts(pidx).pos(1)+1,1:2);
                    if (~isempty(pp))
%                         iy = round(Y0 + sc*(pp(:,2) - objpos.y));
%                         ix = round(X0 + sc*(pp(:,1) - objpos.x));
                        iy = round(Y0 + sc*(pp.y - objpos.y));
                        ix = round(X0 + sc*(pp.x - objpos.x));
                        iy = max(iy,1);
                        iy = min(iy,2*delta);
                        ix = max(ix,1);
                        ix = min(ix,2*delta);
                        iy_hist = ceil(iy/binSize);
                        ix_hist = ceil(ix/binSize);
                        histLoc(iy_hist,ix_hist,pidxLin) = histLoc(iy_hist,ix_hist,pidxLin) + 1;
                        iy_prior = (iy_hist-1)*binSize+1:iy_hist*binSize;
                        ix_prior = (ix_hist-1)*binSize+1:ix_hist*binSize;
                        prior(iy_prior,ix_prior,pidxLin) = prior(iy_prior,ix_prior,pidxLin) + ones(binSize,binSize);
                        nEx(pidxLin) = nEx(pidxLin) + 1;
                    end
                end
            end
        end
        if (~mod(imgidx, 100))
            fprintf(' %d/%d\n',imgidx,length(annolist));
        end
    end
    fprintf(' done\n');
    if (isfield(p,'sigma_gauss'))
        sigma_gauss = p.sigma_gauss;
    else
        sigma_gauss = binSize;
    end
    for pidx = 1:nParts
%         prior(:,:,pidx) = prior(:,:,pidx)/nEx(pidx);
        prior(:,:,pidx) = prior(:,:,pidx)/max(max(prior(:,:,pidx)));
        pr = prior(:,:,pidx);
        pr(pr > 0) = 1;
        filterMask = gausswin(2*sigma_gauss+1)*gausswin(2*sigma_gauss+1)';
        prior(:,:,pidx) = filter2(filterMask, pr);%/nEx(pidx);
        prior(:,:,pidx) = prior(:,:,pidx)/max(max(prior(:,:,pidx)));
    end
%     save(fname, 'prior');
end

labels = {'rankle','rknee','rhip','lhip','lknee','lankle','rwrist','relbow','rshoulder','lshoulder','lelbow','lwrist','neck','tophead'};
tics = [1 size(prior,1)/6:size(prior,1)/6:size(prior,1)];
for pidx = 1:length(parts)
    figure(101);clf;%colormap gray;
    pidxLin = find(p.pidxs == parts(pidx).id);
    if (~isempty(pidxLin))
        pr = prior(:,:,pidxLin);
        imagesc(log(log(pr+1))); hold on; axis equal;
        plot([size(prior,2)/2; size(prior,2)/2], [1; size(prior,1)]', 'k-', 'lineWidth', 2);
        plot([1; size(prior,2)], [size(prior,1)/2; size(prior,1)/2]', 'k-', 'lineWidth', 2);
        title([labels{pidxLin}]);
        set(gca,'XTick',tics,'YTick',tics);
        fImgName = [saveTo '/partDetHist_' labels{pidxLin}];
        print(gcf, '-dpng', [fImgName '.png']);
    end
end
end

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
function [status]=mrLoadRetAnat2SPM(fileRoot)
%function [status]=mrLoadRetAnat2SPM(fileRoot)

global INPLANE;

if (~exist('INPLANE','var'))
    error('INPLANE Structure must exist');
end
if (isempty(INPLANE.anat))
    error('INPLANE.anat must contain data: select view->croppedAnatomy in mrLoadRet to get the anatomies up');
end
if (~exist('fileRoot','var'))
    error('You must enter a file root: e.g. //gwyn/u1/mri/sampleSet/012345/vSPM_test');
end

outDat=INPLANE.anat;

% This is a script to dump out the current INPLANE.anat data in Analyse format so that you can load it in to SPM
% We make some effort to ensure that it conforms to the expected orientation of SPM data
% In the INPLANE structure it's stored as x*y*nSlices
% Now we can't guarantee that we're going to get the orientation right (becasue we're not yet at the stage where we can align funcitonal data to the
% SPM template) but we can try to at least get the nSlices direction right so that the 'coronal' views more or less line up for mrLoadRet and SPM

% Now, just call writeAnatFileSet
fileVector=writeAnatFileSet(outDat,fileRoot);

repo: vistasoft


text:
function [net, gamma, logev] = evidence(net, x, t, num)
%EVIDENCE Re-estimate hyperparameters using evidence approximation.
%
%	Description
%	[NET] = EVIDENCE(NET, X, T) re-estimates the hyperparameters ALPHA
%	and BETA by applying Bayesian re-estimation formulae for NUM
%	iterations. The hyperparameter ALPHA can be a simple scalar
%	associated with an isotropic prior on the weights, or can be a vector
%	in which each component is associated with a group of weights as
%	defined by the INDEX matrix in the NET data structure. These more
%	complex priors can be set up for an MLP using MLPPRIOR. Initial
%	values for the iterative re-estimation are taken from the network
%	data structure NET passed as an input argument, while the return
%	argument NET contains the re-estimated values.
%
%	[NET, GAMMA, LOGEV] = EVIDENCE(NET, X, T, NUM) allows the re-
%	estimation  formula to be applied for NUM cycles in which the re-
%	estimated values for the hyperparameters from each cycle are used to
%	re-evaluate the Hessian matrix for the next cycle.  The return value
%	GAMMA is the number of well-determined parameters and LOGEV is the
%	log of the evidence.
%
%	See also
%	MLPPRIOR, NETGRAD, NETHESS, DEMEV1, DEMARD
%

%	Copyright (c) Ian T Nabney (1996-2001)

errstring = consist(net, '', x, t);
if ~isempty(errstring)
  error(errstring);
end

ndata = size(x, 1);
if nargin == 3
  num = 1;
end

% Extract weights from network
w = netpak(net);

% Evaluate data-dependent contribution to the Hessian matrix.
[h, dh] = nethess(w, net, x, t);
clear h;  % To save memory when Hessian is large
if (~isfield(net, 'beta'))
  local_beta = 1;
end

[evec, evl] = eig(dh);
% Now set the negative eigenvalues to zero.
evl = evl.*(evl > 0);
% safe_evl is used to avoid taking log of zero
safe_evl = evl + eps.*(evl <= 0);

[e, edata, eprior] = neterr(w, net, x, t);

if size(net.alpha) == [1 1]
  % Form vector of eigenvalues
  evl = diag(evl);
  safe_evl = diag(safe_evl);
else
  ngroups = size(net.alpha, 1);
  gams = zeros(1, ngroups);
  logas = zeros(1, ngroups);
  % Reconstruct data hessian with negative eigenvalues set to zero.
  dh = evec*evl*evec';
end

% Do the re-estimation.
for k = 1 : num
  % Re-estimate alpha.
  if size(net.alpha) == [1 1]
    % Evaluate number of well-determined parameters.
    L = evl;
    if isfield(net, 'beta')
      L = net.beta*L;
    end
    gamma = sum(L./(L + net.alpha));
    net.alpha = 0.5*gamma/eprior;
    % Partially evaluate log evidence: only include unmasked weights
    logev = 0.5*length(w)*log(net.alpha);
  else
    hinv = inv(hbayes(net, dh));
    for m = 1 : ngroups
      group_nweights = sum(net.index(:, m));
      gams(m) = group_nweights - ...
	        net.alpha(m)*sum(diag(hinv).*net.index(:,m));
      net.alpha(m) = real(gams(m)/(2*eprior(m)));
      % Weight alphas by number of weights in group
      logas(m) = 0.5*group_nweights*log(net.alpha(m));
    end
    gamma = sum(gams, 2);
    logev = sum(logas);
  end
  % Re-estimate beta.
  if isfield(net, 'beta')
      net.beta = 0.5*(net.nout*ndata - gamma)/edata;
      logev = logev + 0.5*ndata*log(net.beta) - 0.5*ndata*log(2*pi);
      local_beta = net.beta;
  end

  % Evaluate new log evidence
  e = errbayes(net, edata);
  if size(net.alpha) == [1 1]
    logev = logev - e - 0.5*sum(log(local_beta*safe_evl+net.alpha));
  else
    for m = 1:ngroups
      logev = logev - e - ...
	  0.5*sum(log(local_beta*(safe_evl*net.index(:, m))+...
	  net.alpha(m)));
    end
  end
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
params.dictsize = length(wordsbi);
params.numcases = length(labels);

fprintf('CV using dataset l=%d, dictSize=%d, CVNUM=%d\n', ...
    length(allSNumBi), length(wordsbi), params.CVNUM)

% initial = 1.1;
randn('state', 0);
rand('state', 0);

allcounts = [];
allfps = []; allfns = [];
for i=[1:params.CVNUM]
    train_ind = cv_obj.training(i);
    test_ind = cv_obj.test(i);
    assert(0==sum(train_ind == test_ind))

    model = trainfuncp(allSNumBi(train_ind), labels(train_ind), params);
%
    [acc pred softpred] = testfuncp(model, ...
        allSNumBi(test_ind), labels(test_ind), params);

    nblbltst = labels(test_ind);
    fp = sum(nblbltst == 0 & pred == 1);
    fn = sum(nblbltst == 1 & pred == 0);
    allfps = [allfps fp];
    allfns = [allfns fn];
    allcounts = [allcounts acc];
end
allcounts
mean(allcounts)

author: sidaw

repo: nbsvm

save_path: github-repos/MATLAB/sidaw-nbsvm


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  In this script, we perform phase transition analysis
%  of Orthogonal matching pursuit.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all;
clear all;
clc;
rng('default');
% Create the directory for storing results
[status_code,message,message_id] = mkdir('bin');
target_file_path = 'bin/ra_mmv_phase_transition_snr_20db_s_4.mat';
N = 256;
S = 4;
pta = spx.pursuit.PhaseTransitionAnalysis(N);
pta.SNR = 20;

% options for CoSaMP MMV solver
solver_options.RankAwareResidual = true;
P = 2;

% pta.NumTrials = 100;
dict_model = @(M, N) spx.dict.simple.gaussian_dict(M, N);
data_model = @(N, K) spx.data.synthetic.SparseSignalGenerator(N, K, S).gaussian;
recovery_solver = @(Phi, K, y) spx.pursuit.joint.CoSaMP(Phi, K, P, solver_options).solve(y).Z;
pta.run(dict_model, data_model, recovery_solver);
pta.save_results(target_file_path);

author: indigits

repo: sparse-plex


text:
% ------------------------------------------------------------------------
%  Copyright (C)
%  Universitat Politecnica de Catalunya BarcelonaTech (UPC) - Spain
%  University of California Berkeley (UCB) - USA
%
%  Jordi Pont-Tuset <jordi.pont@upc.edu>
%  Pablo Arbelaez <arbelaez@berkeley.edu>
%  June 2014
% ------------------------------------------------------------------------
% This file is part of the MCG package presented in:
%    Arbelaez P, Pont-Tuset J, Barron J, Marques F, Malik J,
%    "Multiscale Combinatorial Grouping,"
%    Computer Vision and Pattern Recognition (CVPR) 2014.
% Please consider citing the paper if you use this code.
% ------------------------------------------------------------------------
function [int_area, int_bbox] = boxes_intersection( bbox1, bbox2 )
    up1    = bbox1(1);
    left1  = bbox1(2);
    down1  = bbox1(3);
    right1 = bbox1(4);

    up2    = bbox2(1);
    left2  = bbox2(2);
    down2  = bbox2(3);
    right2 = bbox2(4);

    int_left  = max(left1,left2);
    int_right = min(right1,right2);
    int_up    = max(up1,up2);
    int_down  = min(down1,down2);

    if (int_left<=int_right) && (int_up<=int_down)
        int_bbox = [int_up, int_left, int_down, int_right];
        int_area = box_area(int_bbox);
    else
        int_bbox = [];
        int_area = 0;
    end
end

author: s-gupta

repo: rcnn-depth

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