text:
function varargout = process_generate_canonical( varargin )
% PROCESS_GENERATE_CANONICAL: Generate SPM canonical surface.
%
% USAGE:     OutputFiles = process_generate_canonical('Run',     sProcess, sInputs)
%         [isOk, errMsg] = process_generate_canonical('Compute', iSubject, iAnatomy=[default], Resolution)
%                          process_generate_canonical('ComputeInteractive', iSubject, iAnatomy=[default], Resolution)

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

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % Description the process
    sProcess.Comment     = 'Generate SPM canonical surfaces';
    sProcess.Category    = 'Custom';
    sProcess.SubGroup    = {'Import', 'Import anatomy'};
    sProcess.Index       = 12;
    sProcess.Description = 'https://github.com/neurodebian/spm12/blob/master/spm_eeg_inv_mesh.m';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'import'};
    sProcess.OutputTypes = {'import'};
    sProcess.nInputs     = 1;
    sProcess.nMinFiles   = 0;
    sProcess.isSeparator = 1;
    % Option: Subject name
    sProcess.options.subjectname.Comment = 'Subject name:';
    sProcess.options.subjectname.Type    = 'subjectname';
    sProcess.options.subjectname.Value   = '';
    % Option: Cortex resolution
    sProcess.options.resolution.Comment = {'5124', '8196', '20484', 'Cortex surface resolution:'};
    sProcess.options.resolution.Type    = 'radio_line';
    sProcess.options.resolution.Value   = 3;
end

%% ===== FORMAT COMMENT =====
function Comment = FormatComment(sProcess) %#ok<DEFNU>
    Comment = sProcess.Comment;
end

%% ===== RUN =====
function OutputFiles = Run(sProcess, sInputs) %#ok<DEFNU>
    OutputFiles = {};
    % Get resolution
    Resolution = sProcess.options.resolution.Value;
    % Get subject name
    SubjectName = file_standardize(sProcess.options.subjectname.Value);
    if isempty(SubjectName)
        bst_report('Error', sProcess, [], 'Subject name is empty.');
        return;
    end
    % Get subject
    [sSubject, iSubject] = bst_get('Subject', SubjectName);
    if isempty(iSubject)
        bst_report('Error', sProcess, [], ['Subject "' SubjectName '" does not exist.']);
        return
    end
    % Call processing function
    [isOk, errMsg] = Compute(iSubject, [], Resolution, 0);
    % Handling errors
    if ~isOk
        bst_report('Error', sProcess, [], errMsg);
    elseif ~isempty(errMsg)
        bst_report('Warning', sProcess, [], errMsg);
    end
    % Return an empty structure
    OutputFiles = {'import'};
end

%% ===== COMPUTE CANONICAL SURFACES =====
function [isOk, errMsg] = Compute(iSubject, iAnatomy, Resolution, isInteractive)
    isOk = 0;
    % Initialize SPM
    [isInstalled, errMsg] = bst_plugin('Install', 'spm12', isInteractive);
    if ~isInstalled
        return;
    end
    bst_plugin('SetProgressLogo', 'spm12');

    % ===== GET SUBJECT =====
    % Get subject
    [sSubject, iSubject] = bst_get('Subject', iSubject);
    if isempty(sSubject)
        errMsg = 'Subject does not exist.';
        return
    end
    % Check if a MRI is available for the subject
    if isempty(sSubject.Anatomy)
        errMsg = ['No MRI available for subject "' SubjectName '".'];
        return
    end
    % Get default MRI if not specified
    if isempty(iAnatomy)
        iAnatomy = sSubject.iAnatomy;
    end

    % ===== VERIFY FIDUCIALS IN MRI =====
    % Load MRI file
    MriFileBst = sSubject.Anatomy(iAnatomy).FileName;
    sMri = in_mri_bst(MriFileBst);
    % If the SCS transformation is not defined: compute MNI transformation to get a default one
    if isempty(sMri) || ~isfield(sMri, 'SCS') || ~isfield(sMri.SCS, 'NAS') || ~isfield(sMri.SCS, 'LPA') || ~isfield(sMri.SCS, 'RPA') || (length(sMri.SCS.NAS)~=3) || (length(sMri.SCS.LPA)~=3) || (length(sMri.SCS.RPA)~=3) || ~isfield(sMri.SCS, 'R') || isempty(sMri.SCS.R) || ~isfield(sMri.SCS, 'T') || isempty(sMri.SCS.T)
        % Issue warning
        errMsg = 'Missing NAS/LPA/RPA: Computing the MNI transformation to get default positions.';
        % Compute MNI normalization
        [sMri, errNorm] = bst_normalize_mni(MriFileBst);
        % Handle errors
        if ~isempty(errNorm)
            errMsg = [errMsg 10 'Error trying to compute the MNI transformation: ' 10 errNorm 10 ...
                'The surfaces will not be properly aligned with the MRI.'];
        end
    end

    % ===== CALL SPM FUNCTIONS =====
    % Empty temporary folder, otherwise it reuses previous files in the folder
    TmpDir = bst_get('BrainstormTmpDir', 0, 'spmcanonical');
    % Save MRI in .nii format
    NiiFile = bst_fullfile(TmpDir, 'spm_canonical.nii');
    out_mri_nii(sMri, NiiFile);
    % Call SPM function
    switch (Resolution)
        case {1,2,3}
            spmMesh = spm_eeg_inv_mesh(NiiFile, Resolution);
        case 4
            spmMesh = ImaGIN_spm_eeg_inv_mesh(NiiFile, Resolution);
        otherwise
            error('Unsupported resolution.');
    end
    % ===== READ OUTPUT SURFACES =====
    % Read transformation from temporary .nii
    niiMri = in_mri_nii(NiiFile, 0, 0, 0);
    % Create surfaces
    sHead   = CreateSurface(sMri, niiMri, export(gifti(spmMesh.tess_scalp),'patch'),  'spm_head');
    sOuter  = CreateSurface(sMri, niiMri, export(gifti(spmMesh.tess_oskull),'patch'), 'spm_outerskull');
    sInner  = CreateSurface(sMri, niiMri, export(gifti(spmMesh.tess_iskull),'patch'), 'spm_innerskull');
    sCortex = CreateSurface(sMri, niiMri, export(gifti(spmMesh.tess_ctx),'patch'),    'spm_cortex');

    % ===== SAVE NEW SURFACES =====
    % Create output filenames
    SurfaceDir    = bst_fileparts(file_fullpath(MriFileBst));
    SpmHeadFile   = file_unique(bst_fullfile(SurfaceDir, sprintf('tess_head_spm_%dV.mat', length(sHead.Vertices))));
    SpmOuterFile  = file_unique(bst_fullfile(SurfaceDir, sprintf('tess_outerskull_spm_%dV.mat', length(sOuter.Vertices))));
    SpmInnerFile  = file_unique(bst_fullfile(SurfaceDir, sprintf('tess_innerskull_spm_%dV.mat', length(sInner.Vertices))));
    SpmCortexFile = file_unique(bst_fullfile(SurfaceDir, sprintf('tess_cortex_spm_%dV.mat', length(sCortex.Vertices))));
    % Save head
    bst_save(SpmHeadFile, sHead, 'v7');
    db_add_surface(iSubject, SpmHeadFile, sHead.Comment);
    % Save outerskull
    bst_save(SpmOuterFile, sOuter, 'v7');
    db_add_surface(iSubject, SpmOuterFile, sOuter.Comment);
    % Save innerskull
    bst_save(SpmInnerFile, sInner, 'v7');
    db_add_surface(iSubject, SpmInnerFile, sInner.Comment);
    % Save cortex
    bst_save(SpmCortexFile, sCortex, 'v7');
    db_add_surface(iSubject, SpmCortexFile, sCortex.Comment);

    % Delete the temporary files
    file_delete(TmpDir, 1, 1);
    % Remove logo
    bst_plugin('SetProgressLogo', []);
    isOk = 1;
end

%% ===== COMPUTE/INTERACTIVE =====
function ComputeInteractive(iSubject, iAnatomy, Resolution) %#ok<DEFNU>
    % Get inputs
    if (nargin < 2) || isempty(iAnatomy)
        iAnatomy = [];
    end
    if (nargin < 3) || isempty(Resolution)
        strRes = java_dialog('question', 'Cortex surface resolution:', 'SPM canonical mesh', [], {'5124', '8196', '20484'}, '20484');
        if isempty(strRes)
            return;
        end
        switch (strRes)
            case '5124',  Resolution = 1;
            case '8196',  Resolution = 2;
            case '20484', Resolution = 3;
        end
    end
    % Open progress bar
    bst_progress('start', 'SPM', 'Generating canonical surfaces...');
    % Compute surfaces
    [isOk, errMsg] = Compute(iSubject, iAnatomy, Resolution, 1);
    % Error handling
    if ~isOk
        bst_error(errMsg, 'SPM canonical surfaces', 0);
    elseif ~isempty(errMsg)
        java_dialog('msgbox', ['Warning: ' errMsg]);
    end
    % Close progress bar
    bst_progress('stop');
end

