text:
vI = cat(3, 255, 102, 51) / 255;
vI = cat(3, 255, 153, 51) / 255;
vI = cat(3, 255, 128, 0) / 255;
vI = cat(3, 204, 127, 102) / 255;

vI = cat(3, 255, 0, 0) / 255;
vI = cat(3, 159, 112, 96) / 255;

% Reds, Yellows, Greens, Cyans, Blues, Magentas
vPhotoshopValues = [25; 75; 35; 65; 45; 55]; %<! Photoshop Neutral

vPhotoshopValues = randi([-200, 300], [6, 1]);
% vPhotoshopValues([1, 2]) = 50;

vPhotoshopValues = [60; 50; 50; 50; 50; 50];

vCoeffValues = (vPhotoshopValues - 50) ./ 50;
mO = ApplyBlackWhiteFilter(vI, vCoeffValues);

round(mO * 255)

% vI = cat(3, 255, 102, 51) / 255;
% round(ConvertRgbToHsl(vI) * 360)
% vI = cat(3, 240, 110, 66) / 255;
% round(ConvertRgbToHsl(vI) * 360)

author: RoyiAvital

repo: StackExchangeCodes


text:
function y = sample_node(CPD, pev)
% SAMPLE_NODE Draw a random sample from P(Xi | x(pi_i), theta_i)  (tabular)
% y = sample_node(CPD, pev)
%
% pev{i} is the value of the i'th parent (if any)

%assert(~any(isemptycell(pev)));

%CPT = CPD_to_CPT(CPD);
%sz = mysize(CPT);
sz = CPD.sizes;
nparents = length(sz)-1;
if nparents > 0
  pvals = cat(1, pev{:});
