%% Trend-Cycle tutorial: dating turning  points and computing cyclical statistics
% Authors:   Filippo Ferroni and  Fabio Canova
% Date:     27/05/2020, revised  15/12/2020

close all; clc; clear all;

addpath ../../cmintools/
addpath ../../bvartools/

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% this program illustrates the use of BB program to  date  turning  points
% and  compute business cycle statistics
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%% exercises  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  exercise 1: Vary the censoring  rules.
%  exercise 2: Use  a  shorter  sample. Do  the  turning  points  coincide?
%  exercise 3: Change  the  thresh  parameter.
%  exercise 4: Change  the  data  set.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%=======================================================================
% LOAD DATA
%=======================================================================
lload=0; % =0 euro data, =1 US data

if lload==0
    % Euro area AWM DATABASE: Quarterly
    [a,b,~] = xlsread('awm19up18');
    % names of variables
    varnames = b(1,2:end);

    % time convention: Q1 = .00 and Q4 =0.75
    time = 1970 : .25 : 2017.75;
    time_start   = find(time==1970.50);
    time_start1  = find(time==1999.50);
    time_end     = find(time==2017.75);
    time_break   = find(time==2007.75);

    % The  CREDIT DATA: Quartely
    [e,f,~]=xlsread('ECB_Credit_gaps.xlsx');
    varnames2= f(1,2:end);

    % real GDP, real consumption, real  investment, GDP defl, HICP,
    % short and long term interest rate, commodity  prices,
    % labor  productivity, total  employment, credit to NFC to  GDP
    yy = [a(:,1) a(:,2) a(:,4) a(:,7) a(:,19) a(:,33) a(:,34) ...
        a(:,35) a(:,42) a(:,29) e(time_start-2:time_end,25)];
    % names: YER PCR ITR YED HICP STN LTN COMPR LPROD LNN

    % data transformations:
    % 1. the log of output, consumption,investment
    ddata(:,1:3) = log(yy(time_start+1 : time_end,1:3));
    % 2. the log/log difference of GDP deflator index and CPI
    %ddata(:,4:5) = log(yy(time_start+1  : time_end,4:5));
    ddata(:,4:5) = diff(log(yy(time_start : time_end,4:5)))*400;
    % 3. the level of short and long term interest rate
    ddata(:,6:7) = yy(time_start +1 : time_end,6:7);
    % 4. the log/log difference of commodity prices
    ddata(:,8) = log(yy(time_start+1  : time_end,8));
    %ddata(:,8) = diff(log(yy(time_start : time_end,8)))*400;
    % 5. the taking the log of labor  productivity and  employment
    ddata(:,9:10) = log(yy(time_start+1 : time_end,9:10));
    %ddata(:,9) = diff(log(yy(time_start+1 : time_end,9)));
    %ddata(:,10) = log(yy(time_start+1 : time_end,10));

    %  data for  credit  to  GDP starts  only  at  time_start1 (1999:25)
    ddata(:,11) = yy(time_start+1:time_end,11);
    endd=length(ddata);

    % pick log output, log consumption, log  investment, log labor
    % productivity, log employment, interest rate, inflation rate
    ddd=ddata(:,[1 2 3 9 10 6 5]);
    T=length(ddd);
    time_data = time(time_start+1:time_end);

else

    % US DATABASE: Quarterly

    % real GDP, urate, real consumption, gdp defl, real  investment,
    % capacity utilization, call rate, 10y goverment  bond  rate
    % names: RGDP, URATE, C, GDPdef, inv, capU, callrate, 10ygbond
    [c,d,~] = xlsread('USdata.xlsx');
    % names of variables
    varnames = d(1,2:end);

    % time convention: Q1 = .00 and Q4 =0.75
    time = 1969.75 : .25 : 2019.50;
    time_start   = find(time==1970.25);
    time_end     = find(time==2019.50);
    time_break   = find(time==2007.75);

    tt=length(c);

    % data transformations:

    % 1. difference log of output, consumption,investment
    yy1(1: tt,[1 2 3]) = log(c(1 : tt,[1 3 5]));
    % 2. difference the log GDP deflator index
    yy1(2: tt,6) = diff(log(c(1: tt,4)))*400;
    % 3. leave  the Urate capU unchanged
    yy1(2: tt,[4 5])= c(2:tt,[2 6]);
    % 4. leave call rate, 10 year unchanged
    yy1(2: tt,7:8)= c(2:tt,7:8);
    % 5. compute detrended C/Y, I/Y ratio
    yy1(1: tt,9) = c(1 : tt,3)./c(1: tt,1);
    yy1(1: tt,10) = c(1 : tt,5)./c(1: tt,1);
    % 6. Compute term  spread
    yy1(1: tt,11) = c(1 : tt,8)-c(1: tt,7);

    ddd=yy1;
    T=length(ddd);
    time_data = time(2:tt);
end

%% Parameters

%frequency
freq   = 'q';        % 'q' for quarterly, 'm' for monthly %
if lload==0
    tstart = [1970 3];
    tend   = [2017 3];
elseif lload==1
    tstart = [1970 3];
    tend   = [2019 3];
end

% cycle parameters
options.turnphase   = 2;
options.phase       = 2;          % censoring rules %
options.cycle       = 5;          % lenght of cycle
options.thresh      = 10.4;       % bypasses phase and cycle restriction if peak to trough is > than thresh

options.nrep     = 1;        % 1 if analyze real data
options.complete = 1;        % if= 1- use complete cycles,if =0 -use incomplete cycles (excess still on complete cycle)

dura=zeros(size(ddd,2),2);  ampl=zeros(size(ddd,2),2);
cumm=zeros(size(ddd,2),2);  excc=zeros(size(ddd,2),2);
durcv=zeros(size(ddd,2),2); ampcv=zeros(size(ddd,2),2);
exccv=zeros(size(ddd,2),2); nott=zeros(size(ddd,2),1);
turn    = time(time_start+1:time_end)';
zz      = time(time_start+1:time_end)';
% tid     = linspace(1970.75,2017.75,T)';

for qq = 1 : size(ddd,2)
    if  lload==0
        % EA
        x = ddd(:,qq);
    else
        % US
        x = ddd(time_start+1:time_end,qq);
        time_data = time(time_start+1:time_end);
    end
    % dating_mbbq
    disp('series')
    disp(qq)

    % Computes turning points and imposes restrictions in one step
    [dt_(qq)] = date_(x, time_data, freq, tstart, tend, options);

    zz      = [zz dt_(qq).st];
    turn    = [turn dt_(qq).trinary];

    dura(qq,:)  = dt_(qq).dura;
    ampl(qq,:)  = dt_(qq).ampl;
    cumm(qq,:)  = dt_(qq).cumm;
    excc(qq,:)  = dt_(qq).excc;
    durcv(qq,:) = dt_(qq).durcv;
    ampcv(qq,:) = dt_(qq).amplcv;
    exccv(qq,:) = dt_(qq).exccv;
    nott(qq,:)  = dt_(qq).notentp;
    %plot(tid,trinary,'linewidth',1)
    %title('Turning  points')

end

pause;

disp(' ')
disp('statistics on average cycle')

disp('duration contractions/duration expansions')
disp(dura)

disp('amplitudes contractions/amplitude expansions')
disp(ampl)

disp('cumulative contractions/cumulative expansions')
disp(cumm)

disp('excess movements percent of triangle area')
disp('contractions/expansions')
disp(excc)

disp('cv of durations contractions/expansions')
disp(durcv)

disp('cv of amplitudes contractions/expansions')
disp(ampcv)

disp('cv of excess movements contractions/expansions')
disp(exccv)

%disp('no of its skipped since no peaks+troughs<=2')
%disp(nott)

%disp('states indicators:contraction=0, expansion=1')
%format  short  g
%round(zz, 2)

disp('concordance index BC phases')
ma=corr(zz,'type','Spearman');
disp(ma(2,3:size(ma,2)))

disp('concordance index turning  points')
qa=corr(turn,'type','Spearman');
disp(qa(2,3:size(qa,2)))

[aa,bb]=size(turn);
tsum=zeros(aa,1);
for kk=1:aa
    ssum=0;
    for  qq=2:bb
        ssum=ssum+turn(kk,qq);
    end
    tsum(kk)=ssum;
end

if lload==0
tid     = linspace(1970.75,2017.75,T)';
else
tid     = linspace(1970.75,2019.50,T-3)';
end
plot(tid,tsum,'linewidth',2)
title('Distribution of Turning  points')

disp('Distribution of  peaks')
dp=[zz(find(tsum>=1),1) tsum(find(tsum>=1),1)];
disp(fix(dp))
disp('Distribution of  throughs')
dt=[zz(find(tsum<=-1),1) tsum(find(tsum<=-1),1)];

disp(fix(dt))

% dates
% peaks    mean 1974.00; 1980.00; 1992.00; 2001.00;  2008:00;  2011.25;
%          mode 1974.50; 1980.00; 1992.00; 2001.25;  2008.00;  2011.50;
% throughs mean 1975.25; 1984.25; 1993.75; 2002.75;  2009:50;  2013.00;
%          mode 1975.00; 1984.25; 1993.75; NaN;      2010.00;  2013.00;

% costructing a  recession indicator
recind=zeros(time_end-time_start,1);

% recession dates
rec1b = find(time==1974.00);
rec1e = find(time==1975.25);

rec2b = find(time==1980.00);
rec2e = find(time==1984.25);

rec3b = find(time==1992.00);
rec3e = find(time==1993.75);

rec4b = find(time==2001.00);
rec4e = find(time==2002.75);

rec5b = find(time==2008.00);
rec5e = find(time==2009.50);

rec6b = find(time==2011.25);
rec6e = find(time==2013.00);

for i=rec1b:rec1e
    recind(i,1)=1;
end

for i=rec2b:rec2e
    recind(i,1)=1;
end
for i=rec3b:rec3e
    recind(i,1)=1;
end

for i=rec4b:rec4e
    recind(i,1)=1;
end

for i=rec5b:rec5e
    recind(i,1)=1;
end

for i=rec6b:rec6e
    recind(i,1)=1;
end

if  lload==0
    save Eurorec recind
else
    save Usarec recind
end
return


text:
function a = one_vs_one(c,hyper)

%=========================================================================
% ONE_VS_ONE one_vs_one object
%=========================================================================
% A=ONE_VS_ONE(C,H) returns an one_vs_one object which trains several
% algorithm C on pairwise problems of class i against class j,
% and is initialized with hyperparameters H. The classifiers are combined
% by outputting the class with the most votes.
%
% Model
%  child=svm            -- classifier to use for each sub-problem
%
% Methods:
%  train, test, get_w
%=========================================================================
% Reference : Multi-class Support Vector Machines
% Author    : Jason Weston , C. Watkins
% Link      : http://citeseer.ist.psu.edu/8884.html
%=========================================================================

  % model
  a.nrofclasses=[];
  if nargin==0
    a.child={svm};
  else
    a.child=c;   %% algorithms to use
    if ~isa(c,'cell')
        a.child={c};
    end;
  end

  %% set to use unsigned output if possible
  for i=1:length(a.child)
    a.child{i}.algorithm.use_signed_output=0;
  end

  p=algorithm('one-vs-one');
  a= class(a,'one_vs_one',p);
  %hyperparams are supplied
  if nargin==2
    eval_hyper;
  end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
% EX_MAXWELL_EIG_MIXED2_SQUARE: solve Maxwell eigenproblem in the unit square, with the second mixed formulation.

% 1) PHYSICAL DATA OF THE PROBLEM
clear problem_data
% Physical domain, defined as NURBS map given in a text file
problem_data.geo_name = 'geo_square.txt';

% Type of boundary conditions
problem_data.nmnn_sides   = [];
problem_data.drchlt_sides = [1 2 3 4];

% Physical parameters
problem_data.c_elec_perm = @(x, y) ones(size(x));
problem_data.c_magn_perm = @(x, y) ones(size(x));

% 2) CHOICE OF THE DISCRETIZATION PARAMETERS
clear method_data
method_data.degree     = [3 3];     % Degree of the bsplines
method_data.regularity = [2 2];     % Regularity of the splines
method_data.nsub       = [10 10];     % Number of subdivisions
method_data.nquad      = [4 4];     % Points for the Gaussian quadrature rule

% 3) CALL TO THE SOLVER
[geometry, msh, space, sp_mul, eigv, eigf] = ...
                       solve_maxwell_eig_mixed2_2d (problem_data, method_data);

% 4) POSTPROCESSING
[eigv, perm] = sort (eigv);

fprintf ('First computed eigenvalues: \n')
disp (eigv(1:6))

figure
sp_plot_solution (eigf(1:space.ndof,perm(9)), space, geometry, [30 30])
title ('8^{th} eigenfunction')

%!demo
%! ex_maxwell_eig_mixed2_square

%!test
%! problem_data.geo_name = 'geo_square.txt';
%! problem_data.nmnn_sides   = [];
%! problem_data.drchlt_sides = [1 2 3 4];
%! problem_data.c_elec_perm = @(x, y) ones(size(x));
%! problem_data.c_magn_perm = @(x, y) ones(size(x));
%! method_data.degree     = [3 3];     % Degree of the bsplines
%! method_data.regularity = [2 2];     % Regularity of the splines
%! method_data.nsub       = [10 10];     % Number of subdivisions
%! method_data.nquad      = [4 4];     % Points for the Gaussian quadrature rule
%! [geometry, msh, space, sp_mul, eigv, eigf] = solve_maxwell_eig_mixed2_2d (problem_data, method_data);
%! [eigv, perm] = sort (eigv);
%! nzeros = numel (find (eigv < 1e-4));
%! assert (msh.nel, 100)
%! assert (space.ndof, 312)
%! assert (sp_mul.ndof, 144)
%! assert (eigv(nzeros+(1:5))/pi^2, [1.00000003326399; 1.00000003326400; 2.00000006652799; 4.00000968384168; 4.00000968384168], 2e-14)

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function jed = ymdf_to_jed_islamic_b ( y, m, d, f )