%% ===== CREATE SURFACE =====
function sTess = CreateSurface(sMri, niiMri, gii, Comment)
    % Transformation to apply to the surface
    ras2vox = inv(niiMri.Header.nifti.vox2ras);
    % Create surfaces structure
    sTess = db_template('SurfaceMat');
    sTess.Comment  = sprintf('%s_%dV', Comment, length(gii.vertices));
    sTess.Vertices = bst_bsxfun(@plus, ras2vox(1:3,1:3)*gii.vertices', ras2vox(1:3,4))';
    sTess.Vertices = bst_bsxfun(@plus, sTess.Vertices, [1 1 1]);
    sTess.Vertices = cs_convert(sMri, 'voxel', 'scs', sTess.Vertices);
    sTess.Faces    = gii.faces(:,[2 1 3]);
    sTess = bst_history('add', sTess, 'spm', ['Canonical surface generated with: ' spm('version')]);
end

author: brainstorm-tools


text:
% Copyright 2014 Jason Heeris, jason.heeris@gmail.com
%
% This file is part of the gammatone toolkit, and is licensed under the 3-clause
% BSD license: https://github.com/detly/gammatone/blob/master/COPYING
function test_gammatonegram()
    % Need:
    %  wave
    %  fs
    %  window_time
    %  hop_time
    %  channels
    %  f_min
    %  f_max

    % Need to mock out:
    %  make_erb_filters output (elide)
    %  centre_freqs (elide)
    %  erb_filterbank (depends on X, SR, N, FMIN)

    % Ensure reproducible tests
    rand('state', [3 1 4 1 5 9 2 7]);

    gammatonegram_inputs = {
        'sawtooth_01', sawtooth(2*pi*10100*[0:22050 - 1]'/22050, 0.5), 22050, 0.025, 0.010, 64, 50; ...
        'sin220_01'  , sin(2*pi*220*[0:4800 - 1]'/48000), 48000, 0.01, 0.01, 64, 50; ...
        'sin220_02'  , sin(2*pi*220*[0:4800 - 1]'/48000), 48000, 0.025, 0.01, 32, 50; ...
        'rand_01'    , rand([1, 4410 - 1]), 44100, 0.02, 0.015, 128, 500; ...
        'rand_02'    , rand([1, 9600 - 1]), 96000, 0.01, 0.005, 256, 20; ...
        'rand_03'    , rand([1, 4800 - 1]), 48000, 0.01, 0.010, 256, 20; ...
    };

    % Mocked intermediate results for unit testing
    gammatonegram_mocks = {};

    % Actual results
    gammatonegram_results = {};

    for tnum=1:size(gammatonegram_inputs)(1)
        [name, wave, fs, twin, thop, chs, fmin] = deal(gammatonegram_inputs{tnum,:});
        res = gammatonegram( ...
                  wave, ...
                  fs, ...
                  twin, ...
                  thop, ...
                  chs, ...
                  fmin, ...
                  0, % fmax is ignored
                  0 % Don't use FFT method
              );

        % This is for mocking the output of the equivalent Python functions
        nwin     = round(twin * fs);
        hopsamps = round(thop * fs);
        f_coefs  = flipud(MakeERBFilters(fs, chs, fmin));
        x_f      = ERBFilterBank(wave, f_coefs);
        x_e      = [x_f .^ 2];
        x_e_cols = size(x_e, 2);
        ncols    = 1 + floor((x_e_cols - nwin) / hopsamps);

        % Mock out the ERB filter functions too
        fcoefs = flipud(MakeERBFilters(fs, chs, fmin));
        erb_fb_output = ERBFilterBank(wave, fcoefs);

        gammatonegram_mocks(tnum, :) = { ...
            erb_fb_output, ...
            x_e_cols ...
        };

        gammatonegram_results(tnum, :) = { ...
            res, ...
            nwin, ...
            hopsamps, ...
            ncols ...
        };

    end;

    results_file = fullfile('..', 'tests', 'data', 'test_gammatonegram_data.mat');
    save(results_file, 'gammatonegram_inputs', 'gammatonegram_mocks', 'gammatonegram_results');
end;

author: detly

repo: gammatone

save_path: github-repos/MATLAB/detly-gammatone


text:
function R=ImportSurferGrd(filename)
% ImportSurferGrid reads an ASCII format surfer .GRD file and returns the
% data matrix in R. NaN values, which are written as 1.7308+38 in surfer
% are replaced with NaN to prevent crashing.
% Input:
%      filename = name of the file to be read, including ".grd" extension
% Output:
% R - matrix containing data (e.g elevation)
%James Ramm 2011
if exist(filename,'file')
grdfile=fopen(filename,'r');    % Open file
else
    error('ImportGrid:NoFile','File does not exist');
end

code=fgetl(grdfile);            % Reads surfer code 'DSAA'
% Grid dimensions (number of nodes, 2nd line of surfer file)
aux=str2num(fgetl(grdfile)); nx=aux(1); % columns
ny=aux(2); % rows
% X limits, 3rd line of .grd file
aux=(fgetl(grdfile)); %xmin=aux(1); xmax=aux(2);
% Y limits, 4th line of .grd file
aux=(fgetl(grdfile)); %ymin=aux(1); ymax=aux(2);
% Z limits, 5th line
aux=(fgetl(grdfile)); %zmin=aux(1); aux(2);
% Reads Z matrix
[R] = fscanf(grdfile, '%f', [nx,ny]);

% transpose matrix as matlab indices are (rows, columns) not (columns, rows).
R = R';
fclose(grdfile);
% Get rid of those crappy surfer nan values and replace with NaN
[I,J]= size(R);

for i = 1:I
    for j=1:J
        if R(i,j)>100000
            R(i,j) = NaN;
        end
    end
end
end

author: Sable

repo: mcbench-benchmarks


text:
function finalStruct = parseKspConfigFile(fileTxt)
    pos1 = strfind(fileTxt,'{');

    pairs = NaN(length(pos1),2);
    for(i=1:length(pos1))
        pos2 = getbracketpairpos(fileTxt,pos1(i));
        pairs(i,:) = [pos1(i), pos2];
    end

    structs = {};
    structsStartEnd = [];
    for(i=1:size(pairs,1))
        pair = pairs(i,:);

        internalPairs = pairs(pairs(:,1) > pair(1) & pairs(:,2) < pair(2),:);

        if(not(isempty(internalPairs)))
            endInd = min(internalPairs(:,1));
        else
            endInd = pair(2)-1;
        end
        startInd = pair(1)+1;
        subText = fileTxt(startInd:endInd);

        priorText = fileTxt(1:pair(1)-1);
        tokens = regexp(priorText,'(\S)+\s*$','tokens');
        if(not(isempty(tokens)))
            className = matlab.lang.makeValidName(tokens{1}{1});
        else
            className = '';
        end

        subTextLines = strtrim(strsplit(subText,'\n'));
        s = struct('ClassName',className);
        for(j=1:length(subTextLines))
            line = subTextLines{j};
            tokens = regexp(line,'(\w+)\s*\=\s*([ \S]+)','tokens');

            if(not(isempty(tokens)))
                token = tokens{1};
                name = matlab.lang.makeValidName(token{1});

                dataStr = attemptStrEval(token{2});
                data = str2double(dataStr);
                if(isnan(data))
                    data = dataStr;
                end

                if(isfield(s,name))
                    if(not(iscell(s.(name))))
                        s.(name) = {s.(name)};
                    end
                    s.(name){end+1} = data;
                else
                    s.(name) = data;
                end
            end
        end

        if(not(isempty(fields(s))))
            structs{end+1} = s; %#ok<AGROW>
            structsStartEnd(end+1,:) = pair; %#ok<AGROW>
        end
    end

    width = structsStartEnd(:,2) - structsStartEnd(:,1);
    [~,I] = sort(width);
    structsStartEnd = structsStartEnd(I,:);
    structs = structs(I);

    for(i=1:length(structsStartEnd))
        structStartEnd = structsStartEnd(i,:);

        greaterStructs = structsStartEnd(structsStartEnd(:,1) < structStartEnd(1) & ...
                                         structsStartEnd(:,2) > structStartEnd(2), :);

        if(not(isempty(greaterStructs)))
            [~,I] = min(greaterStructs(:,2) - greaterStructs(:,1));
            greaterStruct = greaterStructs(I,:);

            [~, index] = ismember(greaterStruct, structsStartEnd, 'rows');

            if(not(isfield(structs{index},structs{i}.ClassName)))
                structs{index}.(structs{i}.ClassName) = structs{i};
            else
                if(not(iscell(structs{index}.(structs{i}.ClassName))))
                    structs{index}.(structs{i}.ClassName) = {structs{index}.(structs{i}.ClassName)};
                end

                structs{index}.(structs{i}.ClassName){end+1} = structs{i};
            end
        end
    end

    finalStruct = structs{end};
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function ind = gpDataIndices(model, dimNo, blockNo)

% GPDATAINDICES Return indices of present data.
% FORMAT
% DESC returns the indices of data which is not missing for a given
% dimension in the GP-LVM.
% ARG model : the model for which the indices are being returned.
% ARG dimNo : the dimension for which the presence of missing data
% is being looked at.
% RETURN ind : indices of training data along that dimension which
% isn't missing.
%
% DESC returns the indices of data which is not missing for a given
% dimension in the GP-LVM and a block number in the PITC approximation.
% ARG model : the model for which the indices are being returned.
% ARG dimNo : the dimension for which the presence of missing data
% is being looked at.
% ARG blockNo : the block number in the PITC approximation for
% which the indices are required.
% RETURN ind : indices of training data along that dimension which
% isn't missing.
%
% SEEALSO : gpCreate
%
% COPYRIGHT : Neil D. Lawrence, 2006

% GP

if nargin > 2
  if ~strcmp(model.approx, 'pitc')
    error('Block number only relevant for pitc approximation.');
  else
    if blockNo == 1
      startVal = 1;
    else
      startVal = model.blockEnd(blockNo-1)+1;
    end
    endVal = model.blockEnd(blockNo);
    if model.isMissingData
      st = min(find(model.indexPresent{dimNo}>=startVal));
      fi = max(find(model.indexPresent{dimNo}<=endVal));
      ind = model.indexPresent{dimNo}(st:fi)';
    else
      ind = startVal:endVal;
    end
  end
else
  if strcmp(model.approx, 'pitc')
    error('Must give block number with PITC approximation');
  else
    if model.isMissingData
      ind = model.indexPresent{dimNo}';
    else
      ind = 1:model.N;
    end
  end
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
% ---------------------------------------------------------
% Copyright (c) 2017, Gurkirt Singh
% This code and is available
% under the terms of MIT License provided in LICENSE.
% Please retain this notice and LICENSE if you use
% this file (or any portion of it) in your project.
% ---------------------------------------------------------
%% This is main script to build tubes and evaluate them %%

function myI01onlineTubes()
close all
data_root = '/home/zhujiagang/realtime-action-detection';
save_root = '/home/zhujiagang/realtime-action-detection/save';
iteration_num_rgb = [120000]; % you can also evaluate on multiple iertations
iteration_num_flow = [120000]; % you can also evaluate on multiple iertations

% add subfolder to matlab paths
addpath(genpath('gentube/'));
addpath(genpath('actionpath/'));
addpath(genpath('eval/'));
addpath(genpath('utils/'));
model_type = 'CONV';

completeList = {...
    {'ucf24','01',{'rgb'},iteration_num_rgb,{'score'}},...
    {'ucf24','01',{'brox'},iteration_num_flow,{'score'}}...
    {'ucf24','01',{'fastOF'},iteration_num_flow,{'score'}}...
    };

alldopts = cell(2,1);
count = 1;
gap=3;

for setind = 1:length(completeList)
    [dataset, listid, imtypes, iteration_nums, costTypes] = enumurateList(completeList{setind});
    for ct = 1:length(costTypes)
        costtype = costTypes{ct};
        for imtind = 1:length(imtypes)
            imgType = imtypes{imtind};
            for iteration = iteration_nums
                for iouthresh=0.1
                    %% generate directory sturcture based on the options
                    dopts = initDatasetOpts(data_root,save_root,dataset,imgType,model_type,listid,iteration,iouthresh,costtype, gap);
                    if exist(dopts.detDir,'dir')
                        alldopts{count} = dopts;
                        count = count+1;
                    end
                end
            end
        end
    end
end

%% For each option type build tubes and evaluate them
for index = 1:count-1
    opts = alldopts{index};
    if exist(opts.detDir,'dir')
        fprintf('Video List %02d :: %s\nAnnotFile :: %s\nImage  Dir :: %s\nDetection Dir:: %s\nActionpath Dir:: %s\nTube Dir:: %s\n',...
            index, opts.vidList, opts.annotFile, opts.imgDir, opts.detDir, opts.actPathDir, opts.tubeDir);
        %% online bbx and prediction scores display given frame level detections
        actionPaths(opts);
    end
end

%% Function to enumrate options
function [dataset,listnum,imtypes,weights,costTypes] = enumurateList(sublist)
dataset = sublist{1}; listnum = sublist{2}; imtypes = sublist{3};
weights = sublist{4};costTypes = sublist{5};

%% Facade function for smoothing tubes and evaluating them
function results = gettubes(dopts)

numActions = length(dopts.actions);
results = zeros(300,6);
counter=1;
class_aps = cell(2,1);

annot = load(dopts.annotFile);
annot = annot.annot;
testvideos = getVideoNames(dopts.vidList);

for alpha = 3
    fprintf('alpha %03d ', alpha);
    % read action paths
    actionpaths = readALLactionPaths(dopts.vidList,dopts.actPathDir,1);
    %% perform temporal trimming
    smoothedtubes = parActionPathSmoother(actionpaths,alpha*ones(numActions,1),numActions);
    fprintf('\n');
end

results(counter:end,:) = [];
result = cell(2,1);
result{2} = class_aps;
result{1} = results;
results = result;

function videos = getVideoNames(split_file)
% -------------------------------------------------------------------------
fid = fopen(split_file,'r');
data = textscan(fid, '%s');
videos  = cell(1);
count = 0;
for i=1:length(data{1})
    filename = cell2mat(data{1}(i,1));
    count = count +1;
    videos{count} = filename;
end

author: gurkirt

repo: realtime-action-detection


text:
function [p,L] = tspsearch(X,m)
%TSPSEARCH Heuristic method for Traveling Salesman Problem (TSP).
%   [P,L] = TSPSEARCH(X,M) gives a tour P of length L. X is either a
%   coordinate matrix of size Nx2 or Nx3 or a symmetric distance matrix.
%   Euclidian distances are used in the coordinate case. M is an integer
%   in the range 1 to N. Default is M = 1.
%
%   METHOD
%   M nearest neighbour tours are generated from randomly selected starting
%   points. Each tour is improved by 2-opt heuristics (pairwise exchange of
%   edges) and the best result is selected.
%
%   EXAMPLES
%
%   X = rand(100,2);
%   [p,L] = tspsearch(X,100);
%   tspplot(p,X)
%
%   % Optimal tour length 1620
%   X = load('hex162.dat');
%   [p,L] = tspsearch(X,10);
%   tspplot(p,X)
%
%   % Optimal tour length 4860
%   X = load('hex486.dat');
%   [p,L] = tspsearch(X);
%   tspplot(p,X)

%   Author: Jonas Lundgren <splinefit@gmail.com> 2012

% Check first argument
[n,dim] = size(X);
if dim == 2 || dim == 3
    % X is a coordinate matrix, compute euclidian distances
    D = distmat(X);
elseif n == dim && min(X(:)) >= 0 && isequal(X,X')
    % X is a distance matrix
    D = X;
else
    mess = 'First argument must be Nx2, Nx3 or symmetric and nonnegative.';
    error('tspsearch:first',mess)
end

% Check second argument
if nargin < 2 || isempty(m)
    m = 1;
elseif ~isscalar(m) || m < 1 || m > n || fix(m) < m
    mess = 'M must be an integer in the range 1 to %d.';
    error('tspsearch:second',mess,n)
end

% Starting points for nearest neighbour tours
s = randperm(n);

Lmin = inf;
for k = 1:m
    % Nearest neighbour tour
	p = greedy(s(k),D);
    % Improve tour by 2-opt heuristics
	[p,L] = exchange2(p,D);
    % Keep best tour
	if L < Lmin
        Lmin = L;
        pmin = p;
	end
end

% Output
p = double(pmin);
L = Lmin;

%--------------------------------------------------------------------------
function D = distmat(X)
%DISTMAT Compute euclidian distance matrix from coordinates

[n,dim] = size(X);
D = zeros(n);
for j = 1:n
    for k = 1:dim
        v = X(:,k) - X(j,k);
        D(:,j) = D(:,j) + v.*v;
    end
end
D = sqrt(D);

%--------------------------------------------------------------------------
function p = greedy(s,D)
%GREEDY Travel to nearest neighbour, starting with node s.

n = size(D,1);
p = zeros(1,n,'uint16');
p(1) = s;

for k = 2:n
    D(s,:) = inf;
    [junk,s] = min(D(:,s)); %#ok
    p(k) = s;
end

%--------------------------------------------------------------------------
function [p,L] = exchange2(p,D)
%EXCHANGE2 Improve tour p by 2-opt heuristics (pairwise exchange of edges).
%   The basic operation is to exchange the edge pair (ab,cd) with the pair
%   (ac,bd). The algoritm examines all possible edge pairs in the tour and
%   applies the best exchange. This procedure continues as long as the
%   tour length decreases. The resulting tour is called 2-optimal.

n = numel(p);
zmin = -1;

% Iterate until the tour is 2-optimal
while zmin < 0

    zmin = 0;
    i = 0;
    b = p(n);

    % Loop over all edge pairs (ab,cd)
    while i < n-2
        a = b;
        i = i+1;
        b = p(i);
        Dab = D(a,b);
        j = i+1;
        d = p(j);
        while j < n
            c = d;
            j = j+1;
            d = p(j);
            % Tour length diff z
            % Note: a == d will occur and give z = 0
            z = (D(a,c) - D(c,d)) + D(b,d) - Dab;
            % Keep best exchange
            if z < zmin
                zmin = z;
                imin = i;
                jmin = j;
            end
        end
    end

    % Apply exchange
    if zmin < 0
        p(imin:jmin-1) = p(jmin-1:-1:imin);
    end

end

% Tour length
q = double(p);
ind = sub2ind([n,n],q,[q(2:n),q(1)]);
L = sum(D(ind));


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

check3Color(img);
[r, c, col] = size(img);
imgOut = zeros(r, c, col);

if(inverse == 0)%forward transform
    norm = zeros(r, c);
    for i=1:3
        norm = norm + img(:,:,i);
    end

    imgOut(:,:,1) = img(:,:,2);

    imgOut(:,:,2) = img(:,:,1) ./ (norm);
    imgOut(:,:,3) = img(:,:,2) ./ (norm);
end

if(inverse == 1)%inverse transform
    Y_over_y = img(:,:,1) ./ img(:,:,3);
    imgOut(:,:,1) = Y_over_y .* img(:,:,2);
    imgOut(:,:,2) = img(:,:,1);
    imgOut(:,:,3) = Y_over_y .* (1.0 - img(:,:,2) - img(:,:,3));
end

imgOut = RemoveSpecials(imgOut);

end

author: banterle


text:
function [Vm]=rigidTransformVertices(V,T,v)

Vm=V;
[I,J,K]=cart2im(V(:,1),V(:,2),V(:,3),v); %Convert to image coordinates
IJK=[I(:) J(:) K(:) ones(size(I(:)))]; %Prepare for mapping
IJK_mapped=(T*IJK')'; %Do mapping
[Vm(:,1),Vm(:,2),Vm(:,3)]=im2cart(IJK_mapped(:,1),IJK_mapped(:,2),IJK_mapped(:,3),v); %Convert mapped image coordinates back to "Cartesian" coordinates

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% 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/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
classdef PlottingVademecumSuperEllipseVsStress < handle

    properties (Access = private)
        phiV
        qOpt
        qMin
        qMax
        qMean
        xi
        rho
        compressedFileName
        vademecum
    end

    properties (Access = private)

    end

    methods (Access = public)

        function obj = PlottingVademecumSuperEllipseVsStress()
            obj.init();
            obj.loadCompressedVademecum();
            obj.compute();
        end

    end

    methods (Access = private)

        function init(obj)
            obj.compressedFileName = 'OptimalSuperEllipses';
        end

        function loadCompressedVademecum(obj)
            d = load(obj.compressedFileName);
            obj.vademecum = d;
        end

        function compute(obj)
            nMx  = size(obj.vademecum.q,1);
            nMy  = size(obj.vademecum.q,2);
            for iMx = 1:nMx
                for iMy = 1:nMy
                    obj.computeDataForPlotter(iMx,iMy)
                    obj.plot(iMx,iMy);
                end
            end
        end

        function computeDataForPlotter(obj,iMx,iMy)
            obj.computeXi(iMx,iMy);
            obj.computeRho(iMx,iMy);
            obj.computePhi(iMx,iMy);
            obj.computeQopt(iMx,iMy);
            obj.computeQminQmax();
            obj.computeQmean();
        end

        function computeXi(obj,iMx,iMy)
            v = obj.vademecum;
            obj.xi = v.xi(iMx,iMy,1);
        end

        function computeRho(obj,iMx,iMy)
            v = obj.vademecum;
            obj.rho = v.rho(iMx,iMy,1);
        end

        function computePhi(obj,iMx,iMy)
            v = obj.vademecum;
            obj.phiV = squeeze(v.phi(iMx,iMy,:));
        end

        function computeQopt(obj,iMx,iMy)
            v = obj.vademecum;
            obj.qOpt = squeeze(v.q(iMx,iMy,:));
        end

        function computeQminQmax(obj)
            s.mxMax = 0.99;
            s.myMax = 0.99;
            s.mxMin = 0.01;
            s.myMin = 0.01;
            s.xi =  obj.xi;
            s.rho = obj.rho;
            sC = SuperellipseExponentBoundsComputer(s);
            obj.qMin = sC.qMin;
            obj.qMax = sC.qMax;
        end

        function computeQmean(obj)
            s.phiV = obj.phiV;
            qMeanC    = SuperEllipseMeanExponentComputer(s);
            obj.qMean = qMeanC.compute(obj.qOpt,obj.xi);
        end

        function plot(obj,iMx,iMy)
            nMx  = size(obj.vademecum.q,1);
            iter = (iMx-1)*nMx + iMy;
            s.phiV   = obj.phiV;
            s.qOpt   = obj.qOpt;
            s.qMin   = obj.qMin;
            s.qMax   = obj.qMax;
            s.qMean  = obj.qMean;
            s.rho    = obj.rho;
            s.xi     = obj.xi;
            p = OptimalSuperEllipseExponentVsGaussianPlotter(s);
            p.plot(iter);
        end

%
%         function isConsistent = checkConsistency(obj,qV,rhoV,xiV,mxV,myV,iMx,iMy,iphi)
%             xi  = xiV(iMx,iMy,iphi);
%             rho = rhoV(iMx,iMy,iphi);
%             q = qV(iMx,iMy,iphi);
%             mx = mxV(iMx);
%             my = myV(iMy);
%             errorMx = norm(mx - SuperEllipseParamsRelator.mx(xi,rho,q))/norm(mx)
%             errorMy = norm(my - SuperEllipseParamsRelator.my(xi,rho,q))/norm(my)
%             isMxEqual = errorMx < 1e-14;
%             isMyEqual = errorMy < 1e-14;
%             isConsistent = isMxEqual && isMyEqual;
%         end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function varargout = drawOrientedBox(box, varargin)
%DRAWORIENTEDBOX Draw centered oriented rectangle.
%
%   Syntax
%   drawOrientedBox(BOX)
%   drawOrientedBox(BOX, 'PropertyName', propertyvalue, ...)
%
%   Description
%   drawOrientedBox(OBOX)
%   Draws an oriented rectangle (or bounding box) on the current axis.
%   OBOX is a 1-by-5 row vector containing box center, dimension (length
%   and width) and orientation (in degrees):
%   OBOX = [CX CY LENGTH WIDTH THETA].
%
%   When OBOX is a N-by-5 array, the N boxes are drawn.
%
%   drawOrientedBox(AX, ...)
%   Specifies the axis to draw to point in. AX should be a handle to a axis
%   object. By default, display on current axis.
%
%   HB = drawOrientedBox(...)
%   Returns a handle to the created graphic object(s). Object style can be
%   modified using syntaw like:
%   set(HB, 'color', 'g', 'linewidth', 2);
%
%   Example
%     % draw an ellipse together with its oriented box
%     elli = [30 40 60 30 20];
%     figure;
%     drawEllipse(elli, 'linewidth', 2, 'color', 'g');
%     hold on
%     box = [30 40 120 60 20];
%     drawOrientedBox(box, 'color', 'k');
%     axis equal;
%
%   See also
%     orientedBox, drawPolygon, drawRect, drawBox, drawCenteredEdge
%

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2011-05-09, using Matlab 7.9.0.529 (R2009b)
% Copyright 2011-2022 INRA - Cepia Software Platform

%% Parses input arguments

% extract handle of axis to draw on
if isAxisHandle(box)
    ax = box;
    box = varargin{1};
    varargin(1) = [];
else
    ax = gca;
end

if length(varargin) > 4 && sum(cellfun(@isnumeric, varargin(1:4))) == 4
    % input given as separate arguments
    cx  = box;
    cy  = varargin{1};
    hl   = varargin{2} / 2;
    hw   = varargin{3} / 2;
    theta   = varargin{4};
    varargin = varargin(5:end);

else
    % input given as packed array
    cx  = box(:,1);
    cy  = box(:,2);
    hl   = box(:,3) / 2;
    hw   = box(:,4) / 2;
    theta = box(:,5);
end

%% Draw each box

% allocate memory for graphical handle
hr = zeros(length(cx), 1);

% iterate on oriented boxes
for i = 1:length(cx)
    % pre-compute angle data
    cot = cosd(theta(i));
    sit = sind(theta(i));

    % x and y shifts
    lc = hl(i) * cot;
    ls = hl(i) * sit;
    wc = hw(i) * cot;
    ws = hw(i) * sit;

    % coordinates of box vertices
    vx = cx(i) + [-lc + ws; lc + ws ; lc - ws ; -lc - ws ; -lc + ws];
    vy = cy(i) + [-ls - wc; ls - wc ; ls + wc ; -ls + wc ; -ls - wc];

    % draw polygons
    hr(i) = plot(ax, vx, vy, varargin{:});
end

%% Format output

if nargout > 0
    varargout = {hr};
end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function [U] = hyperAmee(M, q, Smin, Smax, Imax)
% HYPERAMEE Performs the AMEE algorithm to find q endmembers
%   Performs the Automated Morphological Endmember Extraction (AMEE)
%  algorithm to find q endmembers. If only M is
%  given as input, this function calls hyperHfcVd to estimate the number
%  of endmembers (q) and then hyperPct to reduce dimensionality to (q-1).
%
% Usage
%   [U] = hyperAmee(M, q, Smin)
%   [U] = hyperAmee(M, q, Smin, Smax)
%   [U] = hyperAmee(M, q, Smin, Smax, Imax)
% Inputs
%   M    - 3d matrix of HSI data (m x n x p)
%   q    - Number of endmembers to find
%   Smin - minimum kernel size (Smin x Smin)
%   Smax - maximum kernel size (Smax x Smax)
%   Imax - maximum iterations
% Outputs
%   U - Recovered endmembers (p x q)
%
% References
%   Plaza, Antonio, et al. "Spatial/spectral endmember extraction
% by multidimensional morphological operations." Geoscience and
% Remote Sensing, IEEE Transactions on 40.9 (2002): 2025-2041.

% Error trapping
if nargin < 3 || nargin > 5
    help hyperAmee
    error('Not enough input arguments. See function description.')
elseif nargin == 3
    % i.e. Smax and Imax not given
    Smax = Smin;
    Imax = 1;
elseif nargin == 4
    % i.e. if Imax not given
    if Smin <= Smax
        Imax = round((Smax-Smin)/2)+1;
    else
        help hyperAmee
        error('Smax cannot be less than Smin. See function description.')
    end
elseif nargin == 5
    if Smin == Smax && Imax < 2
        help hyperAmee
        error('Cannot iterate more than once with these S limits.')
    elseif Imax > (Smax-Smin+1)
        Imax = round((Smax-Smin)/2)+1;
    end
end

if ndims(M) ~= 3
    error('Input image must be (m x n x p)');
else
    [h, w, p] = size(M);
end

% Build pixel vectors from M
pixVec = cell(h, w);
for hIter = 1:h
    for wIter = 1:w
        pixVec{hIter,wIter} = squeeze(M(hIter,wIter,:));
    end
end

% Morphological Eccentricity Index Score (MEI)
MEI = zeros(h, w);

% Kernel (structuring element)
for B = round(linspace(Smin, Smax, Imax));
    % Create non-overlapping arrays (to reduce computational load)
    hArray = 1:B:h; hArray(end) = h-B+1;
    wArray = 1:B:w; wArray(end) = w-B+1;

    % Move B though all pixels in M
    for wPixel = wArray;
        for hPixel = hArray;
            % (hPixel,wPixel) defines the top-left pixel of the current kernel
            ker = pixVec(hPixel:hPixel+B-1, wPixel:wPixel+B-1);
            [x, y, mei] = spatialSearch( ker );
            % global (x,y) from local (x,y)
            x = x+hPixel-1;
            y = y+wPixel-1;
            % set MEI value at max pixel location
            MEI(x, y) = mei;
        end
    end
end

% Find q largest MEI values
[tmp,idx] = sort(MEI(:), 'descend');
top       = idx(1:q)';

% Return endmembers
U = cell2mat(pixVec(top));

end % hyperAmee function

function [xMax, yMax, mei] = spatialSearch( ker )
%function [xMax, yMax, mei] = spatialSearch( ker )
% spatialSearch finds the max and min cumulative distances
% between each pixel vector and its neighbors inside a kernel
% of size (B x B) and computes the local MEI.

if ~iscell(ker)
    error('ker is not a cell! See line 57 of hyperAmee.m')
end

B = size(ker,1);

dist = zeros(B);
for i = 1:B;
    for j = 1:B;
        dist(i,j) = Dist(ker{i,j}, ker);
    end
end

% Morphological erosion to find minimum pixel vector in region B
[tmp,rowIdx] = min(dist);
[val,colIdx] = min(tmp);
xMin = rowIdx(colIdx);
yMin = colIdx;

% Morphological dilation to find maximum pixel vector in region B
[tmp,rowIdx] = max(dist);
[val,colIdx] = max(tmp);
xMax = rowIdx(colIdx);
yMax = colIdx;

% MEI computation
mei = Dist(ker{xMin, yMin}, ker{xMax, yMax});

end % spatialSearch function

function dist = Dist(a, C)
%function dist = Dist(a, C)
% Cumulative distance measure between a pixel vector and its neighbor(s).
% Plaza et al. used the Spectral Angle Mapper (SAM) measure.
%
% Inputs
%   a - pixel vector of interest (p x 1)
%   C - all pixels in neighborhood (or a single pixel vector)
% Outpus
%   dist - (cumulative) SAM distance

if ~iscell(C)
    dist = acos(dot(a,C)/(norm(a)*norm(C)));
else
    dist = 0;
    for k = 1:numel(C);
        dist = dist + acos(dot(a,C{k})/(norm(a)*norm(C{k})));
    end
end

end % Dist function

author: davidkun

repo: HyperSpectralToolbox


text:
function W = randInitializeWeights(L_in, L_out)
%RANDINITIALIZEWEIGHTS Randomly initialize the weights of a layer with L_in
%incoming connections and L_out outgoing connections
%   W = RANDINITIALIZEWEIGHTS(L_in, L_out) randomly initializes the weights
%   of a layer with L_in incoming connections and L_out outgoing
%   connections.
%
%   Note that W should be set to a matrix of size(L_out, 1 + L_in) as
%   the first row of W handles the "bias" terms
%

% You need to return the following variables correctly
W = zeros(L_out, 1 + L_in);

% ====================== YOUR CODE HERE ======================
% Instructions: Initialize W randomly so that we break the symmetry while
%               training the neural network.
%
% Note: The first row of W corresponds to the parameters for the bias units
%

epsilon_init = 0.12;
W = rand(L out, 1 + L in) * 2 * epsilon init − epsilon init;

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

end

author: schneems

repo: Octave

save_path: github-repos/MATLAB/schneems-Octave


text:
function outval = AG_evaluation(img)
% OUTVAL = AVG_GRADIENT(IMG)

if nargin == 1
    img = double(img);
    % Get the size of img
    [r,c,b] = size(img);

    dx = 1;
    dy = 1;
    for k = 1 : b
        band = img(:,:,k);
        [dzdx,dzdy] = gradient(band,dx,dy);
        s = sqrt((dzdx .^ 2 + dzdy .^2) ./ 2);
        g(k) = sum(sum(s)) / ((r - 1) * (c - 1));
    end
    outval = mean(g);
else
    error('Wrong number of input!');
end

end

author: Linfeng-Tang

repo: Image-Fusion


text:
function [im_out,alpha_out] = texture_map(V,F,TC,im_in,alpha_in)
  % TEXTURE_MAP  texture map a give mesh (V,F) with a texture (TC,im_in) and save
  % it back to an image (im_out)
  %
  % [im_out,alpha_out] = texture_map(V,F,TC,im_in,alpha_in)
  %
  % Inputs:
  %   V  #V by 2 list of mesh vertex positions
  %   F  #F by 3 list of face indices
  %   TC  #V by 2 list of mesh texture coordinates, these should be between 0
  %     and 1, if they're not they will be normalized
  %   im_in  input texture image
  %   alpha_in  input texture image alpha channel
  % Outputs:
  %   im_out  deformed texture image
  %
  % Example:
  %  % original mesh (V,F) deformed mesh (NV,F)
  %  TC = texture_coords(V);
  %  [oim, oalpha] = texture_map(V,F,TC,im,alpha);
  %

  if ~exist('alpha_in')
    alpha_in = ones(size(im_in,1),size(im_in,2));
  end
  if size(V,2) == 2
    V(:,end+1:3) = 0;
  end

  path_to_TexMapPreview = ...
    '/usr/local/igl/igl_apps/TexMapPreview/TexMapPreview';
  path_to_convert = '/usr/local/bin/convert';
  %path_to_TexMapPreviewConvert = ...
  %  '/usr/local/igl/igl_apps/TexMapPreview/TexMapPreviewConvert.sh';

  prefix = tempprefix();
  png_file_name = [prefix '.png'];
  tga_file_name = [prefix '.tga'];
  obj_file_name = [prefix '.obj'];
  output_tga_file_name = [prefix '-out.tga'];
  output_png_file_name = [prefix '-out.png'];

  % write image to png file
  imwrite(im_in,png_file_name,'Alpha',alpha_in);

  % convert the png image to a .tga image
  [status,result] = ...
    issue( ...
      'export DYLD_LIBRARY_PATH="";', ...
      path_to_convert, ...
      png_file_name, ...
      tga_file_name);
  % error on error
  if(status ~= 0)
      error(result);
  end
  % write mesh to obj file
  writeOBJ(obj_file_name,V,F,TC,F);

  % convert the png image to a .tga image
  [status,result] = ...
    issue( ...
      path_to_TexMapPreview, ...
      obj_file_name, ...
      tga_file_name, ...
      output_tga_file_name);
  % error on error
  if(status ~= 0)
      error(result);
  end

  % convert output tga to png

  % convert the png image to a .tga image
  [status, result] = ...
    issue( ...
      'export DYLD_LIBRARY_PATH="";', ...
      path_to_convert, ...
      output_tga_file_name, ...
      output_png_file_name);
  % error on error
  if(status ~= 0)
      error(result);
  end

  % clean up temporary files
  delete(tga_file_name);
  delete(png_file_name);
  delete(obj_file_name);
  delete(output_tga_file_name);

  if(nargout == 1)
    im_out = imread(output_png_file_name);
  elseif(nargout == 2)
    [im_out,m,alpha_out] = imread(output_png_file_name);
  else
    % clean up png file before error
    delete(output_png_file_name);
    error('Too many outputs');
  end
  % clean up png file
  delete(output_png_file_name);

  function [status, result] = issue(varargin)
  % ISSUE issue a system call with each arguement separated by a space
    command = [];
    for ii = 1:nargin
      command = ['' command ' ' varargin{ii}];
    end
    % show the command
    fprintf('%s\n',command);
    % issue the command
    [status, result] = system( command );
  end

end

author: alecjacobson

repo: gptoolbox


text:
function ok=dfsession(action,fn)
%DFSESSION Clear, load, or save a Distribution Fitting session

%   $Revision: 1.1.6.7 $  $Date: 2004/02/01 22:10:39 $
%   Copyright 2003-2004 The MathWorks, Inc.

% Create a structure with version information
str.ftype = 'Distribution Fitting session';   % type of file
str.version = 1;                              % the most current version
str.allversions = [1];                        % all supported versions
str.properties = {'Color' 'LineStyle' 'LineWidth' 'Marker' 'MarkerSize'};

% Variables we save, and the number required to be in a saved file
varnames = {'ftype' 'version' 'allds' 'dsinfo' 'allfits' ...
            'fitinfo' 'alldists' 'outliers' 'guistate'};
nrequired = 9;

if nargin<2
   fn = '';
end

% Wrap real function with a try/catch block so we can control the legend
oldleg = dfgetset('showlegend');
if isequal(oldleg,'on')
   dfgetset('showlegend','off');
end

ok = true;
%try
   switch(action)
    case 'save'
      ok = savesession(fn,str,varnames);

    case 'load'
      ok = loadsession(fn,str,varnames,nrequired);

    case 'clear'
      ok = clearsession;

    otherwise
      ok = false;
   end
%catch
%end

if isequal(oldleg,'on')
   dfgetset('showlegend','on');
   dfupdatelegend(dfgetset('dffig'));
end

% ----------------------------------------------------------------
function ok=savesession(fn,str,varnames)
%DFSAVESESSION Callback to save a distribution fitting session to a file

% Extract some variables from the input structure
ftype = str.ftype;
version = str.version;
guistate.ftype = dfgetset('ftype');
guistate.dtype = dfgetset('dtype');
guistate.binDlgInfo = dfgetset('binDlgInfo');

% Get all M data set object instances and some properties
dsdb = getdsdb;
allds = find(dsdb);
allds(allds==dsdb) = [];
nds = length(allds);
dsinfo = cell(nds,1);
for j=1:length(allds)
   % Save all datasets
   dj = allds(j);
   dsinfo{j,1} = dj.plot;    % remember this flag separately from object
end

% Get all M fit object instances and some properties
fitdb = getfitdb;
allfits = find(fitdb);
allfits(allfits==fitdb) = [];
nfits = length(allfits);
fitinfo = cell(nfits,1);
for j=1:nfits
   % Save all fit objects separately
   fj = allfits(j);
   fitinfo{j,1} = fj.plot;  % remember this flag separately from object
end

% get the user-defined distributions
alldists = dfgetset('alldistributions');

% Get the outliers (excluded sets)
outdb = getoutlierdb;
outliers = find(outdb);
outliers(outliers==outdb) = [];

% Get file name to use, remember the directory name
olddir = dfgetset('dirname');
filespec = [olddir '*.dfit'];
if isempty(fn)
   [fn,pn] = uiputfile(filespec,'Save Session');
   if isequal(fn,0) || isequal(pn,0)
      ok = false;
      return
   end
   if ~ismember('.',fn)
      fn = [fn '.dfit'];
   end
   dfgetset('dirname',pn);
   fn = [pn fn];
end

% Select a file and save the session variables
le = lasterr;
lasterr('');
try
   save(fn, varnames{:}, '-mat');
catch
   uiwait(errordlg(sprintf('Error saving session file:\n%s', lasterr),...
                   'Save Error','modal'))
   ok = false;
   return
end
lasterr(le);
ok = true;

% ----------------------------------------------------------------
function ok=loadsession(fn,str,varnames,nrequired)
%DFLOADSESSION Callback to load a saved distribution fitting session

import com.mathworks.toolbox.stats.*;

ok = true;
ftype = str.ftype;
version = str.version;
allversions = str.allversions;
properties = str.properties;

% Get file name and load from it, remember the directory name
olddir = dfgetset('dirname');
filespec = [olddir '*.dfit'];

if isempty(fn)
   [fn,pn] = uigetfile(filespec,'Load Session');
   if isequal(fn,0) || isequal(pn,0)
      return
   end
   if ~ismember('.',fn)
      fn = [fn '.cfit'];
   end
   dfgetset('dirname',pn);
   fn = [pn fn];
end

% Clear current session
clearsession;

% Get file contents without adding them to the data bases automatically
dsmgr = DataSetsManager.getDataSetsManager;
fmgr = FitsManager.getFitsManager;
dsmgr.turnOffUDDListener;
fmgr.turnOffUDDListener;
try
   s = load('-mat',fn);
catch
   uiwait(errordlg(sprintf('Error loading session file:\n%s', lasterr),...
                   'Load Error','modal'))
   dsmgr.turnOnUDDListener;
   return
end
dsmgr.turnOnUDDListener;
fmgr.turnOnUDDListener;

for j=1:nrequired
   if ~isfield(s,varnames{j})
      uiwait(errordlg('Not a valid Distribution Fitting session file',...
                   'File Invalid','modal'))
      return
   end
end
if ~isequal(s.ftype,ftype)
   uiwait(errordlg('Not a valid Distribution Fitting session file',...
                   'File Invalid','modal'))
   return
end

if ~ismember(s.version,allversions)
   uiwait(errordlg('Bad version number in Distribution Fitting session file',...
                   'Invalid Version','modal'))
   return
end

% Install the saved distribution definitions
dft = DistributionFitting.getDistributionFitting;
try
   dfsetdistributions(dft,s.alldists);
catch
   uiwait(errordlg(sprintf('Error loading saved session:\n%s',lasterr),...
                   'File Invalid','modal'));
end

% Reset some properties of the gui state
dffig = dfgetset('dffig');
dfsetplottype(dffig, s.guistate.ftype, s.guistate.dtype);
if isfield(s.guistate,'binDlgInfo')
   dfgetset('binDlgInfo',s.guistate.binDlgInfo);
end

hFunctionList = findall(dffig, 'Tag', 'displaylist');
ftypes = getappdata(hFunctionList,'codenames');
value = strmatch(s.guistate.ftype, ftypes, 'exact');
set(hFunctionList, 'Value', value);

isprobplot = isequal(s.guistate.ftype, 'probplot');
if isprobplot
   hDistributionList = findall(dffig, 'Tag', 'typelist');
   dtypes = getappdata(hDistributionList,'okcodenames');
   value = strmatch(s.guistate.dtype, dtypes, 'exact');
   set(hDistributionList, 'Value', value);
end

% Make sure distribution list is accurate
dfupdateppdists(dffig);

% Plot datasets that are flagged for plotting
for j=1:length(s.allds);
   dj = s.allds(j);

   % Now set the plot flag.  Some non-serializable information
   % (line handles, listeners, etc.) will be re-created at this point.
   dj.line = [];
   if s.dsinfo{j}
      dj.plot = 1;
   end
   dsmgr.addDataSet(java(dj), dj.name);
end

% Get exclusion rules (outliers)
omgr = OutliersManager.getOutliersManager;
omgr.init;

% Fix up fit objects
fitdb = getfitdb;
outdb = getoutlierdb;
for j=1:length(s.allfits)
   fj = s.allfits(j);

   % Restore all dataset handles
   dsname = fj.dataset;
   for k=1:length(s.allds)
      if isequal(dsname,s.allds(k).name)
         fj.dshandle = s.allds(k);
         break;
      end
   end

   % Restore all exclusion rule handles
   ername = fj.exclusionrulename;
   if ~isempty(ername)
      erhandle = find(outdb,'name',ername);
      fj.exclusionrule = erhandle;
   end

   % Connect this fit to the fit data base and add to fits manager
   connect(fj,fitdb,'up');

   % Now set the plot flag.  Some non-serializable information
   % (line handles, listeners, etc.) will be re-created at this point.
   fj.line = [];
   if s.fitinfo{j,1}
      fj.plot = 1;
   end
end

% ----------------------------------------------------------------
function ok=clearsession
%DFCLEARSESSION Callback to clear distribution fitting session

ok = true;

% Trigger java listeners to clear all saved java content
import com.mathworks.toolbox.stats.*;
DFToolClearManager.getDFToolClearManager.listenerTrigger;

% Delete all udd fit object instances
fitdb = getfitdb;
ft = down(fitdb);
while(~isempty(ft))
   ftnew = right(ft);
   delete(ft);
   ft = ftnew;
end

% Delete all udd data set object instances
dsdb = getdsdb;
dj = down(dsdb);
while(~isempty(dj))
   if dj.plot
      dj.plot = 0;
   end
   djnew = right(dj);
   delete(dj);
   dj = djnew;
end

% Delete all udd outlier object instances
outdb = getoutlierdb;
outliers = find(outdb);
outliers(outliers==outdb) = [];
delete(outliers);

%init (behind the scenes) analysis and plot

author: zouchuhang

repo: LayoutNet


text:
function theta = lines_imp_angle_2d ( a1, b1, c1, a2, b2, c2 )

%*****************************************************************************80
%
%% LINES_IMP_ANGLE_2D finds the angle between two implicit lines in 2D.
%
%  Discussion:
%
%    The implicit form of a line in 2D is:
%
%      A * X + B * Y + C = 0
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 February 2005
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Adrian Bowyer and John Woodwark,
%    A Programmer's Geometry,
%    Butterworths, 1983.
%
%  Parameters:
%
%    Input, real A1, B1, C1, the implicit parameters of the
%    first line.
%
%    Input, real A2, B2, C2, the implicit parameters of the
%    second line.
%
%    Output, real THETA, the angle between the two lines.
%
  pdotq = a1 * a2 + b1 * b2;
  pnorm = sqrt ( a1 * a1 + b1 * b1 );
  qnorm = sqrt ( a2 * a2 + b2 * b2 );

  theta = r8_acos ( pdotq / ( pnorm * qnorm ) );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ x, w ] = legendre_set_sqrtx_01 ( n )

%*****************************************************************************80
%
%% LEGENDRE_SET_SQRTX_01 sets a Gauss-Legendre rule for SQRT(X) * F(X) on [0,1].
%
%  Discussion:
%
%    The integral:
%
%      Integral ( 0 <= X <= 1 ) SQRT ( X ) * F(X) dX =
%      Integral ( 0 <= Y <= 1 ) 2 * Y**2 * F(Y**2) dY.
%      (using Y = SQRT(X) )
%
%    The quadrature rule:
%
%      Sum ( 1 <= I <= N ) W(I) * F ( X(I) )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 June 2010
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Daniel Zwillinger, editor,
%    CRC Standard Mathematical Tables and Formulae,
%    CRC Press, 30th Edition, 2000, page 696.
%
%  Parameters:
%
%    Input, integer N, the order.
%
%    Output, real X(N), the abscissas.
%
%    Output, real W(N), the weights.
%
  n2 = 2 * n + 1;

  [ x2, w2 ] = legendre_set ( n2 );

  x = zeros ( n, 1 );
  w = zeros ( n, 1 );

  x(1:n,1) = x2(n+2:2*n+1,1).^2;

  w(1:n,1) = 2.0 * w2(n+2:2*n+1,1) .* x(1:n,1);

  return
end


text:
%% Sample execution for Baseline on CFNet-conv3
% hyper-parameters reported in Supp.material for CVPR'17, Table 2 for arXiv version
tracker_par.join.method = 'corrfilt';
tracker_par.net = 'baseline-conv3-on-cfnet-conv3_e100.mat';
tracker_par.net_gray = 'baseline-conv3-on-cfnet-conv3_gray_e70.mat';
tracker_par.scaleStep = 1.034;
tracker_par.scalePenalty = 0.9820;
tracker_par.scaleLR = 0.66;
tracker_par.wInfluence = 0.27;
tracker_par.zLR = 0.008;

[~,~,dist,overlap,~,~,~,~] = run_tracker_evaluation('all', tracker_par);

author: bertinetto

repo: cfnet

save_path: github-repos/MATLAB/bertinetto-cfnet


text:
function [ x, w ] = gen_laguerre_ss_compute ( n, alpha )

%*****************************************************************************80
%
%% GEN_LAGUERRE_SS_COMPUTE computes a generalized Gauss-Laguerre quadrature rule.
%
%  Discussion:
%
%    The integral:
%
%      Integral ( 0 <= X < +oo ) EXP ( - X ) * X^ALPHA * F(X) dX
%
%    The quadrature rule:
%
%      Sum ( 1 <= I <= N ) W(I) * F ( X(I) )
%
%    The integral:
%
%      Integral ( 0 <= X < +oo ) X^ALPHA * F(X) dX
%
%    The quadrature:
%
%      Sum ( 1 <= I <= N ) W(I) * EXP ( X(I) ) * F ( X(I) )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 June 2010
%
%  Author:
%
%    Original FORTRAN77 version by Arthur Stroud, Don Secrest.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Arthur Stroud, Don Secrest,
%    Gaussian Quadrature Formulas,
%    Prentice Hall, 1966,
%    LC: QA299.4G3S7.
%
%  Parameters:
%
%    Input, integer N, the order.
%    N must be at least 1.
%
%    Input, real ALPHA, the exponent of the X factor.
%    Set ALPHA = 0.0 for the simplest rule.
%    ALPHA must be nonnegative.
%
%    Output, real X(N), the abscissas.
%
%    Output, real W(N), the weights.
%
  x = zeros ( n, 1 );
  w = zeros ( n, 1 );
%
%  Set the recursion coefficients.
%
  for i = 1 : n
    b(i) = ( alpha + 2 * i - 1 );
  end

  for i = 1 : n
    c(i) = ( i - 1 ) * ( alpha + i - 1 );
  end

  cc = gamma ( alpha + 1.0 ) * prod ( c(2:n) );

  for i = 1 : n
%
%  Compute an estimate for the root.
%
    if ( i == 1 )

      xval = ( 1.0 + alpha ) * ( 3.0 + 0.92 * alpha ) ...
        / ( 1.0 + 2.4 * n + 1.8 * alpha );

    elseif ( i == 2 )

      xval = xval + ( 15.0 + 6.25 * alpha ) / ( 1.0 + 0.9 * alpha + 2.5 * n );

    else

      r1 = ( 1.0 + 2.55 * ( i - 2 ) ) / ( 1.9 * ( i - 2 ) );

      r2 = 1.26 * ( i - 2 ) * alpha / ( 1.0 + 3.5 * ( i - 2 ) );

      ratio = ( r1 + r2 ) / ( 1.0 + 0.3 * alpha );

      xval = xval + ratio * ( xval - x(i-2) );

    end
%
%  Use iteration to find the root.
%
    [ xval, dp2, p1 ] = gen_laguerre_ss_root ( xval, n, alpha, b, c );
%
%  Set the abscissa and weight.
%
    x(i) = xval;
    w(i) = ( cc / dp2 ) / p1;

  end

  return
end


text:
% new procedure: generate motor seed/map in trialRes space
% find top 100 cells correlating to motor, screen with Rh4+Rh5(+Rh6??) masks

% corr sweep (0.5-0.7) to get top %2 cells, then kmeans, save; rank by stim-lock, save
clear all; close all; clc

%% folder setup
isSaveFig = 1;
isPlotFig = 1;

outputDir = GetOutputDataDir;
saveDir = [];
saveDir{1} = fullfile(outputDir,'motor_map_notlrRes_1_kmeans_030417');
% saveDir{2} = fullfile(outputDir,'motor_map_notlrRes_2_autoclus_030417');
if ~exist(saveDir{1}, 'dir'), mkdir(saveDir{1}), end;
% if ~exist(saveDir{2}, 'dir'), mkdir(saveDir{2}), end;

%% init

hfig = figure;
InitializeAppData(hfig);
ResetDisplayParams(hfig);

setappdata(hfig,'isMotorseed',1);
setappdata(hfig,'isTrialRes',1);

%% run fish
range_fish = GetFishRange;%[1:3,5:18];
M_thres_reg = zeros(3,18);
M_numTopCorr = zeros(1,18);
M_motorseedRegs = cell(1,18);
M_compareMotorCellNumber = zeros(2,18);

for i_fish = range_fish
    ClusterIDs = [11,2];
    [gIX_seed,gIX_seed,M,stim,behavior,M_0] = LoadSingleFishDefault(i_fish,hfig,ClusterIDs);

    MASKs = getappdata(hfig,'MASKs');
    CellXYZ_norm = getappdata(hfig,'CellXYZ_norm');
    absIX = getappdata(hfig,'absIX');

    %% regression (in tRes), thresholding by % of cells (instead of corr thres)

    %     [~,Reg_tRes] = GetTrialAvrLongTrace(hfig,Reg);
    %     [~,M_0_tRes] = GetTrialAvrLongTrace(hfig,M_0);

    Reg = FindClustermeans(gIX_seed,M);
%     Reg_LR = Reg;
    Reg_LR = Reg-repmat(mean(Reg),2,1);

    %     Corr = corr(Reg_tRes_LR',M_0_tRes');
    Corr = corr(Reg_LR',M_0');
    [Corr_sorted,IX_corr] = sort(Corr,2,'descend');

    %     [corr_max,IX] = max(Corr,[],1);
    %     [~,I] = sort(corr_max,'descend');

    M_target_prctcell= [1,2,3];
    for i_numcell = 1%:2,
        nCells_total = size(M_0,1);
        prctcell = M_target_prctcell(i_numcell);
        nCells_target = round(prctcell/100 * nCells_total);

        %% target cell number: counting cells in hindbrain only

        Msk_IDs = 114; % mask for full hindbrain

        % isScreenMskFromAllCells
        cIX = (1:length(absIX))';
        gIX = ones(size(cIX));
        [cIX_hb,gIX_hb] = ScreenCellsWithMasks(Msk_IDs,cIX,gIX,MASKs,CellXYZ_norm,absIX);

        I_hb = ismember(IX_corr,cIX_hb);
        cum_I_hb1 = cumsum(I_hb(1,:));
        cum_I_hb2 = cumsum(I_hb(2,:));
        lastIX{1} = find(cum_I_hb1==nCells_target,1,'first');
        lastIX{2} = find(cum_I_hb2==nCells_target,1,'first');

        % cut off the desired number of cells to display
        for i_lr = 1:2
            cIX = IX_corr(i_lr,1:lastIX{i_lr})';
            gIX = ceil((1:length(cIX))'/(length(cIX)/20));
            M = UpdateIndices_Manual(hfig,cIX,gIX);
            [~,M_tRes] = GetTrialAvrLongTrace(hfig,M);
            setappdata(hfig,'M',M_tRes);

            M_thres_reg(i_lr,i_fish) = Corr_sorted(i_lr,lastIX{i_lr});
            M_compareMotorCellNumber(i_lr,i_fish) = length(cIX);

            %% clustering: trying Autoclus here!
            if i_numcell==1
                %%
                numK = 16; % manual
                gIX = Kmeans_Direct(M_tRes,numK);
            else
                cIX_reg = cIX;%(1:size(M_0,1))';
                [cIX,gIX] = AutoClustering(cIX,gIX,M_0,cIX_reg);
                UpdateIndices_Manual(hfig,cIX,gIX);
            end
            %% Plot anat
            if isPlotFig,

                I = LoadCurrentFishForAnatPlot(hfig,cIX,gIX);
                f{i_lr} = DrawCellsOnAnat(I);
%                 f{i_lr} = figure('Position',[50,100,1400,800]);
%                 % isCentroid,isPlotLines,isPlotBehavior,isPlotRegWithTS
%                 subplot(121)
%                 setappdata(hfig,'isPlotBehavior',1);
%                 setappdata(hfig,'isStimAvr',0);
%                 setappdata(hfig,'isPlotLines',0);
%                 %             UpdateTimeIndex(hfig);
%                 DrawTimeSeries(hfig,cIX,gIX);
%
%                 % right plot
%                 ax = subplot(122);
%                 I = LoadCurrentFishForAnatPlot(hfig,cIX,gIX);
%                 DrawCellsOnAnat(I,ax);

            end
        end

        combineFiguresLR(f{1},f{2});

        % Save plot
        if isSaveFig,
            filename = fullfile(saveDir{i_numcell}, ['Fish',num2str(i_fish),'_motormap_',num2str(prctcell),'%_each']);
            saveas(gcf, filename, 'png');
            close(gcf)
        end

    end
end

%%
% range_fish excludes Fish 4
% M_compareMotorCellNumber(:,4) = NaN;
% figure;bar(M_compareMotorCellNumber')

author: xiuyechen

repo: FishExplorer

path: github-repos/MATLAB/xiuyechen-FishExplorer/FishExplorer-c61392cf0835480d64fc03c15f1992935fdc7106/figure scripts/SensoryMotor/arc code/motormap, tRes versions, comparisons, w Clusters or Autoclus to follow/fig3_motormap_trialRes_lrRes.m


text:
function pcut = addBilinearVariableCuts(p)

if isempty(p.bilinears)
    pcut = p;
else
pcut = emptyNumericalModel;

z = p.bilinears(:,1);
x = p.bilinears(:,2);
y = p.bilinears(:,3);

nn = length(p.c);

still_uncertain = find(abs(p.lb(z)-p.ub(z))>1e-8);
if ~isempty(still_uncertain)

    x_lb = p.lb(x);
    x_ub = p.ub(x);
    y_lb = p.lb(y);
    y_ub = p.ub(y);
    m = length(x);
    one = ones(m,1);
    general_vals =[x_lb.*y_lb one -y_lb -x_lb,x_ub.*y_ub one -y_ub -x_ub,-x_ub.*y_lb -one y_lb x_ub,-x_lb.*y_ub -one y_ub x_lb]';
    general_cols = [one z+1 x+1 y+1 one z+1  x+1 y+1 one z+1 x+1 y+1 one z+1 x+1 y+1]';
    general_row = [1;1;1;1;2;2;2;2;3;3;3;3;4;4;4;4];

    quadratic_row = [1;1;1;2;2 ;2;3; 3; 3;4;4;4;5;5;5;6;6;6];
    quadratic_cols =  [one  z+1 x+1 ,one z+1 x+1 ,one  z+1 x+1, one z+1 x+1 ,one z+1 x+1, one z+1 x+1]';
    x1 = (3*x_ub+x_lb)/4;
    x2 = (x_ub+3*x_lb)/4;
    x3 = (x_ub+x_lb)/2;
    quadratic_vals = [-x_ub.*x_lb -one x_lb+x_ub x_lb.*y_lb one -y_lb-x_lb x_ub.*y_ub one -y_ub-x_ub  x1.*x1 one -x1-x1 x2.*x2 one -x2-x2 x3.*x3 one -x3-x3]';
    m = 1+length(p.c);
    rows = [];
    cols = [];
    vals = [];
    nrow = 0;

    for i =still_uncertain(:)'
        x = p.bilinears(i,2);
        y = p.bilinears(i,3);
        if x~=y
            rows = [rows;general_row+nrow];
            vals = [vals;general_vals(:,i)];
            cols = [cols;general_cols(:,i)];
            nrow = nrow + 4;
        else
            col = quadratic_cols(:,i);
            val = quadratic_vals(:,i);

            rows = [rows;quadratic_row+nrow];
            vals = [vals;val];
            cols = [cols;col];
            nrow = nrow + max(quadratic_row);
        end
    end

    F_temp = sparse(rows,cols,vals,nrow,m);
    keep = find(~isinf(F_temp(:,1)) & ~isnan(F_temp(:,1)));
    F_temp = F_temp(keep,:);
    pcut.F_struc = [F_temp;pcut.F_struc];
    pcut.K.l = pcut.K.l+size(F_temp,1);
end

pcut = mergeNumericalModels(p,pcut);
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
 function [deriv, curv] = trl_dercurv(data, li, curvtype, iblock, nblock)
%function [deriv, curv] = trl_dercurv(data, li, curvtype, iblock, nblock)
%
% evaluate derivatives and curvatures for monoenergetic Poisson transmission
% negative log-likelihoods: f(x) = \sum_i h_i(l), h_i(l) = mi(l) - yi log mi(l)
% mi(l) = bi exp(-l) + ri
%
% in
%	data	{yi, bi, ri}
%				yi is [nb,na].  bi, ri are same or scalar
%				passed to trl_curvature()
%	li	[nb,#view_in_block]
%	curvtype ''		curvature type
%	iblock, nblock		for OS type methods
%
% out
%	deriv	[nb,#]		\dot hi(l)
%	curv	[nb,#]		surrogate curvature for hi at li
%
% Copyright 2004-2-1, Jeff Fessler, The University of Michigan

if nargin < 2, ir_usage, end
if nargin < 3, curvtype = 'pc'; end

yi = data{1};
bi = data{2};
ri = data{3};

if nargin == 5
	ia = iblock:nblock:size(yi,2);
	yi = yi(:,ia);
	if length(bi) > 1
		bi = bi(:,ia);
	end
	if length(ri) > 1
		ri = ri(:,ia);
	end
end

% transmission Poisson likelihood function
ei = exp(-li);
mi = bi .* ei + ri;
deriv = (1 - yi ./ mi) .* (-bi .* ei);

curv = trl_curvature(yi, bi, ri, li, curvtype);

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
% sim_onepulse_shaped.m
% Jamie Near, McGill University 2014.
%
% USAGE:
% out = sim_onepulse_shaped(n,sw,Bfield,linewidth,sys,RF,tp,phCyc,dfdx,G)
%
% DESCRIPTION:
% %This function simulates the effect of a frequency selective or slice
% selective excitation, followed immediately by the acquisition window.
% This is mainly an exercise to see if I can get slice selective
% excitation working.
%
% Note that when simulating a frequency selective pulse, it is okay to
% specify only 8 arguments (no gradient needs to be specified).  If the
% 9th argument, G, is specified and is non-zero, then a slice selective
% pulse is assumed.
%
% INPUTS:
% n         = number of points in fid/spectrum
% sw        = desired spectral width in [Hz]
% Bfield    = main magnetic field strength in [T]
% linewidth = linewidth in [Hz]
% sys       = spin system definition structure
% RF        = RF pulse definition structure (obtain using 'io_loadRFwaveform.m')
% tp        = RF pulse duration in [ms]
% phCyc     = Phase of excitation rf pulse in [degrees].
% dfdx      = if simulating a frequency selective pulse, this argument
%             should be the frequency offset [Hz].  If simulating a slice
%             selective pulse, this argument should be the position offset [cm].
% G         = gradient strength for slice-selective pulse [G/cm];
%
% OUTPUTS:
% out       = simulated spectrum, in FID-A structure format, using pulse-acquire
%             sequence.

function out = sim_onepulse_shaped(n,sw,Bfield,linewidth,sys,RF,tp,phCyc,dfdx,G)

if nargin>9 && G~=0
    simType='g';
else
    simType='f';
end

%Set water to centre
centreFreq=4.65;
for k=1:length(sys)
    sys(k).shifts=sys(k).shifts-centreFreq;
end

%Calculate Hamiltonian matrices and starting density matrix.
[H,d]=sim_Hamiltonian(sys,Bfield);

%BEGIN PULSE SEQUENCE************
if simType=='g'
    d=sim_shapedRF(d,H,RF,tp,90,90+phCyc,dfdx,G);   %slice selective excitation
elseif simType=='f'
    d=sim_shapedRF(d,H,RF,tp,90,90+phCyc,dfdx);     %frequency selective excitation
end
[out,dout]=sim_readout(d,H,n,sw,linewidth,90);      %Readout along y (90 degree phase);
%END PULSE SEQUENCE**************

%Correct the ppm scale:
out.ppm=out.ppm-(4.65-centreFreq);

%Fill in structure header fields:
out.seq='onepulse';
out.te=tp/2;
out.sim='shaped';

%Additional fields for compatibility with FID-A processing tools.
out.sz=size(out.specs);
out.date=date;
out.dims.t=1;
out.dims.coils=0;
out.dims.averages=0;
out.dims.subSpecs=0;
out.dims.extras=0;
out.averages=1;
out.rawAverages=1;
out.subspecs=1;
out.rawSubspecs=1;
out.flags.writtentostruct=1;
out.flags.gotparams=1;
out.flags.leftshifted=0;
out.flags.filtered=0;
out.flags.zeropadded=0;
out.flags.freqcorrected=0;
out.flags.phasecorrected=0;
out.flags.averaged=1;
out.flags.addedrcvrs=1;
out.flags.subtracted=1;
out.flags.writtentotext=0;
out.flags.downsampled=0;
out.flags.isFourSteps=0;

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
function [I J] = idiag(sz, k)
% function [I J] = idiag(sz, k) % OR
% I = itril(sz, k)
%
% Return the subindices [I J] (or linear indices I if single output call)
% in the purpose of extracting the diagonal of the matrix of the size SZ.
% Input k is optional shifting. For k=0, extract from the main
% diagonal. For k>0 -> above the diagonal, k<0 -> below the diagonal
%
% Output is a column and sorted with respect to linear indice
%
% Example:
%
% A = [ 7     5     4
%       4     2     3
%       9     1     9
%       3     5     7 ]
%
% I = idiag(size(A))  % gives [1 6 11]'
% A(I)                % gives [7 2 9]' OR diag(A)
%
% Author: Bruno Luong <brunoluong@yahoo.com>
% Date: 21/March/2009

if isscalar(sz)
    sz = [sz sz];
end
m=sz(1);
n=sz(2);

% Main diagonal by default
if nargin<2
    k=0;
end

% Pay attention to the clipping
I = (1-min(k,0):min(m,n-k)).';
J = I+k;

if nargout<2
    % convert to linear indices
    I = sub2ind([m n], I, J);
end

end % idiag

author: Sable

repo: mcbench-benchmarks


text:
function varargout = ECGBeatFitter(varargin)
%
% ECGBeatFitter(ECG,Phase,ExpParamName),
% Overview: Graphical user interface for ECG approximation with Gaussian kernels.
%
% inputs:
% ECG: a single ECG waveform used for model training
% Phase: the phase corresponding to the ECG waveform
% ExpParamName (optional): The name of the vector containing the exported parameter (default:OptimumParams)
% ExpParamName (optional): The name of the vector containing the exported parameter (default:OptimumParams)
% Title (optional): The title of the plot
%
% outputs:
% OptimizedParameters
%
% Open Source ECG Toolbox, version 2.0, April 2008
% Released under the GNU General Public License
% Copyright (C) 2008  Reza Sameni
% Sharif University of Technology, Tehran, Iran -- GIPSA-LAB, INPG, Grenoble, France
% reza.sameni@gmail.com

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

% ECGBEATFITTER M-file for ECGBeatFitter.fig
%      ECGBEATFITTER, by itself, creates a new ECGBEATFITTER or raises the existing
%      singleton*.
%
%      H = ECGBEATFITTER returns the handle to a new ECGBEATFITTER or the handle to
%      the existing singleton*.
%
%      ECGBEATFITTER('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in ECGBEATFITTER.M with the given input arguments.
%
%      ECGBEATFITTER('Property','Value',...) creates a new ECGBEATFITTER or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before ECGBeatFitter_OpeningFunction gets called.  An
%      unrecognized property name or invalid value makes property
%      application
%      stop.  All inputs are passed to ECGBeatFitter_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help ECGBeatFitter

% Last Modified by GUIDE v2.5 08-Oct-2007 13:42:43

% Begin initialization code - DO NOT EDIT
gui_Singleton = 0;
gui_State = struct('gui_Name',       mfilename, ...
    'gui_Singleton',  gui_Singleton, ...
    'gui_OpeningFcn', @ECGBeatFitter_OpeningFcn, ...
    'gui_OutputFcn',  @ECGBeatFitter_OutputFcn, ...
    'gui_LayoutFcn',  [] , ...
    'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before ECGBeatFitter is made visible.
function ECGBeatFitter_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to ECGBeatFitter (see VARARGIN)

%setappdata(hObject,'SelectedIndeces',0);
mn = varargin{1};
sd = varargin{2};
phase = varargin{3};
handles.mdlError = 'ModelErrorPercentage';
if (nargin>6),
    handles.ExpParamName = varargin{4};
else
    handles.ExpParamName = 'OptimumParams';
end

if (nargin>7),
    set(handles.text3,'String',varargin{5});
end

handles.Linehandles = [];
handles.Initmodhandle = [];
handles.Optmodhandle = [];

handles.SelectedIndeces = [];
handles.ECGmean = mn;
handles.ECGsd = sd;
handles.ECGphase = phase;
handles.OptimizedParameters = [];

axes(handles.axes1);

errorbar(mn,sd/2,'b');
hold on;
plot(mn,'r','linewidth',2);
grid on;
xlabel('Sample index');
ylabel('Amplitude (mV)');
set(handles.axes1,'ButtonDownFcn',@MyButtonDownFcn);
ECGcurve = get(handles.axes1,'Children');
set(ECGcurve,'ButtonDownFcn',@MyButtonDownFcn);
axis tight;
set(handles.figure1,'Pointer','crosshair');

% Choose default command line output for ECGBeatFitter
handles.output = hObject;
%handles.output = 0;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes ECGBeatFitter wait for user response (see UIRESUME)
uiwait(handles.figure1);

% --- Outputs from this function are returned to the command line.
function varargout = ECGBeatFitter_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure

% varargout{1} = handles.output;
varargout{1} = handles.OptimizedParameters;

% --- Executes on button press in Reset.
function Reset_Callback(hObject, eventdata, handles)
% hObject    handle to Reset (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

data = guidata(gcbo);

% delete(data.Initmodhandle);
delete(data.Optmodhandle);
% data.Initmodhandle = [];
data.Optmodhandle = [];

N = length(data.SelectedIndeces);
for i = 1:N
    delete(data.Linehandles(i));
end
data.SelectedIndeces = [];
data.Linehandles = [];

set(data.kernelnumber,'String','');
set(data.opterror,'String','');

guidata(gcbo,data);

% --- Executes on button press in Optimize.
function Optimize_Callback(hObject, eventdata, handles)
% hObject    handle to Optimize (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

data = guidata(gcbo);

% delete(data.Initmodhandle);
delete(data.Optmodhandle);

% localpeaks = zeros(size(data.ECGmean));
I = max(round(data.SelectedIndeces),1);
I = min(I,length(data.ECGmean));

P = length(I);
% localpeaks(I) = 1;

tetai = data.ECGphase(I(1:P));
alphai = 1.2*data.ECGmean(I(1:P));
bi = .04*ones(size(alphai));

options = optimset('TolX',1e-4,'TolFun',1e-4,'MaxIter',100);
InitParams = [alphai bi tetai];

OptParams = nlinfit(data.ECGphase,data.ECGmean,@ECGModel,InitParams,options);
% OptParams = lsqnonlin(@(InitParams) ECGModelError(InitParams,ECGmn,Phasemn,0),InitParams,InitParams-2,InitParams+2,options);
% Model0 = ECGModelError(InitParams,data.ECGmean,data.ECGphase,1);
Model = ECGModelError(OptParams,data.ECGmean,data.ECGphase,1);

ax = axis;
% data.Initmodhandle = plot(Model0,'g');
data.Optmodhandle = plot(Model,'k','LineWidth',2);

L = length(OptParams);

ph = OptParams(2*L/3+1:L);
for i = 1:L/3,
    if(ph(i)<min(data.ECGphase) || ph(i)>max(data.ECGphase))
        ph(i) = mean(data.ECGphase); % place the possible erroneous estimate at the middle of the ECG beat (close to the R-peak)
    end

    pos = find(data.ECGphase>=ph(i),1);

    I = max(pos,1);
    I = min(I,length(data.ECGphase));
    data.SelectedIndeces(i) = I;
    set(data.Linehandles(i),'XData',[I I]);
end
% % %
% % % P = length(I);
% % % % localpeaks(I) = 1;
% % %
% % % tetai = data.ECGphase(I(1:P));

axis(ax);

data.OptimizedParameters = OptParams;

er = 100*mean((Model-data.ECGmean).^2)/mean(data.ECGmean.^2);

set(data.kernelnumber,'String',length(data.SelectedIndeces));
set(data.opterror,'String',round(100*er)/100);

% Modified 25/1/2019 by Reza Sameni:
% assignin('base',handles.mdlError,er);
assignin('caller',handles.mdlError,er);

guidata(gcbo,data);

% --- Executes on button press in ExportData.
function ExportData_Callback(hObject, eventdata, handles)
% hObject    handle to ExportData (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Modified 25/1/2019 by Reza Sameni:
% assignin('base',handles.ExpParamName,handles.OptimizedParameters);
assignin('caller',handles.ExpParamName,handles.OptimizedParameters);

% handles.output = 1; % Function has finished execution
% guidata(hObject, handles);
uiresume(handles.figure1);

function MyButtonDownFcn(hObject, eventdata, handles)
pn = get(gca,'Currentpoint');
% hndl = plot(pn(1,1),pn(1,2),'ro');
ax = axis;
ln = line([pn(1,1),pn(1,1)],[ax(3),ax(4)]);
set(ln,'Color',[0 1 .5]);

data = guidata(gcbo);
data.SelectedIndeces = [data.SelectedIndeces pn(1,1)];
data.Linehandles = [data.Linehandles ln];

guidata(gcbo,data);

% --------------------------------------------------------------------
function Helpmenu_Callback(hObject, eventdata, handles)
% hObject    handle to Helpmenu (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% --------------------------------------------------------------------
function About_1_Callback(hObject, eventdata, handles)
% hObject    handle to About_1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
AboutPath = which('About.htm');
web(AboutPath);

% --------------------------------------------------------------------
function Help_1_Callback(hObject, eventdata, handles)
% hObject    handle to Help_1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

HelpPath = which('Help.htm');
web(HelpPath);

% --- Executes on button press in Exit.
function Exit_Callback(hObject, eventdata, handles)
% hObject    handle to Exit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% assignin('base',handles.ExpParamName,[]);
uiresume(handles.figure1);

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox


text:
function n = p10_n ( )

%*****************************************************************************80
%
%% P10_N returns the number of variables for problem 10.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 October 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer N.  If N is positive, N represents the
%    only legal value for N for this problem.  If N is
%    negative, the absolute value of N is the least legal
%    value of N, but other values are allowable.
%
  n = 2;

  return
end


text:
function stability=step_computeStability(partRate,net,layerID,theConf,fileRoot,partList)
selectPatternRatio=1.0;
patchNumPerPattern=100;
load(sprintf('%s/imdb.mat',fileRoot),'meta');
Name_batch=net.meta.classes.name{1};

if(strcmp(Name_batch,'n01443537'))
    partList=1;
end

convnet=getConvNetPara(net);
net=vl_simplenn_tidy(net);
net=our_vl_simplenn_move(net,'gpu') ;
assert(layerID>0);
stability=computeStability(partRate,net,Name_batch,selectPatternRatio,patchNumPerPattern,partList,theConf,layerID,convnet,meta.dataMean);
disp(stability);
end

repo: interpretableCNN


text:
%A demo of how Adaboost works.
%
%Some base learners have external dependancies. CART requires the stats
%toolbox (classregtree). RBF SVM requires libsvm (which must be above
%bioinformatics toolbox in the path). NeuralNetwork requires netlab.
%
%The demo displays some 2-D data, and shows the development of the decision
%boundary, margins, error rate and margin distribution as a boosted
%ensemble is trained.
%
%The data is split into two classes - red and blue. The decision boundary
%appears as a yellow contour. The margins are shown on a white->black
%gradient, where white is a confident classification as blue and black is a
%confident classification as red.
%
%To add a single learner, use the "Add Learner" button, or "Add 10" to add
%10 learners. Learners are shown in a scrollable list,
%with their adjusted weights beside them. Weight adjustment involves
%normalisation and then scaling so that 1 represents an average learner
%weight. These weights can be adjusted (even made negative), and learners
%can be enabled or disabled individual using the check box. Right-clicking
%on the checkbox disables all other learners. The 'Reset' button enables
%all learners and sets their weights to those computed by Adaboost.
%
%The 'Plot Error' and 'Plot Margins' buttons open new windows to display
%addition information. 'Plot Error' shows the training error (blue) and the
%generalisation error (red) against the size of the ensemble. 'Plot
%Margins' shows the cdf of the margins for the current ensemble. The height
%of the line above 0 indicates the proportion of the data with that margin
%or below.
%
%The two drop-down menus at the top allow a choice of base learner
%algorithms and data generating functions. When either of these are
%changed, the ensemble is completely reset for the new learner type or data
%set.
classdef boosting_demo < handle

  properties (Constant=true)

    %Starting size of the ensemble
    initial_size = 1;
    %Resolution of the decision boundaries and margins. Large numbers make
    %it slower.
    resolution = 200;

  end

  properties

    fig;
    plot_area;
    learner_labels;
    learner_check;
    learner_weights;
    weight_labels;
    learner_count;
    learner_enabled;
    pretend_weights;
    adaboost;
    inputs;
    outputs;
    g_inputs;
    example_weights;
    g_outputs;
    bgc = [0.8 0.8 0.8];
    data_handle;
    data_colors;
    margin_handle;
    boundary_handle;
    margins;
    predictions;
    x_grid;
    y_grid;
    add_learner_button;
    add_10_button;
    enable_all_button;
    error_rate_button;
    error_fig;
    plotting_error;
    margin_distribution_button;
    margin_fig;
    plotting_margin;
    learner_scroll;
    first_learner = 1;
    base_learner;
    base_learners;
    base_learner_names;
    base_learner_popup;
    data_func;
    data_funcs;
    data_func_names;
    data_func_popup;
    data_points;
    running = true;

  end

  methods

    function a = boosting_demo()

      %Number of data points to generate
      a.data_points = 200;
      %To add a new base learner, add the construction function to this
      %list, then add the name to a.base_learner_names. The base learner
      %must have train(inputs, outputs) and outputs = test(inputs)
      %functions.
      base_learner_num = 6;
      a.base_learners = {@()SVM() @()LinearRegression @()CART() @()OnlineNaiveBayes(0) @()NeuralNetwork(10, 0.01, 2, 'linear', 50), @()Stump()};
      a.base_learner_names = {'RBF SVM', 'Linear Regression', 'CART', 'Naive Bayes', 'Neural Network', 'Decision Stump'};
      a.base_learner = a.base_learners{base_learner_num};
      %To add a new data generator, add the function to this list, where c
      %is the number of data points to generate. The function must return
      %[inputs outputs] where inputs is cX2 and outputs is cX1. Then add
      %the name to the second list.
      data_func_num = 3;
      a.data_funcs = {@(c)DataGen.circle(c, 5, 0.5), @(c)DataGen.checkerboard(c, 10, 2, rand * pi), ...
        @(c)DataGen.gaussians(c, [-2 -2], [2 2], [2 2], [2 2])};
      a.data_func_names = {'Circle', 'Checkerboard', 'Gaussians'};
      a.data_func = a.data_funcs{data_func_num};

      %Make the ensemble
      a.adaboost = Adaboost(boosting_demo.initial_size, a.base_learner);
      %Sort out the main window
      screen_size = get(0, 'ScreenSize');
      a.fig = figure('Position', [(screen_size(3) - 800) / 2 (screen_size(4) - 600) / 2 800 600], 'MenuBar', 'none', ...
        'CloseRequestFcn', @(varargin)a.halt(), 'Name', 'Boosting Demo', 'NumberTitle', 'off');

      %Create the error and margin graphs.
      a.error_fig = figure('Name', 'Error Rates', 'NumberTitle', 'off');
      a.margin_fig = figure('Name', 'Margin Distribution', 'NumberTitle', 'off');
      set(a.error_fig, 'visible', 'off');
      set(a.margin_fig, 'visible', 'off');
      figure(a.fig);
      a.plot_area = axes('Units','Pixels');

      hold on;

      colormap('gray');
      %Plot the data
      a.data_handle = scatter([], [], [], '+', 'SizeDataSource', 'a.example_weights', 'XDataSource', 'a.inputs(:, 1)',...
        'YDataSource', 'a.inputs(:, 2)', 'CDataSource', 'a.data_colors', 'HitTest', 'off',...
        'LineWidth', 3);

      %Plot the margin and boundary contours
      min = -5;
      max = 5;
      step = (max - min) / boosting_demo.resolution;
      range = min:step:max;
      [a.x_grid a.y_grid] = meshgrid(range, range);
      a.margins = repmat(range, size(a.x_grid, 2), 1);
      [null a.margin_handle] = contourf(a.x_grid, a.y_grid, a.margins, 'XDataSource', 'a.x_grid', ...
        'YDataSource', 'a.y_grid', 'ZDataSource', 'a.margins', 'LineStyle', 'none');
      [null a.boundary_handle] = contour(a.x_grid, a.y_grid, a.margins > 0, 1, 'y', 'XDataSource', 'a.x_grid', ...
        'YDataSource', 'a.y_grid', 'ZDataSource', 'a.margins > 0', 'LineWidth', 2);

      hold off;

      %Make GUI components
      a.base_learner_popup = uicontrol('Style', 'popupmenu', 'String', a.base_learner_names, ...
        'Callback', @(varargin)a.base_learner_change(), 'Value', base_learner_num);
      a.data_func_popup = uicontrol('Style', 'popupmenu', 'String', a.data_func_names, ...
        'Callback', @(varargin)a.data_func_change(), 'Value', data_func_num);
      a.add_learner_button = uicontrol('String', 'Add Learner', 'Callback', @(varargin)a.add_learner());
      a.add_10_button = uicontrol('String', 'Add 10', 'Callback', @(varargin)a.add_10());
      a.enable_all_button = uicontrol('String', 'Reset', 'Callback', @(varargin)a.enable_all());
      a.error_rate_button = uicontrol('String', 'Plot Error', 'Callback', @(varargin)a.plot_error());
      a.margin_distribution_button = uicontrol('String', 'Plot Margins', ...
        'Callback', @(varargin)a.plot_margins());
      a.learner_scroll = uicontrol('Style', 'Slider', 'min', 0, 'max', 1);

     % listener = handle.listener(a.learner_scroll,'ActionEvent',@(varargin)a.scroll_learners());
      listener = addlistener(a.learner_scroll,'Action',@(varargin)a.scroll_learners());

      setappdata(a.learner_scroll, 'listener', listener);

      %Generate the data
      a.random_data();
      %Train the ensemble (First learner)
      a.train_adaboost();
      %Layout function
      set(a.fig, 'ResizeFcn', @(varargin)a.do_layout());

    end

    %Stop the program, called when the window is closed
    function halt(a)
      a.running = false;
      if (ishandle(a.error_fig))
        close(a.error_fig);
      end
      if (ishandle(a.margin_fig))
        close(a.margin_fig);
      end
      delete(a.fig);
    end

    function do_layout(a)
      %Layout. Do the side bar, and make the data take up the rest of the
      %window
      try
        figure(a.fig);
        pos = get(a.fig, 'Position');
        fig_width = pos(3) - 370;
        fig_height = pos(4) - 60;
        fig_size = min(fig_width, fig_height);
        set(a.plot_area, 'Position', [350 (pos(4) - fig_size - 30) fig_size fig_size]);
        top = pos(4) - 40;
        set(a.base_learner_popup, 'Position', [10 top 150 20]);
        set(a.data_func_popup, 'Position', [170 top 150 20]);
        top = top - 30;
        set(a.error_rate_button, 'Position', [10 top 100 20]);
        set(a.margin_distribution_button, 'Position', [110 top 100 20]);
        set(a.add_10_button, 'Position', [210 top 100 20]);
        top = top - 30;
        set(a.add_learner_button, 'Position', [10 top 150 20]);
        set(a.enable_all_button, 'Position', [170 top 130 20]);
        top = top - 30;
        set(a.learner_scroll, 'Position', [310 20 20 top]);
        %Create new base learner components as necessary
        for i=1:a.adaboost.k_max
          if i > numel(a.learner_labels)
            a.learner_labels(i) = uicontrol('Style','text','String',['Learner ' num2str(i)], ...
              'BackgroundColor',a.bgc,'HorizontalAlignment','left', 'FontSize', 14,...
              'FontWeight', 'bold');
            a.learner_check(i) = uicontrol('Style', 'Checkbox', 'ButtonDownFcn', ...
              @(varargin)a.individual_learner(i), 'Callback', @(varargin)a.enable_learner(i));
            a.learner_weights(i) = uicontrol('Style','Slider','Min',-1,'Max',1);
            a.weight_labels(i) = uicontrol('Style','text','String','0',...
              'FontSize', 14, 'FontWeight', 'bold');
            listener = handle.listener(a.learner_weights(i),'ActionEvent',@(varargin)a.set_weight(i));
            setappdata(a.learner_weights(i), 'listener', listener);
          end
        end
        for i=1:a.first_learner-1
            set(a.learner_labels(i), 'Visible', 'off');
            set(a.learner_check(i), 'Visible', 'off');
            set(a.learner_weights(i), 'Visible', 'off');
            set(a.weight_labels(i), 'Visible', 'off');
        end
        for i=a.first_learner:a.adaboost.k_max
          if (top > 20)
            set(a.learner_labels(i), 'Visible', 'on');
            set(a.learner_check(i), 'Visible', 'on');
            set(a.learner_weights(i), 'Visible', 'on');
            set(a.weight_labels(i), 'Visible', 'on');
            set(a.learner_labels(i), 'Position',[10 top 110 20]);
            set(a.learner_check(i), 'Position',[120 top 20 20]);
            set(a.learner_weights(i), 'Position',[150 top 100 20]);
            set(a.weight_labels(i), 'Position',[260 top 40 20]);
          else
            set(a.learner_labels(i), 'Visible', 'off');
            set(a.learner_check(i), 'Visible', 'off');
            set(a.learner_weights(i), 'Visible', 'off');
            set(a.weight_labels(i), 'Visible', 'off');
          end
          top = top - 30;
        end
      catch exception
        %Buh
      end

    end

    %Called when the base learner type is changed. Reset the ensemble.
    function base_learner_change(a)
      learner_num = get(a.base_learner_popup, 'Value');
      a.base_learner = a.base_learners{learner_num};
      a.new_ensemble();
    end

    %Called when the data generating function is changed. Reset the
    %ensemble and generate new data.
    function data_func_change(a)
      data_num = get(a.data_func_popup, 'Value');
      a.data_func = a.data_funcs{data_num};
      a.random_data();
      a.new_ensemble();
    end

    %Create a new ensemble, and delete any lingering rubbish.
    function new_ensemble(a)
      for i=2:a.learner_count
        delete(a.learner_labels(i));
        delete(a.learner_check(i));
        delete(a.learner_weights(i));
        delete(a.weight_labels(i));
      end
      a.learner_labels = a.learner_labels(1);
      a.learner_check = a.learner_check(1);
      a.learner_weights = a.learner_weights(1);
      a.weight_labels = a.weight_labels(1);
      a.adaboost = Adaboost(boosting_demo.initial_size, a.base_learner);
      a.learner_count = 1;
      a.first_learner = 1;
      set(a.learner_scroll, 'max', 1);
      set(a.learner_scroll, 'Value', 0);
      a.adaboost.k_max = a.learner_count;
      a.adaboost.k_next = a.learner_count;
      a.learner_enabled = true;
      a.train_adaboost();
    end

    %When the learner scrollbar moves.
    function scroll_learners(a)

      a.first_learner = max(ceil(a.learner_count - get(a.learner_scroll, 'Value')), 1);
      a.do_layout();

    end

    %Plot an error against ensemble size graph in a separate window.
    function plot_error(a)

      preds = cell2mat(arrayfun(@(x)a.adaboost.learners{x}.test(a.inputs) .* a.outputs, ...
        1:a.learner_count, 'UniformOutput', 0));
      weighted_preds = preds .* repmat(a.pretend_weights .* a.learner_enabled, numel(a.outputs), 1);
      cumulative_margins = cumsum(weighted_preds, 2);
      g_preds = cell2mat(arrayfun(@(x)a.adaboost.learners{x}.test(a.g_inputs) .* a.g_outputs, ...
        1:a.learner_count, 'UniformOutput', 0));
      g_weighted_preds = g_preds .* repmat(a.pretend_weights .* a.learner_enabled, numel(a.g_outputs), 1);
      g_cumulative_margins = cumsum(g_weighted_preds, 2);
      figure(a.error_fig);
      set(a.error_fig, 'Name', 'Error Rates', 'NumberTitle', 'off');
      clf;
      hold on;
      plot(1:a.learner_count, mean(cumulative_margins < 0), 'LineWidth', 2);
      plot(1:a.learner_count, mean(g_cumulative_margins < 0), 'r', 'LineWidth', 2);
      hold off;
      a.plotting_error = true;

    end

    %Plot the margin distribution in a separate window.
    function plot_margins(a)

      preds = cell2mat(arrayfun(@(x)a.adaboost.learners{x}.test(a.inputs) .* a.outputs, ...
        1:a.learner_count, 'UniformOutput', 0));
      weighted_preds = preds .* repmat(a.pretend_weights .* a.learner_enabled, numel(a.outputs), 1);
      final_margins = sum(weighted_preds, 2);
      range = -1:0.01:1;
      figure(a.margin_fig);
      set(a.margin_fig, 'Name', 'Margin Distribution', 'NumberTitle', 'off');
      clf;
      plot(range, mean(repmat(final_margins, 1, numel(range)) < repmat(range, numel(final_margins), 1)))
      a.plotting_margin = true;

    end

    %Add 10 learners.
    function add_10(a)

      try
        i = 0;
        while (i < 10 && a.running)
          a.add_learner();
          i = i + 1;
        end
      catch
      end

    end

    %Add one learner.
    function add_learner(a)

      a.learner_count = a.learner_count + 1;
      set(a.learner_scroll, 'max', a.learner_count);
      set(a.learner_scroll, 'Value', a.learner_count - a.first_learner);
      a.adaboost.k_max = a.learner_count;
      a.adaboost.k_next = a.learner_count;
      a.learner_enabled(a.learner_count) = true;
      a.adaboost.train(a.inputs, a.outputs);
      a.pretend_weights = a.adaboost.learner_weights;
      new_predictions = a.adaboost.learners{a.learner_count}.test([a.x_grid(:) a.y_grid(:)]);
      a.predictions = [a.predictions new_predictions];
      a.example_weights = abs(a.adaboost.d_weights) * a.data_points * 25;
      a.do_layout();
      a.refresh_weights();
      a.refresh_enabled();
      a.redraw();

    end

    %Train the first learner in Adaboost.
    function train_adaboost(a)

       a.adaboost.train(a.inputs, a.outputs);
       a.learner_count = a.adaboost.k_max;
       a.learner_enabled = ones(1, a.learner_count);
       a.pretend_weights = a.adaboost.learner_weights';
       a.predictions = cell2mat(arrayfun(@(x)a.adaboost.learners{x}.test([a.x_grid(:) a.y_grid(:)]), ...
        1:a.learner_count, 'UniformOutput', 0));
       a.do_layout();
       a.refresh_weights();
       a.refresh_enabled();
       a.redraw();

    end

    %Set the weights back to Adaboost weights.
    function reset_weights(a)
      a.pretend_weights = a.adaboost.learner_weights;
      a.refresh_weights();
    end

    %Enable a base learner
    function enable_learner(a, i)
      a.learner_enabled(i) = ~a.learner_enabled(i);
      set(a.learner_check(i), 'Value', a.learner_enabled(i));
      a.redraw();
    end

    %Enable all base learners.
    function enable_all(a)
      for j=1:a.learner_count
        a.learner_enabled(j) = true;
        set(a.learner_check(j), 'Value', true);
      end
      a.reset_weights();
      a.redraw();
    end

    %Show one single base learner.
    function individual_learner(a, i)

      if strcmp(get(gcf,'SelectionType'), 'normal')
        a.enable_learner(i);
      else
        for j=1:a.learner_count
          if (j ~= i)
            a.learner_enabled(j) = false;
            set(a.learner_check(j), 'Value', false);
          else
            a.learner_enabled(j) = true;
            set(a.learner_check(j), 'Value', true);
          end
        end
        a.redraw();
      end
    end

    %Change the weight of a base learner.
    function set_weight(a, i)
      a.pretend_weights(i) = get(a.learner_weights(i), 'Value');
      a.refresh_weights();
      a.redraw();
    end

    %Redraw the data
    function redraw(a)

      try
        a.margins = reshape(a.predictions * (a.pretend_weights .* a.learner_enabled)', ...
          size(a.x_grid, 1), size(a.x_grid, 2));
        a.margins = a.margins ./ max(abs(a.margins(:)));

        refreshdata(a.margin_handle, 'caller');
        refreshdata(a.boundary_handle, 'caller');
        refreshdata(a.data_handle, 'caller');
        drawnow;
        if (a.plotting_error)
          a.plot_error();
        end
        if (a.plotting_margin)
          a.plot_margins();
        end
      catch
      end

    end

    %Update the base learner weights (normalise etc.)
    function refresh_weights(a)
      a.pretend_weights = a.pretend_weights / sum(abs(a.pretend_weights));
      for i=1:a.learner_count
        set(a.learner_weights(i), 'Value', a.pretend_weights(i));
        set(a.weight_labels(i), 'String', floor(a.data_points * a.pretend_weights(i) * sum(a.learner_enabled)) / a.data_points);
      end
    end

    %Update the checkboxes to show if learners are enabled.
    function refresh_enabled(a)
      for i=1:a.learner_count
        set(a.learner_check(i), 'Value', a.learner_enabled(i));
      end
    end

    %Generate some data from the required function.
    function random_data(a)

      [a.inputs a.outputs] = a.data_func(a.data_points);
      [a.g_inputs a.g_outputs] = a.data_func(a.data_points.*2);
      a.data_colors = [a.outputs == -1 zeros(a.data_points, 1) a.outputs == 1];
      a.example_weights = ones(a.data_points, 1) * a.data_points;

    end

  end

end

author: epfl-lasa