end
switch nparents
 case 0, T = CPD.CPT;
 case 1, T = CPD.CPT(pvals(1), :);
 case 2, T = CPD.CPT(pvals(1), pvals(2), :);
 case 3, T = CPD.CPT(pvals(1), pvals(2), pvals(3), :);
 case 4, T = CPD.CPT(pvals(1), pvals(2), pvals(3), pvals(4), :);
 otherwise,
  psz = sz(1:end-1);
  ssz = sz(end);
  i = subv2ind(psz, pvals(:)');
  T = reshape(CPD.CPT, [prod(psz) ssz]);
  T = T(i,:);
end

if sz(end)==2
  r = rand(1,1);
  if r > T(1)
    y = 2;
  else
    y = 1;
  end
else
  y = sample_discrete(T);
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [b,c,d,x,y,z,qfac] = vox2rasToQform(vox2ras)
% [a,b,c,x,y,z,qfac] = vox2rasToQform(vox2ras)
%
% Converts a vox2ras matrix to NIFTI qform parameters.
% Note: the vox2ras should be 6 DOF. This code mostly just
% follows CH's mriToNiftiQform() in mriio.c.
%
%  hdr.pixdim(1) = qfac;
%  hdr.quatern_b = b;
%  hdr.quatern_c = c;
%  hdr.quatern_d = d;
%  hdr.qoffset_x = x;
%  hdr.qoffset_y = y;
%  hdr.qoffset_z = z;
%  hdr.qform_code = NIFTI_XFORM_SCANNER_ANAT=1;
%

%
% vox2rasToQform.m
%
% Original Author: Doug Greve
%
% Copyright © 2011 The General Hospital Corporation (Boston, MA) "MGH"
%
% Terms and conditions for use, reproduction, distribution and contribution
% are found in the 'FreeSurfer Software License Agreement' contained
% in the file 'LICENSE' found in the FreeSurfer distribution, and here:
%
% https://surfer.nmr.mgh.harvard.edu/fswiki/FreeSurferSoftwareLicense
%
% Reporting: freesurfer@nmr.mgh.harvard.edu
%

a = [];
if(nargin ~= 1)
  fprintf('[a,b,c,x,y,z,qfac] = vox2rasToQform(vox2ras)\n');
  return;
end

x = vox2ras(1,4);
y = vox2ras(2,4);
z = vox2ras(3,4);

d = sqrt(sum(vox2ras(:,1:3).^2));
Mdc = vox2ras(1:3,1:3) ./ repmat(d,[3 1]);
if(det(Mdc) == 0.0)
  fprintf('ERROR: vox2ras determinant is 0\n');
  return;
end

r11 = Mdc(1,1);
r21 = Mdc(2,1);
r31 = Mdc(3,1);
r12 = Mdc(1,2);
r22 = Mdc(2,2);
r32 = Mdc(3,2);
r13 = Mdc(1,3);
r23 = Mdc(2,3);
r33 = Mdc(3,3);

if(det(Mdc) > 0.0) qfac = 1.0;
else
  r13 = -r13;
  r23 = -r23;
  r33 = -r33;
  qfac = -1.0;
end

%  /* following mat44_to_quatern() */

a = r11 + r22 + r33 + 1.0;
if(a > 0.5)
  a = 0.5 * sqrt(a);
  b = 0.25 * (r32-r23) / a;
  c = 0.25 * (r13-r31) / a;
  d = 0.25 * (r21-r12) / a;
else
  xd = 1.0 + r11 - (r22+r33);
  yd = 1.0 + r22 - (r11+r33);
  zd = 1.0 + r33 - (r11+r22);
  if(xd > 1.0)
    b = 0.5 * sqrt(xd);
    c = 0.25 * (r12+r21) / b;
    d = 0.25 * (r13+r31) / b;
    a = 0.25 * (r32-r23) / b;
  elseif( yd > 1.0 )
    c = 0.5 * sqrt(yd);
    b = 0.25 * (r12+r21) / c;
    d = 0.25 * (r23+r32) / c;
    a = 0.25 * (r13-r31) / c;
  else
    d = 0.5 * sqrt(zd);
    b = 0.25 * (r13+r31) / d;
    c = 0.25 * (r23+r32) / d;
    a = 0.25 * (r21-r12) / d;
  end
  if(a < 0.0)
    a = -a;
    b = -b;
    c = -c;
    d = -d;
  end
end

return

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


function [x, y] = trans_equi2fisheye(u, v, R, M, D, fe)

[xcam, ycam, zcam] = trans_equi2cam(u, v, fe);

xr = R(1, 1) * xcam + R(1, 2) * ycam + R(1, 3) * zcam ;
yr = R(2, 1) * xcam + R(2, 2) * ycam + R(2, 3) * zcam ;
zr = R(3, 1) * xcam + R(3, 2) * ycam + R(3, 3) * zcam ;

[x, y] = trans_cam2fisheye(xr, yr, zr, M, D);

end


text:
function b = isnumber(str)

if ~ischar(str)
    b = false;
    return;
end
if isempty(str)
    b = false;
    return;
end

b = all( (str>='0' & str<='9') | (str =='e') | (str=='-' | str=='+' | str=='.' | str == ' '));

author: BUNPC


text:
% This script generates a random deployed wireless network, and
% probabilistic broadcast is simulated.
% In probabilistic broadcast, each node re-transmits the packet according
% to a probability.
% A link on probabilistic broadcast:
% http://lsewww.epfl.ch/Documents/acrobat/SCS03.pdf
% A link on DFS algorithm:
% http://en.wikipedia.org/wiki/Depth-first_search

clear;
rand('state', 5);
numOfNodes = 200;   %   number of nodes
envSize=100;        %   envsizeXenvsize environment
txRange = 15;
global floodProb;   %   the pre-defined probability to broadcast;
global savedTransmission;   % The number of saved re-transmissions;

savedTransmission = 0;
global xLocation;   % Array containing the X-coordinations of wireless nodes;
global yLocation;   % Array containing the Y-coordinations of wireless nodes;
xLocation = rand(numOfNodes,1) * envSize;
yLocation = rand(numOfNodes,1) * envSize;   %x,y coords of nodes

% Generate the adjacent matrix to represent the topology graph of the
% randomly deployed wireless networks;
distMatrix = zeros(numOfNodes,numOfNodes);
for i=1:numOfNodes
   for j=1:numOfNodes
      distMatrix(i,j)=sqrt((xLocation(i)-xLocation(j))^2 + (yLocation(i)-yLocation(j))^2); %distance between node pairs
   end
end
% If the Euclidean distance between two nodes is less than the transmission range, there
% exists a link.
global connMatrix;
connMatrix = ( distMatrix < txRange);  %binary connectivity matrix

% The broadcast will start from a sink node;
sinkNode = 1; % sink node;

% Array visited[] stores the boolean value if the broadcast has reach the node.
global visited;
visited = zeros(1, numOfNodes);

% Show the topology;
figure(1);
plot(xLocation, yLocation, '.');
text(xLocation(sinkNode), yLocation(sinkNode), 'sink');
% title(['p = ' num2str(floodProb)]);
% sink node;

visited(sinkNode) = 1;

% recursive Depth first search is adopted to traverse the whole graph;
DFS(sinkNode);
savedTransmission
return;

author: Sable

repo: mcbench-benchmarks


text:
function aic = gmm_aic(Data,w,Priors,Mu,Sigma,cov_type)
%GMM_AIC Akaike Information criterion
%
%  input ------------------------------------------------------------------
%
%   o Data:     D x N array representing N datapoints of D dimensions.
%
%   o w:        1 x N, set of weights associated with each data point.
%
%   o Priors:   1 x K array representing the prior probabilities of the
%               K GMM components.
%   o Mu:       D x K array representing the centers of the K GMM components.
%
%   o Sigma:    D x D x K array representing the covariance matrices of the
%               K GMM components.
%
%   o cov_type: string, covariance type = 'full','diag' or 'isot'
%
%  output -----------------------------------------------------------------
%
%   o aic   : (1 x 1)
%
%
[D,~]       = size(Data);
K           = length(Priors);
num_param   = K-1 + K * D;

if strcmp(cov_type,'full') == true
    num_param = num_param + K * ( D * ( D - 1)/2 );
elseif strcmp(cov_type,'diag') == true
    num_param = num_param + K * D;
elseif strcmp(cov_type,'isot') == true
    num_param = num_param + K * 1;
else
   error(['no such covariance type: ' cov_type '  only full | diag | isot ']);
end

% compute the loglikelihood of the data given the model
% loglik: (N x 1)
loglik = LogLikelihood_gmm(Data,Priors,Mu,Sigma,w);

aic = - 2 * loglik + 2 * num_param;

end

author: epfl-lasa


text:
function PTKRunTests(tests_or_folder)
    % PTKRunTests. Runs PTK unit tests
    %
    %     tests_to_run - All tests will be run if this argument is omitted
    %
    %
    %
    % Executes all the tests that can be found and reports on the sucess and
    % failures.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2012.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    if (nargin == 0) || isempty(tests_or_folder)
        tests_or_folder = PTKDirectories.GetTestSourceDirectory;
    end

    CoreRunTests(tests_or_folder);
end

author: tomdoel

repo: pulmonarytoolkit


text:
% COINW  Write a Mathematical File from Matlab Values
%
% THIS IS A LOW LEVEL FUNCTION - USE coinWrite() INSTEAD!
%
% coinW uses the COIN Utilities library to open and write the file.
%
%   coinW(prob,path,filetype)
%
%   Input arguments:
%       prob - problem structure (described below)
%       path - full path to the file on your PC.
%       filetype - one of the following types: [mps,lp]
%
%   Problem Structure Fields (all required):
%       H - quadratic objective matrix (sparse, may be empty)
%       f - linear objective vector
%       A - linear inequality and equality constraint matrix (sparse)
%       rl - linear constraint lower bounds
%       ru - linear constraint upper bounds
%       lb - decision variable lower bounds
%       ub - decision variable upper bounds
%       int - integer variable vector (8 bit integer / char)
%       name - problem name (string)
%       sos_type - SOS types (char array)
%       sos_index - SOS indicies (cell array of double vectors)
%       sos_weight - SOS weights (cell array of double vectors)
%       objbias - Objective bias term (optional, negated internally)
%
%
%   Copyright (C) 2011 Jonathan Currie (I2C2)

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
%*******************************************************************
%Initialize the BBCI Toolbox:
%Hand over the directory for your (raw) data, which is in
%the Matlab variable 'MyDataDir' (e.g. 'D:\data' under Windows).

%If you have a tmp subfolder in the same directory,
%Start up the Toolbox as follows:
%"startup_bbci_toolbox('DataDir', MyDataDir);"

%If you want to choose a different folder for temporary files (e.g. 'D:\tmp\')
%Start up the Toolbox as follows:
%"startup_bbci_toolbox('DataDir', MyDataDir, 'TmpDir','/tmp/');"
%*******************************************************************

function startup_bbci_toolbox(varargin)

global BTB

% Find directory of the BBCI Toolbox and add it to the path
BBCI_DIR= fileparts(which(mfilename));
addpath(genpath(BBCI_DIR));
rmpath(genpath(fullfile(BBCI_DIR, '.git')));
% The subfolder 'external' is move to the end of the path, see below.

BBCI_PRIVATE_DIR= fullfile(fileparts(BBCI_DIR), 'bbci_private');

BTB= opt_proplistToStruct(varargin{:});

if ~isfield(BTB, 'TypeChecking'),
  BTB.TypeChecking= 1;
end
BTB= opt_setDefaults(BTB, {'DataDir'   '/home/bbci/data/'   'CHAR'});
if ~exist(BTB.DataDir, 'dir'),
  warning('Optional argument ''DataDir'' should specify an existing folder.');
end

% Guess what the location of other directories could be
BBCI_RAW_DIR= fullfile(BTB.DataDir, 'bbciRaw/');
if ~exist(BBCI_RAW_DIR, 'dir'),
  BBCI_RAW_DIR= BTB.DataDir;
end
BBCI_MAT_DIR= fullfile(BTB.DataDir, 'bbciMat/');
if ~exist(BBCI_MAT_DIR, 'dir'),
  BBCI_MAT_DIR= BTB.DataDir;
end
PYFF_DIR= fullfile(fileparts(BBCI_DIR), 'pyff', 'src');
if ~exist(PYFF_DIR, 'dir'),
  PYFF_DIR= '';
end

props= {'Dir'            BBCI_DIR          'CHAR';
        'DataDir'        ''                'CHAR';
        'RawDir'         BBCI_RAW_DIR      'CHAR';
        'MatDir'         BBCI_MAT_DIR      'CHAR';
        'PrivateDir'     BBCI_PRIVATE_DIR  'CHAR';
        'TmpDir'         ''                'CHAR';
        'PyffDir'        PYFF_DIR          'CHAR';
        'FigDir'         ''                'CHAR';
        'Tp'             struct            'STRUCT';
        'Acq'            struct            'STRUCT';
        'History'        1                 '!BOOL';
        'TypeChecking'   1                 '!BOOL'
       };
[BTB, isdefault]= opt_setDefaults(BTB, props);

% import dependencies if not done yet
if ~exist(fullfile(BTB.Dir, 'external', 'auto_import_performed'), 'file')
    inp = input(sprintf(['\nThis is the first startup of the BBCI toolbox - Welcome!\n\n', ...
        'You can now automatically download and import code from other sources.\n', ...
        'This adds additional functionalities such as ICA and SSD to the BBCI toolbox.\n',...
        'Following startups won''t ask for automatic imports. Manual imports are\n',...
        'possible with ""bbci_import_dependencies()"" at any time.\n\n',...
        'Type <yes> to download all dependencies (approx 10MB): ']), 's');
    if strcmp(inp, 'yes') || strcmp(inp, 'y')
        bbci_import_dependencies();
    end
     % write into file
    fileID = fopen(fullfile(BTB.Dir, 'external', 'auto_import_performed'),'w');
    fprintf(fileID, datestr(now()));
    fclose(fileID);
end

% Put external folders to the end of the path, such that they cannot
% overload functions of the BBCI Toolbox:
IMPORT_PATH= genpath(fullfile(BBCI_DIR, 'external'));
rmpath(IMPORT_PATH);
path(path, IMPORT_PATH);

if exist(BTB.PrivateDir, 'dir'),
  private_folders_to_add= {'utils', 'startup'};
  for kk= 1:length(private_folders_to_add),
    folder= fullfile(BTB.PrivateDir, private_folders_to_add{kk});
    if exist(folder, 'dir'),
      addpath(genpath(folder));
    end
  end
end

if isdefault.TmpDir
  BTB.TmpDir= fullfile(BTB.DataDir, 'tmp');
  if ~exist(BTB.TmpDir, 'dir'),
    fprintf('!! Default TEMP dir not existing at\n  %s\n', BTB.TmpDir);
    fprintf('!! Setting TEMP dir to ''''\n');
    BTB.TmpDir= '';
  end
end

% Information about the test person (Tp)
props= {'Dir'    BTB.TmpDir  'CHAR'
        'Code'   ''          'CHAR'
       };
BTB.Tp= opt_setDefaults(BTB.Tp, props);

% Information about data acquistion
props= {'Prefix'          'a'                     'CHAR'
        'StartLetter'     'a'                     'CHAR'
        'Dir'             BTB.TmpDir              'CHAR'
        'Geometry'        []                      'DOUBLE[1 4]'
        'TriggerFcn'      @bbci_trigger_parport   'FCN'
        'TriggerParam'    {}                      'CELL'
       };
BTB.Acq= opt_setDefaults(BTB.Acq, props);

switch computer
 case 'PCWIN'
  BTB.Acq.IoLib= which('inpout32.dll');
 case 'PCWIN64'
  BTB.Acq.IoLib= which('inpoutx64.dll');
end

if isfield(BTB.Acq, 'IoLib') && isfield(BTB.Acq, 'IoAddr'),
  if isempty(BTB.Acq.TriggerParam) && ...
      isequal(BTB.Acq.TriggerFcn, @bbci_trigger_parport),
    BTB.Acq.TriggerParam= {BTB.Acq.IoLib, BTB.Acq.IoAddr};
  end
else
  fprintf('Parport not installed. Triggers will just be printed.\n');
  BTB.Acq.TriggerFcn= @bbci_trigger_print;
end

evalin('base', 'global BTB');

author: bbci


text:
function [imdepthMin, imdepthMax, imdepthCol, contact, x, y, z] = ...
    getDepthRangeForDisplay(bndinfo, glabels, elabels, dt, v0)
% getMinimumDepth(bndinfo, glabels, v0)
%
% Minimum depth is the depth of the foremost occluder.  Maximum depth is
% the depth of the current region, assuming that it is touching ground at
% its lowest point.  If a region is in the foreground, it's minimum depth
% may be a function, and the minimum depth is equal to the maximum depth.
% Otherwise, the minimum and maximum depths are single (different) values.
%
% log(Z) = log(f*y_c) - log(v_0-v_i)

global DO_DISPLAY;

%% Get depth ordering

ordering = occlusion2depthOrdering(bndinfo, elabels, glabels);

%% Set some parameters

MIN_HV_DIST = 0.01;  % minimum distance between horizon and object contact
SKY_DEPTH = 2*(1 ./ MIN_HV_DIST);

%% Get pixel images for ground, vertical, and sky
labim = glabels(bndinfo.wseg);
vim = labim==2;
gim = labim==1;
sim = labim==3;

[imh, imw] = size(vim);

%% Make sure that horizon is above ground and below sky
scol = sum(sim, 2)>0;
gcol = sum(gim, 2)>0;
v0 = (1-v0) * imh;

try
minv = find(scol);  minv = minv(end-1);
maxv = find(gcol);  maxv = maxv(2);

if v0 < minv || v0 > maxv
    if minv > maxv
        v0 = maxv-1;
    else
        v0 = maxv-imh/20;
    end
end
%disp(num2str([v0 maxv minv]))
catch
end

v0 = 1 - v0/imh;
%disp(num2str(v0))

%figure(1), imagesc(cat(3, vim, gim, sim)), axis image

%% Get possible ground-vertical boundary pixels

stats = regionprops(bndinfo.wseg, 'BoundingBox', 'PixelIdxList');
bbox = vertcat(stats.BoundingBox);
idx = {stats.PixelIdxList};

% bbox = [x1 y1 x2 y2]
bbox = [bbox(:, 1) bbox(:, 2) bbox(:, 1)+bbox(:,3) bbox(:,2)+bbox(:,4)];

vfilt = [ones(3, 1) ; zeros(3, 1)];
gfilt = 1-vfilt;
boundaryim = (imfilter(double(vim), vfilt)==sum(vfilt(:))) & ...
    (imfilter(double(gim), gfilt)==sum(gfilt(:)));
boundaryim(end, :) = vim(end, :);
boundaryim = imdilate(boundaryim, ones(3, 1));

%figure(3), imagesc(boundaryim), axis image

bpts = cell(bndinfo.nseg, 1);
edges = bndinfo.edges.indices;
spLR = bndinfo.edges.spLR;
eim = zeros(bndinfo.imsize);
for k = 1:numel(edges)
    sp1 = spLR(k, 1);
    sp2 = spLR(k, 2);
    bndind = edges{k}(boundaryim(edges{k}));
    if ~isempty(bndind)
        if glabels(sp1)==2 % && glabels(sp2)==1
            bpts{sp1} = [bpts{sp1} ; bndind];
        elseif glabels(sp2)==2 % && glabels(sp1)==1
            bpts{sp2} = [bpts{sp2} ; bndind];
        end
    end
    eim(edges{k}) = 1;
end

%% Get a simple estimate of depth (trace columns upwards)

imdepthCol = zeros(bndinfo.imsize(1:2));
maxg = imh*ones(1, imw);
for y = imh:-1:1
    maxg(gim(y, :)) = y;
    imdepthCol(y, :) = maxg;
end
%imy = repmat((1:imh)', [1 imw]);
% imdepthCol = log(1 ./ max(v0 - (imh-imdepthCol)/imh, MIN_HV_DIST));
% imdepthCol(sim==1) = log(SKY_DEPTH);
for k = 1:numel(idx)
    imdepthCol(k) = (imh-median(imdepthCol(idx{k})))/imh;
end
% if DO_DISPLAY
%     figure(3), imagesc(imdepthCol), axis image
% end

%% Get the depth of each object that has a ground boundary

bx = zeros(numel(bpts), 2);
by = zeros(numel(bpts), 2);
bz = zeros(numel(bpts), 2);
bzmin = zeros(numel(bpts), 2);
bzmax = zeros(numel(bpts), 2);
isvalid = false(numel(bpts), 1);
for k = 1:numel(bpts)
    if numel(bpts{k})>5 || bbox(k, 4)+1>=imh

        [py, px] = getPerimeter(bndinfo, k, 50);
        bx(k, :) = bbox(k, [1 3]);

        cu = [];  cv = [];
        if  1 %bx(k, 2)-bx(k, 1) > imw/10  % only get line for larger objects

            keep = true(size(px));
            for tk = 1:numel(px)
                if ~any(boundaryim(py(tk) + imh*(px(tk)-1)))
                    keep(tk) = false;
                end
            end
            %px = px(keep);  py = py(keep);

            %if k==13, keyboard; end

            try
                [cu, cv] = poly2contacts(dt, px, py, bndinfo.imsize, v0, 1.7, [], 0.1);
            catch
            end
            %plot(cu, cv, 'r*')
            cu = round(cu);  cv = round(cv);
            ind = (boundaryim(cv + (cu-1)*imh)); % | (cv==max(cv));
            cu = cu(ind);  cv = cv(ind);
        end

        [cu, ind] = unique(cu);
        cv = cv(ind);

        if ~isempty(cu)
            if numel(cu) == 1
                ypos = max(py);
                imslope = 0;
            elseif numel(cu)==2
                imslope = (cv(2)-cv(1)) / (cu(2)-cu(1));
                ypos = cv(1) - imslope*cu(1);
            else
                rline = robustfit(cu, cv);
                ypos = rline(1);
                imslope = rline(2);
            end
            by(k, :) = imslope*bx(k, :) + ypos;
            bz(k, :) = (1 ./ max(v0 - (imh-by(k, :))/imh, MIN_HV_DIST));
            isvalid(k) = true;
        end

        contact{k} = [cu(:) cv(:)];

    end
end

bzmin(isvalid, :) = bz(isvalid, :);
bzmax(isvalid, :) = bz(isvalid, :);
unknown = (~isvalid) & (glabels==2);

%% Assign a maximum depth to object of unknown depth and sky
missing = unknown;

bzmax(missing, :) = Inf;

for k = find(missing)'
    bx(k, :) = bbox(k, [1 3]);
    by(k, :) = bbox(k, 4);

    tmpmax = 1 ./ max(v0 - (imh-by(k, :))/imh, MIN_HV_DIST);

    if any(tmpmax < bzmax(k, :))
        bzmax(k, :) =tmpmax;
    end
end

% each region must be in back of things that occlude it (as far back as
% possible without being in front of things that it occludes)
missing2 = missing;
changed = 1;
while changed
    lastbz = bzmax;
    %bzmax(missing2, :) = Inf;
    for k = find(missing2)'
        bzmax(k, :) =  1 ./ max(v0 - (imh-by(k, :))/imh, MIN_HV_DIST);
    end
    %changed = 0;
    for k = find(missing)'
        pos1 = find(ordering==k);

        for k2 = 1:numel(ordering)

            pos2 = find(ordering==k2);
            if pos1 < pos2  && glabels(k2)==2 && ~missing2(k2) ...
                    && (any(spLR(spLR(:, 1)==k, 2)==k2) || any(spLR(spLR(:, 2)==k, 1)==k2))
                % k is directly in front of k2
                imslope = (by(k2, 2)-by(k2, 1))./(bx(k2, 2)-bx(k2, 1));
                ypos = by(k2, 1) - imslope*bx(k2, 1);
                by(k, :) = imslope*bx(k, :) + ypos;
                tmpz = 1 ./ max(v0 - (imh-by(k, :))/imh, MIN_HV_DIST);

                if all(tmpz <= bzmax(k, :))
                    bzmax(k, :) = tmpz; %(1 ./ max(v0 - (imh-by(k, :))/imh, 0.01));
                    changed = 1;
                    missing2(k) = false;
                end

            end
        end
    end
    if all(lastbz==bzmax)
        break;
    end
end

bzmax(isinf(bzmax(:, 1)), :) = 1./ MIN_HV_DIST;

%% Assign a minimum depth to objects of unknown depth;

missing = unknown;

% assign depth as maximum depth of object that occludes it
changed = 1;
while changed
    changed = 0;
    isvalid2 = isvalid;
    for k = find(missing)'
        bx(k, :) = bbox(k, [1 3]);

        bzmin(k, :) = 0; %bzmax(k, :);

        ind = find(spLR(:, 1)==k);
        neighbors = spLR(ind, 2);
        validn = isvalid(neighbors) & (glabels(neighbors)==2);
        neighbors = neighbors(validn);
        ind = ind(validn);

        pos1 = find(ordering==k);

        maxov = 0;
        for tk2 = 1:numel(neighbors)
            k2 = neighbors(tk2);
            pos2 = find(ordering==k2);

            ov = min(bbox([k k2], 3)) - max(bbox([k k2], 1));

            if ~missing(k2) && pos2 < pos1  && glabels(k2)==2 & ...
                    (any(spLR(spLR(:, 1)==k, 2)==k2) || any(spLR(spLR(:, 2)==k, 1)==k2))

                 %(ov > maxov) || (ov > 0.25*bbox(k, 3)-bbox(k,1))
                imslope = (by(k2, 2)-by(k2, 1))./(bx(k2, 2)-bx(k2, 1));
                ypos = by(k2, 1) - imslope*bx(k2, 1);
                by(k, :) = imslope*bx(k, :) + ypos;
                tmpz = 1 ./ max(v0 - (imh-by(k, :))/imh, MIN_HV_DIST);

                if any(tmpz > bzmin(k, :))
                    bzmin(k, :) = tmpz; %(1 ./ max(v0 - (imh-by(k, :))/imh, 0.01));
                    missing(k) = false;
                    isvalid2(k) = true;
                    changed = 1;
                end

            end

        end
    end

    isvalid = isvalid2;

end
bzmin(bzmin(:, 1)==0, :) = bzmax(bzmin(:, 1)==0, :);

% check that is not behind anything that it occludes
changed = 1;
while changed
    changed = 0;
    isvalid2 = isvalid;
    for k = find(missing)'
        bx(k, :) = bbox(k, [1 3]);

        ind = find(spLR(:, 1)==k);
        neighbors = spLR(ind, 2);
        validn = isvalid(neighbors) & (glabels(neighbors)==2);
        neighbors = neighbors(validn);
        ind = ind(validn);

        pos1 = find(ordering==k);

        maxov = 0;
        for tk2 = 1:numel(neighbors)
            k2 = neighbors(tk2);
            pos2 = find(ordering==k2);

            if ~missing(k2) && pos1 < pos2  &&  glabels(k2)==2 && ...
                    (any(spLR(spLR(:, 1)==k, 2)==k2) || any(spLR(spLR(:, 2)==k, 1)==k2))

                imslope = (by(k2, 2)-by(k2, 1))./(bx(k2, 2)-bx(k2, 1));
                ypos = by(k2, 1) - imslope*bx(k2, 1);
                by(k, :) = imslope*bx(k, :) + ypos;
                tmpz = 1 ./ max(v0 - (imh-by(k, :))/imh, MIN_HV_DIST);

                if all(tmpz < bzmin(k, :))
                    bzmin(k, :) = tmpz;
                    missing(k) = false;
                    isvalid2(k) = true;
                    changed = 1;
                end

            end

        end
    end

    isvalid = isvalid2;

end

depthim = zeros([imh imw]);
tmpy = zeros(size(depthim));
tmpz = zeros(size(depthim));
for k = 1:bndinfo.nseg
    [y, x] = ind2sub(bndinfo.imsize, idx{k});
    switch glabels(k)
        case 0
            z = bzmin(k, 1) + (x-bx(k, 1)) ./ ...
                (bx(k, 2)-bx(k, 1)) .* (bzmin(k,2)-bzmin(k,1));
%             tmpy(idx{k}) = (y - (x-bx(k, 1)) ./ ...
%                 (bx(k, 2)-bx(k, 1)) .* (by(k,2)-by(k,1))) * 1.7 ./ ((x-bx(k, 1)) ./ ...
%                 (bx(k, 2)-bx(k, 1)) .* (bzmin(k,2)-bzmin(k,1)));
            tmpy(idx{k}) = -(v0 - (imh-(x-bx(k, 1)) ./ ...
                (bx(k, 2)-bx(k, 1)) .* (by(k,2)-by(k,1)))/imh).*z / 1.4 + 1.7;
            z2 = z;
        case 1
            z = 1 ./ max(v0 - (imh-y)/imh, MIN_HV_DIST);
            tmpy(idx{k}) = 0;
            z2 = z*1.4*1.7;
        case 2
            z = bzmin(k, 1) + (x-bx(k, 1)) ./ ...
                (bx(k, 2)-bx(k, 1)) .* (bzmin(k,2)-bzmin(k,1));
            ty = by(k,1) + (x-bx(k, 1)) ./ (bx(k, 2)-bx(k, 1)) .* (by(k,2)-by(k,1));

%             tmpy(idx{k}) = -(v0 - (imh-(x-bx(k, 1)) ./ ...
%                 (bx(k, 2)-bx(k, 1)) .* (by(k,2)-by(k,1)))/imh).*z / 1.4 + 1.7;

            z2 = 1.4*1.7 ./ max(v0 - (imh-ty)/imh, MIN_HV_DIST);
            %(1.7-y(idx{k})) .* 1.4 ./ (v0 - (imh-ty)./imh+1E-5);
        case 3
            z = SKY_DEPTH;
            tmpy(idx{k}) = 500;
            z2 = (1.7-tmpy(idx{k})) .* 1.4 ./ (v0 - (imh-y)./imh+1E-5);
    end
    depthim(idx{k}) = log(z);
    tmpz(idx{k}) = z2;
end

imdepthMin = depthim;

% x = u.*z./f;

if DO_DISPLAY
    figure(1), imagesc(depthim), axis image, colormap jet
end
%print -f1 -djpeg99 ../tmp/gt_mindepth.jpg

%depthim = zeros([imh imw]);
for k = find(glabels==2)'
    [y, x] = ind2sub(bndinfo.imsize, idx{k});
    zk = bzmax(k, 1) + (x-bx(k, 1)) ./ ...
        (bx(k, 2)-bx(k, 1)) .* (bzmax(k,2)-bzmax(k,1));
    depthim(idx{k}) = log(zk);
end

if DO_DISPLAY
    figure(2), imagesc(depthim), axis image, colormap jet
end
%print -f2 -djpeg99 ../tmp/gt_maxdepth.jpg

imdepthMax = depthim;

x = repmat((1:imw)/imw-0.5, [imh 1]).*tmpz/1.4;
y = 1.7 - (v0-repmat(1-(1:imh)'./imh, [1 imw]))/1.4.*tmpz;
z = tmpz;

%% gets perimeter of region
function [py, px] = getPerimeter(bndinfo, r, npts)

imsize = bndinfo.imsize;
spLR = bndinfo.edges.spLR;
spLR = [spLR ; spLR(:, [2 1])];
eind = find(spLR(:, 1)==r);

failed = false;
order = zeros(1, numel(eind));
order(1) = 1;
for k = 2:numel(order)
    laste = eind(order(k-1));
    ismem = ismember(eind, bndinfo.edges.adjacency{laste});
    next = find(ismem);
    if numel(next)~=1
        failed = true;
        break;
    end
    order(k) = next;
end

if ~failed
    eind = eind(order);
end

edges = cell(numel(eind), 1);
for k = 1:numel(eind)
    if eind(k) <= bndinfo.ne
        edges{k} = bndinfo.edges.indices{eind(k)}(2:end-1);
    else
        edges{k} = bndinfo.edges.indices{eind(k)-bndinfo.ne}(end-1:-1:2);
    end
end
[py, px] = ind2sub(imsize, vertcat(edges{:}));

try
    ind = convhull(double(px), double(py));
catch
    ind = (1:5:numel(px))';
end

if npts > numel(ind) && ~failed
    step = ceil(numel(py) / (npts-numel(ind)));
    ind2 = (1:step:numel(py))';
    ind = unique([ind ; ind2]);
end

py = double(py(ind));
px = double(px(ind));

author: Cloud-CV

repo: object-proposals


text:
function evaluate_tDCF_asvspoof19(CM_SCOREFILE, ASV_SCOREFILE)

% Set t-DCF parameters
cost_model.Pspoof       = 0.05;  % Prior probability of a spoofing attack
cost_model.Ptar         = (1 - cost_model.Pspoof) * 0.99; % Prior probability of target speaker
cost_model.Pnon         = (1 - cost_model.Pspoof) * 0.01; % Prior probability of nontarget speaker
cost_model.Cmiss_asv    = 1;     % Cost of ASV system falsely rejecting target speaker
cost_model.Cfa_asv      = 10;    % Cost of ASV system falsely accepting nontarget speaker
cost_model.Cmiss_cm     = 1;     % Cost of CM system falsely rejecting target speaker
cost_model.Cfa_cm       = 10;    % Cost of CM system falsely accepting spoof

% Load organizer's ASV scores
[asv_attacks, asv_key, asv_score] = textread(ASV_SCOREFILE, '%s %s %f');

% Load CM scores (replace these with the scores of your detectors).
[cm_utt_id, cm_attacks, cm_key, cm_score] = textread(CM_SCOREFILE, '%s %s %s %f');

% Extract target, nontarget and spoof scores from the ASV scores
tar_asv     = asv_score(strcmp(asv_key, 'target'))';
non_asv     = asv_score(strcmp(asv_key, 'nontarget'))';
spoof_asv   = asv_score(strcmp(asv_key, 'spoof'))';

% Extract bona fide (real human) and spoof scores from the CM scores
bona_cm     = cm_score(strcmp(cm_key, 'bonafide'));
spoof_cm    = cm_score(strcmp(cm_key, 'spoof'));

% Fix ASV operating point to EER threshold
[eer_asv, asv_threshold] = compute_eer(tar_asv, non_asv);

% Visualize ASV and CM scores
figure(1);
subplot(121);
[h, xi] = ksdensity(tar_asv);   plot(xi, h, 'gr'); hold on;
[h, xi] = ksdensity(non_asv);   plot(xi, h, 'k');
[h, xi] = ksdensity(spoof_asv); plot(xi, h, 're');
h = plot(asv_threshold, 0, 'mo'); set(h, 'MarkerSize', 10);
legend('Target', 'Nontarget', 'Spoof', 'EER threshold');
xlabel('ASV score'); ylabel('PDF');
title('ASV scores', 'FontSize', 12, 'FontWeight', 'bold');

subplot(122);
[h, xi] = ksdensity(bona_cm);   plot(xi, h, 'gr'); hold on;
[h, xi] = ksdensity(spoof_cm);  plot(xi, h, 're');
legend('Bona fide', 'Spoof');
xlabel('CM score'); ylabel('PDF');
title('CM scores', 'FontSize', 12, 'FontWeight', 'bold');

% Obtain the detection error rates of the ASV system
[Pfa_asv, Pmiss_asv, Pmiss_spoof_asv] = obtain_asv_error_rates(tar_asv, non_asv, spoof_asv, asv_threshold);

% Equal error rate of the countermeasure
[eer_cm, ~] = compute_eer(bona_cm, spoof_cm);

% Compute t-DCF
[tDCF_curve, CM_thresholds] = compute_tDCF(bona_cm, spoof_cm, Pfa_asv, Pmiss_asv, Pmiss_spoof_asv, cost_model, true);

% Minimum normalized t-DCF and the corresponding threshold
[min_tDCF, argmin_tDCF] = min(tDCF_curve);
min_tDCF_threshold = CM_thresholds(argmin_tDCF);

fprintf('ASV SYSTEM\n');
fprintf('     EER               \t= %5.5f %%\t\t Equal error rate (target vs. nontarget discrimination)\n',   100 * eer_asv);
fprintf('     Pfa               \t= %5.5f %%\t\t False acceptance rate of nontargets\n',   100 * Pfa_asv);
fprintf('     Pmiss             \t= %5.5f %%\t\t Miss (false rejection) rate of targets\n',   100 * Pmiss_asv);
fprintf('     1-Pmiss,spoof     \t= %5.5f %%\t Spoof false acceptance rate ("NOT miss spoof trial")\n\n', 100 * (1 - Pmiss_spoof_asv));

fprintf('CM SYSTEM\n');
fprintf('     EER               \t= %f %%\t Equal error rate from CM scores pooled across all attacks. \n\n', 100 * eer_cm);

fprintf('TANDEM\n');
fprintf('     min-tDCF          \t= %f\n',   min_tDCF);

% Plot t-DCF as function of the CM threshold.
figure(2);
plot(CM_thresholds, tDCF_curve); hold on;
h = plot(min_tDCF_threshold, min_tDCF, 'mo'); set(h, 'MarkerSize', 10);
xlabel('CM threshold (operating point)');
ylabel('Norm t-DCF');
title('Normalized t-DCF', 'FontSize', 12, 'FontWeight', 'bold');
h = line([min(CM_thresholds) max(CM_thresholds)], [1 1]); set(h, 'LineStyle', '--', 'Color', 'k');
legend('t-DCF', sprintf('min t-DCF (%5.5f)', min_tDCF), 'Arbitrarily bad CM (Norm t-DCF=1)');
axis([min(CM_thresholds) max(CM_thresholds) 0 1.5]);

% Traditional DET plot using Bosaris toolkit
plot_type = Det_Plot.make_plot_window_from_string('old');
plot_obj = Det_Plot(plot_type, 'Detection error trade-off plot of CM');
plot_obj.set_system(bona_cm, spoof_cm, sprintf('CM system (EER=%2.2f %%)', 100* eer_cm));
plot_obj.plot_steppy_det({'b','LineWidth',2},' ');
plot_obj.display_legend();

end

author: nesl


text:
% Reproduces an experiment similar to the one from figure 2 in "A
% sliding-window kernel RLS algorithm and its application to nonlinear
% channel identification".
%
% MSE performance of sliding-window kernel RLS for two different window
% sizes. Execution time: < 1 minute (Intel Pentium Core2 Duo).
%
% S. Van Vaerenbergh, J. Via, and I. Santamaria, "A Sliding-Window Kernel
% RLS Algorithm and Its Application to Nonlinear Channel Identification,"
% 2006 IEEE International Conference on Acoustics, Speech and Signal
% Processing (ICASSP), May 2006,
%
% This file is part of the Kernel Adaptive Filtering Toolbox for Matlab
% https://github.com/steven2358/kafbox/

clear
close all

%% PARAMETERS
opt.chlen = 4; % channel length, and input embedding
opt.N = 1500; % number of training data
opt.N_switch = 500; % iteration after which the channel switch occurs
opt.N_test = 100; % number of test data
opt.SNR = 20; % SNR in dB

num_sim = 10;

setups{1} = struct('M',75,'c',1E-2,'kerneltype','gauss','kernelpar',5);
setups{2} = struct('M',150,'c',1E-2,'kerneltype','gauss','kernelpar',5);

titles = {'SW-KRLS, M=75','SW-KRLS, M=150'};

%% PREPARE DATA
fprintf('Fig. 2 from "A sliding-window kernel RLS algorithm and its \n');
fprintf('application to nonlinear channel identification".\n');

% perform online learning for each algorithm
fprintf('\n')
N = opt.N;
N_switch = opt.N_switch;
num_alg = length(setups);
MSE = zeros(N,num_alg);
MSE_final = zeros(1,num_alg);
for sim_ind = 1:num_sim
    % generate data
    [X,y,y_ref,X_test,y_test] = generate_channel_switch(opt);

    fprintf('SIM %d:\n',sim_ind)
    for algo_ind = 1:num_alg
        t1 = tic;

        kaf = swkrls(setups{algo_ind});

        fprintf('%9s M=%3d: ',upper(class(kaf)),kaf.M);

        mse = zeros(N,1);
        for i=1:N
            if ~mod(i,floor(N/10)), fprintf('.'); end

            y_est = kaf.evaluate(X_test);
            if i<=N_switch
                mse(i) = mean((y_test(:,1)-y_est).^2);
            else
                mse(i) = mean((y_test(:,2)-y_est).^2);
            end

            kaf.train(X(i,:),y(i));
        end
        MSE(:,algo_ind) = MSE(:,algo_ind) + mse/num_sim;
        mse_final = mean(mse(N-500:N));

        fprintf(' %.2fs. Final MSE=%3.2fdB\n',toc(t1),...
            10*log10(mse_final))
    end
    fprintf('\n');
end
fprintf('\n');

%% OUTPUT

figure; hold all
for algo_ind=1:num_alg
    plot(10*log10(MSE(:,algo_ind)),'LineWidth',1)

    axis([0 N 5*floor(min(10*log10(MSE(:)))/5) 0]);
end
xlabel('iteration')
ylabel('MSE (dB)')
legend(titles)

repo: kafbox


text:
function Xr = pcaRandVec( U, mu, vars, k, n, hypershpere, show )
% Generate random vectors in PCA subspace.
%
% Used to generate random vectors from the subspace spanned by the first k
% principal components.  The points generated come from the gaussian
% distribution from within the subspace. Can optionally generate points on
% the subspace that are also on a hypershpere centered on the origin.  This
% may be useful if the original data points were all from a hypershpere --
% for example they were normalized via imNormalize.  Set the optional
% hypershpere flag to 1 to generate points only on the hypersphere.
%
% USAGE
%  Xr = pcaRandVec( U, mu, vars, k, n, [hypershpere], [show] )
%
% INPUTS
%  U           - returned by pca.m
%  mu          - returned by pca.m
%  vars        - returned by pca.m
%  k           - number of principal coordinates to use
%  n           - number of points to generate
%  hypershpere - [0] generate points on hypersphere (see above)
%  show        - [1] figure to use for display (no display if == 0)
%
% OUTPUTS
%  Xr          - resulting randomly generated vectors
%
% EXAMPLE
%
% See also PCA IMNORMALIZE
%
% Piotr's Computer Vision Matlab Toolbox      Version 2.0
% Copyright 2014 Piotr Dollar.  [pdollar-at-gmail.com]
% Licensed under the Simplified BSD License [see external/bsd.txt]

if( nargin<6 || isempty(hypershpere) ); hypershpere=0; end
if( nargin<7 || isempty(show) ); show=0; end

siz1 = size(mu);   nd = ndims(mu);
sizX = [siz1, n];  D = prod(siz1);

% generate random vectors inside of subspace.
C = diag( vars(1:k).^-1 );
Yr = C * randn(k,n);
Uk = U(:,1:k);
Xr = Uk * Yr;

if(hypershpere)
  % The final point must lie on hypershpere.  Need to scale each element xr
  % of Xr so that after adding it to mu the resulting vector has length D.
  % So need to find k>0 such that (k*xr + mu) has length D.  To find k
  % simply solve the quadratic equation induced by ||(k*xr + mu)||=D,
  % choosing the root such that k>0.  Note that the mean of resulting
  % vector will not necessarily be 0, but the variance will pretty much be
  % 1! Regardless, renomralize at the end.
  muv = mu(:);  muMag = dot(muv,muv);
  for i=1:n
    xr = Xr(:,i);
    rs = roots( [dot(xr,xr), 2 * dot(xr,muv), muMag-D] );
    Xr(:,i) = muv + max(rs)*xr;
  end
  Xr = reshape( Xr, sizX );
  Xr = fevalArrays( Xr, @imNormalize );
else
  % simply add the mean to reshaped Xr
  Xr = reshape( Xr, sizX );
  muRep = repmat(mu, [ones(1,nd), n ] );
  Xr = Xr + muRep;
end

% optionaly show resulting vectors
if(show && (nd==2 || nd==3)); figure(show); montage2(Xr); end

%%% Little test - see if eigenvectors induced by randomly generated vectors
%%% are the same as the original eigenvectors.  Assumes [U,mu,vars] exist.
%   Xr = pcaRandVec( U, mu, vars, 3, 100 );
%   [ Ur, mur, varsr ] = pca( Xr );
%   ind = 3;
%   Uim = reshape( U(:,ind), [ size(mu,1), size(mu,2) ]  );
%   Uimr = reshape( Ur(:,ind), [ size(mu,1), size(mu,2) ]  );
%   if( sum(abs(Uim-Uimr))>sum(abs(Uim+Uimr))) Uimr=Uimr*-1; end; %sign?
%   clf; subplot(3,1,1); im( Uim ); subplot(3,1,2); im( Uimr );
%   subplot(3,1,3); im( Uim - Uimr);
%%%

author: pdollar

repo: toolbox

save_path: github-repos/MATLAB/pdollar-toolbox


text:
function varargout = power(varargin)
%.^	     Pointwise power of a SPHEREFUN.
%
% F.^G returns a SPHEREFUN F to the scalar power G, a scalar F to the
% SPHEREFUN power G, or a SPHEREFUN F to the SPHEREFUN power G.

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

[varargout{1:nargout}] = power@separableApprox(varargin{:});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
%% Creates D-matrix to calculate trends endogenously

function D = TVEregimesDummy(startdate,regimeperiods,T)
% TVEregimesDummy(startdate,regimeperiods,X) creates is a matrix in which
% each row is a selector vector of the active regime

gdates=eval(sprintf(startdate(1:4)))+eval(sprintf(startdate(6)))*0.25:0.25:eval(sprintf(startdate(1:4)))+eval(sprintf(startdate(6)))*0.25-0.25+T/4;
NumberofRegimes=size(regimeperiods,1)+1;

% Create the dummy variable
D=zeros(T,NumberofRegimes);
for it=1:NumberofRegimes-1
    NumberofPeriods=size( regimeperiods(it,:),2)/2;
    for it2=1:NumberofPeriods
        startperiod=eval(sprintf(regimeperiods{it,2*it2-1}(1:4)))+eval(sprintf(regimeperiods{it,2*it2-1}(6)))*0.25;
        endperiod=eval(sprintf(regimeperiods{it,2*it2}(1:4)))+eval(sprintf(regimeperiods{it,2*it2}(6)))*0.25;
        startD=find(gdates==startperiod);
        endD=find(gdates==endperiod);
        D(startD:endD,it+1)=1;
    end
end
if NumberofRegimes>1
    D(:,1)=1-sum(D(:,2:end),2);
else
    D(:,1)=ones(size(D,1),1);
end

end

author: european-central-bank

repo: BEAR-toolbox


text:
function integralHist = integralHistSuperpixels(N)
%comment to be written

N = int16(N);
total_segms = max(max(N));
[height width] = size(N);

integralHist = zeros(height+1,width+1,total_segms);

for sid = 1:total_segms

    superpixelMap = not(N - sid);
    integralHist(:,:,sid) = computeIntegralImage(superpixelMap);
end

end


text:
function [acc_control_learning] = setupAccControlLearning(dataset)
%
% Authors:       Alexander Wischnewski
%
% Description:
%   function used to prepare the variables for the acceleration controller learning algorithms
%
% Input:
%   dataset: Log file with weights used for initialization (if given)

% spread of speed values
acc_control_learning.spread_vx_mps = 10:15:100;
acc_control_learning.spread_ay_req_mps2 = -36:6:36;

% design of basis functions
acc_control_learning.vx_width_mps = 20;
acc_control_learning.ay_width_mps2 = 8;

[VX, AY] = meshgrid(acc_control_learning.spread_vx_mps, acc_control_learning.spread_ay_req_mps2);
acc_control_learning.bf_vx_mps = reshape(VX, numel(VX), 1);
acc_control_learning.bf_ay_req_mps2 = reshape(AY, numel(AY), 1);

% initial weights for lateral acceleration controller learning
% will be updated below if dataset is given
acc_control_learning.w0 = zeros(length(acc_control_learning.bf_vx_mps)+1, 1);

% analyze dataset if available
if(nargin==1)
    % extract weights from log and write to initial weights
    data = load(dataset);
    weights = data.debug_slow.debug_slow_LatAcc_FFweights.Data(end, :);
    acc_control_learning.w0 = weights';
    % calculate learned feedforward part for the end of the log
    visz_vx_mps = 0:5:100;
    visz_ay_req_mps2 = -40:0.5:40;
    [Visz_X, Visz_Y] = meshgrid(visz_vx_mps, visz_ay_req_mps2);
    Visz_Z = zeros(length(visz_ay_req_mps2), length(visz_vx_mps));
    for i = 1:1:length(visz_vx_mps)
        for j = 1:1:length(visz_ay_req_mps2)
            % apply basis functions to each sample point
            dist = ((acc_control_learning.bf_vx_mps - visz_vx_mps(i))./acc_control_learning.vx_width_mps).^2 + ...
                           ((acc_control_learning.bf_ay_req_mps2 - visz_ay_req_mps2(j))./acc_control_learning.ay_width_mps2).^2;
            Visz_Z(j, i) = weights*[exp(-dist); visz_ay_req_mps2(j)];
        end
    end
    figure;
    surface(Visz_X, Visz_Y, Visz_Z);
    xlabel('Velocity in mps');
    ylabel('Lateral acc. target in mps2');
    zlabel('Corrective steering angle in rad');
end

author: TUMFTM


text:
function output = ...
  EstimateInitialF0AndDeviation(probability_structure, trajectory, option)
% estimate F0 initial value from instantaneous frequency map, residual map
% and given trajectory
%
% output = ...
%   EstimateInitialF0AndDeviation(probability_structure, trajectory)
%
% output = ...
%   EstimateInitialF0AndDeviation(probability_structure, trajectory, option)
%
% Arguments
%
%   probability_structure: structure variable, following fields are used
%     inst_freq_map: matrix, time-frequency map of instantaneous frequency
%     residual_map: matrix, time-frequency map of residuals
%     center_frequency_list: vector, list of filter center frequency
%     channels_in_octave: scalar, number of channels in one octave
%  trajectory: vector: sequence of F0 trajectory in channel ID number
%  option: structure variable with the following fields
%    residual_margin: scalar, safeguard for residual selection
%    frequency_margin: scalar, safeguard for limiting instantaneous
%      frequency
%
% Return value
%
%  output: structure variable with the following fields
%    f0_initial: vector, initial estimate of F0 sequence (Hz)
%    sd_initial: vector, initial estimate of relative deviation
%    temporal_positions: vector, frame center location (s)
%    residual_margin: scalar, actual value used
%    frequency_margin: scalar, actual value used
%    sd_on_trajectory: vector, residual value on give trajectory
%    elapsed_time: scalar, time needed to process this function

% Copyright 2016 Google Inc. All Rights Reserved
% Author: hidekik@google.com (Hideki Kawahara)

start_tic = tic;
residual_margin = 1.5; % default value
frequency_margin = 1.15; % default value
narginchk(2, 3);
if nargin == 3
  if isfield(option, 'residual_margin')
    residual_margin = option.residual_margin;
  end;
  if isfield(option, 'frequency_margin')
    frequency_margin = option.frequency_margin;
  end;
end;
inst_freq_map = probability_structure.inst_freq_map;
residual_map = probability_structure.residual_map;
fc_list = probability_structure.center_frequency_list;
n_channels = size(inst_freq_map,1);
n_frames = size(inst_freq_map, 2);
f0_initial = zeros(n_frames, 1);
sd_initial = zeros(n_frames, 1);
sample_f0_floor = ...
  probability_structure.center_frequency_list(min(trajectory));
sample_f0_ceiling = ...
  probability_structure.center_frequency_list(max(trajectory));
sd_on_trajectory = zeros(n_frames, 1); % for debug
n_in_octave = probability_structure.channels_in_octave;
scan_base = -ceil(n_in_octave / 2 + 2):ceil(n_in_octave / 4);
for ii = 1:n_frames
  base_index = trajectory(ii) + scan_base;
  base_index = base_index(base_index > 0 & base_index <= n_channels);
  reference_residual = residual_map(trajectory(ii), ii);
  sd_on_trajectory(ii) = sqrt(reference_residual);
  base_index = base_index(residual_map(base_index, ii) < ...
    reference_residual * residual_margin);
  reference_if = fc_list(trajectory(ii));
  tmp_residual = residual_map(base_index, ii);
  tmp_if = inst_freq_map(base_index, ii);
  tmp_if = max(reference_if / frequency_margin, ...
    min(reference_if * frequency_margin, tmp_if));
  if length(base_index) > 1
    best_weight = GetBestMixingWeight(sqrt(tmp_residual));
  else
    best_weight = 1;
  end;
  f0_initial(ii) = sum(best_weight .* tmp_if);
  sd_initial(ii) = sqrt(sum(best_weight.^2 .* tmp_residual));
end;
output = ...
  struct('f0_initial', max(sample_f0_floor, ...
                           min(sample_f0_ceiling, f0_initial)), ...
         'sd_initial', sd_initial, ...
         'temporal_positions', probability_structure.temporal_positions, ...
         'sample_f0_floor', sample_f0_floor, ...
         'sample_f0_ceiling', sample_f0_ceiling, ...
         'residual_margin', residual_margin, ...
         'frequency_margin', frequency_margin, ...
         'sd_on_trajectory', sd_on_trajectory, ...
         'elapsed_time', toc(start_tic));
end

author: google


text:
function [Y, lbls, Ytest, lblstest] = lvmLoadData(dataset, seedVal)

% LVMLOADDATA Load a latent variable model dataset.
% FORMAT
% DESC loads a data set for a latent variable modelling problem.
% ARG dataset : the name of the data set to be loaded. Currently
% the possible names are 'robotWireless',
% 'robotTwoLoops', 'robotTraces', 'robotTracesTest', 'cmu35gplvm',
% 'cmu35Taylor', 'cmu35walkJog', 'vowels', 'stick', 'brendan',
% 'digits', 'twos', 'oil', 'oilTest', 'oilValid', 'oil100',
% 'swissRoll', 'missa'.
% RETURN Y : the training data loaded in.
% RETURN lbls : a set of labels for the data (if there are no
% labels it is empty).
% RETURN Ytest : the test data loaded in. If no test set is
% available it is empty.
% RETURN lblstest : a set of labels for the test data (if there are
% no labels it is empty).
%
% SEEALSO : mapLoadData, datasetsDirectory
%
% COPYRIGHT : Neil D. Lawrence, 2004, 2005, 2006, 2008, 2009

% DATASETS

  if nargin > 1
    randn('seed', seedVal)
    rand('seed', seedVal)
  end

  % get directory

  baseDir = datasetsDirectory;
  dirSep = filesep;
  lbls = [];
  lblstest = [];
  switch dataset
   case 'YaleSubset6_1'
          load([baseDir 'YaleSubset6.mat']);
          lbls = [height;width];
   case 'missa'
       load([baseDir 'miss-americaHD.mat']);
       lbls = [288 ;360];
   case 'movielens'
    try
      load([baseDir 'movielens.mat']);
    catch
      [void, errid] = lasterr;
      if strcmp(errid, 'MATLAB:load:couldNotReadFile');
        fileName = [baseDir dirSep 'movielens' dirSep 'large' dirSep 'ratings.dat'];
        [users, films, ratings, timeStamp] = textread(fileName, '%n::%n::%n::%n');
        ind = randperm(size(users, 1));
        users = users(ind, :);
        films = films(ind, :);
        ratings = ratings(ind, :);
        numUsers = max(users);
        numFilms = max(films);

        numRatings = size(users, 1);
        numTrainRatings = ceil(0.8*numRatings);
        Y = spalloc(numFilms, numUsers, numTrainRatings);
        Ytest = spalloc(numFilms, numUsers, numRatings-numTrainRatings);
        indTrain = sub2ind(size(Y), films(1:numTrainRatings), users(1:numTrainRatings));
        indTest = sub2ind(size(Ytest), films(numTrainRatings+1:numRatings), users(numTrainRatings+1:numRatings));
        Y(indTrain) = ratings(1:numTrainRatings);
        Ytest(indTest) = ratings(numTrainRatings+1:numRatings);
        save([baseDir 'movielens.mat'], 'Y', 'lbls', 'Ytest', 'lblstest');
      else
        error(lasterr);
      end
    end

   case {'movielensSmall1', 'movielensSmall2', 'movielensSmall3', 'movielensSmall4', 'movielensSmall5'}

    partNo = str2num(dataset(end));
    uTrain = load([baseDir dirSep 'movielens' dirSep 'small' dirSep 'u' num2str(partNo) '.base']);
    numUsers = max(uTrain(:, 1));
    numFilms = max(uTrain(:, 2));
    numRatings = size(uTrain, 1);
    Y = spalloc(numFilms, numUsers, numRatings);

    for i = 1:size(uTrain, 1);
      Y(uTrain(i, 2), uTrain(i, 1)) = uTrain(i, 3);
    end
    meanY = mean(Y(find(Y)));
    Y(find(Y)) = (Y(find(Y))-meanY);
    uTest = load([baseDir dirSep 'movielens' dirSep 'small' dirSep 'u' num2str(partNo) '.test']);
    numTestRatings = size(uTest, 1);
    Ytest = spalloc(numFilms, numUsers, numTestRatings);
    for i = 1:size(uTest, 1);
      Ytest(uTest(i, 2), uTest(i, 1)) = uTest(i, 3);
    end
    Ytest(find(Ytest)) = (Ytest(find(Ytest))-meanY);
    stdY = sqrt(var(Y(find(Y))));
    Y = Y/stdY;
    Ytest = Ytest/stdY;

   case 'crabs'
    fid = fopen([baseDir 'crabs.dat']);
    if fid == -1
      error(['No such file name ' fileName])
    end
    readLine = fgets(fid);
    counter = 0;
    data = [];
    while readLine ~= -1
      readLine = fgets(fid);
      counter = counter + 1;
      parts = stringSplit(readLine, ' ');
      if parts{1} == 'B' && parts{2} == 'F'
        lbls(counter, :) = [0 0 0 1];
      elseif parts{1} == 'B' && parts{2} == 'M'
        lbls(counter, :) = [0 0 1 0];
      elseif parts{1} == 'O' && parts{2} == 'F'
        lbls(counter, :) = [0 1 0 0];
      else
        lbls(counter, :) = [1 0 0 0];
      end
      counter2 = 0;
      for i = 3:size(parts, 2)
        if ~isempty(str2num(parts{i}))
          counter2 = counter2 + 1;
          data(counter, counter2) = str2num(parts{i});
        end
      end
    end
    fclose(fid);
    data = data(:, 2:end);
    Y = data./repmat(sum(data, 2), 1, size(data, 2));

   case 'spellman'
    try
      load([baseDir 'spellman.mat']);

    catch
      [void, errid] = lasterr;
      if strcmp(errid, 'MATLAB:load:couldNotReadFile');
        fileName = ([baseDir 'combined.csv']);
        [geneName, ...
         y1, y2, y3, y4, y5, y6, y7, y8, y9, ...
         y10, y11, y12, y13, y14, y15, y16, y17, y18, ...
         y19, y20, y21, y22, y23, y24, y25, y26, y27, ...
         y28, y29, y30, y31, y32, y33, y34, y35, y36, ...
         y37, y38, y39, y40, y41, y42, y43, y44, y45, ...
         y46, y47, y48, y49, y50, y51, y52, y53, y54, ...
         y55, y56, y57, y58, y59, y60, y61, y62, y63, ...
         y64, y65, y66, y67, y68, y69, y70, y71, y72, ...
         y73, y74, y75, y76, y77] = ...
            textread(fileName, ['%q ' ...
                            '%f %f %f %f %f %f %f %f %f' ...
                            '%f %f %f %f %f %f %f %f %f' ...
                            '%f %f %f %f %f %f %f %f %f' ...
                            '%f %f %f %f %f %f %f %f %f' ...
                            '%f %f %f %f %f %f %f %f %f' ...
                            '%f %f %f %f %f %f %f %f %f' ...
                            '%f %f %f %f %f %f %f %f %f' ...
                            '%f %f %f %f %f %f %f %f %f' ...
                            '%f %f %f %f %f'], ...
                     'headerlines', 1, 'whitespace', '', 'delimiter', ',');

        Y = [y1 y2 y3 y4 y5 y6 y7 y8 y9 ...
             y10, y11, y12, y13, y14, y15, y16, y17, y18, ...
             y19, y20, y21, y22, y23, y24, y25, y26, y27, ...
             y28, y29, y30, y31, y32, y33, y34, y35, y36, ...
             y37, y38, y39, y40, y41, y42, y43, y44, y45, ...
             y46, y47, y48, y49, y50, y51, y52, y53, y54, ...
             y55, y56, y57, y58, y59, y60, y61, y62, y63, ...
             y64, y65, y66, y67, y68, y69, y70, y71, y72, ...
             y73, y74, y75, y76, y77]';
        lbls = [repmat([1 0 0 0 0], 4, 1);
                repmat([0 1 0 0 0], 18,1);
                repmat([0 0 1 0 0], 24,1);
                repmat([0 0 0 1 0], 17,1);
                repmat([0 0 0 0 1], 14,1)];
        save([baseDir 'spellman.mat'], 'Y', 'lbls');
      else
        error(lasterr);
      end
    end
   case 'robotWireless'
    Ydat = parseWirelessData([baseDir 'uw-floor.txt']);
    Ydat = (Ydat + 85)/15;
    lbls = 'connect';
    Y = Ydat(1:215, :);
    Ytest = Ydat(216:end, :);

   case 'robotTwoLoops'
    Y = csvread([baseDir 'TwoLoops.slam'], 1, 0);
    Y = Y(1:floor(end/2), 4:end);
    Y(find(Y==-100))=-NaN;
    Y = (Y + 85)/15;
    lbls = 'connect';

   case 'robotTraces'
    Y = csvread([baseDir 'Trace-3rdFloor-01.uwar.slam'], 1, 0);
    Y = Y(1:floor(end/2), 4:end);
    Y(:, [3 4 38]) = []; % Remove columns of missing data.
    Y(find(Y==-100))=NaN;
    Y = (Y + 85)/15;
    lbls = 'connect';

   case 'robotTracesTest'
    Y = csvread([baseDir 'Trace-3rdFloor-01.uwar.slam'], 1, 0);
    Y = Y(ceil(end/2):end, 4:end);
    Y(:, [3 4 38]) = []; % Remove columns of missing data.
    Y(find(Y==-100))=NaN;
    Y = (Y + 85)/15;
    lbls = 'connect';

   case 'osr'
    t = linspace(-pi, pi, 256);
    window = (cos(t) + 1)/2;
    [x, fs] = wavread([baseDir 'osr' dirSep 'OSR_uk_000_0020_8k.wav']);
    [r, theta] = stft(x, fs, 256, window);
    Y = [r(:, 1:32) theta(:, 1:32)];

   case 'cmu49BalanceArm'
    [Y, lbls, Ytest, lblstest] = lvmLoadData('cmu49Balance');
    % Extract left arm only
    ind = [41:47 49:50];
    Y = Y(:, ind);
    Ytest = Ytest(:, ind);

   case 'cmu49Balance'

    try
      load([baseDir 'cmu49Balance.mat']);
    catch
      [void, errid] = lasterr;
      if strcmp(errid, 'MATLAB:load:couldNotReadFile');
        sampleEvery = 32;
        skel = acclaimReadSkel([baseDir 'mocap' dirSep 'cmu' dirSep '49' dirSep '49.asf']);
        examples = {'18', '19'};
        testExamples = {'20'};
        % Label differently for each sequence
        exlbls = eye(length(examples));
        testexlbls = eye(length(testExamples));
        totLength = 0;
        totTestLength = 0;
        for i = 1:length(examples)
          [tmpchan, skel] = acclaimLoadChannels([baseDir 'mocap' dirSep 'cmu' dirSep '49' dirSep '49_' ...
                              examples{i} '.amc'], skel);
          tY{i} = tmpchan(1:sampleEvery:end, :);
          tlbls{i} = repmat(exlbls(i, :), size(tY{i}, 1), 1);
          totLength = totLength + size(tY{i}, 1);
        end
        Y = zeros(totLength, size(tY{1}, 2));
        lbls = zeros(totLength, size(tlbls{1}, 2));
        endInd = 0;
        for i = 1:length(tY)
          startInd = endInd + 1;
          endInd = endInd + size(tY{i}, 1);
          Y(startInd:endInd, :) = tY{i};
          lbls(startInd:endInd, :) = tlbls{i};
        end
        for i = 1:length(testExamples)
          [tmpchan, skel] = acclaimLoadChannels([baseDir 'mocap' dirSep 'cmu' dirSep '49' dirSep '49_' ...
                              testExamples{i} '.amc'], skel);
          tYtest{i} = tmpchan(1:sampleEvery:end, :);
          tlblstest{i} = repmat(testexlbls(i, :), size(tYtest{i}, 1), 1);
          totTestLength = totTestLength + size(tYtest{i}, 1);
        end
        Ytest = zeros(totTestLength, size(tYtest{1}, 2));
        lblstest = zeros(totTestLength, size(tlblstest{1}, 2));
        endInd = 0;
        for i = 1:length(tYtest)
          startInd = endInd + 1;
          endInd = endInd + size(tYtest{i}, 1);
          Ytest(startInd:endInd, :) = tYtest{i};
          lblstest(startInd:endInd, :) = tlblstest{i};
        end
        save([baseDir 'cmu49Balance.mat'], 'Y', 'lbls', 'Ytest', 'lblstest');
      else
        error(lasterr);
      end
    end

   case 'cmu07WalkFeet'
    [Y, lbls, Ytest, lblstest] = lvmLoadData('cmu07Walk');
    % Extract left arm only
    ind = [7:12 14:19];
    Y = Y(:, ind);
    Ytest = Ytest(:, ind);

   case 'cmu07Walk'

    try
      load([baseDir 'cmu07Walk.mat']);
    catch
      [void, errid] = lasterr;
      if strcmp(errid, 'MATLAB:load:couldNotReadFile');
        sampleEvery = 4;
        skel = acclaimReadSkel([baseDir 'mocap' dirSep 'cmu' dirSep '07' dirSep '07.asf']);
        examples = {'01', '02'};
        testExamples = {'03'};
        % Label differently for each sequence
        exlbls = eye(length(examples));
        testexlbls = eye(length(testExamples));
        totLength = 0;
        totTestLength = 0;
        for i = 1:length(examples)
          [tmpchan, skel] = acclaimLoadChannels([baseDir 'mocap' dirSep 'cmu' dirSep '07' dirSep '07_' ...
                              examples{i} '.amc'], skel);
          tY{i} = tmpchan(1:sampleEvery:end, :);
          tlbls{i} = repmat(exlbls(i, :), size(tY{i}, 1), 1);
          totLength = totLength + size(tY{i}, 1);
        end
        Y = zeros(totLength, size(tY{1}, 2));
        lbls = zeros(totLength, size(tlbls{1}, 2));
        endInd = 0;
        for i = 1:length(tY)
          startInd = endInd + 1;
          endInd = endInd + size(tY{i}, 1);
          Y(startInd:endInd, :) = tY{i};
          lbls(startInd:endInd, :) = tlbls{i};
        end
        for i = 1:length(testExamples)
          [tmpchan, skel] = acclaimLoadChannels([baseDir 'mocap' dirSep 'cmu' dirSep '07' dirSep '07_' ...
                              testExamples{i} '.amc'], skel);
          tYtest{i} = tmpchan(1:sampleEvery:end, :);
          tlblstest{i} = repmat(testexlbls(i, :), size(tYtest{i}, 1), 1);
          totTestLength = totTestLength + size(tYtest{i}, 1);
        end
        Ytest = zeros(totTestLength, size(tYtest{1}, 2));
        lblstest = zeros(totTestLength, size(tlblstest{1}, 2));
        endInd = 0;
        for i = 1:length(tYtest)
          startInd = endInd + 1;
          endInd = endInd + size(tYtest{i}, 1);
          Ytest(startInd:endInd, :) = tYtest{i};
          lblstest(startInd:endInd, :) = tlblstest{i};
        end
        save([baseDir 'cmu07Walk.mat'], 'Y', 'lbls', 'Ytest', 'lblstest');
      else
        error(lasterr);
      end
    end

   case 'cmu35gplvm'
    [Y, lbls, Ytest, lblstest] = lvmLoadData('cmu35WalkJog');
    skel = acclaimReadSkel([baseDir 'mocap' dirSep 'cmu' dirSep '35' dirSep '35.asf']);
    [tmpchan, skel] = acclaimLoadChannels([baseDir 'mocap' dirSep 'cmu' dirSep '35' dirSep '35_01.amc'], skel);

    Ytest = Ytest(find(lblstest(:, 2)), :);
    lblstest = lblstest(find(lblstest(:, 2)), 2);

    %left indices
    xyzInd = [2];
    xyzDiffInd = [1 3];
    rotInd = [4 6];
    rotDiffInd = [5];
    generalInd = [7:38 41:47 49:50 53:59 61:62];

    jointAngles  = asin(sin(pi*Y(:, generalInd)/180));
    jointAnglesTest  = asin(sin(pi*Ytest(:, generalInd)/180));

    endInd = [];
    for i = 1:size(lbls, 2)
      endInd = [endInd max(find(lbls(:, i)))];
    end
    catJointAngles = [];
    xyzDiff = [];
    catSinCos = [];
    startInd = 1;
    for i = 1:length(endInd)
      ind1 = startInd:endInd(i)-1;
      ind2 = startInd+1:endInd(i);
      catJointAngles = [catJointAngles; ...
                        jointAngles(ind2, :)];
      xyzDiff = [xyzDiff;
                 Y(ind1, xyzDiffInd) - Y(ind2, xyzDiffInd) ...
                 Y(ind2, xyzInd)];
      catSinCos = [catSinCos; ...
                   sin(pi*Y(ind2, rotInd)/180) ...
                   sin(pi*Y(ind1, rotDiffInd)/180)-sin(pi*Y(ind2, rotDiffInd)/180) ...
                   cos(pi*Y(ind2, rotInd)/180) ...
                   cos(pi*Y(ind1, rotDiffInd)/180)-cos(pi*Y(ind2, rotDiffInd)/180)];
      startInd = endInd(i)+1;
    end
    Y = [catJointAngles xyzDiff catSinCos];
    lbls = [];

    endInd = [];
    for i = 1:size(lblstest, 2)
      endInd = [endInd max(find(lblstest(:, i)))];
    end
    catJointAnglesTest = [];
    xyzDiffTest = [];
    catSinCosTest = [];
    startInd = 1;
    for i = 1:length(endInd)
      ind1 = startInd:endInd(i)-1;
      ind2 = startInd+1:endInd(i);
      catJointAnglesTest = [catJointAnglesTest; ...
                          jointAnglesTest(ind2, :)];
      xyzDiffTest = [xyzDiffTest;
                     Ytest(ind1, xyzDiffInd) - Ytest(ind2, xyzDiffInd) ...
                     Ytest(ind2, xyzInd)];
      catSinCosTest = [catSinCosTest; ...
                       sin(pi*Ytest(ind2, rotInd)/180) ...
                       sin(pi*Ytest(ind1, rotDiffInd)/180)-sin(pi*Ytest(ind2, rotDiffInd)/180) ...
                       cos(pi*Ytest(ind2, rotInd)/180) ...
                       cos(pi*Ytest(ind1, rotDiffInd)/180)-cos(pi*Ytest(ind2, rotDiffInd)/180)];
      startInd = endInd(i)+1;
    end
    Ytest = [catJointAnglesTest xyzDiffTest catSinCosTest];
    lblstest = [];

   case 'cmu35Taylor'
    % The CMU 35 data set as Graham Taylor et al. have it in their
    % NIPS 2006 paper.
    [rawY, lbls, rawYtest, lblstest] = lvmLoadData('cmu35WalkJog');
    rawYtest = rawYtest(45:end, :);
    lblstest = lblstest(45:end, 2);
    skel = acclaimReadSkel([baseDir 'mocap' dirSep 'cmu' dirSep '35' dirSep '35.asf']);
    [tmpchan, skel] = acclaimLoadChannels([baseDir 'mocap' dirSep 'cmu' dirSep '35' dirSep '35_01.amc'], skel);
    posInd = [1 2 3];
    rotInd = [4 5 6];
    generalInd = [7:38 41:47 49:50 53:59 61:62];

    rawY = [rawY(:, posInd) sind(rawY(:, rotInd)) cosd(rawY(:, rotInd)) asind(sind(rawY(:, generalInd)))];
    rawYtest = [rawYtest(:, posInd) sind(rawYtest(:, rotInd)) cosd(rawYtest(:, rotInd)) asind(sind(rawYtest(:, generalInd)))];

    endInd = [];
    for i = 1:size(lbls, 2)
      endInd = [endInd max(find(lbls(:, i)))];
    end
    Y = [];
    startInd = 1;
    for i = 1:length(endInd)
      ind1 = startInd:endInd(i)-1;
      ind2 = startInd+1:endInd(i);
      Y = [Y; rawY(ind1, :) rawY(ind2, :)];
      startInd = endInd(i)+1;
    end
    lbls = [];

    endInd = [];
    for i = 1:size(lblstest, 2)
      endInd = [endInd max(find(lblstest(:, i)))];
    end
    Ytest = [];
    startInd = 1;
    for i = 1:length(endInd)
      ind1 = startInd:endInd(i)-1;
      ind2 = startInd+1:endInd(i);
      Ytest = [Ytest; rawYtest(ind1, :) rawYtest(ind2, :)];
      startInd = endInd(i)+1;
    end
    lblstest = [];

   case {'dur', 'cmp'}
    [Y, void, lbls] = synthLoadData(dataset);

   case {'cmpdur'}
    [Y1, void, lbls] = synthLoadData('cmp');
    [Y2, void, lbls] = synthLoadData('dur');
    Y = [Y1, Y2];

   case {'grid_vowels', ...
         'grid_consonants', ...
         'grid_approximants', ...
         'grid_uplosives', ...
         'grid_vplosives', ...
         'grid_plosives', ...
         'grid_nasals', ...
         'grid_ufricatives', ...
         'grid_vfricatives', ...
         'grid_fricatives', ...
         'grid_affricatives', ...
         'grid_silence', ...
         'grid_other', ...
         'grid_all'}
    try
      load([baseDir dataset '.mat']);
    catch

      fullList = {'iy2', 'iy4', 'iy5', 'iy6', ...     % Vowels (iy, ih, eh, ey, ae, aa, aw, ah, ao, ow, uw, ax, ia)
                'ih1', 'ih3', 'ih5', ...
                'eh1', 'eh2', 'eh4', 'eh5', 'eh6', ...
                'ey1', 'ey4', 'ey5', ...
                'ae3', ...
                'aa4', ...
                'aw6', ...
                'ay2', 'ay3', 'ay4', 'ay5', ...
                'ah5', ...
                'ao5', ...
                'ow4', 'ow5', ...
                'uw2', 'uw4', 'uw5', 'uw6', ...
                'ax4', 'ax6', ...
                'ia5', ...
                'w2', 'w3', 'w4', 'w5', ... % approximants -- kind of between a vowel and consonant (w, y, r, l)
                'y4', ...
                'r2', 'r5', ...
                'l1', 'l2', 'l4', 'l6', ...
                'b1', 'b2', 'b3', 'b4', ... % Consonants  - plosives, voiced (b, d, g)
                'd2', 'd4', ...
                'g2', 'g6', ...
                'p1', 'p4', 'p6', ... % Consonants - plosives, unvoiced (p, t, k)
                't1', 't2', 't3', 't4', 't5',...
                'k4', 'k5', ...
                'm4', ... % Consonants - nasals (m, n)
                'n1', 'n2', 'n3', 'n4', 'n5', 'n6', ...
                'f4', 'f5', ... % Consonants - fricatives, unvoiced (f, s, th, sh)
                's1', 's4', 's5', 's6', ...
                'th5', ...
                'v4', 'v5', ... % Consonant - fricatives, voiced (v, z, dh, zh)
                'z4', 'z5', 'z6', ...
                'dh3', ...
                ...%(but consonsant can also be grouped strong/weak:   strong = s, sh, z, zh;  weak=f,v,th,dh)
                'ch4', ...  % Consonant - affricative  (basically a blend between a plosive and a fricative) (ch  =  t + sh, jh = d + zh)
                'jh4', ...
                'sil', ... %silence
                'sp'}; % ask Jon

      vowelInd = 1:33;
      consonantInd = 34:84;
      approximantInd = 34:44;
      vplosivesInd = 45:52;
      uplosivesInd = 53:62;
      nasalsInd = 63:69;
      ufricativesInd = 70:76;
      vfricativesInd = 77:82;
      affricativesInd = 83:84;
      silenceInd = 85;
      otherInd = 86;
      constanants = {};
      all = {};
      switch dataset(6:end)
       case 'vowels'
        phoneList = fullList(vowelInd);
       case 'consonants'
        phoneList = fullList(consonantInd);
       case 'approximants'
        phoneList = fullList(approximantInd);
       case 'uplosives'
        phoneList = fullList(uplosivesInd);
       case 'vplosives'
        phoneList = fullList(vplosivesInd);
       case 'plosives'
        phoneList = fullList([uplosivesInd vplosivesInd]);
       case 'nasals'
        phoneList = fullList(nasalsInd);
       case 'ufricatives'
        phoneList = fullList(ufricativesInd);
       case 'vfricatives'
        phoneList = fullList(vfricativesInd);
       case 'fricatives'
        phoneList = fullList([ufricativesInd vfricativesInd]);
       case 'affricatives'
        phoneList = fullList(affricativesInd);
       case 'silence'
        phoneList = fullList(silenceInd);
       case 'other'
        phoneList = fullList(otherInd);
       case 'all'
        phoneList = fullList;
      end

      motherDir = '/local/data/datasets/synthesis/models/';
      fileLoc = '/ver3/cmp/monophoneC.mmf';
      numId = 34;
      i = 1;
      idStr = ['id' num2str(i)];
      fileName = [motherDir idStr fileLoc];
      [meanValCell, varValCell] = htkLoadMmf(phoneList, fileName);
      meanVals = [meanValCell{:, :}];
      varVals = [varValCell{:, :}];
      meanVals = zeros(numId, size(meanVals, 2));
      varVals = zeros(numId, size(varVals, 2));

      for i = 1:numId
        idStr = ['id' num2str(i)];
        fileName = [motherDir idStr fileLoc];
        [meanValCell, varValCell] = htkLoadMmf(phoneList, fileName);
        meanVals(i, :) = [meanValCell{:, :}];
        varVals(i, :) = [varValCell{:, :}];
      end
      meanVals = meanVals./sqrt(varVals);
      varVals = 0.5*log(varVals);
      meanMean = mean(meanVals);
      meanLogStd = mean(varVals);
      for i = 1:numId
        meanVals(i, :) = meanVals(i, :) - meanMean;
        varVals(i, :) = varVals(i, :) - meanLogStd;
      end
      varMeanVals = var(meanVals(:));
      varVarVals = var(varVals(:));
      Y = [meanVals/sqrt(varMeanVals) varVals/sqrt(varVarVals)];
      lbls = [1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 0, 1, 0, 1, 0, ...
              0, 0, 0, 0, 0, 1, 1, 1, 0, 1, 0, 1, 0, 0]';
      lbls = [lbls ~lbls];
      save([baseDir dataset '.mat'], 'Y', 'lbls', 'meanMean', 'meanLogStd', 'varMeanVals', 'varVarVals');
    end
    if size(lbls, 2) == 1
      lbls = [lbls ~lbls];
      save([baseDir dataset '.mat'], 'Y', 'lbls', 'meanMean', 'meanLogStd', 'varMeanVals', 'varVarVals');
    end
    retLbls = cell(2,1);
    retLbls{1} = lbls;
    %/~retLbls{2} = {'JB', 'Martin Cooke', 'Stuart Cunningham', 'Sue Harding', ...
%                  'Jonny Laidler', 'Stuart Wriggley', 'Gillian', 'Mike Stannett', ...
%                  'John Karn', 'UG', 'Anna', 'UG', ...
%                  'UG', 'UG', 'Female', 'UG', ...
%                  'UG', 'UG',  'UG', 'UG', ...
%                  'UG', 'UG', 'UG', 'UG', ...
%                  'UG', 'James Carmichael', 'Rob Mill', 'Matt Gibson', ...
%                  'Sarah Simpson', 'Vinny Wan', 'Lucy', 'NDL', ...
%                  'UG', 'UG'};
                  %~/
                  retLbls{2} = {'JB', 'MPC', 'Stuart Cunningham', 'Sue Harding', ...
                                'Jonny Laidler', 'Stuart Wriggley', 'South Yorkshire G', 'MS', ...
                                'JK', 'UG 1', 'South Yorkshire A', 'UG 2', ...
                                'UG 3', 'UG 4', 'UG 5', 'UG 6', ...
                                'UG 7', 'UG 8',  'UG 9', 'UG 10', ...
                                'UG 11', 'UG 12', 'UG 13', 'UG14', ...
                                'UG', 'West Indies', 'Rob Mill', 'Scottish', ...
                                'Sarah Simpson', 'Vinny Wan', 'South Yorkshire L', 'NDL', ...
                                'UG 15', 'UG 16'};
    lbls = retLbls;

   case 'cmu35WalkJog'
    try
      load([baseDir 'cmu35WalkJog.mat']);
    catch
      [void, errid] = lasterr;
      if strcmp(errid, 'MATLAB:load:couldNotReadFile');
        skel = acclaimReadSkel([baseDir 'mocap' dirSep 'cmu' dirSep '35' dirSep '35.asf']);
        examples = ...
            {'01', '02', '03', '04', '05', '06', '07', '08', '09', '10', ...
             '11', '12', '13', '14', '15', '16', '17', '19', '20', ...
             '21', '22', '23', '24', '25', '26', '28', '30', '31', '32', '33', '34'};
        testExamples = {'18', '29'};
        % Label differently for each sequence
        exlbls = eye(31);
        testexlbls = eye(2);
        totLength = 0;
        totTestLength = 0;
        for i = 1:length(examples)
          [tmpchan, skel] = acclaimLoadChannels([baseDir 'mocap' dirSep 'cmu' dirSep '35' dirSep '35_' ...
                              examples{i} '.amc'], skel);
          tY{i} = tmpchan(1:4:end, :);
          tlbls{i} = repmat(exlbls(i, :), size(tY{i}, 1), 1);
          totLength = totLength + size(tY{i}, 1);
        end
        Y = zeros(totLength, size(tY{1}, 2));
        lbls = zeros(totLength, size(tlbls{1}, 2));
        endInd = 0;
        for i = 1:length(tY)
          startInd = endInd + 1;
          endInd = endInd + size(tY{i}, 1);
          Y(startInd:endInd, :) = tY{i};
          lbls(startInd:endInd, :) = tlbls{i};
        end
        for i = 1:length(testExamples)
          [tmpchan, skel] = acclaimLoadChannels([baseDir 'mocap' dirSep 'cmu' dirSep '35' dirSep '35_' ...
                              testExamples{i} '.amc'], skel);
          tYtest{i} = tmpchan(1:4:end, :);
          tlblstest{i} = repmat(testexlbls(i, :), size(tYtest{i}, 1), 1);
          totTestLength = totTestLength + size(tYtest{i}, 1);
        end
        Ytest = zeros(totTestLength, size(tYtest{1}, 2));
        lblstest = zeros(totTestLength, size(tlblstest{1}, 2));
        endInd = 0;
        for i = 1:length(tYtest)
          startInd = endInd + 1;
          endInd = endInd + size(tYtest{i}, 1);
          Ytest(startInd:endInd, :) = tYtest{i};
          lblstest(startInd:endInd, :) = tlblstest{i};
        end
        save([baseDir 'cmu35WalkJog.mat'], 'Y', 'lbls', 'Ytest', 'lblstest');
      else
        error(lasterr);
      end
    end

   case 'vowels'
    load([baseDir 'jon_vowel_data']);
    Y = [a_raw; ae_raw; ao_raw; ...
         e_raw; i_raw; ibar_raw; ...
         o_raw; schwa_raw; u_raw];
    Y(:, [13 26]) = [];
    lbls = [];
    for i = 1:9
      lbl = zeros(1, 9);
      lbl(i) = 1;
      lbls = [lbls; repmat(lbl, size(a_raw, 1), 1)];
    end

   case 'stick'
    Y = mocapLoadTextData([baseDir 'run1']);
    Y = Y(1:4:end, :);
    lbls = 'connect';

   case 'brendan'
    load([baseDir 'frey_rawface.mat']);
    Y = double(ff)';

   case 'usps'

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

    load([baseDir 'usps_train.mat']);
    [ALL_T, sortIndices] = sort(ALL_T);
    ALL_DATA = ALL_DATA(sortIndices(:), :);
    Y = ALL_DATA;
    lbls = zeros(size(ALL_T,1),10);
    for digit = 0:9;
      ind1 = min(find(ALL_T==digit));
      ind2 = max(find(ALL_T==digit));
      lbls(ind1:ind2,digit+1) = ones(size(ind1:ind2,2),1);
    end

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

    load([baseDir 'usps_test.mat']);
    [ALL_T, sortIndices] = sort(ALL_T);
    ALL_DATA = ALL_DATA(sortIndices(:), :);
    Ytest = ALL_DATA;
    lblstest = zeros(size(ALL_T,1),10);
    for digit = 0:9;
      ind1 = min(find(ALL_T==digit));
      ind2 = max(find(ALL_T==digit));
      lblstest(ind1:ind2,digit+1) = ones(size(ind1:ind2,2),1);
    end

   case 'digits'

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

    load([baseDir 'usps_train.mat']);
    % Extract 600 of digits 0 to 4
    [ALL_T, sortIndices] = sort(ALL_T);
    ALL_DATA = ALL_DATA(sortIndices(:), :);
    Y = [];
    lbls = [];
    numEachDigit = 600;
    for digit = 0:4;
      firstDigit = min(find(ALL_T==digit));
      Y = [Y; ALL_DATA(firstDigit:firstDigit+numEachDigit-1, :)];
      lbl = zeros(1, 5);
      lbl(digit+1) = 1;
      lbls = [lbls; repmat(lbl, numEachDigit, 1)];
    end

   case 'twos'
    % load data
    load([baseDir 'twos']);
    Y = 2*a-1;

   case 'oil'
    load([baseDir '3Class.mat']);
    Y = DataTrn;
    lbls = DataTrnLbls;

   case 'oilTest'
    load([baseDir '3Class.mat']);
    Y = DataTst;
    lbls = DataTstLbls;

   case 'oilValid'
    load([baseDir '3Class.mat']);
    Y = DataVdn;
    lbls = DataVdnLbls;

   case 'oil100'
    randn('seed', 1e5);
    rand('seed', 1e5);
    load([baseDir '3Class.mat']);
    Y = DataTrn;
    lbls = DataTrnLbls;
    indices = randperm(size(Y, 1));
    indices = indices(1:100);
    Y = Y(indices, :);
    lbls = lbls(indices, :);

   case 'swissRoll'
    load([baseDir 'swiss_roll_data']);
    Y = X_data(:, 1:1000)';

   case 'swissRollFull'
    load([baseDir 'swiss_roll_data']);
    Y = X_data(:, 1:3000)';

    %/~
   case 'fourWalks'
    try
      load([baseDir 'fourWalks.mat']);
    catch
      [void, errid] = lasterr;
      if strcmp(errid, 'MATLAB:load:couldNotReadFile');

        subj = {'35', '10', '12', '16'};
        examples = {'02', '04', '01', '15'};
        exlbls = eye(4);
        startVals = [55, 222, 22, 62];
        endVals = [338, 499, 328, 342];
        diffInd = [1 3 5];
        totLength = 0;
        totTestLength = 0;
        for i = 1:length(examples)
          skel = acclaimReadSkel([baseDir 'mocap' dirSep 'cmu' dirSep subj{i} dirSep subj{i} ...
                              '.asf']);
          [tmpchan, skel] = acclaimLoadChannels([baseDir 'mocap' dirSep 'cmu' dirSep subj{i} dirSep subj{i} '_' ...
                              examples{i} '.amc'], skel);
          tY{i} = tmpchan(startVals(i):4:endVals(i), :);
          % Make given indices differences.
          tY{i}(:, diffInd) = tmpchan(startVals(i)+1:4:endVals(i)+1, diffInd)...
              - tY{i}(:, diffInd);
          tlbls{i} = repmat(exlbls(i, :), size(tY{i}, 1), 1);
          totLength = totLength + size(tY{i}, 1);
        end
        Y = zeros(totLength, size(tY{1}, 2));
        lbls = zeros(totLength, size(tlbls{1}, 2));
        endInd = 0;
        for i = 1:length(tY)
          startInd = endInd + 1;
          endInd = endInd + size(tY{i}, 1);
          Y(startInd:endInd, :) = tY{i};
          lbls(startInd:endInd, :) = tlbls{i};
        end
        Ytest = [];
        lblstest = [];
        subj = {'35', '12', '16', '12', '07', '07', '08', '08'};
        examples = {'03', '02', '21', '03', '01', '02', '01', '02'};
        testexlbls = eye(8);
        startVals = ones(1, 8);
        endVals = 200*ones(1, 8);
        totLength = 0;
        totTestLength = 0;
        for i = 1:length(examples)
          skel = acclaimReadSkel([baseDir 'mocap' dirSep 'cmu' dirSep subj{i} dirSep subj{i} ...
                              '.asf']);
          [tmpchan, skel] = acclaimLoadChannels([baseDir 'mocap' dirSep 'cmu' dirSep subj{i} dirSep subj{i} '_' ...
                              examples{i} '.amc'], skel);
          tYtest{i} = tmpchan(startVals(i):4:endVals(i), :);
          tYtest{i}(:, diffInd) = tmpchan(startVals(i)+1:4:endVals(i)+1, diffInd)...
              - tYtest{i}(:, diffInd);
          tlblstest{i} = repmat(testexlbls(i, :), size(tYtest{i}, 1), 1);
          totLength = totLength + size(tYtest{i}, 1);
        end
        Ytest = zeros(totTestLength, size(tYtest{1}, 2));
        lblstest = zeros(totTestLength, size(tlblstest{1}, 2));
        endInd = 0;
        for i = 1:length(tYtest)
          startInd = endInd + 1;
          endInd = endInd + size(tYtest{i}, 1);
          Ytest(startInd:endInd, :) = tYtest{i};
          lblstest(startInd:endInd, :) = tlblstest{i};
        end
        save([baseDir 'fourWalks.mat'], 'Y', 'lbls', 'Ytest', 'lblstest');
      else
        error(lasterr);
      end
    end

   case 'walkSitJog'
    Y = mocapLoadTextData([baseDir 'walkSitJog']);
    Y = Y(1:4:end, :)*200;

   case 'robotAtrium'
    Y = parseWirelessData([baseDir 'uw-atrium.txt']);
    Y = (Y + 85)/15;
    Y = Y(1:281, :);
    lbls = 'connect'

    Ytest = parseWirelessData([baseDir 'uw-atrium.txt']);
    Ytest = (Ytest + 85)/15;
    Ytest = Ytest(282:end, :);

   case 'accWalkSitJog'
    try
      load([baseDir 'walkSitJog.mat']);
    catch
      [void, errid] = lasterr;
      if strcmp(errid, 'MATLAB:load:couldNotReadFile')
        skel = acclaimReadSkel([baseDir 'mocap' dirSep 'cmu' dirSep '86' dirSep '86.asf']);
        [channels, skel] = acclaimLoadChannels([baseDir 'mocap' dirSep 'cmu' dirSep '86' dirSep '86_10.amc'], skel);
        % Remove root z & x pos channels.
        channels(:, [1 3]) = zeros(size(channels, 1), 2);
        Y = channels(1:4:end, :);
        save([baseDir 'walkSitJog.mat'], 'Y')
      else
        error(lasterr);
      end
    end
   case 'horse'
    load([baseDir 'horse.dat']);
    horse([133, 309], :) = [];
    ordinalAtr = [1 2 7 8 9 10 11 12 13 14 15 17 18 21];
    horse(:, ordinalAtr) = horse(:, ordinalAtr) - 1;

    % Set labels as outcome --- 1 lived, 2 died, 3 was euthanized
    lbls = zeros(size(horse, 1), 3);
    for i = 1:size(horse, 1)
      lbls(i, horse(i, 23)) = 1;
    end

    horse(:, 23:28) = [];
    horse(:, 3) = [];
    Y = horse;
    % Normalise Gaussian variables.
    for i = [3 4 5 15 18 19 21];
      va = var(Y(find(~isnan(Y(:, i))), i));
      Y(:, i) = Y(:, i)/va;
    end

   case 'pitch'
    load('-ascii', [baseDir 'bballpitch.mat']);
    Y = bballpitch(:,8:end);

   case 'ratemaps'
    % Fix seeds
    Y = loadRateMap([baseDir 'ratemaps'])';
    Y = Y(1:10000, :);

   case 'cepstral'
    Y = load([baseDir 'cepvecs']);
    Y = Y(1:5000, :);

   case 'vowels90'
    load([baseDir 'jon_vowel_data']);
    Y = [a_raw(1:30:end, :); ae_raw(1:30:end, :); ao_raw(1:30:end, :); ...
         e_raw(1:30:end, :); i_raw(1:30:end, :); ibar_raw(1:30:end, :); ...
         o_raw(1:30:end, :); schwa_raw(1:30:end, :); u_raw(1:30:end, :)];
    lbls = [];
    for i = 1:9
      lbl = zeros(1, 9);
      lbl(i) = 1;
      lbls = [lbls; repmat(lbl, 10, 1)];
    end
   case 'vowels3'
    load([baseDir 'jon_vowel_data']);
    Y = [a_raw(1:3:end, :); ae_raw(1:3:end, :); ao_raw(1:3:end, :); ...
         e_raw(1:3:end, :); i_raw(1:3:end, :); ibar_raw(1:3:end, :); ...
         o_raw(1:3:end, :); schwa_raw(1:3:end, :); u_raw(1:3:end, : ...
                                                      )];
    lbls = [];
    for i = 1:9
      lbl = zeros(1, 9);
      lbl(i) = 1;
      lbls = [lbls; repmat(lbl, 100, 1)];
    end
   case 'breakdance'
    Y = mocapLoadTextData([baseDir 'EricCamper04'], 0);
    Y = Y(1:4:end, :);

   case 'dance'
    Y = mocapLoadTextData([baseDir 'dancemotion'], 0);
    Y = Y(1:4:end, :);

   case 'danceReduced'
    Y = mocapLoadTextData([baseDir 'dancemotion'], 0);
    Y = Y(1:16:end, :);
    %~/

   otherwise
    error('Unknown data set requested.')

  end
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
% OP_CURLU_CURLV_3D: assemble the matrix A = [a(i,j)], a(i,j) = (coeff curl u_j, curl v_i), with vector-valued curl.
%
%   mat = op_curlu_curlv_3d (spu, spv, msh, epsilon);
%   [rows, cols, values] = op_curlu_curlv_3d (spu, spv, msh, epsilon);
%
% INPUT:
%
%   spu:   structure representing the space of trial functions (see sp_vector/sp_evaluate_col)
%   spv:   structure representing the space of test functions  (see sp_vector/sp_evaluate_col)
%   msh:   structure containing the domain partition and the quadrature rule (see msh_cartesian/msh_evaluate_col)
%   coeff: physical parameter
%
% OUTPUT:
%
%   mat:    assembled matrix
%   rows:   row indices of the nonzero entries
%   cols:   column indices of the nonzero entries
%   values: values of the nonzero entries
%
% Copyright (C) 2009, 2010 Carlo de Falco
% Copyright (C) 2011, 2017 Rafael Vazquez
%
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.

%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program.  If not, see <http://www.gnu.org/licenses/>.

function varargout = op_curlu_curlv_3d (spu, spv, msh, coeff)

  rows = zeros (msh.nel * spu.nsh_max * spv.nsh_max, 1);
  cols = zeros (msh.nel * spu.nsh_max * spv.nsh_max, 1);
  values = zeros (msh.nel * spu.nsh_max * spv.nsh_max, 1);

  jacdet_weights = msh.jacdet .* msh.quad_weights .* coeff;

  ncounter = 0;
  for iel = 1:msh.nel
    if (all (msh.jacdet(:,iel)))
      curlu_iel = reshape (spu.shape_function_curls(:, :, :, iel), 3, msh.nqn, 1, spu.nsh_max);
      curlv_iel = reshape (spv.shape_function_curls(:, :, :, iel), 3, msh.nqn, spv.nsh_max, 1);

      jacdet_iel = reshape (jacdet_weights(:,iel), [1,msh.nqn,1,1]);

      jacdet_curlu = bsxfun (@times, jacdet_iel, curlu_iel);
      tmp1 = sum (bsxfun (@times, jacdet_curlu, curlv_iel), 1);
      elementary_values = reshape (sum (tmp1, 2), spv.nsh_max, spu.nsh_max);

      [rows_loc, cols_loc] = ndgrid (spv.connectivity(:,iel), spu.connectivity(:,iel));
      indices = rows_loc & cols_loc;
      rows(ncounter+(1:spu.nsh(iel)*spv.nsh(iel))) = rows_loc(indices);
      cols(ncounter+(1:spu.nsh(iel)*spv.nsh(iel))) = cols_loc(indices);
      values(ncounter+(1:spu.nsh(iel)*spv.nsh(iel))) = elementary_values(indices);
      ncounter = ncounter + spu.nsh(iel)*spv.nsh(iel);
    else
      warning ('geopdes:jacdet_zero_at_quad_node', 'op_curlu_curlv_3d: singular map in element number %d', iel)
    end
  end

  if (nargout == 1)
    varargout{1} = sparse (rows(1:ncounter), cols(1:ncounter), ...
                           values(1:ncounter), spv.ndof, spu.ndof);
  elseif (nargout == 3 || nargout == 0)
    varargout{1} = rows(1:ncounter);
    varargout{2} = cols(1:ncounter);
    varargout{3} = values(1:ncounter);
  else
    error ('op_curlu_curlv_3d: wrong number of output arguments')
  end

end

%% COPY OF THE FIRST VERSION OF THE FUNCTION (MORE UNDERSTANDABLE)
%
% function mat = op_curlu_curlv_3d (spu, spv, msh, coeff)
%
%   mat = spalloc(spv.ndof, spu.ndof, 1);
%   for iel = 1:msh.nel
%     if (all (msh.jacdet(:,iel)))
%       mat_loc = zeros (spv.nsh(iel), spu.nsh(iel));
%       for idof = 1:spv.nsh(iel)
%         ishc = reshape(spv.shape_function_curls(:,:,idof,iel), 3, []);
%         for jdof = 1:spu.nsh(iel)
%           jshc = reshape (spu.shape_function_curls(:,:,jdof,iel), 3, []);
% % The cycle on the quadrature points is vectorized
%           %for inode = 1:msh.nqn
%             mat_loc(idof, jdof) = mat_loc(idof, jdof) + ...
%               sum (msh.jacdet(:,iel) .* msh.quad_weights(:, iel) .*...
%               sum (ishc .* jshc, 1).' .* coeff (:, iel));
%           %end
%         end
%       end
%       mat(spv.connectivity(:,iel), spu.connectivity(:,iel)) = ...
%         mat(spv.connectivity(:,iel), spu.connectivity(:,iel)) + mat_loc;
%     else
%       warning ('geopdes:jacdet_zero_at_quad_node', 'op_curlu_curlv_3d: singular map in element number %d', iel)
%     end
%   end
%
% end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function [state] = tapas_sampler_mc3g(data, model, inference, state)
%% Sampler for population mcmc (mc3) with generalized temperature.
%
% aponteeduardo@gmail.com
% copyright (C) 2016
%

ollh = state.llh{1};
olpp = state.llh{2};

% Get the joint
%for i = 2:nn - 1
%    lpp = lpp + sum(state.llh{i}, 1);
%end

T = model.graph{1}.htheta.T;
[nv, nc] = size(T);
v = reshape([1:nv * nc], nv, nc);

% Doesn't need to do anything
if nc == 1
    return
end

c = 0;
for l = 1 : inference.mc3it
    s1 = ceil(rand() * nv);
    s2 = ceil(rand() * (nc - 1));
    % Realinate
    [s1, s2i] = ind2sub([nv, nc], v(s1, s2));
    [s1, s2j] = ind2sub([nv, nc], v(s1, s2 + 1));

    % Now we need to adjust the probability
    plpp1 = model.graph{2}.llh_sn(state.graph{2}, ...
        state.graph{3}, model.graph{2}.htheta, s1, s2, s2j);
    plpp2 = model.graph{2}.llh_sn(state.graph{2}, ...
        state.graph{3}, model.graph{2}.htheta, s1, s2 + 1, s2i);

    % Compute the ratio
    p = exp((ollh(s1, s2i) .* T(s1, s2 + 1) + ...
        ollh(s1, s2j) .* T(s1, s2) + ...
        plpp1 + plpp2) - ...
        (ollh(s1, s2i) .* T(s1, s2) + ...
        ollh(s1, s2j) .* T(s1, s2 + 1) + ...
        olpp(s1, s2) + olpp(s1, s2 + 1) ));
    if rand() < p
        c = c + 1;
        v(s1, [s2, s2 + 1]) = v(s1, [s2 + 1, s2]);
        olpp(s1, s2) = plpp1;
        olpp(s1, s2 + 1) = plpp2;
    end
end
%:fprintf(1, 'Swap rate %0.5f\n', c/inference.mc3it);
state.llh{1} = state.llh{1}(v);
state.llh{2} = olpp;
state.graph{2}.y = state.graph{2}.y(v);

end

author: translationalneuromodeling

repo: tapas


text:
%************************************************************************
%                       Bezier_circle                                   *
%                                                                       *
%           This program for practising CAO - EMMC11                    *
%           Professor :     Pierre BECKERS                              *
%           (Aerospace Laboratory, University of Liege, Belgium)        *
%           Student   :     NGUYEN QUOC DUAN                            *
%                          ( ngqduan@yahoo.com )                        *
%                                                                       *
%************************************************************************

%*** Problem ************************************************************
% Given 6 points (pi = 3.1416) :                                        *
%   Pi = [rcos(-pi/3+i*pi/3)    rsin(-pi/3+i*pi/3)] for  i = 0,1,...5   *
% Defined the Bezier curve starting in P0, ending in P5 and passing     *
% through points P1, P2, P3,P4. Explain the method and show all         *
% hypotheses.                                                           *
% ***********************************************************************

clear all; close all; clc;

% INPUT DATA

r=1;            % the length of the radius
pi=4*atan(1);   % the value of pi
P=[];           % coordinates matrix of 6 initial points
for i=1:6
    P(i,:)=[r*cos(-pi/3+(i-1)*pi/3) r*sin(-pi/3+(i-1)*pi/3)];
end

% COMPUTE THE COEFFICIENTS MATRICES

n=5;                % the polynomial degree of the Bezier curve
coefx=zeros(3,3);
coefy=zeros(3,3);
t=[0:0.2:1];        % matrix of the values of parameter at 6 given points
for k=1:3
    for j=1:3
        i=j-1;
        coefx(k,j)=C(i,n)*(B_para(i,n,t(k))-B_para(n-i,n,t(k)));
        coefy(k,j)=C(i,n)*(B_para(i,n,t(k))+B_para(n-i,n,t(k)));
    end
end

% COMPUTE THE COORDINATES OF BEZIER POINTS

% Initial values
xB=zeros(6,1);          % x-coordinates of 6 Bezier points
yB=zeros(6,1);          % y-coordinates of 6 Bezier points
% Solve the two equations systems
    xB(1:3)=coefx\P(1:3,1);  % x-coordinates of the first 3 Bezier points
    yB(1:3)=coefy\P(1:3,2);  % y-coordinates of the first 3 Bezier points
% Apply the symmetry of 6 given points
for i=4:6
    xB(i)=-xB(6-i+1);
    yB(i)=yB(6-i+1);
end

% OUTPUT THE SOLUTIONS
number=[0:5]';
disp('6 Bezier points coordinates : ' )
[number xB yB]

% PLOT FIGURE
figure(1);
h=gcf;
set(h,'name','Bezier control points & Bezier curve');
set(h,'NumberTitle','off');
axis equal;
title('Bezier control points & Bezier curve');

% Conecting 6 Bezier points by straight segments
plot(xB,yB,'r-')
hold on;

% Locate 4 Bezier control points
text(xB(2),yB(2),'Q1');
text(xB(3),yB(3),'Q2');
text(xB(4),yB(4),'Q3');
text(xB(5),yB(5),'Q4');
hold on;
% plot Bezier curve
precision=0.01;
for k=0:1/precision
    tp=k*precision;
    xp(k+1)=0;
    yp(k+1)=0;
    for j=1:6
        i=j-1;
        xpj=xB(j)*C(i,n)*B_para(i,n,tp);
        ypj=yB(j)*C(i,n)*B_para(i,n,tp);
        xp(k+1)=xpj+xp(k+1);
        yp(k+1)=ypj+yp(k+1);
    end
end
plot(xp,yp,'g-')
hold on;

% Locate 6 given points
text(P(1,1),P(1,2),'Po');
text(P(2,1),P(2,2),'P1');
text(P(3,1),P(3,2),'P2');
text(P(4,1),P(4,2),'P3');
text(P(5,1),P(5,2),'P4');
text(P(6,1),P(6,2),'P5');


text:
function [U,D] = unique_images(dir_name)
  % UNIQUE_IMAGES Generate a list of unique and duplicate images in a
  % directory.
  %
  % [U,D] = unique_images(dir_name)
  %
  % Inputs:
  %   dir_name  path to directory
  % Outputs:
  %   U list of unique images/image paths
  %   D list of duplicate iamges, so that [U,D] are all images
  %

  assert(isdir(dir_name));
  files = imdir(dir_name);
  U = {};
  ucount = 0;
  D = {};
  dcount = 0;
  valid = true(size(files));

  for ii = 1:numel(files)
    fii = [dir_name '/' files(ii).name];
    imii = imread(fii);
    for jj = (ii+1):numel(files)
      imjj = imread([dir_name '/' files(jj).name]);
      if all(size(imii) == size(imjj)) && all(imii(:) == imjj(:))
        valid(jj) = false;
      end
    end
  end
  U = files(valid);
  D = files(~valid);

end

author: alecjacobson

repo: gptoolbox


text:
% Copyright (C) 1993-2014, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB 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 Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

%%begin

% A quadrotor is a simple flying vehicle with 4 propellers with their thrust
% axes pointing upwards.  This demo shows an animation of a flying quadrotor
% with a nested control system.

sl_quadrotor

% Running the simulation shows the vehicle takes off and flying in a cirle.
sim('sl_quadrotor');

author: Allopart

repo: rbpf-gmapping


text:
function [confmatrix] = cfmatrix2(actual, predict, classlist, per, printout)
% CFMATRIX2 calculates the confusion matrix for any prediction
% algorithm ( prediction algorithm generates a list of classes to which
% each test feature vector is assigned );
%
% Outputs: confusion matrix
%
%                 Actual Classes
%                   p       n
%              ___|_____|______|
%    Predicted  p'|     |      |
%      Classes  n'|     |      |
%
%           Also the TP, FP, FN and TN are output for each class based on
%           http://en.wikipedia.org/wiki/Confusion_matrix
%           The Precision, Sensitivity and Specificity for each class have
%           also been added in this update along with the overall accuracy
%           of the model ( ModelAccuracy ).
%
%
% Further description of the outputs:
%
% True Postive [TP] = Condition Present + Positive result
% False Positive [FP] = Condition absent + Positive result [Type I error]
% False (invalid) Negative [FN] = Condition present + Negative result [Type II error]
% True (accurate) Negative [TN] = Condition absent + Negative result
% Precision(class) = TP(class) / ( TP(class) + FP(class) )
% Sensitivity(class) = Recall(class) = TruePositiveRate(class)
% = TP(class) / ( TP(class) + FN(class) )
% Specificity ( mostly used in 2 class problems )=
% TrueNegativeRate(class)
% = TN(class) / ( TN(class) + FP(class) )
%
% Inputs:
%
% 1. actual / 2. predict
% The inputs provided are the 'actual' classes vector
% and the 'predict'ed classes vector. The actual classes are the classes
% to which the input feature vectors belong. The predicted classes are the
% class to which the input feature vectors are predicted to belong to,
% based on a prediction algorithm.
% The length of actual class vector and the predicted class vector need to
% be the same. If they are not the same, an error message is displayed.
% 3. classlist
% The third input provides the list of all the classes {p,n,...} for which
% the classification is being done. All classes are numbers.
% 4. per = 1/0 (default = 0)
% This parameter when set to 1 provides the values in the confusion matrix
% as percentages. The default provides the values in numbers.
% 5. printout = 1/0 ( default = 1 )
% This parameter when set to 1 provides output on the matlab terminal and
% can be used to suppress output by setting to 0. ( default = 1 ). Assuming
% 'printout' of output use case would be more common and at the same time
% provided option to suppress output when the number of classes can be very
% large.
%
% Example:
% >> a = [ 1 2 3 1 2 3 1 1 2 3 2 1 1 2 3];
% >> b = [ 1 2 3 1 2 3 1 1 1 2 2 1 2 1 3];
% >> Cf = cfmatrix2(a, b, [1 2 3], 0, 1);
% is equivalent to
% >> Cf = cfmatrix2(a, b);
% The values of classlist(unique from actual), per(0), printout(1) are set
% to the respective defaults.
%
%
% [Avinash Uppuluri: avinash_uv@yahoo.com: Last modified: 03/28/2012]
%
% Changes added for 03/28/2012 upload
% a. Pre-initialize confmatrix
% b. Simplified logic making the code more readable and faster;
%    (based on comments from an interviewer who reviewed the code)
% c. Provide input variable 'printout' as an option to suppress output to
%    screen ( output to display is still the default (printout = 1)
%    assuming that will be the more common use case ).
% d. Added Precision(class), Sensitivity(class), Specificity(class) and
%    the overall accuracy of model calculations.

% If classlist not entered: make classlist equal to all
% unique elements of actual
if (nargin < 2)
   error('Not enough input arguments. Need atleast two vectors as input');
elseif (nargin == 2)
    classlist = unique(actual); % default values from actual
    per = 0;
    printout = 1;
elseif (nargin == 3)
    per = 0; % default is numbers and input 1 or higher for percentage
    printout = 1;
elseif (nargin == 4)
    printout = 1; % default is silent output ( 0 ); one or higher printsout
elseif (nargin > 5)
   error('Too many input arguments.');
end

if (length(actual) ~= length(predict))
    error('First two inputs need to be vectors with equal size.');
elseif ((size(actual,1) ~= 1) && (size(actual,2) ~= 1))
    error('First input needs to be a vector and not a matrix');
elseif ((size(predict,1) ~= 1) && (size(predict,2) ~= 1))
    error('Second input needs to be a vector and not a matrix');
end
format short g;
n_class = length(classlist);
confmatrix = zeros(n_class);
line_two = '----------';
line_three = '_________|';

for i = 1:n_class
    for j = 1:n_class
        m = (predict == classlist(i) ...
           & actual  == classlist(j));
        confmatrix(i,j) = sum(m);
    end
    line_two = strcat(line_two,'---',num2str(classlist(i)),'-----');
    line_three = strcat(line_three,'__________');
end

% True Postive [TP] = Condition Present + Positive result
% False Positive [FP] = Condition absent + Positive result [Type I error]
% False (invalid) Negative [FN] = Condition present + Negative result [Type II error]
% True (accurate) Negative [TN] = Condition absent + Negative result
% Precision(class) = TP(class) / ( TP(class) + FP(class) )
% Sensitivity(class) = Recall(class) = TruePositiveRate(class)
% = TP(class) / ( TP(class) + FN(class) )
% Specificity ( mostly used in 2 class problems )=
% TrueNegativeRate(class)
% = TN(class) / ( TN(class) + FP(class) )

TPFPFNTN    = zeros(4, n_class);
Precision   = zeros(1, n_class);
Sensitivity = zeros(1, n_class);
Specificity = zeros(1, n_class);

temps1 = sprintf('    TP  ');
temps2 = sprintf('    FP  ');
temps3 = sprintf('    FN  ');
temps4 = sprintf('    TN  ');
temps5 = sprintf('Preci.  ');
temps6 = sprintf('Sensi.  ');
temps7 = sprintf('Speci.  ');

for i = 1:n_class
    TPFPFNTN(1, i) = confmatrix(i,i); % TP
    temps1 = strcat(temps1,sprintf(' |   %2.2f    ',TPFPFNTN(1, i)));
    TPFPFNTN(2, i) = sum(confmatrix(i,:))-confmatrix(i,i); % FP
    temps2 = strcat(temps2,sprintf(' |   %2.2f    ',TPFPFNTN(2, i) ));
    TPFPFNTN(3, i) = sum(confmatrix(:,i))-confmatrix(i,i); % FN
    temps3 = strcat(temps3,sprintf(' |   %2.2f    ',TPFPFNTN(3, i) ));
    TPFPFNTN(4, i) = sum(confmatrix(:)) - sum(confmatrix(i,:)) -...
        sum(confmatrix(:,i)) + confmatrix(i,i); % TN
    temps4 = strcat(temps4,sprintf(' |   %2.2f    ',TPFPFNTN(4, i) ));
    % Precision(class) = TP(class) / ( TP(class) + FP(class) )
    Precision(i)   = TPFPFNTN(1, i) / sum(confmatrix(i,:));
    temps5 = strcat(temps5,sprintf(' |   %1.2f    ',Precision(i) ));
    % Sensitivity(class) = Recall(class) = TruePositiveRate(class)
    % = TP(class) / ( TP(class) + FN(class) )
    Sensitivity(i) = TPFPFNTN(1, i) / sum(confmatrix(:,i));
    temps6 = strcat(temps6,sprintf(' |   %1.2f    ',Sensitivity(i) ));
    % Specificity ( mostly used in 2 class problems )=
    % TrueNegativeRate(class)
    % = TN(class) / ( TN(class) + FP(class) )
    Specificity(i) = TPFPFNTN(4, i) / ( TPFPFNTN(4, i) + TPFPFNTN(2, i) );
    temps7 = strcat(temps7,sprintf(' |   %1.2f    ',Specificity(i) ));
end

ModelAccuracy = sum(diag(confmatrix))/sum(confmatrix(:));
temps8 = sprintf('Model Accuracy is %1.2f ',ModelAccuracy);

if (per > 0) % ( if > 0 implies true; < 0 implies false )
    confmatrix = (confmatrix ./ length(actual)).*100;
end

if ( printout > 0 ) % ( if > 0 printout; < 0 no printout )
    disp('------------------------------------------');
    disp('             Actual Classes');
    disp(line_two);
    disp('Predicted|                     ');
    disp('  Classes|                     ');
    disp(line_three);

    for i = 1:n_class
        temps = sprintf('       %d             ',i);
        for j = 1:n_class
            temps = strcat(temps,sprintf(' |    %2.1f    ',confmatrix(i,j)));
        end
        disp(temps);
        clear temps
    end
    disp('------------------------------------------');

    disp('------------------------------------------');
    disp('             Actual Classes');
    disp(line_two);
    disp(temps1); disp(temps2); disp(temps3); disp(temps4);
    disp(temps5); disp(temps6); disp(temps7);
    disp('------------------------------------------');
    disp(temps8);
    disp('------------------------------------------');
end
clear temps1 temps2 temps3 temps4 temps5 temps6 temps7 temps8

author: Sable