%*****************************************************************************80
%
%% YMDF_TO_JED_ISLAMIC_B converts an Islamic B YMDF date to a JED.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Edward Richards,
%    Algorithm E,
%    Mapping Time, The Calendar and Its History,
%    Oxford, 1999, pages 323-324.
%
%  Parameters:
%
%    Input, integer Y, M, D, real F, the YMDF date.
%
%    Output, real JED, the corresponding Julian Ephemeris Date.
%

%
%  Check the date.
%
  [ y, m, d, ierror ] = ymd_check_islamic ( y, m, d );

  if ( ierror ~= 0 )
    jed = -1.0;
    return
  end
%
%  Convert the calendar date to a computational date.
%
  y_prime = y + 5519 - floor ( ( 12 - m ) / 12 );
  m_prime = mod ( m + 11, 12 );
  d_prime = d - 1;
%
%  Convert the computational date to a JED.
%
  j1 = floor ( ( 10631 * y_prime + 14 ) / 30 );

  j2 = floor ( ( 2951 * m_prime + 51 ) / 100 );

  jed = j1 + j2 + d_prime - 7664 - 0.5 + f;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [X, R] = imstack2vectors(S, MASK)
%IMSTACK2VECTORS Extracts vectors from an image stack.
%   [X, R] = imstack2vectors(S, MASK) extracts vectors from S, which
%   is an M-by-N-by-n stack array of n registered images of size
%   M-by-N each (see Fig. 11.24). The extracted vectors are arranged
%   as the rows of array X. Input MASK is an M-by-N logical or
%   numeric image with nonzero values (1s if it is a logical array)
%   in the locations where elements of S are to be used in forming X
%   and 0s in locations to be ignored. The number of row vectors in X
%   is equal to the number of nonzero elements of MASK. If MASK is
%   omitted, all M*N locations are used in forming X.  A simple way
%   to obtain MASK interactively is to use function roipoly. Finally,
%   R is an array whose rows are the 2-D coordinates containing the
%   region locations in MASK from which the vectors in S were
%   extracted to form X.

%   Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
%   Digital Image Processing Using MATLAB, Prentice-Hall, 2004
%   $Revision: 1.6 $  $Date: 2003/11/21 14:37:21 $

% Preliminaries.
[M, N, n] = size(S);
if nargin == 1
   MASK = true(M, N);
else
   MASK = MASK ~= 0;
end

% Find the set of locations where the vectors will be kept before
% MASK is changed later in the program.
[I, J] = find(MASK);
R = [I, J];

% Now find X.

% First reshape S into X by turning each set of n values along the third
% dimension of S so that it becomes a row of X. The order is from top to
% bottom along the first column, the second column, and so on.
Q = M*N;
X = reshape(S, Q, n);

% Now reshape MASK so that it corresponds to the right locations
% vertically along the elements of X.
MASK = reshape(MASK, Q, 1);

% Keep the rows of X at locations where MASK is not 0.
X = X(MASK, :);

repo: weiyanmin


text:
function y = ft_preproc_slidingrange(dat, width, varargin)

% FT_PREPROC_SLIDINGRANGE computes the range of the data in a sliding time
% window of the width specified. Width should be an odd number (since the
% window needs to be centered on an individual sample).
%
% Use as
%   y = ft_preproc_slidingrange(dat, width, ...)
%
% Optional key-value pair arguments are:
%   'normalize', whether to normalize the range of the data with the square
%                root of the window size

% Copyright (C) 2012, Donders Centre for Cognitive Neuroimaging, Nijmegen, NL
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

normalize = ft_getopt(varargin, 'normalize', false);

% preprocessing fails on channels that contain NaN
if any(isnan(dat(:)))
  ft_warning('FieldTrip:dataContainsNaN', 'data contains NaN values');
end

if mod(width+1, 2)
  ft_error('width should be an odd number');
end

% compute half width
h = (width-1)/2;

n = size(dat,2);
minval = zeros(size(dat));
maxval = zeros(size(dat));

for i=1:n
  begsample = i-h;
  endsample = i+h;
  if begsample<1
    begsample = 1;
  end
  if endsample>n
    endsample=n;
  end
  minval(:,i) = min(dat(:,begsample:endsample),[],2);
  maxval(:,i) = max(dat(:,begsample:endsample),[],2);
end

y = maxval - minval;

if istrue(normalize)
  y = y ./ sqrt(width);
end

author: spm


text:
function varargout = icatb_read_gzip_nii(fileIn, varargin)
%% Read data from gzip nifti file (*.nii.gz) directly using GZIPInputstream.
%
% Inputs:
%
% 1. fileIn - Nifti gzip file name
% 2. varargin - Variable arguments passed in pairs
%
%   a. read_hdr_only - Optional argument. By default, both header and data information is returned. If set to 1, only necessary fields in header
%   are returned (slices and timeNo variables are ignored).
%   b. slices - Axial slice/slices. You could enter row vector of slices. By
%   default, all slices are loaded.
%   c. timepoints - Timepoints of interest. By default, all timepoints are
%   returned.
%   d.mask - Mask of interest. Specify full file path or use boolean mask
%   or indices.
%   e. use_spm -  Gzip files are un-zipped and spm volume functions are used to read nifti
%   files
%   f. buffer_size - Default buffer size is set to 2^14 bytes (memory efficient) when reading
%   using GZIPInputstream. Max buffer size you could specify is 2^31-1.
%
% Outputs:
% [hdr, V] - if read_hdr_only == 1. hdr is header and V is spm volume
% structure.
% [data, hdr, V] if read_hdr_only == 0. Data is 4D array (x,y,z,t) or 2D
% array (x*y*z, t) if mask is specified, hdr is header and V os spm volume structure.
%
%

NIFTI_GZ = 0;

GZIPINFO.isLargeFile = 1;
GZIPINFO.buffer_size = 2^14;

%% Parse inputs
read_hdr_only = 0;
buffer_size = 2^14;
isLargeFile = 0;
try
    buffer_size = GZIPINFO.buffer_size;
catch
end

try
    isLargeFile = GZIPINFO.isLargeFile;
catch
end

isSPMRead = NIFTI_GZ;
if (isempty(isSPMRead))
    isSPMRead = 0;
end
slices = [];
timeNo = [];

for n = 1:2:length(varargin)
    if (strcmpi(varargin{n}, 'read_hdr_only'))
        read_hdr_only = varargin{n + 1};
    elseif (strcmpi(varargin{n}, 'slices'))
        slices = varargin{n + 1};
    elseif (strcmpi(varargin{n}, 'timepoints') || strcmpi(varargin{n}, 'time'))
        timeNo = varargin{n + 1};
    elseif (strcmpi(varargin{n}, 'mask'))
        mask = varargin{n + 1};
    elseif (strcmpi(varargin{n}, 'use_spm'))
        isSPMRead = varargin{n + 1};
    elseif (strcmpi(varargin{n}, 'buffer_size'))
        buffer_size = varargin{n + 1};
    end
end

if (nargout > 3)
    error('Max number of output arguments can be returned is 3');
end

if (read_hdr_only)
    if (nargout > 2)
        error('Max number of output arguments can be returned is 2 with read header only');
    end
end

if (~usejava('jvm'))
    error('Requires jvm to read gzip nifti (*.nii.gz) file');
end

jar_p = fileparts(which('icatb_read_gzip_nii.m'));

if (~isdeployed)
    javaaddpath(fullfile(jar_p, 'icatb_gz.jar'));
end

import icatb_gz.*;

if (~isSPMRead)
    try
        zclass = icatb_gz.read_gzip();
    catch
        isSPMRead = 1;
    end
end

nBytesToRead = 348;
doSwapBytes = 0;

%% Use GZIPInputstream
if (read_hdr_only)
    fid  = java.io.FileInputStream(fileIn);
    zid  = java.util.zip.GZIPInputStream(fid);
    if (nargout == 1)
        zid.skip(40);
        if (~isSPMRead)
            inByteArray = zclass.read(zid, 16);
        else
            inByteArray = readInitialBytes(zid, 16);
        end
    else
        if (~isSPMRead)
            inByteArray = zclass.read(zid, nBytesToRead);
        else
            inByteArray = readInitialBytes(zid, nBytesToRead);
        end
    end
    try
        zid.close();
        fid.close();
    catch
    end
    hdr = read_hdr(inByteArray);
    if (nargout > 1)
        gzV = getVol(fileIn, hdr);
    end

else

    if (~isSPMRead)
        [inByteArray, hdr, doSwapBytes] = readGzip(fileIn, buffer_size, isLargeFile);
        gzV = getVol(fileIn, hdr);
    else
        % Uncompress file and read header
        [hdr, orig_gzfn, gzV, XYZ] = getSpmVol(fileIn, timeNo);
    end

end

if (read_hdr_only)

    varargout{1} = hdr;
    if (nargout == 2)
        varargout{2} = gzV;
    end

    return;

end

dims = hdr.dim(2:5);

if (~exist('slices', 'var') || isempty(slices))
    slices = 1:dims(3);
end

if (~exist('timeNo', 'var') || isempty(timeNo))
    timeNo = 1:dims(4);
end

byteA = double(hdr.bitpix / 8);

convertDataTo2D = 1;

%% Mask
if (exist('mask', 'var') && ~isempty(mask))
    if (ischar(mask))
        mask = spm_read_vols(spm_vol(mask));
        mask(isfinite(mask)==0) = 0;
    end
    if (numel(mask) ~= prod(dims(1:3)))
        maskF = zeros(dims(1:3));
        maskF(mask) = 1;
    else
        maskF = mask;
    end
    maskF = reshape(maskF, dims(1:3));
else
    convertDataTo2D = 0;
    maskF = (ones(dims(1:3)) ~= 0);
end

mask_inds = find(maskF(:, :, slices) ~= 0);

if (~convertDataTo2D)
    data = zeros(dims(1), dims(2), length(slices), length(timeNo));
else
    data = zeros(length(mask_inds), length(timeNo));
end

for t = 1:length(timeNo)
    endT = 0;
    for z = 1:length(slices)
        startT = endT  + 1;

        if (~isSPMRead)

            if (convertDataTo2D)
                mInds = find(maskF(:, :, slices(z)) ~= 0);
            end

            position =  double(byteA*((timeNo(t)-1)*prod(dims(1:3)) + (slices(z)-1)*prod(dims(1:2)))) + double(hdr.vox_offset);

            tmp = inByteArray(position + 1: position + dims(1)*dims(2)*byteA);

            if (~convertDataTo2D)
                if (~doSwapBytes)
                    tmp = double(typecast(tmp, hdr.precision));
                else
                    tmp = double(swapbytes(typecast(tmp, hdr.precision)));
                end
                tmpDat = reshape(tmp, [dims(1), dims(2)]);
            else
                if (~doSwapBytes)
                    tmpDat = double(typecast(tmp, hdr.precision));
                else
                    tmpDat = double(swapbytes(typecast(tmp, hdr.precision)));
                end
                tmpDat = tmpDat(mInds);
            end

            endT = endT + length(tmpDat);

        else

            mInds = find(maskF(:, :, slices(z)) ~= 0);
            tmpa = squeeze(XYZ(1, :, :, slices(z)));
            tmpb = squeeze(XYZ(2, :, :, slices(z)));
            tmpa = tmpa(mInds);
            tmpb = tmpb(mInds);
            tmpDat = [];
            if (~isempty(mInds))
                tmpDat = spm_sample_vol(gzV(t), tmpa, tmpb, slices(z)*ones(size(mInds)), 0);
            end
            endT = endT + length(tmpDat);
            if (~convertDataTo2D)
                tmpDat = reshape(tmpDat, dims(1), dims(2));
            end

        end

        clear tmp;

        if (~convertDataTo2D)
            data(:, :, z, t) = tmpDat;
        else
            data(startT:endT, t) = tmpDat;
        end

    end

end

if (~isSPMRead)
    if (hdr.scl_slope ~=0)
        data = data*hdr.scl_slope + hdr.scl_inter;
    end
end

data (isfinite(data) == 0) = 0;

varargout{1} = data;
varargout{2} = hdr;

if (nargout == 3)
    varargout{3} = gzV(1);
end

if (exist('orig_gzfn', 'var'))
    doCleanUpFiles(orig_gzfn);
end

function [hdr, doSwapBytes] = read_hdr(inByteArray)
%% Read header

doSwapBytes = 0;
corr_types = {'uint8', 'int16', 'int32', 'single', 'double', 'int8', 'uint16', 'uint32'};
codes   = [2, 4, 8, 16, 64, 256, 512, 768];

% necessary fields
if (length(inByteArray) == 16)
    dim = typecast(inByteArray, 'int16');
    if (dim(1) < 1 || dim(1) > 7)
        doSwapBytes = 1;
        dim = swapbytes(typecast(inByteArray, 'int16'));
    end
    hdr.dim = dim(:)';
    return;
end

chk = typecast(inByteArray(41:56), 'int16');
if (chk(1) < 1 || chk(1) > 7)
    doSwapBytes = 1;
end

fieldsIn  = {   'dim',    'int16',  (41:56);
    'datatype',           'int16',  (71:72);
    'bitpix',             'int16',  (73:74);
    'vox_offset',         'single', (109:112);
    'pixdim',             'single', (77:108);
    'scl_slope',          'single', (113:116);
    'scl_inter',          'single', (117:120);
    'qform_code',         'int16',  (253:254);
    'sform_code',         'int16',  (255:256);
    'quatern_b',          'single', (257:260);
    'quatern_c',          'single', (261:264);
    'quatern_d',          'single', (265:268);
    'qoffset_x',          'single', (269:272);
    'qoffset_y',          'single', (273:276);
    'qoffset_z',          'single', (277:280);
    'sizeof_hdr',         'int32',  (1:4);
    'data_type',          'char',   (5:14);
    'db_name',            'char',   (15:32);
    'extents',            'int32',  (33:36);
    'session_error',      'int16',  (37:38);
    'regular',            'char',   (39:39);
    'dim_info',           'uint8',   (40:40);
    'intent_p1',          'single', (57:60);
    'intent_p2',          'single', (61:64);
    'intent_p3',          'single', (65:68);
    'intent_code',        'int16',  (69:70);
    'slice_start',        'int16',  (75:76);
    'slice_end',          'int16',  (121:122);
    'slice_code',         'uint8',   (123:123);
    'xyzt_units',         'uint8',   (124:124);
    'cal_max',            'single', (125:128);
    'cal_min',            'single', (129:132);
    'slice_duration',     'single', (133:136);
    'toffset',            'single', (137:140);
    'glmax',              'single', (141:144);
    'glmin',              'single', (145:148);
    'descrip',            'char',   (149:228);
    'aux_file',           'char',   (229:252);
    'srow_x',             'single', (281:296);
    'srow_y',             'single', (297:312);
    'srow_z',             'single', (313:328);
    'intent_name',        'char',   (329:344);
    'magic',              'char',   (345:348)};

for nF = 1:size(fieldsIn, 1)
    if (strcmpi(fieldsIn{nF, 2}, 'char'))
        tmp = char(inByteArray(fieldsIn{nF, 3}));
        if (size(tmp, 1) > 1)
            tmp = tmp';
        end
    else
        if (~doSwapBytes)
            tmp = typecast(inByteArray(fieldsIn{nF, 3}), fieldsIn{nF, 2});
        else
            tmp = swapbytes(typecast(inByteArray(fieldsIn{nF, 3}), fieldsIn{nF, 2}));
        end
    end
    if (isnumeric(tmp))
        tmp = double(tmp(:)');
    end
    hdr.(fieldsIn{nF, 1}) = tmp;
end

hdr.dim = hdr.dim(:)';
hdr.pixdim = hdr.pixdim(:)';

hdr.precision = [];
chk = find(codes == hdr.datatype);
if (~isempty(chk))
    hdr.precision = corr_types{chk};
end

if isempty(hdr.precision)
    error('Unknown data type');
end

function [inByteArray, hdr, doSwapBytes] = readGzip(fileIn, buffer_size, isLargeFile)
%% Read large data in chunks
%

import icatb_gz.*;
zclass = icatb_gz.read_gzip();
hdrSize = 348;

tmpFileInfo = dir(fileIn);
tmpGBytes = tmpFileInfo(1).bytes/1024/1024/1024;

% if (tmpGBytes > 0.5)
%     isLargeFile = 1;
% end

if (isLargeFile)

    fid  = java.io.FileInputStream(fileIn);
    zid  = java.util.zip.GZIPInputStream(fid);

    try

        hdrBytes = zclass.read(zid, hdrSize);

        [hdr, doSwapBytes] = read_hdr(typecast(hdrBytes, 'uint8'));

        voxel_offset = hdr.vox_offset;
        bytesA = hdr.bitpix/8;
        byte_array_size = voxel_offset + prod(hdr.dim(2:5))*bytesA;
        inByteArray = zeros(byte_array_size, 1, 'int8');
        inByteArray(1:hdrSize) = hdrBytes;
        zid.close();
        fid.close();

        fid  = java.io.FileInputStream(fileIn);
        zid  = java.util.zip.GZIPInputStream(fid);
        zid.skip(voxel_offset);
        count = voxel_offset;
        while 1

            count = count + 1;

            tmp = zclass.read(zid, buffer_size);

            len = length(tmp);

            if (len == 0)
                break;
            end

            endCount = count + len - 1;
            inByteArray(count:endCount) = tmp;
            count = endCount;

        end

        inByteArray = typecast(inByteArray, 'uint8');

        zid.close();
        fid.close();

    catch

        zid.close();
        fid.close();

    end

else

    inByteArray = zclass.read(fileIn, buffer_size);
    inByteArray = typecast(inByteArray, 'uint8');
    [hdr, doSwapBytes] = read_hdr(inByteArray);

end

function V = getVol(fname, hdr)
%% Create spm volume structure from header

fname_i = fname;
[pathstr, fn, extn] = fileparts(fname);

extn = strrep(lower(extn), '.gz', '');
fname = fullfile(pathstr, [fn, extn]);

mat = decode_qform0(hdr);

V   = struct('fname', fname,...
    'dim',   hdr.dim(2:4),...
    'dt',   [hdr.datatype, 0],...
    'pinfo', [hdr.scl_slope hdr.scl_inter hdr.vox_offset]',...
    'mat',   mat,...
    'n',     1,...
    'descrip', 'NIFTI');

h = nifti;
h.mat = V.mat;
h.mat0 = V.mat;
fp  = fopen(fname_i, 'r', 'native');
[~,~,mach] = fopen(fp);
fclose(fp);
dat = file_array(V.fname, hdr.dim(2:5),[hdr.datatype, strfind(mach,'be')], hdr.vox_offset, hdr.scl_slope, hdr.scl_inter);
h.dat = dat;
h.diminfo.slice = 3;
h.diminfo.slice_time.code = hdr.slice_code;
h.diminfo.slice_time.start = hdr.slice_start+1;
h.diminfo.slice_time.end = hdr.slice_end+1;
h.timing.toffset = hdr.toffset;
h.timing.tspace = hdr.pixdim(5);
matcodes = {'UNKNOWN', 'Scanner', 'Aligned', 'Talairach', 'MNI152'};
h.mat_intent = matcodes{hdr.sform_code+1};
h.mat0_intent = matcodes{hdr.qform_code+1};
V.private = h;

function M = decode_qform0(hdr)
% Decode qform info from NIFTI-1 headers.
% _______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

%
% $Id: decode_qform0.m 3131 2009-05-18 15:54:10Z guillaume $

dim    = double(hdr.dim);
pixdim = double(hdr.pixdim);
if ~isfield(hdr,'magic') || hdr.qform_code <= 0,
    flp = spm_flip_analyze_images;
    %disp('------------------------------------------------------');
    %disp('The images are in a form whereby it is not possible to');
    %disp('tell the left and right sides of the brain apart.');
    %if flp,
    %    disp('They are assumed to be stored left-handed.');
    %else
    %    disp('They are assumed to be stored right-handed.');
    %end;
    %disp('------------------------------------------------------');

    %R     = eye(4);
    n      = min(dim(1),3);
    vox    = [pixdim(2:(n+1)) ones(1,3-n)];

    if ~isfield(hdr,'origin') || ~any(hdr.origin(1:3)),
        origin = (dim(2:4)+1)/2;
    else
        origin = double(hdr.origin(1:3));
    end;
    off     = -vox.*origin;
    M       = [vox(1) 0 0 off(1) ; 0 vox(2) 0 off(2) ; 0 0 vox(3) off(3) ; 0 0 0 1];

    % Stuff for default orientations
    if flp, M = diag([-1 1 1 1])*M; end;
else

    % Rotations from quaternions
    R = Q2M(double([hdr.quatern_b hdr.quatern_c hdr.quatern_d]));

    % Translations
    T = [eye(4,3) double([hdr.qoffset_x hdr.qoffset_y hdr.qoffset_z 1]')];

    % Zooms.  Note that flips are derived from the first
    % element of pixdim, which is normally unused.
    n = min(dim(1),3);
    Z = [pixdim(2:(n+1)) ones(1,4-n)];
    Z(Z<0) = 1;
    if pixdim(1)<0, Z(3) = -Z(3); end;
    Z = diag(Z);

    M = T*R*Z;

    % Convert from first voxel at [1,1,1]
    % to first voxel at [0,0,0]
    M = M * [eye(4,3) [-1 -1 -1 1]'];
end;
return;

function M = Q2M(Q)
% Generate a rotation matrix from a quaternion xi+yj+zk+w,
% where Q = [x y z], and w = 1-x^2-y^2-z^2.
% See: http://skal.planet-d.net/demo/matrixfaq.htm
% _______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

%
% $Id: Q2M.m 1143 2008-02-07 19:33:33Z spm $

Q = Q(1:3); % Assume rigid body
w = sqrt(1 - sum(Q.^2));
x = Q(1); y = Q(2); z = Q(3);
if w<1e-7,
    w = 1/sqrt(x*x+y*y+z*z);
    x = x*w;
    y = y*w;
    z = z*w;
    w = 0;
end;
xx = x*x; yy = y*y; zz = z*z; ww = w*w;
xy = x*y; xz = x*z; xw = x*w;
yz = y*z; yw = y*w; zw = z*w;
M = [...
    (xx-yy-zz+ww)      2*(xy-zw)      2*(xz+yw) 0
    2*(xy+zw) (-xx+yy-zz+ww)      2*(yz-xw) 0
    2*(xz-yw)      2*(yz+xw) (-xx-yy+zz+ww) 0
    0              0              0  1];
return;

function varargout = getSpmVol(fileIn, timeNo)
%% Unzip file and read info
%

gzfn = gunzip (fileIn, tempdir);
orig_gzfn = char(gzfn);

if (nargout > 2)
    gzfn = rename_4d_file(orig_gzfn);
    if (~exist('timeNo', 'var') || isempty(timeNo))
        timeNo = (1:size(gzfn, 1));
    end
    gzV = spm_vol(deblank(gzfn(timeNo, :)));
    hdr = gzV(1).private.hdr;

    hdr.dim = double(hdr.dim);

else

    hdr = icatb_read_hdr(deblank(orig_gzfn));
    hdr.dim = double(hdr.dime.dim);

end

varargout{1} = hdr;
varargout{2} = orig_gzfn;

if (nargout >= 3)
    varargout{3} = gzV;
end

if (nargout == 4)
    XYZ = icatb_get_voxel_coords(gzV(1).dim(1:3));
    XYZ = reshape(XYZ, [3, gzV(1).dim(1:3)]);
    varargout{4} = XYZ;
end

function doCleanUpFiles(orig_gzfn)

try
    for nF = 1:size(orig_gzfn, 1)
        delete(deblank(orig_gzfn(nF, :)));
    end
catch
end

function inByteArray = readInitialBytes(zid, nBytesToRead)
%% Use read method in matlab
%

if (isfinite(nBytesToRead))
    inByteArray = zeros(nBytesToRead, 1, 'uint8');
else
    inByteArray = [];
    inByteArray = uint8(inByteArray);
end

count = 0;
tmp = 0;
while ((tmp ~= -1) && (count < nBytesToRead))
    count = count + 1;
    tmp = uint8(read(zid));
    inByteArray(count) = tmp;
end

function [hdr,otherendian] = icatb_read_hdr(fname)
% Read (SPM customised) Analyze header
% fname       - .hdr filename
% hdr         - structure containing Analyze header
% otherendian - byte swapping necessary flag
%_______________________________________________________________________
% Copyright (C) 2005 Wellcome Department of Imaging Neuroscience

% John Ashburner
% $Id: spm_read_hdr.m 112 2005-05-04 18:20:52Z john $

otherendian = 0;
fname = deblank(fname);
if (strcmpi(fname(end-2:end), '.gz'))
    hdr = icatb_read_gzip_nii(fname, 'read_hdr_only', 1);
    return;
end

fid         = fopen(fname,'r','native');
%otherendian = 0;
if (fid > 0)
    dime = read_dime(fid);
    if dime.dim(1)<0 | dime.dim(1)>15, % Appears to be other-endian
        % Re-open other-endian
        fclose(fid);
        if spm_platform('bigend'), fid = fopen(fname,'r','ieee-le');
        else,                      fid = fopen(fname,'r','ieee-be'); end;
        otherendian = 1;
        dime = read_dime(fid);
    end;
    %hk       = read_hk(fid);
    %hist     = read_hist(fid);
    %hdr.hk   = hk;
    hdr.dime = dime;
    %hdr.hist = hist;

    % SPM specific bit - unused
    %if hdr.hk.sizeof_hdr > 348,
    %	spmf = read_spmf(fid,dime.dim(5));
    %	if ~isempty(spmf),
    %		hdr.spmf = spmf;
    %	end;
    %end;

    fclose(fid);
else,
    hdr = [];
    otherendian = NaN;
    %error(['Problem opening header file (' fopen(fid) ').']);
end;
return;
%_______________________________________________________________________
%_______________________________________________________________________
function hk = read_hk(fid)
% read (struct) header_key
%-----------------------------------------------------------------------
fseek(fid,0,'bof');
hk.sizeof_hdr 		= fread(fid,1,'int32');
hk.data_type  		= mysetstr(fread(fid,10,'uchar'))';
hk.db_name    		= mysetstr(fread(fid,18,'uchar'))';
hk.extents    		= fread(fid,1,'int32');
hk.session_error	= fread(fid,1,'int16');
hk.regular			= mysetstr(fread(fid,1,'uchar'))';
hk.hkey_un0			= mysetstr(fread(fid,1,'uchar'))';
if isempty(hk.hkey_un0), error(['Problem reading "hk" of header file (' fopen(fid) ').']); end;
return;
%_______________________________________________________________________
%_______________________________________________________________________
function dime = read_dime(fid)
% read (struct) image_dimension
%-----------------------------------------------------------------------
fseek(fid,40,'bof');
dime.dim		= fread(fid,8,'int16')';
%dime.vox_units	= mysetstr(fread(fid,4,'uchar'))';
%dime.cal_units	= mysetstr(fread(fid,8,'uchar'))';
%dime.unused1	= fread(fid,1,'int16');
%dime.datatype	= fread(fid,1,'int16');
%dime.bitpix		= fread(fid,1,'int16');
%dime.dim_un0	= fread(fid,1,'int16');
%dime.pixdim		= fread(fid,8,'float')';
%dime.vox_offset	= fread(fid,1,'float');
%dime.funused1	= fread(fid,1,'float');
%dime.funused2	= fread(fid,1,'float');
%dime.funused3	= fread(fid,1,'float');
%dime.cal_max	= fread(fid,1,'float');
%dime.cal_min	= fread(fid,1,'float');
%dime.compressed	= fread(fid,1,'int32');
%dime.verified	= fread(fid,1,'int32');
dime.glmax		= fread(fid,1,'int32');
dime.glmin		= fread(fid,1,'int32');
if isempty(dime.glmin), error(['Problem reading "dime" of header file (' fopen(fid) ').']); end;
return;
%_______________________________________________________________________
%_______________________________________________________________________
function hist = read_hist(fid)
% read (struct) data_history
%-----------------------------------------------------------------------
fseek(fid,148,'bof');
hist.descrip	= mysetstr(fread(fid,80,'uchar'))';
hist.aux_file	= mysetstr(fread(fid,24,'uchar'))';
hist.orient		= fread(fid,1,'uchar');
hist.origin		= fread(fid,5,'int16')';
hist.generated	= mysetstr(fread(fid,10,'uchar'))';
hist.scannum	= mysetstr(fread(fid,10,'uchar'))';
hist.patient_id	= mysetstr(fread(fid,10,'uchar'))';
hist.exp_date	= mysetstr(fread(fid,10,'uchar'))';
hist.exp_time	= mysetstr(fread(fid,10,'uchar'))';
hist.hist_un0	= mysetstr(fread(fid,3,'uchar'))';
hist.views		= fread(fid,1,'int32');
hist.vols_added	= fread(fid,1,'int32');
hist.start_field= fread(fid,1,'int32');
hist.field_skip	= fread(fid,1,'int32');
hist.omax		= fread(fid,1,'int32');
hist.omin		= fread(fid,1,'int32');
hist.smax		= fread(fid,1,'int32');
hist.smin		= fread(fid,1,'int32');
if isempty(hist.smin), error(['Problem reading "hist" of header file (' fopen(fid) ').']); end;
return;
%_______________________________________________________________________
%_______________________________________________________________________
function spmf = read_spmf(fid,n)
% Read SPM specific fields
% This bit may be used in the future for extending the Analyze header.

fseek(fid,348,'bof');
mgc = fread(fid,1,'int32');    % Magic number
if mgc ~= 20020417, spmf = []; return; end;

for j=1:n,
    spmf(j).mat    = fread(fid,16,'double'); % Orientation information
    spmf(j).unused = fread(fid,384,'uchar'); % Extra unused stuff
    if length(spmf(j).unused)<384,
        error(['Problem reading "spmf" of header file (' fopen(fid) ').']);
    end;
    spmf(j).mat = reshape(spmf(j).mat,[4 4]);
end;
return;
%_______________________________________________________________________
%_______________________________________________________________________
function out = mysetstr(in)
tmp = find(in == 0);
tmp = min([min(tmp) length(in)]);
out = setstr([in(1:tmp)' zeros(1,length(in)-(tmp))])';
return;
%_______________________________________________________________________
%_______________________________________________________________________

function files = rename_4d_file(files, fileNumber)
%For nifti and analyze files rename the files by adding a number at the
%end

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

if ~isempty(files)

    % Do pattern match to check IMG or NII files
    files = cellstr(files);

    % MATCH IMG OR NII
    checkNII = regexpi(files, '\.nii$|\.gz$');
    checkIMG = regexpi(files, '\.img$');

    good_inds1 = icatb_good_cells(checkNII);
    good_inds2 = icatb_good_cells(checkIMG);

    % Get good cells
    good_inds = (good_inds1 | good_inds2);
    good_inds = find(good_inds ~= 0);

    clear good_inds1 good_inds2 checkNII checkIMG;

    % If IMG or NII files exist check the headers of the image files
    if ~isempty(good_inds)
        imFiles = files(good_inds);
        if (~isempty(fileNumber)) && (length(fileNumber) == 1) && (fileNumber == 1)
            imFiles = strcat(imFiles, ',1');
            files(good_inds) = imFiles;
        else
            % Loop over img or nii files
            for nIm = 1:length(imFiles)
                currentFile = imFiles{nIm};

                try
                    numFiles = getDims(currentFile);
                catch
                    numFiles = 1;
                end

                if (isempty(fileNumber))
                    tempFileNum = (1:numFiles);
                else
                    tempFileNum = fileNumber;
                end

                tempFileNum(tempFileNum > numFiles) = [];

                if ~isempty(tempFileNum)
                    tempFiles = [repmat(currentFile, length(tempFileNum), 1), repmat(',', length(tempFileNum), 1), numberToString(tempFileNum(:))];
                else
                    tempFiles = '';
                end
                files{good_inds(nIm)} = tempFiles;
            end
            % End loop over img or nii files
            files = cellstr(char(files));
            ind = icatb_good_cells(files);
            files = files(ind);
        end
    end
    % End for checking IMG or NII files

    if ~isempty(files)
        files = char(files);
    end

end

function str = numberToString(nums)
% number to string. pad space after the number

try
    str = arrayfun(@num2str, nums, 'uniformoutput', false);
catch
    str = cell(length(nums), 1);
    for n = 1:length(nums)
        str{n} = num2str(nums(n));
    end
end

str = char(str);

function tp = getDims(currentFile)

if (~strcmpi(currentFile(end-2:end),'.gz'))
    ni = nifti(currentFile);
    tp = ni.dat.dim(4);
else
    ni = icatb_read_hdr(currentFile);
    tp = ni.dim(5);
end

function ind = icatb_good_cells(mycell)
%% Find good cells

if ~iscell(mycell)
    mycell = {mycell};
end

ind = cellfun('isempty', mycell);

% Good cells
ind = (ind == 0);

% Convert ind to row vector
ind = ind(:)';

function xyz = icatb_get_voxel_coords(dims)
%% Get voxel coords
%

xdim = dims(1);
ydim = dims(2);
zdim = dims(3);

[xords, yords] = ndgrid(1:xdim, 1:ydim);
xords = xords(:)'; yords = yords(:)';
CrPl    = 1:zdim;
zords   = CrPl(:)*ones(1,xdim*ydim);
xyz   = [repmat(xords,1,numel(CrPl)); ...
    repmat(yords,1,numel(CrPl)); ...
    reshape(zords',1,[])];


text:
% Generating VOCopts varible of VOC evalution toolkit for HollywoodHeads dataset
%
% This code is based on VOCinit.m from VOC2012 devkit.
% Copyright Tuan-Hung VU - tuanhungvu@gmail.com

clear VOCopts

VOCopts.dataset='HollywoodHeads';

% change this path to a writable directory for your results
curPath = pwd;

% change this path to point to your copy of the HollywoodHeads dataset
VOCopts.datadir= fullfile(curPath, 'data');

VOCopts.resdir=fullfile(curPath, 'results', VOCopts.dataset, 'VOCEval');
if ~exist(VOCopts.resdir, 'dir')
    mkdir(VOCopts.resdir);
end

% change this path to a writable local directory for the annotation cache
VOCopts.localdir=fullfile(curPath, 'results', VOCopts.dataset, 'VOCEval', 'anns');
if ~exist(VOCopts.localdir, 'dir')
    mkdir(VOCopts.localdir);
end

% initialize the training set
VOCopts.trainset='train'; % use train for development

% initialize the test set
VOCopts.testset='test'; % use validation data for development test set

% initialize main paths
VOCopts.annopath=fullfile(VOCopts.datadir, VOCopts.dataset, 'Annotations', '%s.xml');
VOCopts.imgpath=fullfile(VOCopts.datadir, VOCopts.dataset, 'JPEGImages', '%s.jpeg');
VOCopts.imgsetpath=fullfile(VOCopts.datadir, VOCopts.dataset, 'Splits', '%s.txt');

VOCopts.detrespath=fullfile(VOCopts.resdir, ['%s_det_' VOCopts.testset '_%s.txt']);

VOCopts.classes={'head'};
VOCopts.nclasses=length(VOCopts.classes);

% overlap threshold
VOCopts.minoverlap=0.5;

% annotation cache for evaluation
VOCopts.annocachepath=fullfile(VOCopts.localdir, '%s_anno.mat');

author: aosokin


text:
% See how well partial Kalman filter updates work

seed = 0;
rand('state', seed);
randn('state', seed);
nlandmarks = 6;
T = 12;

[A,B,C,Q,R,Qbig,Rbig,init_x,init_V,robot_block,landmark_block,...
	  true_landmark_pos, true_robot_pos, true_data_assoc, ...
	  obs_rel_pos, ctrl_signal] = mk_linear_slam(...
	      'nlandmarks', nlandmarks, 'T', T, 'ctrl', 'leftright', 'data-assoc', 'cycle');

% exact
[xe, Ve] = kalman_filter(obs_rel_pos, A, C, Qbig, Rbig, init_x, init_V, ...
				     'model', true_data_assoc, 'u', ctrl_signal, 'B', B);

% approx
%k = nlandmarks-1; % exact
k = 3;
ndx = {};
for t=1:T
  landmarks = unique(true_data_assoc(t:-1:max(t-k,1)));
  tmp = [landmark_block(:, landmarks) robot_block'];
  ndx{t} = tmp(:);
end

[xa, Va] = kalman_filter(obs_rel_pos, A, C, Qbig, Rbig, init_x, init_V, ...
			 'model', true_data_assoc, 'u', ctrl_signal, 'B', B, ...
		       'ndx', ndx);

nrows = 10;
stepsize = T/(2*nrows);
ts = 1:stepsize:T;

if 1 % plot

clim = [0 max(max(Va(:,:,end)))];

figure(2)
if 0
  imagesc(Ve(1:2:end,1:2:end, T))
  clim = get(gca,'clim');
else
  i = 1;
  for t=ts(:)'
    subplot(nrows,2,i)
    i = i + 1;
    imagesc(Ve(1:2:end,1:2:end, t))
    set(gca, 'clim', clim)
    colorbar
  end
end
suptitle('exact')

figure(3)
if 0
  imagesc(Va(1:2:end,1:2:end, T))
  set(gca,'clim', clim)
else
  i = 1;
  for t=ts(:)'
    subplot(nrows,2,i)
    i = i+1;
    imagesc(Va(1:2:end,1:2:end, t))
    set(gca, 'clim', clim)
    colorbar
  end
end
suptitle('approx')

figure(4)
i = 1;
for t=ts(:)'
  subplot(nrows,2,i)
  i = i+1;
  Vd = Va(1:2:end,1:2:end, t) - Ve(1:2:end,1:2:end,t);
  imagesc(Vd)
  set(gca, 'clim', clim)
  colorbar
end
suptitle('diff')

end % all plot

for t=1:T
  %err(t)=rms(xa(:,t), xe(:,t));
  err(t)=rms(xa(1:end-2,t), xe(1:end-2,t)); % exclude robot
end
figure(5);plot(err)
title('rms mean pos')

for t=1:T
  i = 1:2*nlandmarks;
  denom = Ve(i,i,t) + (Ve(i,i,t)==0);
  Vd =(Va(i,i,t)-Ve(i,i,t)) ./ denom;
  Verr(t) = max(Vd(:));
end
figure(6); plot(Verr)
title('max relative Verr')

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function ray = rayTracer(ray,ord,rMax)
%+========================================================================+
%|                                                                        |
%|           OPENRAY - LIBRARY FOR TRI-DIMENSIONAL RAY TRACING            |
%|           openRay is part of the GYPSILAB toolbox for Matlab           |
%|                                                                        |
%| COPYRIGHT : Matthieu Aussal (c) 2017-2018.                             |
%| PROPERTY  : Centre de Mathematiques Appliquees, Ecole polytechnique,   |
%| route de Saclay, 91128 Palaiseau, France. All rights reserved.         |
%| LICENCE   : This program is free software, distributed in the hope that|
%| it will be useful, but WITHOUT ANY WARRANTY. Natively, you can use,    |
%| redistribute and/or modify it under the terms of the GNU General Public|
%| License, as published by the Free Software Foundation (version 3 or    |
%| later,  http://www.gnu.org/licenses). For private use, dual licencing  |
%| is available, please contact us to activate a "pay for remove" option. |
%| CONTACT   : matthieu.aussal@polytechnique.edu                          |
%| WEBSITE   : www.cmap.polytechnique.fr/~aussal/gypsilab                 |
%|                                                                        |
%| Please acknowledge the gypsilab toolbox in programs or publications in |
%| which you use it.                                                      |
%|________________________________________________________________________|
%|   '&`   |                                                              |
%|    #    |   FILE       : rayTracer.m                                   |
%|    #    |   VERSION    : 0.41                                          |
%|   _#_   |   AUTHOR(S)  : Matthieu Aussal                               |
%|  ( # )  |   CREATION   : 14.03.2017                                    |
%|  / 0 \  |   LAST MODIF : 01.04.2018                                    |
%| ( === ) |   SYNOPSIS   : Ray tracer with octree speedup                |
%|  `---'  |                                                              |
%+========================================================================+

% Infos
disp('====> RAY TRACING <====')
tps = time();

% Mesh dimension
Nelt = length(ray.msh);

% Prepare sphere
tree = rayTreeInit(ray);

% Available ray
ind = find( (ray.dst<rMax) & (sum(ray.dir.^2,2)~=0) );

% Current order
n = length(ray.pos)-1;

% Infos
disp([' ~~> Tree with ',num2str(length(tree)),' stage(s) - Elapsed time is ', ...
                num2str(time()-tps),' seconds.'])

% Iterative loop
while ~isempty(ind) && (n < ord)
    % Infos
    tps = time();

    % Full repartition for smallest meshes
    if (Nelt < 50)
        Iray = cell(Nelt,1);
        for el = 1:Nelt
            Iray{el} = ind;
        end

    % Hierarchical repartition
    else
        Iray = rayTree(ray,tree,ind);
    end

    % Intersection between mesh and ray
    ray = rayCollision(ray,Iray);

    % Update available ray
    ind = find( (ray.dst<rMax) & (sum(ray.dir.^2,2)~=0) );

    % Order incrementation
    n = n + 1;

    % Infos
    disp([' + Step ',num2str(n), ' - Elapsed time is ', ...
                num2str(time()-tps),' seconds.'])
end
end

function tps = time()
tps = clock;
tps = tps(4)*3600 + tps(5)*60 + tps(6);
end

%     % Representation graphique
%     figure
%     hold on
%     for el = 1:Nelt
%         subElt = ray.msh.sub(el);
%         subRay = ray.sub(Iray{el});
%         plot(subElt,1)
%         plot(subRay)
%         pause
%     end
%     hold off

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [output] = mf(input)
% A wrapper of MF (Multi-deviation Fusion method)
%
% -------------------------------------------------------------------------
% @article{fu2016mf,
%   title={A fusion-based enhancing method for weakly illuminated images},
%   author={Fu, Xueyang and Zeng, Delu and Huang, Yue and Liao, Yinghao and Ding, Xinghao and Paisley, John},
%   journal={Signal Processing},
%   volume={129},
%   pages={82--96},
%   year={2016},
%   publisher={Elsevier}
% }
% -------------------------------------------------------------------------
%
% OpenCE https://github.com/baidut/OpenCE
% Author: Zhenqiang Ying
% 2017-5-12

if nargin == 0
   I = imload; % imread('yellowlily.jpg'); %
   J = mf(I);
   ezFig(I,J);
   return;
end

Require MF

if ~isfloat(input), input = im2double(input); end
if isfloat(input), input = (input)*255; end

output = multi_fusion(input);
end

author: baidut

repo: BIMEF

save_path: github-repos/MATLAB/baidut-BIMEF


text:
function [g1, g2] = simXsimKernDiagGradient(simKern1, simKern2, t, covDiag)

% SIMXSIMKERNDIAGGRADIENT Gradient for the diagonal between two SIM kernels.
% FORMAT
% DESC computes cross gradient of parameters in the diagonal of a cross
% kernel between two sim kernels for the multiple output kernel.
% ARG simKern1 : the kernel structure associated with the first SIM
% kernel.
% ARG simKern2 : the kernel structure associated with the second SIM
% kernel.
% ARG t : inputs for which kernel is to be computed.
% ARG covDiag : gradient of the objective function with respect to
% the elements of the cross kernel matrix.
% RETURN g1 : gradient of the parameters of the first kernel, for
% ordering see simKernExtractParam.
% RETURN g2 : gradient of the parameters of the second kernel, for
% ordering see simKernExtractParam.
%
% SEEALSO : simXsimKernDiagCompute.m, simKernParamInit, simKernExtractParam
%
% COPYRIGHT : Mauricio A. Alvarez, 2010

% KERN

if size(t, 2) > 1
  error('Input can only have one column');
end

if simKern1.inverseWidth ~= simKern2.inverseWidth
  error('Kernels cannot be cross combined if they have different inverse widths.')
end
% The normalisation in the SIM kernel arises from the use of the normalised
% version of the RBF kernel. Hence, both have to be normalised or not.
if ~isfield(simKern1, 'isNormalised')
    isSim1Normalised = false;
else
    isSim1Normalised = simKern1.isNormalised;
end
if ~isfield(simKern2, 'isNormalised')
    isSim2Normalised = false;
else
    isSim2Normalised = simKern2.isNormalised;
end
if isSim1Normalised ~= isSim2Normalised
    error('Both SIM kernels have to be either normalised or not.');
end

sigma = sqrt(2/simKern1.inverseWidth);

if (~simKern1.isStationary || ~ simKern2.isStationary)
    [h1, dh1_dD1, dh1_dD2, dh1_dsigma] = ...
        simXsimComputeDiagH(t, simKern1.decay, simKern2.decay, simKern1.delay, simKern2.delay, sigma);
else
    [h1, dh1_dD1, dh1_dD2, dh1_dsigma] = ...
        simXsimComputeDiagHStat(t, simKern1.decay, simKern2.decay, simKern1.delay, simKern2.delay, sigma);
end

% Avoid making the expensive call twice unless really necessary
if ((simKern1.decay == simKern2.decay) && (simKern1.delay == simKern2.delay)),
  h2 = h1;
  dh2_dD2 = dh1_dD1;
  dh2_dD1 = dh1_dD2;
  dh2_dsigma = dh1_dsigma;
elseif (~simKern1.isStationary) || (~simKern2.isStationary)
  [h2, dh2_dD2, dh2_dD1, dh2_dsigma] = ...
      simXsimComputeDiagH(t, simKern2.decay, simKern1.decay, simKern2.delay, simKern1.delay, sigma);
else
    [h2, dh2_dD2, dh2_dD1, dh2_dsigma] = ...
      simXsimComputeDiagHStat(t, simKern2.decay, simKern1.decay, simKern2.delay, simKern1.delay, sigma);
end

dK_dsigma = dh1_dsigma + dh2_dsigma;

if isfield(simKern1, 'isVarS') && (simKern1.isVarS)
    K = 0.5 * (h1 + h2);
    if ~isSim1Normalised
        K = sqrt(pi) * K;
        dk_dD1 = (sum(covDiag.*dh1_dD1) + sum(covDiag.*dh2_dD1))*0.5*sqrt(pi)*sigma;
        dk_dD2 = (sum(covDiag.*dh1_dD2) + sum(covDiag.*dh2_dD2))*0.5*sqrt(pi)*sigma;
        dk_dsigma = sum(covDiag.*(dK_dsigma*0.5*sqrt(pi)*sigma + K));
    else
        dk_dD1 = (sum(covDiag.*dh1_dD1) + sum(covDiag.*dh2_dD1))*0.5;
        dk_dD2 = (sum(covDiag.*dh1_dD2) + sum(covDiag.*dh2_dD2))*0.5;
        dk_dsigma = 0.5 * sum(covDiag.*dK_dsigma);
    end
    dk_dinvWidth = -0.5*sqrt(2)/(simKern1.inverseWidth* ...
        sqrt(simKern1.inverseWidth))*dk_dsigma;
    % only pass the gradient with respect to the inverse width to one
    % of the gradient vectors ... otherwise it is counted twice.
    g1 = real([dk_dD1 dk_dinvWidth]);
    g2 = real([dk_dD2 0]);
else
    if isfield(simKern1, 'isNegativeS') && (simKern1.isNegativeS == true)
        C1 = simKern1.sensitivity;
        C2 = simKern2.sensitivity;
    else
        C1 = sqrt(simKern1.variance);
        C2 = sqrt(simKern2.variance);
    end
    K = 0.5 * (h1 + h2);
    var2 = C1*C2;
    if ~isSim1Normalised
        K = sqrt(pi) * K;
        dk_dD1 = (sum(covDiag.*dh1_dD1) + sum(covDiag.*dh2_dD1))*0.5*sqrt(pi)*sigma*var2;
        dk_dD2 = (sum(covDiag.*dh1_dD2) + sum(covDiag.*dh2_dD2))*0.5*sqrt(pi)*sigma*var2;
        dk_dsigma = sum(covDiag.*(dK_dsigma*0.5*sqrt(pi)*sigma + K))*var2;
        dk_dC1 = sigma * C2 * sum(covDiag.*K);
        dk_dC2 = sigma * C1 * sum(covDiag.*K);
    else
        dk_dD1 = (sum(covDiag.*dh1_dD1) + sum(covDiag.*dh2_dD1))*0.5*var2;
        dk_dD2 = (sum(covDiag.*dh1_dD2) + sum(covDiag.*dh2_dD2))*0.5*var2;
        dk_dsigma = 0.5 * var2 * sum(covDiag.*dK_dsigma);
        dk_dC1 = C2 * sum(covDiag.*K);
        dk_dC2 = C1 * sum(covDiag.*K);
    end
    if isfield(simKern1, 'isNegativeS') && simKern1.isNegativeS
        dk_dSim1Variance = dk_dC1;
    else
        dk_dSim1Variance = dk_dC1*0.5/C1;
    end
    if isfield(simKern2, 'isNegativeS') && simKern2.isNegativeS
        dk_dSim2Variance = dk_dC2;
    else
        dk_dSim2Variance = dk_dC2*0.5/C2;
    end
    dk_dinvWidth = -0.5*sqrt(2)/(simKern1.inverseWidth* ...
        sqrt(simKern1.inverseWidth))*dk_dsigma;
    % only pass the gradient with respect to the inverse width to one
    % of the gradient vectors ... otherwise it is counted twice.
    g1 = real([dk_dD1 dk_dinvWidth dk_dSim1Variance]);
    g2 = real([dk_dD2 0 dk_dSim2Variance]);
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function map = tab10(N)
% Qualitative colormap from MatPlotLib, for plots using the line ColorOrder.
% In MatPlotLib 2 it is named VEGA10, for MatPlotLib 3 was renamed TAB10.
%
% Copyright (c) 2017-2020 Stephen Cobeldick
%
%%% Syntax:
%  map = tab10
%  map = tab10(N)
%
% For MatPlotLib 2.0 improved colormaps were created for plot lines of
% categorical data. The new colormaps are introduced here:
% <http://matplotlib.org/2.0.0rc2/users/dflt_style_changes.html>
% VEGA10/TAB10 is the default Line Color Order for MatPlotLib 2 and 3.
%
% MATLAB axes ColorOrder (note that this is NOT the axes COLORMAP):
% <https://www.mathworks.com/help/matlab/creating_plots/defining-the-color-of-lines-for-plotting.html>
% <https://www.mathworks.com/help/matlab/graphics_transition/why-are-plot-lines-different-colors.html>
%
%% Examples %%
%
%%% PLOT using matrices:
% N = 10;
% axes('ColorOrder',tab10(N),'NextPlot','replacechildren')
% X = linspace(0,pi*3,1000);
% Y = bsxfun(@(x,n)sqrt(n)*sin(x+2*n*pi/N), X(:), 1:N);
% plot(X,Y, 'linewidth',4)
%
%%% PLOT in a loop:
% N = 10;
% set(0,'DefaultAxesColorOrder',tab10(N))
% X = linspace(0,pi*3,1000);
% Y = bsxfun(@(x,n)sqrt(n)*sin(x+2*n*pi/N), X(:), 1:N);
% for n = 1:N
%     plot(X(:),Y(:,n), 'linewidth',4);
%     hold all
% end
%
%%% LINE using matrices:
% N = 10;
% set(0,'DefaultAxesColorOrder',tab10(N))
% X = linspace(0,pi*3,1000);
% Y = bsxfun(@(x,n)sqrt(n)*cos(x+2*n*pi/N), X(:), 1:N);
% line(X(:),Y)
%
%% Input and Output Arguments %%
%
%%% Inputs (*=default):
% N = NumericScalar, N>=0, an integer to define the colormap length.
%   = *[], use the length of the current figure's colormap (see COLORMAP).
%
%%% Outputs:
% map = NumericMatrix, size Nx3, a colormap of RGB values between 0 and 1.
%
% See also TAB20 TAB20B TAB20C SET VIRIDIS TWILIGHT LINES COLORMAP PARULA

if nargin<1 || isnumeric(N)&&isequal(N,[])
	N = size(get(gcf,'colormap'),1);
else
	assert(isscalar(N)&&isfinite(N)&&isreal(N),...
		'SC:tab10:N:NotRealFiniteScalarNumeric',...
		'First argument <N> must be a real finite numeric scalar.')
end
%
hex = ['#1f77b4';'#ff7f0e';'#2ca02c';'#d62728';'#9467bd';'#8c564b';'#e377c2';'#7f7f7f';'#bcbd22';'#17becf'];
raw = sscanf(hex.','#%2x%2x%2x',[3,Inf]).';
%
map = raw(1+mod(0:N-1,size(raw,1)),:) / 255;
%
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%tab10

author: kavli-ntnu


text:
function plotData(D)

data = D.data;
Jumps = D.Jumps;
P = D.P;

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                           States                              %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
figH = figure(1); clf;
set(figH,'Name','States','NumberTitle','off')

subplot(4,1,1);hold on;
plot(data.time,data.state.th);
ylabel('Angle (rad)')
xlabel('Time (s)')
dottedLine(Jumps,axis);

subplot(4,1,2); hold on;
plot(data.time,data.state.dth);
ylabel('Rate (rad)')
xlabel('Time (s)')
dottedLine(Jumps,axis);

subplot(4,1,3); hold on
plot(data.time,data.state.x,'r');
plot(data.time,data.state.y,'b');
legend('x','y');
ylabel('position (m)')
xlabel('Time (s)')
dottedLine(Jumps,axis);

subplot(4,1,4); hold on
plot(data.time,data.state.dx,'r');
plot(data.time,data.state.dy,'b');
legend('dx','dy');
ylabel('velocity (m/s)')
xlabel('Time (s)')
dottedLine(Jumps,axis);

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                                Traces                                   %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

L = P.L;
LL = P.FullLength;
x0 = data.state.x;
y0 = data.state.y;
x1 = x0 - L*sin(data.state.th);
y1 = y0 + L*cos(data.state.th);
x2 = x0 - LL*sin(data.state.th);
y2 = y0 + LL*cos(data.state.th);
bnd = [min(min(x0),min(x2)),max(max(x0),max(x2))];

figH = figure(3); clf; hold on;
set(figH,'Name','Traces','NumberTitle','off')

plot(bnd,[0,0],'k-','LineWidth',1);
plot(x0,y0,'r','LineWidth',3)
plot(x1,y1,'m','LineWidth',3)
plot(x2,y2,'b','LineWidth',3)
legend('ground','O','CoM','Tip');
axis 'equal';
xlabel('Horizontal Position (m)');
ylabel('Vertical Position (m)');

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                           Forces and Energy                             %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
figH = figure(2); clf;
set(figH,'Name','Energy','NumberTitle','off')

subplot(3,1,1); hold on;
plot(data.time,data.contact.h,'r');
plot(data.time,data.contact.v,'b');
legend('Horizontal','Vertical')
ylabel('Contact Force (N)')
xlabel('Time (s)')
dottedLine(Jumps,axis);

subplot(3,1,2); hold on;
E_pot = data.energy.potential;
E_kin = data.energy.kinetic;
plot(data.time,E_pot,'b');
plot(data.time,E_kin,'r');
plot(data.time,E_kin + E_pot,'k','LineWidth',2);
legend('Potential','Kinetic','Total','Location','NorthWest')
ylabel('Energy (J)')
xlabel('Time (s)')
dottedLine(Jumps,axis);

subplot(3,1,3); hold on;
plot(data.time,E_kin + E_pot,'k','LineWidth',2);
ylabel('Total Energy (J)')
xlabel('Time (s)')
dottedLine(Jumps,axis);

figH = figure(4); clf; hold on
set(figH,'Name','Contacts','NumberTitle','off')

th = data.state.th*180/pi;
plot(th,data.contact.h,'r');
plot(th,data.contact.v,'b');
legend('Horizontal','Vertical')
ylabel('Contact Force (N)')
xlabel('Angle (deg)')

end

%%%% SUB FUNCTIONS %%%%

function dottedLine(time,AXIS)

for i=1:length(time)
    %Plots a dotted line between phases
    plot(time(i)*[1;1],[AXIS(3);AXIS(4)],'k:','LineWidth',1);
end

end

author: MatthewPeterKelly

repo: dscTutorials


text:
%LUT  Performs a look-up table transform of an array
%
%     dst = cv.LUT(src, lut)
%
% ## Input
% * __src__ input array of 8-bit elements (`uint8` or `int8`).
% * __lut__ look-up table of 256 elements; in case of multi-channel input
%   array `src`, the table should either have a single channel (in this case
%   the same table is used for all channels) or the same number of channels as
%   in the input array.
%
% ## Output
% * __dst__ output array of the same size and number of channels as `src`, and
%   the same depth as `lut`.
%
% The function cv.LUT fills the output array with values from the look-up
% table. Indices of the entries are taken from the input array. That is, the
% function processes each element of src as follows:
%
%     dst(I) = lut(src(I) + d)
%
% where:
%
%     d = { 0   if src has uint8 depth
%         { 128 if src has  int8 depth
%
% See also: cv.convertScaleAbs, cv.convertTo, cv.applyColorMap, intlut
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
% The COBRAToolbox: testEFlux.m
% Aerobic and Anaerobic.txt are obtained from GEO (GSM1126445 and GSM1010240, respectively)
%
% Purpose:
%     - Tests the eFlux function
% Author:
%     - Original file: Thomas Pfau

currentDir = pwd;
% initialize the test
fileDir = fileparts(which('testEFlux.m'));
cd(fileDir);

% requires access to GEO to download data.
% Matlabs tolerances and precision seem incompatible for this function.
solverPkgs = prepareTest('needsLP',true,'toolboxes',{'bioinformatics_toolbox'},'needsWebAddress','https://www.ncbi.nlm.nih.gov/geo/query/','excludeSolvers',{'matlab','pdco'});

model = getDistributedModel('ecoli_core_model.mat');
model = removeGenesFromModel(model,'s0001','keepReactions',true);
geoaddress = 'https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?form=text&acc=%s&view=full';
anaerobic = geosoftread(urlread(sprintf(geoaddress,'GSM1010240')));
aerobic = geosoftread(urlread(sprintf(geoaddress,'GSM1126445')));
expressionControl.preprocessed = false;
expressionControl.value = cell2mat(aerobic.Data(:,2));
expressionControl.target = cellfun(@(x) x(1:5),aerobic.Data(:,1),'Uniform',false);
expressionCondition.preprocessed = false;
expressionCondition.value = cell2mat(anaerobic.Data(:,2));
expressionCondition.target = cellfun(@(x) x(1:5),anaerobic.Data(:,1),'Uniform',false);
% There is a potential cycle between reaction 44 and 89 in the model.
% an additional cycle exists for reactions 42,43,68,92, which we will also
% ignore.
% we will ignore those reactions for the assertions.
toAssert = true(numel(model.rxns),1);
toAssert([42,43,68,92,44,89]) = false;

% return to original directory
for k = 1:numel(solverPkgs.LP)
    changeCobraSolver(solverPkgs.LP{k},'LP');
    fprintf('Testing EFlux with %s ...\n',solverPkgs.LP{k})
    [fChangeOrig,~,solContOrig,solCondOrig] = eFlux(model,expressionControl,expressionCondition);
    %anaerobic growth less than aerobic.
    assert(fChangeOrig < 1);
    oxpos = ismember(model.rxns,'EX_o2(e)');
    % oxygen uptake in aerobic > than anaerobic (i.e. the value is smaller, exchange reaction).
    assert(solContOrig.x(oxpos) < solCondOrig.x(oxpos))
    % also test the soft bounds
    [fChangeSoft,~,solContSoft,solCondSoft] = eFlux(model,expressionControl,expressionCondition,'softBounds',true);
    %since this is punished with a value of 1, this is exactly the same as
    %the original solution.
    assert(abs(fChangeSoft-fChangeOrig) < 1e-4);

    assert(all(abs(solContOrig.x(toAssert)-solContSoft.x(toAssert)) < 1e-4));
    assert(all(abs(solCondOrig.x(toAssert)-solCondSoft.x(toAssert)) < 1e-4));
    % Different use of gpr rules
    [fChange,~,solCont,solCond] = eFlux(model,expressionControl,expressionCondition,'minSum',false);
    % results are still similar.
    assert(fChange < 1);
    oxpos = ismember(model.rxns,'EX_o2(e)');
    % oxygen uptake in aerobic > than anaerobic (i.e. the value is smaller, exchange reaction).
    assert(solCont.x(oxpos) < solCond.x(oxpos))

    %test this with a 1% noise.
    [fChangeNoise,stderr,solContNoise,solCondNoise] = eFlux(model,expressionControl,expressionCondition,'testNoise',true,'noiseStd',expressionControl.value*0.01);
    assert(abs(fChangeNoise-fChangeOrig) < 1e-4);
    assert(stderr > 0) % this should could be wrong at some point, but it is so unlikely....
end
cd(currentDir)

author: opencobra

repo: cobratoolbox


text:
function F = conj(F)
%CONJ   Complex conjugate of a CHEBFUN.
%   CONJ(F) is the complex conjugate of F.
%
% See also REAL, IMAG.

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

% Loop over the columns of F:
for j = 1:numel(F)

    % Conjugate the point values:
    F(j).pointValues = conj(F(j).pointValues);

    % Conjugate the FUNs:
    for k = 1:numel(F(j).funs)
        F(j).funs{k} = conj(F(j).funs{k});
    end
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function sift_arr = sp_normalize_sift(sift_arr)
% normalize SIFT descriptors (after Lowe)
% URL: http://www.cs.illinois.edu/homes/slazebni/research/SpatialPyramid.zip

% find indices of descriptors to be normalized (those whose norm is larger than 1)
tmp = sqrt(sum(sift_arr.^2, 2));
normalize_ind = find(tmp > 1);

sift_arr_norm = sift_arr(normalize_ind,:);
sift_arr_norm = sift_arr_norm ./ repmat(tmp(normalize_ind,:), [1 size(sift_arr,2)]);

% suppress large gradients
sift_arr_norm(find(sift_arr_norm > 0.2)) = 0.2;

% finally, renormalize to unit length
tmp = sqrt(sum(sift_arr_norm.^2, 2));
sift_arr_norm = sift_arr_norm ./ repmat(tmp, [1 size(sift_arr,2)]);

sift_arr(normalize_ind,:) = sift_arr_norm;

author: adikhosla

repo: feature-extraction


text:
load('table_T.mat')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Params for UAV
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% UAV == Zagi
%     p_drone.mass = 1.56;
%     p_drone.Jx   = 0.1147;
%     p_drone.Jy   = 0.0576;
%     p_drone.Jz   = 0.1712;
%     p_drone.Jxz  = 0.0015;

%%% UAV == Aerosonde
p_drone.mass = 13.5;
p_drone.Jx   = 0.8244;
p_drone.Jy   = 1.135;
p_drone.Jz   = 1.759;
p_drone.Jxz  = .1204;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute gammas
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

p_drone.gamma  = p_drone.Jx*p_drone.Jz - p_drone.Jxz^2;
p_drone.gamma1 = p_drone.Jxz*(p_drone.Jx-p_drone.Jy+p_drone.Jz)/p_drone.gamma;
p_drone.gamma2 = (p_drone.Jz*(p_drone.Jz-p_drone.Jy)+p_drone.Jxz^2)/p_drone.gamma;
p_drone.gamma3 = p_drone.Jz/p_drone.gamma;
p_drone.gamma4 = p_drone.Jxz/p_drone.gamma;
p_drone.gamma5 = (p_drone.Jz-p_drone.Jx)/p_drone.Jy;
p_drone.gamma6 = p_drone.Jxz/p_drone.Jy;
p_drone.gamma7 = ((p_drone.Jx-p_drone.Jy)*p_drone.Jx+p_drone.Jxz^2)/p_drone.gamma;
p_drone.gamma8 = p_drone.Jx/p_drone.gamma;

% Aerodynamic coefficients
p_drone.S_wing          = 0.55;     % Surface of the plane
p_drone.b              = 2.8956;
p_drone.c              = 0.18994;
p_drone.S_prop         = 0.2027;    % Surface of the propeller
p_drone.k_motor        = 80;        % k motor
p_drone.k_TP           = 0;
p_drone.k_omega        = 0;
p_drone.e              = 0.9;

p_drone.kb = @(Va) 0.5*p_drone.b/Va;
p_drone.kc = @(Va) 0.5*p_drone.c/Va;

% Longitudinal dynamics:
%   lift and drag (F_L, F_D)
p_drone.CL0            = 0.28;     % C_L_0
p_drone.CL_alpha       = 3.45;     % C_L_alpha
p_drone.CL_q           = 0.0;      % C_L_q
p_drone.CL_de          = -0.36;    % C_L_delta_e
p_drone.CD0            = 0.03;     % C_D_0
p_drone.CD_alpha       = 0.30;     % C_D_alpha
p_drone.CD_p           = 0.0437;   % C_D_p
p_drone.CD_q           = 0.0;      % C_D_q
p_drone.CD_de          = 0.0;      % C_D_delta_e

%   moment around y (m)
p_drone.Cm0            = -0.02338; % C_m_0
p_drone.Cm_alpha       = -0.38;    % C_m_alpha
p_drone.Cm_q           = -3.6;     % C_m_q
p_drone.Cm_de          = -0.5;     % C_m_delta_e

% Lateral dynamics
%   F_y
p_drone.CY0            = 0.0;
p_drone.CY_beta        = -0.98;
p_drone.CY_p           = 0.0;
p_drone.CY_r           = 0.0;      % C_Y_r
p_drone.CY_da          = 0.0;      % C_Y_delta_a
p_drone.CY_dr          = -0.17;    % C_Y_delta_r

%   moment around x (l) and z (n)
p_drone.Cl0            = 0.0;      % C_l_0
p_drone.Cl_beta        = -0.12;    % Cl_beta
p_drone.Cl_p           = -0.26;    % C_l_p
p_drone.Cl_r           = 0.14;     % C_l_r
p_drone.Cl_da          = 0.08;     % C_l_delta_aileron
p_drone.Cl_dr          = 0.105;    % C_l_delta_rudder
p_drone.Cn0            = 0.0;      % C_n_0
p_drone.Cn_beta        = 0.25;     % C_n_beta
p_drone.Cn_p           = 0.022;    % C_n_p
p_drone.Cn_r           = -0.35;    % C_n_r
p_drone.Cn_da          = 0.06;     % C_n_delta_a
p_drone.Cn_dr          = -0.032;   % C_n_delta_r

p_drone.C_prop         = 1.0;      % propeller proportional coeff fx=Sprop*Cprop*dp
p_drone.M              = 50;
p_drone.epsilon        = 0.1592;
p_drone.alpha0         = 0.4712;

p_drone.Cp0     = p_drone.gamma3 * p_drone.Cl0     + p_drone.gamma4 * p_drone.Cn0;
p_drone.Cp_beta = p_drone.gamma3 * p_drone.Cl_beta + p_drone.gamma4 * p_drone.Cn_beta;
p_drone.Cp_p    = p_drone.gamma3 * p_drone.Cl_p    + p_drone.gamma4 * p_drone.Cn_p;
p_drone.Cp_r    = p_drone.gamma3 * p_drone.Cl_r    + p_drone.gamma4 * p_drone.Cn_r;
p_drone.Cp_da   = p_drone.gamma3 * p_drone.Cl_da   + p_drone.gamma4 * p_drone.Cn_da;
p_drone.Cp_dr   = p_drone.gamma3 * p_drone.Cl_dr   + p_drone.gamma4 * p_drone.Cn_dr;
p_drone.Cr0     = p_drone.gamma4 * p_drone.Cl0     + p_drone.gamma8 * p_drone.Cn0;
p_drone.Cr_b    = p_drone.gamma4 * p_drone.Cl_beta + p_drone.gamma8 * p_drone.Cn_beta;
p_drone.Cr_p    = p_drone.gamma4 * p_drone.Cl_p    + p_drone.gamma8 * p_drone.Cn_p;
p_drone.Cr_r    = p_drone.gamma4 * p_drone.Cl_r    + p_drone.gamma8 * p_drone.Cn_r;
p_drone.Cr_da   = p_drone.gamma4 * p_drone.Cl_da   + p_drone.gamma8 * p_drone.Cn_da;
p_drone.Cr_dr   = p_drone.gamma4 * p_drone.Cl_dr   + p_drone.gamma8 * p_drone.Cn_dr;

% first cut at initial conditions
% p_drone.pn_T    = xyu_T(1);  % initial North position
% p_drone.pe_T    = xyu_T(2);  % initial East position
% p_drone.pd_T    = xyu_T(3);  % initial Down position (negative altitude)
% p_drone.vx_T    = xyu_T(4); % initial velocity along body x-axis
% p_drone.vy_T    = xyu_T(5);  % initial velocity along body y-axis
% p_drone.vz_T    = xyu_T(6);  % initial velocity along body z-axis
% p_drone.phi_T   = xyu_T(7);  % initial roll angle
% p_drone.theta_T = xyu_T(8);  % initial pitch angle
% p_drone.psi_T   = xyu_T(9);  % initial yaw angle
% p_drone.p_T     = xyu_T(10);  % initial body frame roll rate
% p_drone.q_T     = xyu_T(11);  % initial body frame pitch rate
% p_drone.r_T     = xyu_T(12);  % initial body frame yaw rate

% p_drone.Va_T = xyu_T(13);
% p_drone.alpha_T = xyu_T(14);
% p_drone.beta_T = xyu_T(15);

% p_drone.de_T = xyu_T(16);
% p_drone.da_T = xyu_T(17);
% p_drone.dr_T = xyu_T(18);
% p_drone.dt_T = xyu_T(19);

% Values for the Aerosonde UAV
p_drone.Va0 = 35;        % m/s (~85 mph)

% Trim conditions used for static PID tuning
Va_Td = 35;
gamma_Td = deg2rad(0);
R_Td = Inf;

Td = [Va_Td, gamma_Td, R_Td];

% Compute PID control
% p_drone.C = compute_pid_controller(Td, P)
[Va_idx, gamma_idx, R_idx]  = Td2idx(Td);
xyu_T = table_T(:,Va_idx, gamma_idx, R_idx);
p_drone.phi_T   = 0;  % initial roll angle
p_drone.theta_T = 0;  % initial pitch angle
p_drone.psi_T   = 0;  % initial yaw angle
p_drone.Va_T = 35;
[a,TF] = compute_tf_model(xyu_T, p_drone, p_physics);
% PVa_T = xyu_T(13);  % Extract trimmed state

% Altitude state machine parameters
p_drone.altitude_take_off_zone = 40;
p_drone.h_theta_hold_zone      = 10;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Tuning PID loops

% Lateral dynamics
% Roll_aileron-hold
ephi_max = deg2rad(15);
p_drone.kp_phi = p_drone.da_max / ephi_max;
ksi_phi = KSI_OPT;
wn_phi = sqrt(p_drone.kp_phi * a.phi2);
p_drone.kd_phi = (2*ksi_phi*wn_phi - a.phi1) / a.phi2;

p_drone.tau_phi = 0.05; % [seconds] <=> 1 autopilot loop

% Course_roll-hold
W_chi_phi = 10; % Bandwith separation factor
% Between 5 and 10. Safety corresponds to 10.
wn_chi = wn_phi / W_chi_phi; %[rad/sec]
ksi_chi = KSI_OPT;
p_drone.kp_chi = 2*ksi_chi*wn_chi*p_drone.Va_T/p_physics.gravity;
p_drone.ki_chi = wn_chi^2*p_drone.Va_T/p_physics.gravity;

% Sideslip_rudder-hold
%TODO

% Longitudinal dynamics
% Pitch_elevator-hold
etheta_max  = deg2rad(10);
p_drone.kp_theta  = p_drone.de_max / etheta_max * sign(a.theta3);
wn_theta = sqrt(a.theta2 + p_drone.kp_theta*a.theta3);
ksi_theta = KSI_OPT;
if isfinite(a.theta3)
    p_drone.kd_theta  = (2*ksi_theta*wn_theta - a.theta1) / a.theta3;
else
    p_drone.kd_theta  = 0;
    disp("kd_theta is null!");
end
p_drone.tau_theta = 0.1; % [seconds] <=> 10 autopilot loops

% Altitude_pitch-hold
% Scale factor between theta_c and theta
K_theta_DC = p_drone.kp_theta*a.theta3/(a.theta2+p_drone.kp_theta*a.theta3);
W_h_theta = 10; % Between 5 and 15. Safety corresponds to 15.
wn_h = wn_theta / W_h_theta; % [rad/sec]
ksi_h = 1;
p_drone.kp_h  = 2*ksi_h*wn_h/(K_theta_DC*p_drone.Va_T);
p_drone.ki_h  = wn_h^2 / (K_theta_DC*p_drone.Va_T);

% Airspeed_thrust-hold
wn_va_dt = 10; % [rad/sec] Need to be tuned.
ksi_va_dt = 1;
p_drone.kp_va_dt = (2*ksi_va_dt*wn_va_dt - a.va1)/a.va2;
p_drone.ki_va_dt = wn_va_dt^2/a.va2;

% Airspeed_pitch-hold
W_va_theta = 7; % [rad/sec] 10 by default, can be lower.
wn_va_theta = wn_theta / W_va_theta;
ksi_va_theta = KSI_OPT;
p_drone.kp_va_theta = (a.va1 - 2*ksi_va_theta*wn_va_theta)/(K_theta_DC*p_physics.gravity);
p_drone.ki_va_theta = -wn_va_theta^2/(K_theta_DC*p_physics.gravity);

% first cut at initial conditions
p_drone.pn_T    = xyu_T(1);  % initial North position
p_drone.pe_T    = xyu_T(2);  % initial East position
p_drone.pd_T    = xyu_T(3)-200;  % initial Down position (negative altitude)
p_drone.vx_T    = xyu_T(4); % initial velocity along body x-axis
p_drone.vy_T    = xyu_T(5);  % initial velocity along body y-axis
p_drone.vz_T    = xyu_T(6);  % initial velocity along body z-axis
p_drone.phi_T   = xyu_T(7);  % initial roll angle
p_drone.theta_T = xyu_T(8);  % initial pitch angle
p_drone.psi_T   = xyu_T(9);  % initial yaw angle
p_drone.p_T     = xyu_T(10);  % initial body frame roll rate
p_drone.q_T     = xyu_T(11);  % initial body frame pitch rate
p_drone.r_T     = xyu_T(12);  % initial body frame yaw rate

p_drone.Va_T    = xyu_T(13);
p_drone.alpha_T = xyu_T(14);
p_drone.beta_T  = xyu_T(15);

p_drone.de_T    = xyu_T(16);
p_drone.da_T    = xyu_T(17);
p_drone.dr_T    = xyu_T(18);
p_drone.dt_T    = xyu_T(19);

run('param_sensors')

author: lis-epfl

repo: swarmlab

save_path: github-repos/MATLAB/lis-epfl-swarmlab


text:
function linplus_test385 ( )

%*****************************************************************************80
%
%% TEST385 tests R8GE_PLU.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 March 2009
%
%  Author:
%
%    John Burkardt
%
  m = 5;
  n = 4;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST385\n' );
  fprintf ( 1, '  For a matrix in general storage,\n' );
  fprintf ( 1, '  R8GE_PLU returns the PLU factors of a matrix.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix rows M =    %d\n', m );
  fprintf ( 1, '  Matrix columns N = %d\n', n );

  a = r8ge_random ( m, n, seed );

  r8ge_print ( m, n, a, '  Matrix A:' );
%
%  Compute the PLU factors.
%
  [ p, l, u ] = r8ge_plu ( m, n, a );

  r8ge_print ( m, m, p, '  Factor P:' );

  r8ge_print ( m, m, l, '  Factor L:' );

  r8ge_print ( m, n, u, '  Factor U:' );

  plu(1:m,1:n) = p(1:m,1:m) * ( l(1:m,1:m) * u(1:m,1:n) );

  r8ge_print ( m, n, plu, '  Product P*L*U:');

  return
end


text:
function sR=icassoExp(sR)
%function sR=icassoExp(sR)
%
%PURPOSE
%
%To prepare Icasso result structure for exploratory analysis, i.e.,
%to compute (dis)similarity matrix, clustering, and projection.
%
%EXAMPLES OF BASIC USAGE
%
%First we produce an Icasso result structure...
%
%   load megdata
%   sR=icassoEst('randinit',megdata,15,'lastEig',20);
%
%The following command performs the default Icasso clustering
%procedure:
%
%   sR=icassoExp(sR);
%
%The next step would be to return results launch the visualizations
%of the results: See icassoShow, icassoViz, and icassoResult.
%
%You can customize the Icasso procedure by using this script
%as a model but changing the optional input parameters. See also
%icassoCluster and icassoProjection.
%
%DETAILS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%icassoExp performs two functions icassoCluster and icassoProjection:
%
%Step 1
%
%Cluster the estimates
%
%Details:
% 1. Compute similarities (S) between the estimates
% 2. Store the similarity matrix S in field sR.cluster.similarity.
% 3. Record the similarity function into sR.cluster.simfcn, in
%    this case, 'abscorr'. Other possibility is 'power'. You can
%    also specify an explicit similarity matrix
% 4. Compute dissimilarities D from S using function
%    sim2dis (simply D=1-S). If you want something else, you can give another
%    function name
% 5. Compute hierarchical clustering according to the readily computed
%    dissimilarities using average-linkage strategy. Other
%    possibilities: 'CL' (comptele link) 'SL' (single-link)
% 6. Store the clustering into field sR.cluster.partition and the
%    strategy in sR.cluster.strategy
% 7. Compute a clustering solution validity index, R-index up to L
%    clusters (according to the partition and dissimilarities D)
%    ('rdim' sets  L equal to the (reduced) data dimension)
% 9. Store the index into field sR.cluster.index.R
%
%   sR=icassoCluster(sR,'strategy','AL','simfcn','abscorr','s2d','sim2dis','L','rdim');
%
%Step 2
%
%Compute the coordinates for correlation graph using Curvilinear
%Component Analysis. Other possibilities: 'mmds' (principal
%coordinates), 'sammon' (Sammon's projection). See icassoProjection
%for details.
%
%Note that the dissimilarities are slightly differently scaled than in clustering:
%
%Details:
%1. Compute similarity-to-dissimilarity transformation that is
%   found to be good for the proximity preserving projection:
%   D=sqrt(1-sR.cluster.similarity) (i.e., call function sqrtsim2dis)
%2. compute CCA on this dissimilarity matrix using 75 epochs and default
%   parameters (see icassoProjection).
%
%  sR=icassoProjection(sR,'cca','s2d','sqrtsim2dis','epochs',75);

%COPYRIGHT NOTICE
%This function is a part of Icasso software library
%Copyright (C) 2003-2005 Johan Himberg
%
%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 any later version.
%
%This program is distributed in the hope that it will be useful,
%but WITHOUT ANY WARRANTY; without even the implied warranty of
%MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%GNU General Public License for more details.
%
%You should have received a copy of the GNU General Public License
%along with this program; if not, write to the Free Software
%Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.

% ver 1.21 040305

sR=icassoCluster(sR,'strategy','AL','simfcn','abscorr','s2d','sim2dis','L','rdim');
sR=icassoProjection(sR,'cca','s2d','sqrtsim2dis','epochs',75);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%CLASSC Convert classifier to normalized classifier (yielding confidences)
%
%  V = CLASSC(W)
%  V = W*CLASSC
%  D = CLASSC(A*W)
%  D = A*W*CLASSC
%  D = CLASSC(A,W)
%
% INPUT
%  W  Trained or untrained classifier
%  A  Dataset
%
% OUTPUT
%  V  Normalized classifier producing confidences instead of
%     densities or distances (after training if W is untrained)
%
% DESCRIPTION
% The trained or untrained classifier W may yield densities or unnormalised
% confidences. The latter holds for two-class discriminants like FISHERC
% and SVC as well as for neural networks. Such classifiers use or should
% use CNORMC to convert distances to confidences. In multi-class problems
% as well as in combining schemes they do not produce normalises
% confidences. These outcomes, like the density outcomes of classifiers
% liek QDC, LDC and PARZENC, can be converted by CLASSC into confidences:
% the sum of the outcomes will be one for every object.
%
% In case W is a one-dimensional mapping, it is converted into a two-class
% classifier, provided that during the construction a class label was
% supplied. If not, the mapping cannot be converted and an error is
% generated.
%
% CLASSC lists the outcomes on the screen in case no output argument is
% supplied. Also true and estimated labels are supplied.
%
% SEE ALSO (PRTools Guide)
% MAPPINGS, DATASETS, CNORMC, LABELD

% Copyright: R.P.W. Duin, r.p.w.duin@37steps.com
% Faculty EWI, Delft University of Technology
% P.O. Box 5031, 2600 GA Delft, The Netherlands

% $Id: classc.m,v 1.5 2010/02/23 15:21:54 duin Exp $

function w = classc(w,flag)

		if nargin < 2, flag = 0; end % flag forces non=combiner behavior avoiding recursion

	if (nargin == 0)

		% Untrained mapping.
		w = prmapping('classc','combiner',flag);

	elseif (ismapping(w))

		% If mapping is stacked or parallel, recurse over the individual
		% sub-mappings and call CLASSC for each of them.
		if ((isstacked(w)) | (isparallel(w))) & (flag == 0)
			v = cell(1,length(w.data));
			for j = 1:length(w.data)
				if ismapping(w.data{j}) % the parallel combiner may have nonmapping data
					v{j} = feval(mfilename,w.data{j});
				else
					v{j} = w.data{j};
				end
			end
			w = setdata(w,v);
			w = feval(mfilename,w,1);  % and here CLASSC is called for the combiner avoiding recursion
		else
			conv = get(w,'out_conv');
			if (conv < 1)
				% Set the "normalization" bit in the mapping's output conversion flag
				w = set(w,'out_conv',conv+2);
			end;
		end
	elseif (isdataset(w))
		if ismapping(flag)
			if nargout == 1
				w = feval(mfilename,w*flag);
			else
				feval(mfilename,w*flag);
				clear w;
			end
			return
		end
		w = w*normm;
		w = w*costm;
		if nargout == 0 % list outcomes on the screen
			ww = +w;
			ss = repmat('-',1,9*size(ww,2));
			fprintf('\n True   Estimated        Class \nLabels    Labels      Confidences\n');
			fprintf('------------------%s\n',ss);
			nlab = getnlab(w);
			[wmax,K] = max(ww,[],2);
			lablist = getlablist(w);
			if ~isempty(lablist) & ~ischar(lablist)
				nlab = lablist(nlab);
				K = lablist(K);
			end
			for j=1:size(ww,1)
				if (nlab(j) ~= K(j))
					fprintf(' %3.0f    ->%3.0f    ',nlab(j),K(j));
				else
					fprintf(' %3.0f      %3.0f    ',nlab(j),K(j));
				end
				fprintf('  %7.4f',ww(j,:));
				fprintf('\n');
			end
			lablist = getlablist(w);
			if ischar(lablist)
				fprintf('\n');
				for j=1:size(lablist,1)
					fprintf('   %i %s\n',j,lablist(j,:));
				end
			end

			clear w;
		end
	else
		error('input should be mapping or dataset');
	end

return

author: marianux

repo: ecg-kit

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


text:
function lambda = rutis4_eigenvalues ( n )

%*****************************************************************************80
%
%% RUTIS4_EIGENVALUES returns the eigenvalues of the RUTIS4 matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Output, real LAMBDA(N,1), the eigenvalues.
%
  lambda = zeros ( n, 1 );

  for i = 1 : n
    angle = ( i ) * pi / ( 2 * ( n + 1 ) );
    lambda(i,1) = 64.0 * ( cos ( angle ) )^6;
  end

  return
end


text:
% labels01 = (labelsPM+1)/2; % maps -1->0, +1->1
% labelsPM = (2*labels01)-1; % maps 0,1 -> -1,1


text:
classdef CompositeResultsPrinter < ResultsPrinter

    properties (Access = protected)
        hasGaussData
        printers
        simulationStr
    end

    methods (Access = public)

        function  printResults(obj,iter,fileID)
            obj.setSimulationStrToPrinters()
            for iprinter = 1:numel(obj.printers)
                p = obj.printers{iprinter};
                p.printResults(iter,fileID);
            end
        end

    end

    methods (Access = protected)

        function init(obj,d)
            obj.init@ResultsPrinter(d);
            obj.createPrinters(d);
            obj.obtainHasGaussData();
        end

        function createHeadPrinter(obj,d,dh)
            for iprinter = 1:numel(obj.printers)
                p = obj.printers{iprinter};
                p.createHeadPrinter(d,dh);
                if p.getHasGaussData()
                    h = p.getHeadPrinter();
                    obj.headPrinter = h;
                    return
                else
                    h = p.getHeadPrinter();
                    obj.headPrinter = h;
                end
            end
        end

        function storeFieldsToPrint(obj,d)
            for iprinter = 1:numel(obj.printers)
                p = obj.printers{iprinter};
                p.storeFieldsToPrint(d);
            end
        end

    end

    methods (Access = private)

        function obtainHasGaussData(obj)
            obj.hasGaussData = false;
            for iprinter = 1:numel(obj.printers)
                p = obj.printers{iprinter};
                hG = p.getHasGaussData();
                obj.hasGaussData = hG || obj.hasGaussData;
            end
        end

        function setSimulationStrToPrinters(obj)
            s = obj.simulationStr;
            for iprinter = 1:numel(obj.printers)
                p = obj.printers{iprinter};
                p.setSimulationStr(s);
            end
        end
    end

    methods (Access = protected, Abstract)
        createPrinters(obj,d)
    end

end


text:
function val = viewGetRetinotopy(vw,param,varargin)
% Get data from various view structures
%
% This function is wrapped by viewGet. It should not be called by anything
% else other than viewGet.
%
% This function retrieves information from the view that relates to a
% specific component of the application.
%
% We assume that input comes to us already fixed and does not need to be
% formatted again.

if notDefined('vw'), vw = getCurView; end
if notDefined('param'), error('No parameter defined'); end

mrGlobals;
val = [];

switch param

    case 'framestouse'
        % Return a vector of time frames in the current or specified
        % scan to be used for coranal (block) analyses
        %   frames = viewGet(vw,'frames to use');
        %   scan = 1; frames = viewGet(vw,'frames to use',scan);
        if isempty(varargin) || isempty(varargin{1})
            scan = viewGet(vw, 'CurScan');
        else
            scan = varargin{1};
        end
        dt         = viewGet(vw, 'dtStruct');
        blockParms = dtGet(dt,'bparms',scan);
        if checkfields(blockParms, 'framesToUse')
            val = blockParms.framesToUse;
        else
            val = 1: viewGet(vw,'nFrames',scan);
        end

    case 'rmfile'
        % Return the path to the currently loaded retinotopy model.
        %   rmFile = viewGet(vw, 'retinotopy model file');
        if checkfields(vw, 'rm', 'retinotopyModelFile')
            val = vw.rm.retinotopyModelFile;
        else
            val = [];
        end
    case 'rmmodel'
        % Return the currently loaded retinotopy model struct.
        %   rm = viewGet(vw, 'retinotopy model');
        if checkfields(vw, 'rm', 'retinotopyModels')
            val = vw.rm.retinotopyModels;
        else
            val = [];
        end

    case 'rmcurrent'
        % Return the currently selected retinotopy model struct. Note that
        % there may be multiple models loaded.
        %   rm = viewGet(vw, 'rm current model');
        if checkfields(vw, 'rm', 'retinotopyModels')
            val = vw.rm.retinotopyModels{ viewGet(vw, 'rmModelNum') };
        else
            val = [];
        end
    case 'rmmodelnames'
        % Return the description of currently loaded retinotopy models.
        %   models = viewGet(vw, 'rm model names');
        %   models = viewGet(vw, 'retinotopy model names');
        models = viewGet(vw, 'Retinotopy Model');
        val = cell(1,numel(models));
        for n = 1:numel(models)
            val{n} = rmGet(models{n},'description');
        end
    case 'rmparams'
        % Return the retinotopy model parameters.
        %   params = viewGet(vw, 'Retinotopy Parameters');
        if checkfields(vw,'rm')
            val = vw.rm.retinotopyParams;
        end
    case 'rmstimparams'
        % Return the retinotopy model stimulus parameters. This is a subset
        % of the retinopy model parameters.
        %   stimParams = viewGet(vw, 'RM Stimulus Parameters');
        if checkfields(vw,'rm','retinotopyParams','stim')
            val = vw.rm.retinotopyParams.stim;
        end
    case 'rmmodelnum'
        % Return the retinotopy model number that is currently selected.
        % (There may be more than one model loaded.)
        %   modelNum = viewGet(vw, 'Retinotopy Model Number');
        if checkfields(vw, 'rm', 'modelNum') && ...
                ~isempty(vw.rm.modelNum)
            val = vw.rm.modelNum;
        else
            val = rmSelectModelNum(vw);
        end
    case 'rmhrf'
        % Return the hrf struct for the current retinopy model. This struct
        % contains a descriptor (such as 'two gammas (SPM style)') and the
        % parameters associated with this function.
        %   rmhrf = viewGet(vw, 'Retinotopy model HRF');
        if checkfields(vw,'rm','retinotopyParams','stim')
            val1 = vw.rm.retinotopyParams.stim.hrfType;
            switch(lower(val1))
                case {'one gamma (boynton style)','o','one gamma' 'b' 'boynton'}
                    val2 = vw.rm.retinotopyParams.stim.hrfParams{1};
                case {'two gammas (spm style)' 't' 'two gammas' 'spm'}
                    val2 = vw.rm.retinotopyParams.stim.hrfParams{2};
                case {'impulse' 'no hrf' 'none'}
                    val2 = [];
                otherwise
                    val2 = [];
            end
            val = {val1 val2};
        end

    otherwise
        error('Unknown viewGet parameter');

end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% Gnufft_test.m
% Test the Gnufft object (vs exact Gdsft)

%% create Gnufft class object
if 1 || ~isvar('A'), printm 'setup Gnufft_test'
	im plc 3 3
	if 1 % 2d
		N = [32 30];
		J = [6 7];
	%	N = [32 1]; J = [5 1];
		omega = linspace(0, 10*2*pi, 201)'; % crude spiral:
		omega = pi*[cos(omega) sin(omega)].*omega(:,[1 1])/max(omega);
		if im, im subplot 1, plot(omega(:,1), omega(:,2), '.'), end
	else % 3d
		N = [16 12 14];
		J = [6 7 5];
		tmp = linspace(0, 10*2*pi, 201)';
		omega = pi*[cos(tmp) sin(tmp)].*tmp(:,[1 1])/max(tmp);
		omega(:,3) = pi * tmp / max(tmp); % spiral cone
		if im, im subplot 1, plot3(omega(:,1), omega(:,2), omega(:,3), '.'), end
	end

	K = 2*N;
	if N(2) == 1
		args = {omega(:,1), N(1), J(1), K(1)}; omega(:,2) = 0;
	else
		args = {omega, N, J, K};
	end

%	mask = true(N); mask(1,1) = false;
% todo: mask!

	args = {args{:}, 'table', 2^10, 'minmax:kb'}; % test with table

	A = Gnufft(args);

	Ad = Gdsft(omega, N);
end

%% test save/load
if 0
	save('/tmp/A.mat', 'A')
return
end

%% test data
if 1 || ~isvar('x'), printm 'setup data'
	if length(N) == 2
		x = zeros(N);
		x(5:25,10:25) = 1;
		x(15:20,15:20) = 2;
		x(15,5) = 2;
		if N(2) == 1
			x = x(:,5,1);
		end
	else
		rng(0)
		x = rand(N);
	end
	im(2, x, '\x')

	yd = Ad * x;

	if length(N) == 2
		n1 = ([0:N(1)-1]/N(1) - 0.5)*2*pi;
		if N(2) == 1 % 1D case
			yd_g = interp1(omega(:,1), yd, n1);
		else
			n2 = ([0:N(2)-1]/N(2) - 0.5)*2*pi;
			[nn1 nn2] = ndgrid(n1, n2);
			yd_g = griddata(omega(:,1), omega(:,2), yd,  nn1, nn2);
			yd_g(isnan(yd_g)) = 0;
		end
		im(3, abs(yd_g), '$|\y_d|$'), cbar
	end
end

%% build gram
if 1, printm 'Gnufft gram'
	wi = [1:size(omega,1)]';
	T = build_gram(A, wi);
	y2 = T * x(:);
	y1 = A' * (wi .* (A * x(:)));
	max_percent_diff y1 y2
%	equivs(y1, y2)
prompt
end

%% compare forward
if 1, printm 'forward'
	yn = A * [x(:) x(:)]; % test with two
	yn = yn(:,1);

	if length(N) == 2
		if N(2) == 1 % 1D case
			yn_g = interp1(omega(:,1), yn, n1);
		else
			yn_g = griddata(omega(:,1), omega(:,2), yn,  nn1, nn2);
			yn_g(isnan(yn_g)) = 0;
		end

		im(4, abs(yn_g), '$|\y_g|$'), cbar
		im(5, abs(yd_g - yn_g), '$|\y_g - \y_d|$'), cbar
	end
	max_percent_diff yd yn
end

%% compare adjoint
if 1, printm 'adjoint'
	yb = ones(size(omega,1), 1);
	xd = Ad' * yb;
	xd = iembed(Ad, xd);
	xn = A' * [yb yb]; % test with two
	xn = xn(:,1);
	xn = iembed(A, xn);

	im(7, fftshift(abs(xd)), '|back dtft|'), cbar
	im(8, fftshift(abs(xn)), '|back nufft|'), cbar
	im(9, fftshift(abs(xn-xd)), '|back err|'), cbar
	max_percent_diff xd xn
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function C = hcp_lattice(I,J,K)
  % HCP_LATTICE generate the center locations of a hexagonal close-packed
  % lattice of unit-radius spheres.
  %
  % Inputs:
  %   I  #I list of x-coordinate indices
  %   J  #I list of y-coordinate indices
  %   K  #I list of z-coordinate indices
  % Outputs:
  %   C  #I by 3 list of center locations
  %
  % Example:
  %   [I,J,K] = meshgrid(0:9,0:9,0:9);
  %   C = hcp_lattice(I(:),J(:),K(:));
  %   [SV,SF] = lloyd_sphere(300);
  %   [V,F] = repmesh(SV,SF,C);
  i = I(:);
  j = J(:);
  k = K(:);
  % https://en.wikipedia.org/wiki/Close-packing_of_equal_spheres#Simple_hcp_lattice
  C = [2*i + mod(j+k,2), sqrt(3)*(j+1/3*mod(k,2)), 2*sqrt(6)/3*k];
end

author: alecjacobson

repo: gptoolbox


text:
function [node,elem,T] = squarequadmesh(square,h)
%% SQUAREMESH uniform mesh of a square
%
% [node,elem] = squarequadmesh([x0,x1,y0,y1],h) generates a uniform mesh of the
% rectangle [x0,x1]*[y0,y1] with mesh size h.
%
% Example
%
%   [node,elem] = squarequadmesh([0,1,0,2],0.5);
%   showmesh(node,elem);
%   findnode(node);
%   findquadelem(node,elem);
%
% See also: squaregradmeshquad, cubehexmesh
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

%% Generate nodes
x0 = square(1); x1 = square(2);
y0 = square(3); y1 = square(4);
% [x,y] = meshgrid(x0:h:x1,y1:-h:y0);
[x,y] = ndgrid(x0:h:x1,y0:h:y1);
node = [x(:),y(:)];

%% Generate elements
ni = size(x,1); % number of rows
nj = size(x,2);
N = size(node,1);
nodeidx = reshape(1:N,ni,nj);
t2nidxMap = nodeidx(1:ni-1,1:nj-1);
k = t2nidxMap(:);
elem = [k k+ni k+ni+1 k+1];
NT = size(elem,1);
% 4 k+1 --- k+ni+1 3
%    |        |
% 1  k  ---  k+ni  2

%% Check if further structure is needed
if nargout<=2
    return;
end

%% Generate edges
Ne1 = (ni-1)*nj;
Ne2 = ni*(nj-1);
edge = zeros(Ne1+Ne2,2,'int32');
% vertical edges
e2nidxMap = nodeidx(1:ni-1,1:nj);
k = e2nidxMap(:);
edge(1:Ne1,:) = [k k+1];
% horizontal edges
e2nidxMap = nodeidx(1:ni,1:nj-1);
k = e2nidxMap(:);
edge(Ne1+(1:Ne2),:) = [k k+ni];

%% Generate index map between element and edges
elem2edge = zeros(NT,4,'int32');
edge2elem = zeros(Ne1+Ne2,2,'int32');
k = 1:NT;
veidx = k; % veidx = sub2ind([ni-1,nj],i,j);
heidx = Ne1 + (1:Ne2); % linear indexing
heidx = reshape(heidx,ni,nj-1); % matrix indexing
heidx = heidx(1:ni-1,:); % without top edges
elem2edge(k,1) = heidx(:);
elem2edge(k,3) = heidx(:)+1;
elem2edge(k,2) = veidx + (ni-1);
elem2edge(k,4) = veidx;
% vertical edges
edge2elem(veidx,2) = k;   % right element
edge2elem(veidx+ni-1,1) = k; % left element
% horizontal edges
edge2elem(heidx,1) = k;
edge2elem(heidx+1,2) = k;

%% Generate boundary edges
% for boundary edges, set the another element idx to zero
idx1 = edge2elem(:,1) == 0; % left and top
idx2 = edge2elem(:,2) == 0;
bdEdge = [edge(idx1,[2 1]); edge(idx2,:)];
% switch edge index on the left and top to have a consistent orientation

%% Auxstructure
T = struct('edge',edge,'elem2edge',elem2edge,'edge2elem',edge2elem,'bdEdge',bdEdge);
% neighbor and bdElem can be easily get from 2-D index of all elements

repo: Aero-matlab


text:
function hypersphere_integrals_test ( )

%*****************************************************************************80
%
%% HYPERSPHERE_INTEGRALS_TEST tests the HYPERSPHERE_INTEGRALS library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 January 2014
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'HYPERSPHERE_INTEGRALS_TEST\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the HYPERSPHERE_INTEGRALS library.\n' );

  hypersphere_integrals_test01 ( );
  hypersphere_integrals_test02 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'HYPERSPHERE_INTEGRALS_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end

author: johannesgerer

repo: jburkardt-m