text:
function [alg] = learnDisturbanceQuantile(input_signal, alg)

%__________________________________________________________________________
%% Documentation
%
% Authors:      Alexander Wischnewski (alexander.wischnewski@tum.de)
%
% Start Date:   21.11.2019
% Last Update:  12.03.2020
%
% Description:
%   learns a quantile estimate for the given input signal using a batch algorithm (Wischnewski, 2020)
%
% Inputs:
%   input_signal    Input signal used for quantile estimation
%   alg             algorithm intermediate results and parameters
%
% Algorithm intermediate results and parameters
%   alg.r_est               Estimate for the q-quantile
%   alg.q_target            Target percentage of the quantile
%   alg.lambda              Learning rate
%   alg.N                   Samples for one batch
%   alg.count_lower         Count of samples in this batch lower than the quantile estimate
%   alg.count               Count of samples in this batch
%__________________________________________________________________________
%% Algorithm

% if the input signal is smaller than the current quantile estimate,
% increase counter of signals lower than quantile estimate. This is used to estimate the
% percentage of the data covered by the estimate.
if(input_signal < alg.r_est)
    alg.count_lower = alg.count_lower + 1;
end
alg.count = alg.count + 1;
%
% if the batch size is reached, update the estimate
if(alg.count >= alg.N)
    alg.r_est = alg.r_est - alg.lambda*(alg.count_lower/alg.count - alg.q_target);
    alg.count_lower = 0;
    alg.count = 0;
end

author: TUMFTM


text:
function colornames_deltaE(palette,rgb)
% Create a figure comparing the color difference (deltaE) calculations used in COLORNAMES.
%
% (c) 2014-2019 Stephen Cobeldick
%
%%% Syntax:
%  colornames_deltaE(palette,rgb)
%
% Creates a figure showing the supplied colormap as horizontal color bands,
% overlaid with columns of the closest named colors from the selected
% palette. Each column shows one color difference (deltaE) calculation.
%
% Note: Requires the function COLORNAMES and its associated MAT file (FEX 48155).
%
% For more information on color difference concepts and formulae:
% https://en.wikipedia.org/wiki/Color_difference
% http://www.colorwiki.com/wiki/Delta_E:_The_Color_Difference
%
%% Examples %%
%
% colornames_deltaE('html4',jet(18))
%
% colornames_deltaE('x11',summer(18))
%
% colornames_deltaE('matlab',jet(18))
%
%% Input and Output Arguments %%
%
%%% Inputs (all inputs are optional):
% palette = CharRowVector, the name of a palette supported by COLORNAMES.
% rgb     = Numeric Array, size Nx3, each row is an RGB triple (0<=rgb<=1).
%
%%% Outputs:
% none
%
% See also COLORNAMES COLORNAMES_CUBE COLORNAMES_VIEW MAXDISTCOLOR COLORMAP

% Get list of palettes and color difference metrics:
[fnm,~,dtE] = colornames();
%
% Define default values:
if nargin<2
	N = 16;
	map = get(0,'defaultFigureColormap');
	rgb = interp1(linspace(1,N,size(map,1)),map,1:N);
end
if nargin<1
	palette = fnm{1+rem(round(now*1e7),numel(fnm))};
end
%
persistent fgh axh txh
%
% Text lightness threshold:
thr = 0.54;
%
if isempty(fgh)||~ishghandle(fgh)
	fgh = figure('HandleVisibility','callback', 'IntegerHandle','off',...
		'NumberTitle','off', 'Name',mfilename, 'Color','white', 'Toolbar','none');
	axh = axes('Parent',fgh, 'Visible','off', 'XTick',[], 'YTick',[],...
		'Units','normalized', 'Position',[0,0,1,1]);
else
	try %#ok<TRYNC>
		cla(axh)
	end
end
%
assert(ischar(palette)&&isrow(palette),'First input <palette> must be a 1xN char.')
idz = strcmpi(palette,fnm);
assert(any(idz),'Palette ''%s'' is not supported. Call COLORNAMES() to list all palettes.',palette)
set(fgh,'Name',sprintf('%s (palette = ''%s'')',mfilename,fnm{idz}))
%
assert(ismatrix(rgb)&&size(rgb,2)==3&&isreal(rgb)&&all(0<=rgb(:)&rgb(:)<=1),...
	'Second input argument can be a colormap of RGB values (size Nx3).')
colormap(axh,rgb);
%
N = size(rgb,1);
x = [0;0;1;1];
y = [0;1;1;0];
X = repmat(x,1,N);
Y = bsxfun(@plus,y,N-1:-1:0);%0:N-1);
patch(X,Y,1:N, 'Parent',axh, 'EdgeColor','none', 'FaceColor','flat', 'CDataMapping','direct');
%
dEn = numel(dtE);
[cnc,RGB] = cellfun(@(t)colornames(palette,rgb,t),dtE, 'uni',false);
BAW = cellfun(@(c)(c*[0.298936;0.587043;0.114021])<thr,RGB, 'uni',false);
%
tmp = @(s,n) text((2*n-1)*ones(1,N)/(2*dEn), mean(Y.',2), zeros(1,N),...
	s, 'Parent',axh, 'HorizontalAlignment','center');
txh = cellfun(tmp, cnc, num2cell(1:dEn), 'uni',false);
tmp = @(h,c,b) set(h(:), {'BackgroundColor'},num2cell(c,2), {'Color'},num2cell(b(:,[1,1,1]),2));
cellfun(tmp, txh, RGB, BAW)
%
set(axh,'YLim',[0,N+1]);
text((1:2:2*dEn)/(2*dEn), N+ones(1,dEn)/2, zeros(1,dEn), dtE(:),...
	'Parent',axh, 'HorizontalAlignment','center', 'Color','black');
%
drawnow()
%
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%colornames_deltaE

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
%TRAINED_MAPPING Define trained mapping
%
%   W = TRAINED_MAPPING(A,DATA,DIM)
%
% INPUT
%   A      - Dataset used for training
%   DATA   - Data (cell araay or structure) to be stored in the data-field
%            of the mapping in order to transfer it to the execution part
%   DIM    - Dimensionality of output space.
%
% OUTPUT
%   W      - Mapping
%
% DESCRIPTION
% This is a simplified version of the definition of a trained mapping. It
% calls PRMAPPING and derives all needed information from the dataset A used
% for training the mapping. In DATA everything should be stored needed for
% the execution of the mapping, either in a structure or by a cell array.
%
% SEE ALSO (PRTools Guide)
% MAPPINGS, PRMAPPING, TRAINED_CLASSIFIER, DEFINE_MAPPING, MAPPING_TASK

% Copyright: Robert P.W. Duin, prtools@rduin.nl

function w = trained_mapping(varargin)

  [a,data,out_size] = setdefaults(varargin,[],[],0);
  fname = callername;
  mapname = getname(feval(fname));
  if isdataset(a)
    if out_size == 0, out_size = getsize(a,3); end
    w = prmapping(fname,'trained',data,getlablist(a),size(a,2),out_size);
  else
    if out_size == 0, out_size = 1; end
    w = prmapping(fname,'trained',data,[],size(a,2),out_size);
  end
  w = setname(w,mapname);

return

%CALLERNAME
%
%	NAME = CALLERNAME
%
% Returns the name the calling function

function name = callername

[ss ,i] = dbstack;
if length(ss) < 3
	name = [];
else
	name = ss(3).name;
end

author: marianux

repo: ecg-kit

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


text:
%kkbitnand 'Output = NOT(Input 1) OR NOT(Input 2 or Constant)'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros kbitnand.pane file
%
% Parameters:
% InputFile: i1 'Input 1', required: 'First input data object'
% OutputFile: o 'Output', required: 'Resulting output data object'
% InputFile: i2 'Input 2', optional: 'Second input data object'
% Integer: real 'Constant', default: 1: 'Constant value'
%
% Example: o = kkbitnand({i1, i2}, {'i1','';'o','';'i2','';'real',1})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% kbitnand - Output = NOT(Input 1)  OR  NOT(Input 2 or Constant)
%
%  DESCRIPTION
% The bitwise NAND operator performs a bitwise NAND between each point in
% \fBInput 1" and either \fBInput 2\fP or the \fBConstant\fP value, which
% ever is specified by the user.
%
%  "Data Type - Single Input" 5
% .cI $DATAMANIP/repos/shared/man/sections/value_type_1input
%
%  "Data Type - Two Input Objects" 5
% .cI $DATAMANIP/repos/shared/man/sections/value_type_2input
% The bitwise operators will not operate on float or complex data types.
%
%  "Map Data" 5
% The bitwise operations have not been written to be fully polymorphic yet.
% They does not
% check for map data, and will therefore always operate on the value data,
% even if a map exists.  This will most likely corrupt indexing into the map.
% In the case of a single input, it is recommended to use the "Copy to
% Value" (kcptoval) segment operator to temporarily move the map data into
% the value segment, run the bitwise operation on the data, then move it
% back to the map with the "Copy from Value" (kcpfromval) operator.
% When operating with two input objects, where at least one of the objects
% contains map data, data should be mapped prior to running the bitwise
% operator.  The "Map Data" operator (kmapdata) can be used to perform
% the mapping operation.
%
%  "Validity Mask - Single Input" 5
% .cI $DATAMANIP/repos/shared/man/sections/mask_1input
%
%  "Validity Mask - Two Input Objects" 5
% Masking has not been implemented yet for the bitwise operators.  Therefore,
% only the mask from the first input object will be transferred to the output.
%
%  "Input Objects of Different Sizes" 5
% .cI $DATAMANIP/repos/shared/man/sections/resize_2input
% The value used to pad the data when the input files are not the same size is
% zero.
%
%  "Explicit Location and Time Data - Single Input" 5
% .cI $DATAMANIP/repos/shared/man/sections/loc_and_time_1input
%
%  "Explicit Location and Time Data - Two Input Objects" 5
% The bitwise operations have not been extended to understand location and
% time data.  Therefore, only location and time data present in the first input
% object will be transferred to the output.
%
% Executing \fBNAND" runs the program "kbitwise\fP with the -nand flag.
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
% DATAMANIP::kbitwise
%
%  RESTRICTIONS
%
% THIS ROUTINE is still under development, and has not been modified to fully
% support the polymorphic data model.  See paragraphs above for discussions
% concerning the polymorphic data segments.
%
% All processing is currently performed as UNSIGNED long, so the operations
% may not work properly on negative values.
%
% The bitwise operators will not operate on float or complex data types.
%
% In the future, the padding value used for this operation when the input
% images are of different sizes, will be modified from binary 0 to binary 1.
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

function varargout = kkbitnand(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kkbitnand(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i1', '__input';'o', '__output';'i2', '__input';'real', 1};
maxval={0,0,1,0};
minval={0,0,1,0};
istoggle=[0,0,1,1];
was_set=istoggle * 0;
paramtype={'InputFile','OutputFile','InputFile','Integer'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'kbitwise"  -nand'],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
classdef RegistriesClass < handle

    properties
        userfuncdir
        igroup
        isubj
        isess
        irun
        funcReg
        config
        filename;
    end

    methods

        % ----------------------------------------------------------------
        function obj = RegistriesClass(mode)
            %
            % Syntax:
            %   reg = RegistriesClass();
            %   reg = RegistriesClass('normal');
            %   reg = RegistriesClass('empty');
            %   reg = RegistriesClass('reload');
            %
            % Class containing an array of registries named 'group',
            % 'subj', and 'run'
            %
            %
            % Added mode argument to allow 'empty' option to avoid loading
            % registry but make the methods AND config options available.
            % This is useful for deleting saved Registry.mat when unit testing.
            if nargin==0
                mode = 'normal';
            end

            obj.igroup = 1;
            obj.isubj = 2;
            obj.isess = 3;
            obj.irun = 4;
            obj.funcReg = FuncRegClass().empty();
            obj.filename = '';

            % Get the parameter items from config file relevant to this class
            obj.userfuncdir = FindUserFuncDir(obj);
            if isempty(obj.userfuncdir)
                return
            end
            if strcmp(mode, 'empty')
                return;
            end
            if strcmp(mode, 'reset')
                obj.DeleteSaved();
                return;
            end

            % Check if saved registry exists. If so load that and exit
            if ~strcmp(mode, 'reload')
                if exist([obj.userfuncdir{1}, 'Registry.mat'], 'file')
                    obj.filename = [obj.userfuncdir{1}, 'Registry.mat'];
                    r = load(obj.filename, 'reg');
                    if isa(r.reg, 'RegistriesClass') && ~isempty(r.reg)
                        obj.Copy(r.reg);
                        if obj.IsValid()
                            return;
                        end
                    end
                end
            end

            obj.Load();

            % Save registry for next time
            obj.Save();
        end

        % ----------------------------------------------------------------------------------
        function Load(obj, type)
            if nargin==1
                type = 'all';
            end

            switch(type)
                case 'all'
                    obj.funcReg(obj.igroup) = FuncRegClass('group');
                    obj.funcReg(obj.isubj) = FuncRegClass('subj');
                    obj.funcReg(obj.isess) = FuncRegClass('sess',2);
                    obj.funcReg(obj.irun) = FuncRegClass('run');
                case 'group'
                    obj.funcReg(obj.igroup) = FuncRegClass('group');
                case 'subj'
                    obj.funcReg(obj.isubj) = FuncRegClass('subj');
                case 'sess'
                    obj.funcReg(obj.isess) = FuncRegClass('sess');
                case 'run'
                    obj.funcReg(obj.irun) = FuncRegClass('run');
            end
        end

        % ----------------------------------------------------------------------------------
        function Save(obj)
            reg = obj;
            save([obj.userfuncdir{1}, 'Registry.mat'], 'reg');
        end

        % ----------------------------------------------------------------------------------
        function Reload(obj, type)
            if nargin==1
                type = 'all';
            end

            % Delete saved resistry file
            obj.DeleteSaved();

            % Reload and resave the registry
            obj.Load(type);
            obj.Save();
        end

        % -------------------------------------------------------------------------------
        function err = AddEntry(obj, funcname)
            switch(funcname(1:5))
                case 'hmrG_'
                    itype = obj.igroup;
                case 'hmrS_'
                    itype = obj.isubj;
                case 'hmrE_'
                    itype = obj.isess;
                case 'hmrR_'
                    itype = obj.irun;
                otherwise
                    err = 1;
                    return;
            end
            err = obj.funcReg(itype).AddEntry(funcname);
        end

        % -------------------------------------------------------------------------------
        function err = ReloadEntry(obj, funcname)
            switch(funcname(1:5))
                case 'hmrG_'
                    itype = obj.igroup;
                case 'hmrS_'
                    itype = obj.isubj;
                case 'hmrE_'
                    itype = obj.isess;
                case 'hmrR_'
                    itype = obj.irun;
                otherwise
                    err = 1;
                    return;
            end
            err = obj.funcReg(itype).ReloadEntry(funcname);
        end

        % ----------------------------------------------------------------------------------
        function DeleteSaved(obj)
            if exist([obj.userfuncdir{1}, 'Registry.mat'], 'file')
                delete([obj.userfuncdir{1}, 'Registry.mat']);
            end
        end

        % ----------------------------------------------------------------------------------
        function Copy(obj, obj2)
            obj.igroup = obj2.igroup;
            obj.isubj = obj2.isubj;
            obj.isess = obj2.isess;
            obj.irun = obj2.irun;
            for ii=1:length(obj2.funcReg)
                obj.funcReg(ii) = FuncRegClass(obj2.funcReg(ii));
            end
        end

        % ----------------------------------------------------------------
        function idx = IdxGroup(obj)
            idx=[];
            if isempty(obj)
                return;
            end
            idx = obj.igroup;
        end

        % ----------------------------------------------------------------
        function idx = IdxSubj(obj)
            idx=[];
            if isempty(obj)
                return;
            end
            idx = obj.isubj;
        end

        % ----------------------------------------------------------------
        function idx = IdxRun(obj)
            idx=[];
            if isempty(obj)
                return;
            end
            idx = obj.irun;
        end

        % ----------------------------------------------------------------
        function idx = IdxSess(obj)
            idx=[];
            if isempty(obj)
                return;
            end
            idx = obj.isess;
        end

        % ----------------------------------------------------------------------------------
        function b = IsEmpty(obj)
            b = true;
            if isempty(obj)
                return;
            end
            if isempty(obj.funcReg)
                return;
            end
            if obj.funcReg(obj.irun).IsEmpty()
                return;
            end
            b = false;
        end

        % ----------------------------------------------------------------------------------
        function usagename = GetUsageName(obj, fcall)
            usagename = '';
            if isempty(obj)
                return;
            end
            if nargin<2
                return;
            end
            for ii=1:length(obj.funcReg)
                usagename = obj.funcReg(ii).GetUsageName(fcall);
                if ~isempty(usagename)
                    break;
                end
            end
        end

        % ----------------------------------------------------------------------------------
        function usagenames = GetUsageNames(obj, funcname)
            usagenames = {};
            if isempty(obj)
                return;
            end
            if nargin<2
                return;
            end
            for ii=1:length(obj.funcReg)
                usagenames = obj.funcReg(ii).GetUsageNames(funcname);
                if ~isempty(usagenames)
                    break;
                end
            end
        end

        % ----------------------------------------------------------------------------------
        function fcallstr = GetFuncCallStrDecoded(obj, key, usagename)
            fcallstr = '';
            if isempty(obj)
                return;
            end
            if nargin<3
                return;
            end
            for ii=1:length(obj.funcReg)
                fcallstr = obj.funcReg(ii).GetFuncCallStrDecoded(key, usagename);
                if ~isempty(fcallstr)
                    break;
                end
            end
        end

        % ----------------------------------------------------------------------------------
        function fcall = GetFuncCallDecoded(obj, key, usagename)
            fcall = FuncCallClass().empty();
            if isempty(obj)
                return;
            end
            if nargin<3
                return;
            end
            for ii=1:length(obj.funcReg)
                fcall = obj.funcReg(ii).GetFuncCallDecoded(key, usagename);
                if ~isempty(fcall)
                    break;
                end
            end
        end

        % ----------------------------------------------------------------------------------
        function fcall = FindClosestMatch(obj, fcall0)
            fcall = FuncCallClass().empty();
            if isempty(obj)
                return;
            end
            if nargin<2
                return;
            end
            for ii=1:length(obj.funcReg)
                fcall = obj.funcReg(ii).FindClosestMatch(fcall0);
                if ~isempty(fcall)
                    break;
                end
            end
        end

        % ----------------------------------------------------------------------------------
        function fname = GetSavedRegistryPath(obj)
            fname = obj.filename;
        end

        % ----------------------------------------------------------------------------------
        function n = GetNumUsages(obj, funcname)
            n = [];
            if isempty(obj)
                return;
            end
            if nargin<2
                return;
            end
            for ii = 1:length(obj.funcReg)
                n = obj.funcReg(ii).GetNumUsages(funcname);
                if n>0
                    break;
                end
            end
        end

        % ----------------------------------------------------------------------------------
        function Import(obj, funcpath)
            [~, fname, ext] = fileparts(funcpath);
            if exist([obj.userfuncdir{1}, fname, ext], 'file') == 2
                q = MenuBox('Function already exists in Registry folder. Do you want to replece it',{'Yes','No'});
                if q==2
                    return;
                end
            end
            copyfile(funcpath, obj.userfuncdir{1});
            if strncmp(fname, 'hmrG_', 5)
                type = 'group';
            elseif strncmp(fname, 'hmrS_', 5)
                type = 'subj';
            elseif strncmp(fname, 'hmrE_', 5)
                type = 'sess';
            elseif strncmp(fname, 'hmrR_', 5)
                type = 'run';
            else
                type = 'all';
            end
            obj.Reload(type);
        end

        % ----------------------------------------------------------------------------------
        function b = IsValid(obj)
            b = false;
            regfile = dir([obj.userfuncdir{1}, 'Registry.mat']);
            for ii = 1:length(obj.funcReg)
                try
                    if obj.funcReg(ii).DateLastModified() > datetime(regfile.date,'locale','system')
                        return;
                    end
                catch
                    % pass
                end
            end
            b = true;
        end

        % ----------------------------------------------------------------------------------
        function entry = GetEntryByName(obj, entry_name)
            entry = [];
            for i=1:size(obj.funcReg, 2)
                for j=1:size(obj.funcReg(i).entries, 2)
                    if strcmp(obj.funcReg(i).entries(j).name, entry_name) || strcmp(obj.funcReg(i).entries(j).uiname, entry_name)
                        entry = obj.funcReg(i).entries(j);
                    end
                end
            end
        end

    end

end

author: BUNPC


%% Example script to visualize the aircraft simulation data
% Add the path of the aircraft_3d_animation function
addpath('../src/');
% path of the *.mat file containing the 3d model information
model_info_file = '../3d_models/saab_gripen_3d_model.mat';
% Load the simulation data
% load('scissors_maneuver.mat')
% load('breakaway_maneuver.mat')
% load('split_s_maneuver.mat')
load('departure.mat')
% define the reproduction speed factor
speedx = 1;
% Do you want to save the animation in a mp4 file? (0.No, 1.Yes)
isave_movie = 0;
% Movie file name
movie_file_name = '';

% -------------------------------------------------------------------------
% The frame sample time shall be higher than 0.02 seconds to be able to
% update the figure (CPU/GPU constraints)
frame_sample_time = max(0.02, tout(2)-tout(1));
% Resample the time vector to modify the reproduction speed
t_new   = tout(1):frame_sample_time*(speedx):tout(end);
% Resample the recorded data
act     = interp1(tout, act, t_new','linear');
stick   = interp1(tout, stick, t_new','linear');
y_new   = interp1(tout, yout, t_new','linear');
% We have to be careful with angles with ranges
y_new(:, 7)  = atan2(interp1(tout, sin(yout(:, 7)), t_new','linear'), interp1(tout, cos(yout(:, 7)), t_new','linear')) * 180 / pi;
y_new(:, 8)  = atan2(interp1(tout, sin(yout(:, 8)), t_new','linear'), interp1(tout, cos(yout(:, 8)), t_new','linear')) * 180 / pi;
y_new(:, 9)  = atan2(interp1(tout, sin(yout(:, 9)), t_new','linear'), interp1(tout, cos(yout(:, 9)), t_new','linear')) * 180 / pi;
% Assign the data
heading_deg           =  y_new(:, 7);
pitch_deg             =  y_new(:, 8);
bank_deg              =  y_new(:, 9);
roll_command          = -stick(:, 2);
pitch_command         = -stick(:, 1);
angle_of_attack_deg   =  y_new(:, 2) * 180 / pi;
angle_of_sideslip_deg =  y_new(:, 3) * 180 / pi;
fligh_path_angle_deg  =  y_new(:, 22) * 180 / pi;
mach                  =  y_new(:, 21);
altitude_ft           = -y_new(:, 12);
nz_g                  =  y_new(:, 19);
% Flight control surfaces
le     = act(:, 9);
dr     = act(:, 8);
df1    = act(:, 6);
df2    = act(:, 5);
df3    = act(:, 4);
df4    = act(:, 3);
dfp    = 0.5 * (act(:, 1) + act(:, 2));
% Control array assignation
% (modify the order according to your particular 3D model)
controls_deflection_deg = [dfp(:), dfp(:), le(:), le(:), dr(:), 0.5*(df1(:)+df2(:)), 0.5*(df3(:)+df4(:))];

%% Run aircraft_3d_animation function
% -------------------------------------------------------------------------
aircraft_3d_animation(model_info_file,...
    heading_deg, ...            Heading angle [deg]
    pitch_deg, ...              Pitch angle [deg]
    bank_deg, ...               Roll angle [deg]
    roll_command, ...           Roll  stick command [-1,+1] [-1 -> left,            +1 -> right]
    pitch_command, ...          Pitch stick command [-1,+1] [-1 -> full-back stick, +1 -> full-fwd stick]
    angle_of_attack_deg, ...    AoA [deg]
    angle_of_sideslip_deg, ...  AoS [deg]
    fligh_path_angle_deg, ...   Flight path angle [deg]
    mach, ...                   Mach number
    altitude_ft, ...            Altitude [ft]
    nz_g,  ...                  Vertical load factor [g]
    controls_deflection_deg, ...Flight control deflection (each column is a control surface)
    frame_sample_time, ...      Sample time [sec]
    speedx, ...                 Reproduction speed
    isave_movie, ...            Save the movie? 0-1
    movie_file_name);           % Movie file name


text:
 function varargout = repout(in)
%function varargout = repout(in)
% extension of "deal" to take a single input to multiple outputs, e.g.,
% [a b c] = repout(in)

if nargin < 1, ir_usage, end
for ii=1:nargout
	varargout{ii} = in;
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
  function fmap = mri_field_map_cyl()
%| synthesize field map associated with a cylinder of uniform susceptibility
%| within an extended object also having uniform susceptibility
%| cylinder is aligned with static field

% todo: use
%| lucas:93:dii
%| aksit:07:tpm
%| yoder:04:msw

fail 'todo'


text:
function display(X)
%DISPLAY   Display a CHEBFUN2V.
%   DISPLAY(F) outputs important information about the CHEBFUN2V F to the
%   command window, including its domain of definition, length (number of pivots
%   used to represent it), and a summary of its structure.
%
%   It is called automatically when the semicolon is not used at the end of a
%   statement that results in a CHEBFUN2V.

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

if ( isequal(get(0, 'FormatSpacing'), 'compact') )
	disp([inputname(1), ' =']);
	disp(X);
else
	disp(' ');
	disp([inputname(1), ' =']);
	disp(' ');
    disp(X);
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
 function [ hHead,hTail,hPlate ] = arrow3d( orgin,r,h1,h2,dir)
%ARROW3D Summary of this function goes here
%  Detailed explanation goes here

%
% orgin=[0 0 0.2]';
% r=1;
% h=90;
% dir='y';

n=25;
% p=0.15;

hTail=Cylinder(orgin,r/2,h1,dir,n,'open');
hold on

if dir=='x'
    org=orgin+[h1; 0; 0];
elseif dir=='y'
    org=orgin+[0; h1; 0];
elseif dir=='z'
    org=orgin+[0; 0; h1];
end

xlabel('x')
ylabel('y')
zlabel('z')

[hHead hPlate]=Cone( org,r,h2,dir,n,'closed' );
set(hTail,'EdgeAlpha',0)

author: Sable

repo: mcbench-benchmarks


text:
function cr = slcorrectrate_blks(scores, blocks, clabels, qlabels, op, varargin)
%SLCORRECTRATE_BLKS Computes the correct rate based on blockwise scores
%
% $ Syntax $
%   - cr = slcorrectrate_blks(scores, blocks, clabels, qlabels, op, ...)
%
% $ Arguments $
%   - scores:       the score matrix
%   - blocks:       the cell array of block limits
%   - clabels:      the class labels of reference samples
%   - qlabels:      the query labels of query samples
%   - op:           the score attribute
%   - cr:           the computed classification correct rate
%
% $ Remarks $
%   - An extension of slcorrectrate to support blockwise scores
%
% $ History $
%   - Created by Dahua Lin, on Aug 9th, 2006
%   - Modified by Dahua Lin on Aug 16th, 2006
%     - Based on new slclassify to support multiple schemes
%

%% parse and verify input arguments

if nargin < 5
    raise_lackinput('slcorrectrate_blks', 5);
end

%% Make decision

n = length(qlabels);
decisions = slclassify_blks(scores, n, blocks, clabels, op, varargin{:});

%% Evaluate correct rate

qlabels = qlabels(:)';
cr = sum(decisions == qlabels) / length(qlabels);

author: lmthang


text:
function value = he_triple_product_integral ( i, j, k )

%*****************************************************************************80
%
%% HE_TRIPLE_PRODUCT_INTEGRAL: integral of He(i,x)*He(j,x)*He(k,x)*e^(-x^2/2).
%
%  Discussion:
%
%    VALUE = integral ( -oo < x < +oo ) He(i,x)*He(j,x)*He(k,x) exp(-x^2/2) dx
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 March 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Dongbin Xiu,
%    Numerical Methods for Stochastic Computations: A Spectral Method Approach,
%    Princeton, 2010,
%    ISBN13: 978-0-691-14212-8,
%    LC: QA274.23.X58.
%
%  Parameters:
%
%    Input, integer I, J, K, the polynomial indices.
%
%    Output, real VALUE, the value of the integral.
%
  s = floor ( ( i + j + k ) / 2 );

  if ( s < i || s < j || s < k )
    value = 0.0;
  elseif ( mod ( i + j + k, 2 ) ~= 0 )
    value = 0.0;
  else
    value = r8_factorial ( i ) / r8_factorial ( s - i ) ...
          * r8_factorial ( j ) / r8_factorial ( s - j ) ...
          * r8_factorial ( k ) / r8_factorial ( s - k );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [expected_overlaps, evaluated_lengths, practical_difference] = estimate_expected_overlap(tracker, experiment, sequences, varargin)
% estimate_expected_overlap Estimates expected average overlap for
% different sequence lengths
%
% This function estimates the expected average overlap for sequence of a
% given lengths based on the data
%
% Input:
% - tracker (struct): A valid tracker descriptor.
% - experiment (struct): A valid experiment descriptor.
% - sequences (cell): An array of valid sequence descriptors.
% - varargin[Lengths] (vector): A vector of sequence lengths for which the
% overlap should be evaluated.
% - varargin[Weights] (vector): A vector of per-sequence weigths that indicate
% how much does each sequence contributes to the estimate.
% - varargin[Tags] (cell): A set of tags for which to perform
% calculation. If not set then only 'all' is used.
%
% Output:
% - expected_overlaps (vector): Expected overlaps for corresponding
% lengths.
% - evaluated_lengths (vector): A filtered array of lengths (removed duplicates).
% - practical_difference (vector): An estimate of the practical difference for
% corresponding expected overlap.

lengths = [];
weights = ones(numel(sequences), 1);
tags = {'all'};

for j=1:2:length(varargin)
    switch lower(varargin{j})
        case 'lengths', lengths = varargin{j+1};
        case 'weights', weights = varargin{j+1};
        case 'tags', tags = varargin{j+1};
        otherwise, error(['unrecognized argument ' varargin{j}]);
    end
end

context.failures = {};
context.overlaps = {};
context.practical = {};
context.sources = [];
context = iterate(experiment, tracker, sequences, 'iterator', @collect_segments, 'context', context);
failures = context.failures;
segments = context.overlaps;
practical = context.practical;
occurences = hist(context.sources, max(context.sources));
% if isempty(context.practical)
%   occurences = 0;
% else
%   occurences = hist(context.sources, max(context.sources));
% end

if isempty(lengths)
    maxlen = max(cellfun(@(x) numel(x), segments, 'UniformOutput', true));
    lengths = 1:maxlen;
end

% sort and remove duplicates
lengths = unique(lengths);

skipping = experiment.parameters.skip_initialize;

fragments_count = sum(cellfun(@(x) numel(x) + 1, failures, 'UniformOutput', true));
fragments_length = max(lengths);

tag_count = numel(tags);

if isempty(segments)
    expected_overlaps = zeros(0, tag_count);
    practical_difference = zeros(0, tag_count);
    evaluated_lengths = [];
    return;
end

expected_overlaps = zeros(numel(lengths), tag_count);
practical_difference = zeros(numel(lengths), tag_count);

for l = 1:tag_count

    sequence_weights = weights(context.sources(:));
    frequency = occurences(context.sources(:));
    sequence_weights = sequence_weights(:) ./ frequency(:);

    tag = tags{l};

    fragments = nan(fragments_count, fragments_length);
    fpractical = nan(fragments_count, fragments_length);
    fweights = nan(fragments_count, 1);
    f = 1;
    for i = 1:numel(segments)
        % calculate number of failures and their positions in the trajectory
        F = numel(failures{i});
        if F > 0
            % add first part of the trajectory to the fragment list
            points = failures{i}' + skipping;
            points = [1, points(points <= numel(segments{i}))];

            for j = 1:numel(points)-1;
                o = segments{i}(points(j):points(j+1)); o(isnan(o)) = 0;
                fragments(f, :) = 0;
                fragments(f, 1:min(numel(o), fragments_length)) = o;

                o = practical{i}(points(j):points(j+1)); o(isnan(o)) = 0;
                fpractical(f, :) = 0;
                fpractical(f, 1:min(numel(o), fragments_length)) = o;

                w = numel(sequence_query_tag(sequences{context.sources(i)}, tag, points(j):(points(j+1)))) ...
                    / (points(j+1) - points(j) + 1);

                fweights(f) = sequence_weights(i) * w;

                f = f + 1;
            end;

            o = segments{i}(points(end):end); o(isnan(o)) = 0;
            fragments(f, 1:min(numel(o), fragments_length)) = o;
            o = practical{i}(points(end):end); o(isnan(o)) = 0;
            fpractical(f, 1:min(numel(o), fragments_length)) = o;

            w = numel(sequence_query_tag(sequences{context.sources(i)}, tag, points(end):length(segments{i}))) ...
                / (sequences{context.sources(i)}.length - points(end) + 1);

            fweights(f) = sequence_weights(i) * w;

            f = f + 1;
        else
        % process also last part of the trajectory - segment without failure
            if numel(segments{i}) >= fragments_length
                % tracker did not fail on this sequence and it is longer than
                % observed interval
                fragments(f, :) = segments{i}(1:fragments_length);
                fpractical(f, :) = practical{i}(1:fragments_length);

                w = numel(sequence_query_tag(sequences{context.sources(i)}, tag, 1:fragments_length)) ...
                    / fragments_length;
            else
                fragments(f, 1:numel(segments{i})) = segments{i};
                fpractical(f, 1:numel(practical{i})) = practical{i};

                w = numel(sequence_query_tag(sequences{context.sources(i)}, tag)) ...
                    / sequences{context.sources(i)}.length;
            end

            fweights(f) = sequence_weights(i) * w;
            f = f + 1;
        end
    end

    for e = 1:size(expected_overlaps, 1)
        len = lengths(e);
        % do not calculate for Ns == 1: overlap on first frame is always NaN
        if len == 1
            expected_overlaps(e, l) = 1;
            continue;
        end

        usable = ~isnan(fragments(:, len));

        if ~any(usable)
            continue;
        end;

        % for each len get a single number - average overlap
        expected_overlaps(e, l) = sum(mean(fragments(usable, 2:len), 2) .* fweights(usable)) ./ sum(fweights(usable));
        practical_difference(e, l) = sum(mean(fpractical(usable, 2:len), 2) .* fweights(usable)) ./ sum(fweights(usable));

    end

end

evaluated_lengths = lengths;

end

function context = collect_segments(event, context)

switch (event.type)

    case 'sequence_enter'

        sequence = event.sequence;

        sequence_directory = fullfile(event.tracker.directory, event.experiment.name, ...
            sequence.name);

        switch event.experiment.type
            case {'supervised', 'realtime'}

                for i = 1:event.experiment.parameters.repetitions

                    result_file = fullfile(sequence_directory, sprintf('%s_%03d.txt', event.sequence.name, i));

                    if i == 4 && is_deterministic(sequence, 3, sequence_directory)
                        break;
                    end;

                    if ~exist(result_file, 'file')
                        continue;
                    end;

                    trajectory = read_trajectory(result_file);

                    [~, frames] = estimate_accuracy(trajectory, sequence);

                    [~, failures] = estimate_failures(trajectory, sequence);

                    practical = sequence_get_frame_value(sequence, 'practical');

                    context.failures{end+1} = failures(failures <= sequence.length);
                    context.overlaps{end+1} = frames;
                    context.sources(end+1) = event.sequence_index;

                    if isempty(practical)
                        context.practical{end+1} = zeros(sequence.length, 1);
                    else
                        context.practical{end+1} = practical;
                    end

                end;

            otherwise, error(['unsupported type ', event.experiment.type]);
        end

end;

end

author: votchallenge

repo: toolkit-legacy


text:
function [im,stamp] = fmfreadframe(fid,loc,nr,nc,bytes_per_chunk,data_format)

if fseek(fid,loc,'bof') == 0,
  [im,stamp] = fmf_read_frame(fid,nr,nc,bytes_per_chunk,data_format);
else
  im = -1;
  stamp = -1;
end

author: kristinbranson

repo: JAABA


text:
% eeg_load_difdata_script - Script to load ERP difference waves
% Load data (680 rows x 124 cols)
eegpath = 'D:\data_emse\ptsdpet\link14hz\'; cd(eegpath);
ext = '_link14hz';
groups = {'c' 'p'};
cond = {'tac_oac'};
subs = [1:10];

% Define sample rate & epoch parameters
sample_rate = 2.5;
epoch_start = -200;
epoch_end = 1500;
points = 681;
timeArray = meshgrid(epoch_start:sample_rate:epoch_end,1)';
timeNonZero = find(timeArray);
timeZero = find(timeArray == 0);

for g = 1:length(groups),
    for i=1:length(subs),
        for c=1:length(cond),

            % Load data
            file = sprintf('%s%s%02d%s%s.txt',eegpath,groups{g},subs(i),cond{c},ext);
            data = sprintf(  '%s%02d%s%s',            groups{g},subs(i),cond{c},ext);
            if ~exist(data,'var')
                fprintf('loading %s\n',file);
                load(file);

                % Interpolate the Zero value, only useful for ascii from
                % Neuroscan 3x or 4x AVG files
                volt = eval(data);
                InterpZero = interp1( timeArray(timeNonZero), volt, 0, 'cubic' );
                volt = [volt(1:timeZero-1,:); InterpZero; volt(timeZero:end,:)];
                eval(strcat(data,' = volt;'));

            else
                fprintf('Already loaded %s\n',file);
                % Save translated data file

                %dat = eval(data); dat = dat';
                %file = sprintf('%s%s%02d%s%s.dat',eegpath,groups{g},subs(i),cond{c},ext);
                %eeg_write_ascii(file,dat,'%12.4f');
            end
        end
    end
end

clear g i c file data timeNonZero volt InterpZero

author: PatternRecognition

repo: OpenBMI


text:
function [C, C_index, C_dist] = DC(S_prefix, S_number, S_profix, r, k)
%
% Discord Candidate
% Authors: Dragomir Yankov, Eamonn Keogh and Umaa Rebbapragada
% Original paper :
% Disk Aware Discord Discovery: Finding Unusual Time Series
% in Terabyte Sized Datasets
% Website : http://www.cs.ucr.edu/~eamonn/selected_publications.htm
% e-mail : {dyankov,eamonn}@cs.ucr.edu, urebbapr@cs.tufts.edu
% Programmer: Wei-Chih Lai
%
% Inputs
%   S_prefix: a string as prefix for path of files
%   S_number: 1x1 integer means number of files
%   S_profix: a string as profix for path of files
%   r: 1x1 integer means threshold
%   k: 1x1 integer means number of top outliers
%
% Outputs
%   C: kxn vector
%      outlier instances
%   C_index: kx1 vector
%      index of outliers in origin datasets
%   C_dist: kx1 vector
%      distance of outliers
%
    assert(k > 0);
    [C, C_index] = DC_Selection(S_prefix, S_number, S_profix, r);
    [C, C_dist] = DC_Refinement(S_prefix, S_number, S_profix, C, C_index, r);
    [~, index] = sort(C_dist, 'descend');
    index = index(1:min(size(C_dist, 1), k), 1);
    C = C(index, :);
    C_index = C_index(index, :);
    C_dist = C_dist(index, :);

author: dsmi-lab-ntust

repo: AnomalyDetectionToolbox


text:
%  Author: Relja Arandjelovic (relja@relja.info)

% This file contains a few examples of how to test our NetVLAD on standard object/image retrieval benchmarks, refer to README.md for setup instructions, and to our project page for all relevant information (e.g. our paper): http://www.di.ens.fr/willow/research/netvlad/

%  The code samples use the GPU by default, if you want to use the CPU instead (quite slow!), add `'useGPU', false` to the affected function call (serialAllFeats)

% Set the MATLAB paths
setup;

% ---------- Use/test our networks

% Load our network
netID= 'vd16_pitts30k_conv5_3_vlad_preL2_intra_white';
paths= localPaths();
load( sprintf('%s%s.mat', paths.ourCNNs, netID), 'net' );
net= relja_simplenn_tidy(net);

% --- Oxford Buildings / Paris Buildings

dsetName= 'ox5k'; % change to `'paris'` for Paris Buildings

% Define if you want to use the query ROI information.
% The original evaluation procedure for Oxford / Paris DOES use ROIs,
% so we really recommend to keep this as is, but we also include the setting
% where they are not used in order to be able to compare to some recent works.
% Please continue using ROIs as this makes it easier to compare numbers
% accross different papers, and it is the proper evaluation procedure for
% these datasets.
% Note that, for ease of implementation, if you are using ROIs the code will
% crop the query images and save them into a new folder
% `images_crop_<amount_of_crop>` inside the dataset root
% (i.e. makes sure this is writable). The crop is made by extending the ROI
% by half of the receptive field size, as this ensures that the centres of
% all 'features' % (i.e. conv responses) are inside the ROI; this is a setting
% similar to the original test procedure where a SIFT feature is kept iff its
% centre is inside the ROI. Note that this is ROI extension is only valid
% for single-scale extraction of features, if you want multiple scales then
% the implementation needs to be different
%
% To use the ROI: specify the receptive field size (can be 0 for a tight crop)
% Not to use the ROI (not recommended): specify a negative number

useROI= true;
if useROI
    lastConvLayer= find(ismember(relja_layerTypes(net), 'custom'),1)-1; % Relies on the fact that NetVLAD, Max and Avg pooling are implemented as a custom layer and are the first custom layer in the network. Change if you use another network which has other custom layers before
    netBottom= net;
    netBottom.layers= netBottom.layers(1:lastConvLayer);
    info= vl_simplenn_display(netBottom);
    clear netBottom;
    recFieldSize= info.receptiveFieldSize(:, end);
    assert(recFieldSize(1) == recFieldSize(2));% we are assuming square receptive fields, otherwise dbVGG needs to change to account for non-square
    recFieldSize= recFieldSize(1);
    strMode= 'crop';
else
    recFieldSize= -1;
    strMode= 'full';
end

dbTest= dbVGG(dsetName, recFieldSize);

% Compute db image representations (images have different resolutions so batchSize is constrained to 1)
dbFeatFn= sprintf('%s%s_%s_db.bin', paths.outPrefix, netID, dbTest.name);
serialAllFeats(net, dbTest.dbPath, dbTest.dbImageFns, dbFeatFn, 'batchSize', 1);

if useROI
    qFeatFn = sprintf('%s%s_%s_q_%d.bin', paths.outPrefix, netID, dbTest.name, recFieldSize);
    serialAllFeats(net, dbTest.qPath, dbTest.qImageFns, qFeatFn, 'batchSize', 1);
end

% Load the image representations
dbFeat= fread( fopen(dbFeatFn, 'rb'), inf, 'float32=>single');
dbFeat= reshape(dbFeat, [], dbTest.numImages);
nDims= size(dbFeat, 1);
if useROI
    qFeat= fread( fopen(qFeatFn, 'rb'), [nDims, dbTest.numImages], 'float32=>single');
    assert(size(qFeat,2)==dbTest.numQueries);
else
    qFeat= dbFeat(:, dbTest.queryIDs);
end

% Measure recall@N
mAP= relja_retrievalMAP(dbTest, struct('db', dbFeat, 'qs', qFeat), true);
relja_display('Performance on %s (%s), %d-D: %.2f', dbTest.name, strMode, size(dbFeat,1), mAP*100);

% Try 256-D
D= 256;
mAPsmall= relja_retrievalMAP(dbTest, struct( ...
    'db', relja_l2normalize_col(dbFeat(1:D,:)), ...
    'qs', relja_l2normalize_col(qFeat(1:D,:)) ...
    ), true);
relja_display('Performance on %s (%s), %d-D: %.2f', dbTest.name, strMode, D, mAPsmall*100);

% --- Holidays

useRotated= false; % for rotated Holidays
dbTest= dbHolidays(useRotated); % it will automatically downscale images to (1024x768) pixels, as per the original testing procedure (see the Holidays website)

% Set the output filename for the database (query images are a subset)
dbFeatFn= sprintf('%s%s_%s_db.bin', paths.outPrefix, netID, dbTest.name);

% Compute db image representations (images have different resolutions so batchSize is constrained to 1)
serialAllFeats(net, dbTest.dbPath, dbTest.dbImageFns, dbFeatFn, 'batchSize', 1);

% Load the image representations
dbFeat= fread( fopen(dbFeatFn, 'rb'), inf, 'float32=>single');
dbFeat= reshape(dbFeat, [], dbTest.numImages);
nDims= size(dbFeat, 1);
qFeat= dbFeat(:, dbTest.queryIDs);

mAP= relja_retrievalMAP(dbTest, struct('db', dbFeat, 'qs', qFeat), true);
relja_display('Performance on %s, %d-D: %.2f', dbTest.name, size(dbFeat,1), mAP*100);

% Try 256-D
D= 256;
mAPsmall= relja_retrievalMAP(dbTest, struct( ...
    'db', relja_l2normalize_col(dbFeat(1:D,:)), ...
    'qs', relja_l2normalize_col(qFeat(1:D,:)) ...
    ), true);
relja_display('Performance on %s, %d-D: %.2f', dbTest.name, D, mAPsmall*100);

author: Relja

repo: netvlad

save_path: github-repos/MATLAB/Relja-netvlad


text:
function [dimMismatch, isTranspose] = VerifyTiedLayers(tiedLayers)

for j=1:length(tiedLayers)
    dim(:,j) = tiedLayers{j}.dim;
end
dimMismatch = [];
isTranspose = [];
for j=2:length(tiedLayers)
    if sum(abs(dim(:,j)-dim(:,1)))
        if sum(abs(dim(end:-1:1,j)-dim(:,1)))
            dimMismatch(j) = 1;
            fprintf('Error: dimension mismatch between layers that share the same weight matrix\n');
        else
            isTranspose(j) = 1;
        end
    else
        isTranspose(j) = 0;
    end
end

end

author: singaxiong

repo: SignalGraph


text:
function Xhat = gsp_jft(G,X)
%GSP_JFT Compute the Joint time-vertex Fourier Transform
%   Usage:  Xhat = gsp_jft(G,X)
%
%   Input parameters:
%         G          : Time-Vertex graph structure
%         X          : Time-Vertex signal
%   Output parameters:
%         Xhat       : Joint Time-Vertex Fourier Transform of X
%
%
%   'gsp_jft(G,X)' computes the joint time-vertex Fourier transform of the time-vertex
%    signal $X$ with respect to the Fourier basis of the graph G: U_G and the DFT basis U_T.
%
%   .. X_hat = U_G' * X * conj(U_T)
%
%
%   To compute the Fourier basis of a graph G, you can use the function::
%
%           G = gsp_compute_fourier_basis(G);
%
%   Example:::
%
%           N = 50; T=200;
%           G = gsp_sensor(N);
%           G = gsp_jtv_graph(G,T);
%           G = gsp_compute_fourier_basis(G);
%           X = sin((1:N)'*(1:T)*pi/(4*N));
%           Xhat = gsp_jft(G,X);
%           imagesc(fftshift(abs(Xhat),2));
%
%   See also: gsp_ijft, gsp_compute_fourier_basis
%

% Author : Francesco Grassi
% Date   : September 2016

if isempty(G.jtv.NFFT)
    NFFT = size(X,2);
else
    NFFT = G.jtv.NFFT;
end

normalize = 1/sqrt(NFFT);

switch G.jtv.transform
    case 'dft'
        Xhat = fft(gsp_gft(G,X),NFFT,2)*normalize;
    case 'dct'
        if size(X,3)>1
            Xhat = ipermute(dct(permute(gsp_gft(G,X),[2 1 3]),NFFT),[2 1 3]);
        else
            Xhat = dct(gsp_gft(G,X).',NFFT).';
        end

    otherwise
        error('Unknown transform');
end

end

repo: gspbox


text:
function ompspeedtest
%OMPSPEEDTEST Test the speed of the OMP functions.
%  OMPSPEEDTEST invokes the three operation modes of OMP and compares
%  their speeds. The function automatically selects the number of signals
%  for the test based on the speed of the system.
%
%  To run the test, type OMPSPEEDTEST from the Matlab prompt.
%
%  See also OMPDEMO.

%  Ron Rubinstein
%  Computer Science Department
%  Technion, Haifa 32000 Israel
%  ronrubin@cs
%
%  August 2009

% random dictionary %

n = 512;
L = 1000;
T = 8;

D = randn(n,L);
D = D*diag(1./sqrt(sum(D.*D)));    % normalize the dictionary

% select signal number according to computer speed %

x = randn(n,20);
tic; omp(D,x,[],T,'messages',-1); t=toc;
signum = ceil(20/(t/20));     % approximately 20 seconds of OMP-Cholesky

% generate random signals %

X = randn(n,signum);

% run OMP  %

printf('\nRunning OMP-Cholesky...');
tic; omp(D,X,[],T,'messages',4); t1=toc;

printf('\nRunning Batch-OMP...');
tic; omp(D,X,D'*D,T,'messages',1); t2=toc;

printf('\nRunning Batch-OMP with D''*X specified...');
tic; omp(D'*X,D'*D,T,'messages',1); t3=toc;

% display summary  %

printf('\n\nSpeed summary for %d signals, dictionary size %d x %d:\n', signum, n, L);
printf('Call syntax        Algorithm               Total time');
printf('--------------------------------------------------------');
printf('OMP(D,X,[],T)      OMP-Cholesky            %5.2f seconds', t1);
printf('OMP(D,X,G,T)       Batch-OMP               %5.2f seconds', t2);
printf('OMP(DtX,G,T)       Batch-OMP with D''*X     %5.2f seconds\n', t3);

author: tiepvupsu

repo: DICTOL

save_path: github-repos/MATLAB/tiepvupsu-DICTOL


text:
function SupervisedClassification_Callback(hObject, eventdata, handles)

% This function uses a convolutional neural network, trained in
% "TrainSupervisedClassifier_Callback.m", to classify USVs.

[FileName,PathName] = uigetfile(fullfile(handles.data.squeakfolder,'Clustering Models','*.mat'),'Select Network');
net = load([PathName FileName],'ClassifyNet','wind','noverlap','nfft','imageSize');

if exist('net','var') ~= 1
    errordlg('Network not be found. Is this file a trained CNN?')
    return
end

if exist(handles.data.settings.detectionfolder,'dir')==0
    errordlg('Please Select Detection Folder')
    uiwait
    load_detectionFolder_Callback(hObject, eventdata, handles)
    handles = guidata(hObject);  % Get newest version of handles
end

options.imageSize = [128, 128, 1];
[ClusteringData, Class, options.freqRange, options.maxDuration, options.spectrogram] = CreateClusteringData(handles, 'scale_duration', true, 'fixed_frequency', true);

% Resize the images to match the input image size
images = zeros([options.imageSize, size(ClusteringData, 1)]);
for i = 1:size(ClusteringData, 1)
    images(:,:,:,i) = imresize(ClusteringData.Spectrogram{i}, options.imageSize(1:2));
end
% wind=options.spectrogram.windowsize;
% noverlap=options.spectrogram.overlap;
% nfft=options.spectrogram.nfft;
imageSize=options.imageSize;

h = waitbar(0,'Initializing');
for j=1:height(ClusteringData);
    waitbar(j/height(ClusteringData), h, ['Classifying Image ' num2str(j) ' of ' num2str(height(ClusteringData))]);
    X = images(:,:,:,j) ./ 256;
    [Cl, sc] = classify(net.ClassifyNet, X);
    clustAssign(j,1)=Cl;
end
close(h);
clusterName=unique(clustAssign);
saveChoice =  questdlg('Update files with new classifications?','Save Classification','Yes','No','Yes');
switch saveChoice
    case 'Yes'
        UpdateCluster(ClusteringData, clustAssign, clusterName, zeros(1,height(ClusteringData)));
        update_folders(hObject, eventdata, handles);
        if isfield(handles,'current_detection_file')
            loadcalls_Callback(hObject, eventdata, handles, true)
        end
    case 'No'
        return
end
end

author: DrCoffey

repo: DeepSqueak

save_path: github-repos/MATLAB/DrCoffey-DeepSqueak


text:
function [Ff, res] = loadPROC(root, mname, datExp, npl, plane_um, pix_um, f0)

root    = fullfile(root, mname,datExp);
fdir    = dir(root);
blockstring = fdir(3).name;
root    = fullfile(root,blockstring);

clear Fc

ik = 0;
med0 = [];
igpl = [1:npl];
statall = [];
redcell = [];

for iplane = igpl
    fname = sprintf('F_%s_%s_plane%d_Nk*_proc.mat', mname, datExp, iplane);

    fs = dir((fullfile(root, fname)));
    fname = fs(1).name;
    load(fullfile(root, fname))

    ds = (ceil((npl+1)/2) - iplane)/npl;
    %
    Fcell = dat.F.Fcell;

    [NN NT] = size(Fcell{1});

    iscell  = logical(dat.cl.iscell);
    statall = cat(1, statall, dat.stat(iscell)');
    redcell = cat(1, redcell, dat.cl.redcell(iscell));

    for j = 1:length(Fcell)
        F = Fcell{j}(iscell, :);
        F = register_F(F', ds)';
        Fc{iplane,j} = F;
    end
    %
    icells = find(iscell);
    for ip = 1:length(icells)
        ik = ik+1;
        med0(ik,:) = [pix_um*dat.stat(icells(ip)).med plane_um * iplane];
    end
end
%%

Ls = cellfun(@(x) size(x,2), Fc);
Lmax = min(Ls(igpl, :), [], 1);
%
Ff  = [];
for iplane = igpl
    F0 = [];
    for j = [1:size(Fc,2)]
        F0 = cat(2, F0, Fc{iplane, j}(:, 1:Lmax(j)));
    end
    Ff = cat(1, Ff, F0);
end
Ff = Ff';

med = med0;
if npl>=9
    goodcell = remove_doubles2(Ff, med);
    %
    Ff      = Ff(:, goodcell);
    med     = med(goodcell, :);
    statall = statall(goodcell);
    redcell = redcell(goodcell);
end
%
nCum = [0 cumsum(Lmax)];
%
[NT NN] = size(Ff);

Ff = single(Ff);

%%
tic
[NT NN] = size(Ff);
parfor i = 1:NN
    Ff(:,i) = Ff(:,i) - ordfilt2(Ff(:,i), 31 * ceil(f0/3), true(201 * ceil(f0/3),1), 'symmetric');
    Ff(:,i) = Ff(:,i) - median(Ff(:,i));
end

res.med     = med;
res.statall = statall;
res.redcell = redcell;
%%

res.nCum = nCum;

author: cortex-lab


text:
function obj = updateTernaryColorbar(obj,colorbarIndex)

% title: ...[DONE]
% titleside: ...[DONE]
% titlefont: ...[DONE]
% thickness: ...[DONE]
% thicknessmode: ...[DONE]
% len: ...[DONE]
% lenmode: ...[DONE]
% x: ...[DONE]
% y: ...[DONE]
% autotick: ...[DONE]
% nticks: ...[DONE]
% ticks: ...[DONE]
% showticklabels: ...[DONE]
% tick0: ...[DONE]
% dtick: ...[DONE]
% ticklen: ...[DONE]
% tickwidth: ...[DONE]
% tickcolor: ...[DONE]
% tickangle: ...[NOT SUPPORTED IN MATLAB]
% tickfont: ...[DONE]
% exponentformat: ...[DONE]
% showexponent: ...[NOT SUPPORTED IN MATLAB]
% xanchor: ...[DONE]
% yanchor: ...[DONE]
% bgcolor: ...[DONE]
% outlinecolor: ...[DONE]
% outlinewidth: ...[DONE]
% borderwidth: ...[NOT SUPPORTED IN MATLAB]
% bordercolor: ...[NOT SUPPORTED IN MATLAB]
% xpad: ...[DONE]
% ypad: ...[DONE]

%-FIGURE STRUCTURE-%
figureData = get(obj.State.Figure.Handle);

%-PLOT DATA STRUCTURE- %
try
    colorbarData = get(obj.State.Colorbar(colorbarIndex).Handle);
catch
    disp('could not extract ColorBar data');
end

%-STANDARDIZE UNITS-%
colorbarunits = colorbarData.Units;
set(obj.State.Colorbar(colorbarIndex).Handle,'Units','normalized');

%-------------------------------------------------------------------------%

%-colorbar position-%
colorbar.xanchor = 'left';
colorbar.yanchor = 'bottom';
colorbar.x = colorbarData.Position(1)*1.025;
colorbar.y = colorbarData.Position(2);

%-------------------------------------------------------------------------%

%-exponent format-%
colorbar.exponentformat = obj.PlotlyDefaults.ExponentFormat;

%-------------------------------------------------------------------------%

% get colorbar title and labels
colorbarTitle = colorbarData.Label;

if isHG2
    colorbarTitleData = get(colorbarTitle);
    colorbarYLabel = colorbarTitle;
    colorbarYLabelData = get(colorbarTitle);
    colorbarXLabelData.String = [];
else
    colorbarTitleData = get(colorbarTitle);
    colorbarXLabel = colorbarData.XLabel;
    colorbarXLabelData = get(colorbarXLabel);
    colorbarYLabel = colorbarData.YLabel;
    colorbarYLabelData = get(colorbarYLabel);
end

%-colorbar title-%
if ~isempty(colorbarTitleData.String)
    colorbar.title = parseString(colorbarTitleData.String,colorbarTitleData.Interpreter);

elseif ~isempty(colorbarXLabelData.String)
    colorbar.title = parseString(colorbarXLabelData.String,colorbarXLabelData.Interpreter);

elseif ~isempty(colorbarYLabelData.String)
    colorbar.title = parseString(colorbarYLabelData.String,colorbarYLabelData.Interpreter);
end

%-------------------------------------------------------------------------%

%-STANDARDIZE UNITS-%
titleUnits = colorbarTitleData.Units;
titleFontUnits = colorbarTitleData.FontUnits;
yLabelUnits = colorbarYLabelData.Units;
yLabelFontUnits = colorbarYLabelData.FontUnits;
set(colorbarTitle,'Units','data');
set(colorbarYLabel,'Units','data');
set(colorbarYLabel,'FontUnits','points');

if ~isHG2
    xLabelUnits = colorbarXLabelData.Units;
    xLabelFontUnits = colorbarXLabelData.FontUnits;
    set(colorbarTitle,'FontUnits','points');
    set(colorbarXLabel,'Units','data');
    set(colorbarXLabel,'FontUnits','points');
end

if ~isempty(colorbarTitleData.String)
    if colorbarTitleData.Rotation == 90
        colorbar.titleside = 'right';
    else
        colorbar.titleside = 'top';
    end

    colorbar.titlefont.family = matlab2plotlyfont(colorbarTitleData.FontName);
    col = 255*colorbarTitleData.Color;
    colorbar.titlefont.color = sprintf('rgb(%f,%f,%f)', col);
    colorbar.titlefont.size = 1.20 * colorbarTitleData.FontSize;

elseif ~isempty(colorbarXLabelData.String)
    colorbar.titleside = 'right';
    colorbar.titlefont.family = matlab2plotlyfont(colorbarXLabelData.FontName);
    col = 255*colorbarXLabelData.Color;
    colorbar.titlefont.color = sprintf('rgb(%f,%f,%f)', col);
    colorbar.titlefont.size = 1.20 * colorbarXLabelData.FontSize;

elseif ~isempty(colorbarYLabelData.String)
    colorbar.titleside = 'bottom';
    colorbar.titlefont.family = matlab2plotlyfont(colorbarYLabelData.FontName);
    col = 255*colorbarYLabelData.Color;
    colorbar.titlefont.color = sprintf('rgb(%f,%f,%f)', col);
    colorbar.titlefont.size = 1.20 * colorbarYLabelData.FontSize;
end

%-REVERT UNITS-%
set(colorbarTitle,'Units',titleUnits);
set(colorbarTitle,'FontUnits',titleFontUnits);
set(colorbarYLabel,'Units',yLabelUnits);
set(colorbarYLabel,'FontUnits',yLabelFontUnits);

if ~isHG2
    set(colorbarXLabel,'Units',xLabelUnits);
    set(colorbarXLabel,'FontUnits',xLabelFontUnits);
end

%-------------------------------------------------------------------------%

%-some colorbar settings-%

lineWidth = colorbarData.LineWidth*obj.PlotlyDefaults.AxisLineIncreaseFactor;
tickLength = min(obj.PlotlyDefaults.MaxTickLength,...
    max(colorbarData.TickLength(1)*colorbarData.Position(3)*obj.layout.width,...
    colorbarData.TickLength(1)*colorbarData.Position(4)*obj.layout.height));

colorbar.thicknessmode = 'fraction';
colorbar.thickness = colorbarData.Position(3);
colorbar.tickwidth = lineWidth;
colorbar.ticklen = tickLength;

colorbar.lenmode = 'fraction';
colorbar.len = colorbarData.Position(4)*1.025;
colorbar.outlinewidth = lineWidth;

%-------------------------------------------------------------------------%

% orientation vertical check
orientVert = colorbar.len > colorbar.thickness;

%-------------------------------------------------------------------------%

%-coloration-%

if isHG2
    col = 255*colorbarData.Color;
else
    if orientVert
        col = 255*colorbarData.YColor;
    else
        col = 255*colorbarData.XColor;
    end
end

colorbarColor = sprintf('rgb(%f,%f,%f)', col);

colorbar.outlinecolor = colorbarColor;
colorbar.tickcolor = colorbarColor;
colorbar.tickfont.color = colorbarColor;

%-------------------------------------------------------------------------%

%-axis tickfont-%
colorbar.tickfont.size = colorbarData.FontSize;
colorbar.tickfont.family = matlab2plotlyfont(colorbarData.FontName);

%-------------------------------------------------------------------------%

%-colorbar pad-%
colorbar.xpad = obj.PlotlyDefaults.MarginPad;
colorbar.ypad = obj.PlotlyDefaults.MarginPad;

%-------------------------------------------------------------------------%

%-set ticklabels-%
nticks = length(colorbarData.Ticks);

if isHG2
    if isempty(colorbarData.Ticks)
        %-hide tick labels-%
        colorbar.ticks = '';
        colorbar.showticklabels = false;

    else

        %-tick direction-%
        switch colorbarData.TickDirection
            case 'in'
                colorbar.ticks = 'inside';
            case 'out'
                colorbar.ticks = 'outside';
        end

        if strcmp(colorbarData.TickLabelsMode,'auto')
            colorbar.autotick = true;
            colorbar.nticks = length(colorbarData.Ticks) + 1; %nticks = max ticks (so + 1)

        else

            %-show tick labels-%
            if isempty(colorbarData.TickLabels)
                colorbar.showticklabels = false;

            else
                colorbar.autotick = false;
                colorbar.tickvals = colorbarData.Ticks;
                colorbar.ticktext = colorbarData.TickLabels;
            end
        end
    end

else
    if orientVert
        if isempty(colorbarData.YTick)
            %-show tick labels-%
            colorbar.ticks = '';
            colorbar.showticklabels = false;
        else
            %-tick direction-%
            switch colorbarData.TickDir
                case 'in'
                    colorbar.ticks = 'inside';
                case 'out'
                    colorbar.ticks = 'outside';
            end

            if strcmp(colorbarData.YTickLabelMode,'auto')
                %-autotick-%
                colorbar.autotick = true;
                %-numticks-%
                colorbar.nticks = length(colorbarData.YTick) + 1; %nticks = max ticks (so + 1)
            else
                %-show tick labels-%
                if isempty(colorbarData.YTickLabel)
                    colorbar.showticklabels = false;
                else
                    %-autotick-%
                    colorbar.autotick = false;
                    %-tick0-%
                    colorbar.tick0 = str2double(colorbarData.YTickLabel(1,:));
                    %-dtick-%
                    colorbar.dtick = str2double(colorbarData.YTickLabel(2,:)) - str2double(colorbarData.YTickLabel(1,:));
                end
            end
        end
    else
        if isempty(colorbarData.XTick)
            %-show tick labels-%
            colorbar.ticks = '';
            colorbar.showticklabels = false;
        else
            %-tick direction-%
            switch colorbarData.TickDir
                case 'in'
                    colorbar.ticks = 'inside';
                case 'out'
                    colorbar.ticks = 'outside';
            end

            if strcmp(colorbarData.XTickLabelMode,'auto')
                %-autotick-%
                colorbar.autotick = true;
                %-numticks-%
                colorbar.nticks = length(colorbarData.XTick) + 1;
            else
                %-show tick labels-%
                if isempty(colorbarData.XTickLabel)
                    colorbar.showticklabels = false;
                else
                    %-autotick-%
                    colorbar.autotick = false;
                    %-tick0-%
                    colorbar.tick0 = str2double(colorbarData.XTickLabel(1,:));
                    %-dtick-%
                    colorbar.dtick = str2double(colorbarData.XTickLabel(2,:)) - str2double(colorbarData.XTickLabel(1,:));
                end
            end
        end
    end
end

%-------------------------------------------------------------------------%

%-colorbar bg-color-%
if ~isHG2

    if ~ischar(colorbarData.Color)
        col = 255*colorbarData.Color;
    else
        col = 255*figureData.Color;
    end

    obj.layout.plot_bgcolor = sprintf('rgb(%f,%f,%f)', col);
end

%-------------------------------------------------------------------------%

%-ASSOCIATED DATA-%
if isfield(colorbarData.UserData,'dataref')
    colorbarDataIndex = colorbarData.UserData.dataref;
else
    colorbarDataIndex = findColorbarData(obj,colorbarIndex);
end

if (nticks ~= 0)
    colorIndex = linspace(0, 1, nticks);
    colorData = linspace(0, 1, nticks-1);
    m = 1;

    for n = 1:nticks-1
        col = 1-colorData(n);
        colorscale{m} = {colorIndex(n), sprintf('rgb(%f,%f,%f)', 255*[col, col, col])};
        colorscale{m+1} = {colorIndex(n+1), sprintf('rgb(%f,%f,%f)', 255*[col, col, col])};
        m = 2*n+1;
    end
    obj.data{colorbarDataIndex}.marker.color = colorbarData.Ticks;
else
    colorscale = {{0, 'rgb(255,255,255)'}, {1, 'rgb(0,0,0)'}}
end

obj.data{colorbarDataIndex}.marker.colorscale = colorscale;
obj.data{colorbarDataIndex}.marker.colorbar = colorbar;
obj.data{colorbarDataIndex}.showscale = true;

%-------------------------------------------------------------------------%

%-REVERT UNITS-%
set(obj.State.Colorbar(colorbarIndex).Handle,'Units',colorbarunits);

end

author: plotly


text:
function [lng] = spm_lg_gamma(p,b)
% Log of generalised gamma function
% FORMAT [lng] = spm_lg_gamma(p,b)
%
% p       - dimension parameter
% b       - degrees of freedom type parameter
%__________________________________________________________________________
%
% References:
% * Bayesian Inference in Statistical Analysis, Box & Tiao, 1992, p. 427.
% * Aspects of Multivariate Statistical Theory, R.J. Muirhead, p. 62.
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Will Penny
% $Id: spm_lg_gamma.m 2784 2009-02-24 19:11:20Z guillaume $

if b <= (p-1)/2
   warning('Parameter out of range');
   lng = NaN;
   return
end

lng = (p*(p-1)/2) * gammaln(0.5);
for alpha = 1:p
  lng = lng + gammaln(b+0.5*(alpha-1));
end

author: spm


text:
function blas2_test06 ( )

%*****************************************************************************80
%
%% BLAS2_TEST06 tests DTRMV.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 April 2014
%
%  Author:
%
%    John Burkardt
%
  m = 5;
  lda = m;
  n = m;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'BLAS2_TEST06\n' );
  fprintf ( 1, '  For a triangular matrix A,\n' );
  fprintf ( 1, '  DTRMV computes y := A * x or y := A'' * x\n' );

  for test = 1 : 2

    uplo = 'U';

    if ( test == 1 )
      trans = 'N';
    else
      trans = 'T';
    end

    diag = 'N';

    for j = 1 : n
      for i = 1 : j
        a(i,j) = i + j;
      end
      for i = j + 1 : m
        a(i,j) = 0.0;
      end
    end

    incx = 1;
    for i = 1 : n
      x(i) = i;
    end

    x = dtrmv ( uplo, trans, diag, n, a, lda, x, incx );

    if ( trans == 'N' )
      r8vec_print ( n, x, '  Result y = A * x' );
    else
      r8vec_print ( n, x, '  Result y = A'' * x' );
    end

  end

  return
end


text:
function varargout=sumann(varargin)
%
% report=sumann(recName,annName,stopTime,qrsAnnotationsOnly)
%
%    Wrapper to WFDB SUMANN:
%         http://www.physionet.org/physiotools/wag/sumann-1.htm
%
% Reads a WFDB annotation file and summarize its contents.
%
% Ouput Parameters:
%
% report
%       String with the contaning summary of the contents, including the
%       number of annotations of each type as well the duration and number of
%      episodes of each rhythm and signal quality.
%
%Input Parameters:
% recName
%       String specifying the WFDB record file.
%
% annName
%       String specifying the reference WFDB annotation file.
%
% stopTime (Optional)
%       String specifying the stop time in WFDB format (default is end of
%       record).
%
% qrsAnnotationsOnly (Optional)
%       1x1 Boolean. If true, summarize QRS annotation only (default = 0).
%
%
% Written by Ikaro Silva, 2013
% Last Modified: -
% Version 1.0
% Since 0.9.0
%
% %Example (this will generate a /mitdb/100.qrs file at your directory):
%
% report=sumann('mitdb/100','atr');
%
%
%
% See also RDANN, MXM, WFDBTIME, BXB

%endOfHelp

javaWfdbExec=getWfdbClass('sumann');

%Set default pararamter values
inputs={'recName','annName','stopTime','qrsAnnotationsOnly'};
recName=[];
annName=[];
stopTime=[];
qrsAnnotationsOnly=0;
for n=1:nargin
    if(~isempty(varargin{n}))
        eval([inputs{n} '=varargin{n};']);
    end
end

wfdb_argument={'-r',recName,'-a',annName};

if(~isempty(stopTime))
     wfdb_argument{end+1}='-t';
    wfdb_argument{end+1}=stopTime;
end
if(qrsAnnotationsOnly)
     wfdb_argument{end+1}='-q';
end

report=javaWfdbExec.execToStringList(wfdb_argument);
if(nargout>0)
   varargout{1}=report;
end

author: ikarosilva

repo: wfdb-app-toolbox


text:
% [Cs, Es] = curve_ext_multi(Tx, log2(fs), nc, lambda, clwin)
%
% Consecutively extract maximum energy / minimum curvature curves
% from a synchrosqueezing representation.  As curves are extracted,
% their associated energies are zeroed out in the synsq
% representation.  Curve extraction is then again performed on
% the remaining representaton.
%
% For more details, see help curve_ext.
%
%
%
% Input:
%  Tx, fs, lambda - same as input to curve_ext() (lambda default = 1e3)
%  nc - Number of curves to extract
%  clwin - frequency clearing window; after curve extraction, a window
%    of frequencies [Cs(:,i)-clwin:Cs(:,i)+clwin] is removed from the
%    original representation. (default = 2)
%
% Output:
%  Cs - [N x nc] matrix of curve indices (N==size(Tx,2))
%  Es - [nc x 1] vector of associated (logarithmic) energies
%
%---------------------------------------------------------------------------------
%    Synchrosqueezing Toolbox
%    Authors: Eugene Brevdo (http://www.math.princeton.edu/~ebrevdo/)
%---------------------------------------------------------------------------------
function [Cs, Es] = curve_ext_multi(Tx, lfs, nc, lambda, clwin)
if nargin<5, clwin = 4; end
if nargin<4, lambda = 1e3; end

[na, N] = size(Tx);

Cs = zeros(N, nc);
Es = zeros(nc, 1);

for ni=1:nc
    [Cs(:,ni), Es(ni)] = curve_ext(Tx, lfs, lambda);

    % Remove this curve from the representation
    % Max/min frequencies for each time step
    for cli=0:clwin
        fb = max(1, Cs(:,ni)-cli);
        fe = min(na, Cs(:,ni)+cli);
        Tx([0:N-1]'*na+fb) = sqrt(eps);
        if cli>0, Tx([0:N-1]'*na+fe) = sqrt(eps); end
    end
end

author: ebrevdo

repo: synchrosqueezing


text:
%NEURC Automatic neural network classifier
%
% 	W = NEURC(A,UNITS)
% 	W = A*NEURC([],UNITS)
% 	W = A*NEURC(UNITS)
%
% INPUT
%   A      Dataset
%   UNITS  Number of units
%          Default: 0.2 x size smallest class in A.
%
% OUTPUT
%   W      Trained feed-forward neural network mapping
%
% DESCRIPTION
% Automatically trained feed-forward neural network classifier with UNITS
% units in a single hidden layer. Training, by LMNC, is stopped when the
% performance on an artificially generated tuning set of 1000 samples per
% class (based on k-nearest neighbour interpolation) does not improve anymore.
%
% NEURC always tries three random initialisations, with fixed random seeds,
% and returns the best result according to the tuning set. This is done in
% order to obtain a reproducable result.
%
% If UNITS is NaN it is optimised by REGOPTC. This may take a long
% computing time and is often not significantly better than the default.
%
% Uses the Mathworks' neural network toolbox.
%
% SEE ALSO (PRTools Guide)
% MAPPINGS, DATASETS, LMNC, BPXNC, GENDATK, REGOPTC

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

% $Id: neurc.m,v 1.9 2008/07/03 09:11:44 duin Exp $

function argout = neurc (varargin)

  checktoolbox('nnet');

	mapname = 'AutoNeuralNet';
  argin = shiftargin(varargin,'integer');
  argin = setdefaults(argin,[],[]);

  if mapping_task(argin,'definition')
    argout = define_mapping(argin,'untrained',mapname);

  elseif mapping_task(argin,'training')			% Train a mapping.

    [a,units] = deal(argin{:});
    n_attempts = 3;			% Try three different random initialisations.
    [m,k] = size(a);
    if isempty(units)
      cs = classsizes(a);
      units = ceil(0.2*min(cs));
    end

  	if isnan(units) % optimize complexity parameter: number of neurons
			defs = {[]};
			parmin_max = [1,30];
			w = regoptc(a,mfilename,{units},defs,[1],parmin_max,testc([],'soft'),0);
			return
  	end

		islabtype(a,'crisp');
		isvaldfile(a,1,2); % at least 1 object per class, 2 classes
		a = testdatasize(a);
		a = setprior(a,getprior(a,0));

		% train a network.
		% Reproducability: always use same seeds.
		randstate = randreset(1); opt_err = inf; opt_mapping = [];

		% Try a number of random initialisations.
		s = sprintf('%i neural network initializations: ',n_attempts);
		prwaitbar(n_attempts,s);
		for attempt = 1:n_attempts
			prwaitbar(n_attempts,attempt,[s int2str(attempt)]);
			prwarning(4,'training with initialisation %d of %d',attempt,n_attempts);
			t = gendatk(a,1000,2,1); 			% Create tuning set based on training set.
			w = lmnc(a,units,inf,[],t);		% Find LMNC mapping.
			e = t*w*testc;								% Calculate classification error.
			if (e < opt_err)
				% If this is the best of the three repetitions, store it.
				opt_mapping = w; opt_err = e;
			end
		end
		prwaitbar(0);
    randreset(randstate); % return original state

		% Output is best network found.
		argout = setname(opt_mapping,mapname);

  else % Evaluation

    [a,w] = deal(argin{1:2});
		nodatafile(a);
		data = getdata(w);

		if (length(data) > 1)

	    % "Old" neural network - network is second parameter: unpack.
  	  data = getdata(w); weights = data{1};
    	pars = data{2}; numlayers = length(pars);

	    output = a;                       % Output of first layer: dataset.
  	  for j = 1:numlayers-1
    	  % Number of inputs (n_in) and outputs (n_out) of neurons in layer J.
      	n_in = pars(j); n_out = pars(j+1);

	      % Calculate output of layer J+1. Note that WEIGHTS contains both
  	    % weights (multiplied by previous layer's OUTPUT) and biases
    	  % (multiplied by ONES).

      	this_weights = reshape(weights(1:(n_in+1)*n_out),n_in+1,n_out);
	      output = sigm([output,ones(m,1)]*this_weights);

  	    % Remove weights of this layer.
    	  weights(1:(n_in+1)*n_out) = [];
	    end
		else
			% "New" neural network: unpack and simulate using the toolbox.
			net = data{1};
			output = sim(net,+a')';
		end;

		% 2-class case, therefore 1 output: 2nd output is 1-1st output.
		if (size(output,2) == 1)
			output = [output (1-output)];
		end

		% Output is mapped dataset.
		argout = setdat(a,output,w);

	end

return

author: marianux

repo: ecg-kit


text:
function [varargout]=subHex(varargin)

% function [E,V,C,CV]=subHex(E,V,n,splitMethod)
% ------------------------------------------------------------------------
%
% Kevin Mattheus Moerman
% gibbon.toolbox@gmail.com
%
%
% 2017/08/31 Added function description to top of function
% ------------------------------------------------------------------------

%% Parse input

switch nargin
    case 2
        E=varargin{1};
        V=varargin{2};
        n=1;
        splitMethod=1;
    case 3
        E=varargin{1};
        V=varargin{2};
        n=varargin{3};
        splitMethod=1;
    case 4
        E=varargin{1};
        V=varargin{2};
        n=varargin{3};
        splitMethod=varargin{4};
end

C=(1:1:size(E,1))'; %Element colors or indices
CV=[];

%%
if n>0
    for qIter=1:1:n
        switch splitMethod
            case 1
                %% Mid edge sets
                edgeMat=[E(:,[1 2]); E(:,[2 3]);  E(:,[3 4]); E(:,[4 1]);... %top
                    E(:,[5 6]); E(:,[6 7]);  E(:,[7 8]); E(:,[8 5]);... %bottom
                    E(:,[1 5]); E(:,[2 6]);  E(:,[3 7]); E(:,[4 8])];   %top-bottom edges

                E_sort=sort(edgeMat,2); %Sorted edges matrix
                [~,ind1,~]=unique(E_sort,'rows');
                edgeMat=edgeMat(ind1,:);

                numPoints = size(V,1);
                numEdges = size(edgeMat,1);

                % Get indices of the four edges associated with each face
                A = sparse(edgeMat(:,1),edgeMat(:,2),(1:numEdges)+numPoints,numPoints,numPoints,numEdges);
                A = max(A,A'); %Copy symmetric

                %Indices for A matrix
                indA_12=E(:,1)+(E(:,2)-1)*numPoints;
                indA_23=E(:,2)+(E(:,3)-1)*numPoints;
                indA_34=E(:,3)+(E(:,4)-1)*numPoints;
                indA_41=E(:,4)+(E(:,1)-1)*numPoints;

                indA_56=E(:,5)+(E(:,6)-1)*numPoints;
                indA_67=E(:,6)+(E(:,7)-1)*numPoints;
                indA_78=E(:,7)+(E(:,8)-1)*numPoints;
                indA_85=E(:,8)+(E(:,5)-1)*numPoints;

                indA_15=E(:,1)+(E(:,5)-1)*numPoints;
                indA_26=E(:,2)+(E(:,6)-1)*numPoints;
                indA_37=E(:,3)+(E(:,7)-1)*numPoints;
                indA_48=E(:,4)+(E(:,8)-1)*numPoints;

                %Get indices for vertex array
                indV_12=full(A(indA_12));
                indV_23=full(A(indA_23));
                indV_34=full(A(indA_34));
                indV_41=full(A(indA_41));

                indV_56=full(A(indA_56));
                indV_67=full(A(indA_67));
                indV_78=full(A(indA_78));
                indV_85=full(A(indA_85));

                indV_15=full(A(indA_15));
                indV_26=full(A(indA_26));
                indV_37=full(A(indA_37));
                indV_48=full(A(indA_48));

                %% Mid element
                indV_mid=(1:1:size(E,1))'+numPoints+size(edgeMat,1);

                %% Mid face

                %Element faces matrix
                F =[E(:,[4 3 2 1]);... %top
                    E(:,[5 6 7 8]);... %bottom
                    E(:,[1 2 6 5]);... %side 1
                    E(:,[3 4 8 7]);... %side 2
                    E(:,[2 3 7 6]);... %front
                    E(:,[5 8 4 1]);]; %back

                F_sort=sort(F,2); %Sorted edges matrix
                [~,ind1,ind2]=unique(F_sort,'rows');
                F=F(ind1,:);

                indV_midFace=(1:1:size(F_sort,1))';
                indOffset=numPoints+size(edgeMat,1)+size(E,1);

                indV_midFace4321=ind2(indV_midFace((1-1)*size(E,1)+(1:size(E,1))))+indOffset;
                indV_midFace5678=ind2(indV_midFace((2-1)*size(E,1)+(1:size(E,1))))+indOffset;
                indV_midFace1265=ind2(indV_midFace((3-1)*size(E,1)+(1:size(E,1))))+indOffset;
                indV_midFace3487=ind2(indV_midFace((4-1)*size(E,1)+(1:size(E,1))))+indOffset;
                indV_midFace2376=ind2(indV_midFace((5-1)*size(E,1)+(1:size(E,1))))+indOffset;
                indV_midFace5841=ind2(indV_midFace((6-1)*size(E,1)+(1:size(E,1))))+indOffset;

                %% Create element array
                Es=[E(:,1) indV_12 indV_midFace4321 indV_41 indV_15 indV_midFace1265 indV_mid indV_midFace5841;...%Corner hex 1
                    indV_12 E(:,2) indV_23 indV_midFace4321 indV_midFace1265 indV_26 indV_midFace2376 indV_mid;... %Corner hex 2
                    indV_midFace4321 indV_23 E(:,3) indV_34 indV_mid indV_midFace2376 indV_37 indV_midFace3487;... %Corner hex 3
                    indV_41 indV_midFace4321 indV_34 E(:,4) indV_midFace5841 indV_mid indV_midFace3487 indV_48;... %Corner hex 4
                    indV_15 indV_midFace1265 indV_mid indV_midFace5841 E(:,5) indV_56 indV_midFace5678 indV_85 ;...%Corner hex 5
                    indV_midFace1265 indV_26 indV_midFace2376 indV_mid  indV_56 E(:,6) indV_67 indV_midFace5678;... %Corner hex 6
                    indV_mid indV_midFace2376 indV_37 indV_midFace3487 indV_midFace5678 indV_67 E(:,7) indV_78 ;... %Corner hex 7
                    indV_midFace5841 indV_mid indV_midFace3487 indV_48 indV_85 indV_midFace5678 indV_78 E(:,8) ;... %Corner hex 8
                    ];

                %% Create vertex array
                Vn=0.5*(V(edgeMat(:,1),:)+V(edgeMat(:,2),:)); %new mid-edge points
                Vm=patchCentre(E,V); %Mid element nodes

                Vf=zeros(size(F,1),3);
                for q=1:1:size(V,2)
                    X=V(:,q);
                    Vf(:,q)=mean(X(F),2);
                end
                Vs = [V; Vn; Vm; Vf]; %Join point sets

                CV=[0*ones(size(V,1),1); 1*ones(size(Vn,1),1); 2*ones(size(Vm,1),1); 3*ones(size(Vf,1),1);];

            case 2

                %% Mid edge sets
                edgeMat=[E(:,[1 5]); E(:,[2 6]);  E(:,[3 7]); E(:,[4 8])];   %top-bottom edges

                E_sort=sort(edgeMat,2); %Sorted edges matrix
                [~,ind1,~]=unique(E_sort,'rows');
                edgeMat=edgeMat(ind1,:);

                numPoints = size(V,1);
                numEdges = size(edgeMat,1);

                % Get indices of the four edges associated with each face
                A = sparse(edgeMat(:,1),edgeMat(:,2),(1:numEdges)+numPoints,numPoints,numPoints,numEdges);
                A = max(A,A'); %Copy symmetric

                %Indices for A matrix
                indA_15=E(:,1)+(E(:,5)-1)*numPoints;
                indA_26=E(:,2)+(E(:,6)-1)*numPoints;
                indA_37=E(:,3)+(E(:,7)-1)*numPoints;
                indA_48=E(:,4)+(E(:,8)-1)*numPoints;

                %Get indices for vertex array
                indV_15=full(A(indA_15));
                indV_26=full(A(indA_26));
                indV_37=full(A(indA_37));
                indV_48=full(A(indA_48));

                %% Create element array
                Es=[E(:,1:4) indV_15 indV_26 indV_37 indV_48;...
                    indV_15 indV_26 indV_37 indV_48 E(:,5:8)];

                %% Create vertex array
                Vn=0.5*(V(edgeMat(:,1),:)+V(edgeMat(:,2),:)); %new mid-edge points
                Vs = [V; Vn;]; %Join point sets
                CV=[0*ones(size(V,1),1); 1*ones(size(Vn,1),1);];

            case 3

                %% Mid edge sets
                edgeMat=[E(:,[1 2]); E(:,[5 6]);  E(:,[7 8]); E(:,[3 4])]; %left-right edges

                E_sort=sort(edgeMat,2); %Sorted edges matrix
                [~,ind1,~]=unique(E_sort,'rows');
                edgeMat=edgeMat(ind1,:);

                numPoints = size(V,1);
                numEdges = size(edgeMat,1);

                % Get indices of the four edges associated with each face
                A = sparse(edgeMat(:,1),edgeMat(:,2),(1:numEdges)+numPoints,numPoints,numPoints,numEdges);
                A = max(A,A'); %Copy symmetric

                %Indices for A matrix
                indA_12=E(:,1)+(E(:,2)-1)*numPoints;
                indA_56=E(:,5)+(E(:,6)-1)*numPoints;
                indA_78=E(:,7)+(E(:,8)-1)*numPoints;
                indA_34=E(:,3)+(E(:,4)-1)*numPoints;

                %Get indices for vertex array
                indV_12=full(A(indA_12));
                indV_56=full(A(indA_56));
                indV_78=full(A(indA_78));
                indV_34=full(A(indA_34));

                %% Create element array
                Es=[E(:,1) indV_12 indV_34 E(:,4) E(:,5) indV_56 indV_78 E(:,8);...
                    indV_12 E(:,2) E(:,3) indV_34 indV_56  E(:,6) E(:,7) indV_78];

                %% Create vertex array
                Vn=0.5*(V(edgeMat(:,1),:)+V(edgeMat(:,2),:)); %new mid-edge points
                Vs = [V; Vn;]; %Join point sets
                CV=[0*ones(size(V,1),1); 1*ones(size(Vn,1),1);];

            case 4

                %% Mid edge sets
                edgeMat=[E(:,[4 1]); E(:,[8 5]);  E(:,[6 7]); E(:,[2 3])]; %front-back edges

                E_sort=sort(edgeMat,2); %Sorted edges matrix
                [~,ind1,~]=unique(E_sort,'rows');
                edgeMat=edgeMat(ind1,:);

                numPoints = size(V,1);
                numEdges = size(edgeMat,1);

                % Get indices of the four edges associated with each face
                A = sparse(edgeMat(:,1),edgeMat(:,2),(1:numEdges)+numPoints,numPoints,numPoints,numEdges);
                A = max(A,A'); %Copy symmetric

                %Indices for A matrix
                indA_41=E(:,4)+(E(:,1)-1)*numPoints;
                indA_85=E(:,8)+(E(:,5)-1)*numPoints;
                indA_67=E(:,6)+(E(:,7)-1)*numPoints;
                indA_23=E(:,2)+(E(:,3)-1)*numPoints;

                %Get indices for vertex array
                indV_41=full(A(indA_41));
                indV_85=full(A(indA_85));
                indV_67=full(A(indA_67));
                indV_23=full(A(indA_23));

                %% Create element array
                Es=[E(:,[1 2]) indV_23 indV_41 E(:,[5 6]) indV_67 indV_85;...
                    indV_41 indV_23 E(:,[3 4]) indV_85 indV_67 E(:,[7 8])];

                %% Create vertex array
                Vn=0.5*(V(edgeMat(:,1),:)+V(edgeMat(:,2),:)); %new mid-edge points
                Vs = [V; Vn;]; %Join point sets
                CV=[0*ones(size(V,1),1); 1*ones(size(Vn,1),1);];

        end

        %% Override input
        C=repmat(C,[size(Es,1)/size(E,1),1]);
        E=Es;
        V=Vs;
    end
end

varargout{1}=E;
varargout{2}=V;
varargout{3}=C;
varargout{4}=CV;

%%
% _*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:
function isgn = i4row_compare ( m, n, a, i, j )

%*****************************************************************************80
%
%% I4ROW_COMPARE compares two rows of an I4ROW.
%
%  Example:
%
%    Input:
%
%    M = 3, N = 4, I = 2, J = 3
%
%    A = (
%      1  2  3  4
%      5  6  7  8
%      9 10 11 12 )
%
%    Output:
%
%      ISGN = -1
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 April 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns.
%
%    Input, integer A(M,N), an array of M rows of vectors of length N.
%
%    Input, integer I, J, the rows to be compared.
%    I and J must be between 1 and M.
%
%    Output, integer ISGN, the results of the comparison:
%    -1, row I < row J,
%     0, row I = row J,
%    +1, row J < row I.
%

%
%  Check that I and J are legal.
%
  if ( i < 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4ROW_COMPARE - Fatal error!\n' );
    fprintf ( 1, '  Row index I is less than 1.\n' );
    fprintf ( 1, '  I = %d\n', i );
    error ( 'I4ROW_COMPARE - Fatal error!' );
  elseif ( m < i )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4ROW_COMPARE - Fatal error!\n' );
    fprintf ( 1, '  Row index I is out of bounds.\n' );
    fprintf ( 1, '  I = %d\n', i );
    fprintf ( 1, '  Maximum legal value is M = %d\n', m );
    error ( 'I4ROW_COMPARE - Fatal error!' );
  end

  if ( j < 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4ROW_COMPARE - Fatal error!\n' );
    fprintf ( 1, '  Row index J is less than 1.\n' );
    fprintf ( 1, '  J = %d\n', j );
    error ( 'I4ROW_COMPARE - Fatal error!' );
  elseif ( m < j )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4ROW_COMPARE - Fatal error!\n' );
    fprintf ( 1, '  Row index J is out of bounds.\n' );
    fprintf ( 1, '  J = %d\n', j );
    fprintf ( 1, '  Maximum legal value is M = %d\n', m );
    error ( 'I4ROW_COMPARE - Fatal error!' );
  end

  isgn = 0;

  if ( i == j )
    return
  end

  k = 1;

  while ( k <= n )

    if ( a(i,k) < a(j,k) )
      isgn = -1;
      return
    elseif ( a(j,k) < a(i,k) )
      isgn = +1;
      return
    end

    k = k + 1;

  end

  return
end


text:
function [ DCM ] = tapas_dcm_euler_make_indices(DCM)
% [ DCM ] = tapas_dcm_euler_make_indices(DCM)
%
% Get the indices for the Euler integration
%
%   Input:
%   	DCM         - model structure
%
%   Output:
%       DCM         - model structure with indices
%

% ----------------------------------------------------------------------
%
% Authors: Stefan Fraessle (stefanf@biomed.ee.ethz.ch), Ekaterina I. Lomakina
%
% Copyright (C) 2016-2018 Translational Neuromodeling Unit
%                         Institute for Biomedical Engineering
%                         University of Zurich & ETH Zurich
%
% This file is part of the TAPAS rDCM Toolbox, which is released under the
% terms of the GNU General Public License (GPL), version 3.0 or later. You
% can redistribute and/or modify the code under the terms of the GPL. For
% further see COPYING or <http://www.gnu.org/licenses/>.
%
% Please note that this toolbox is in an early stage of development. Changes
% are likely to occur in future releases.
%
% ----------------------------------------------------------------------

% length of input
L = size(DCM.U.u,1);

% number of regions
nr = size(DCM.a, 1);

% check whether there is a delay specified
if ~isfield('DCM','delay')
    DCM.delay = ones(nr,1);
end

% create index array
Indices = repmat(1 : L, 1, 1);

% get the indices that coincide wiht the data timepoints
idx = zeros(1, L);
idx(DCM.delay:floor(DCM.Y.dt/DCM.U.dt):end) = 1;

% asign those timings
Indices = Indices(idx>0);
DCM.M.idx = Indices;

end

author: translationalneuromodeling

repo: tapas


text:
function MatingPool = MatingSelection(PopObj,N,div)
% The mating selection of PESA-II

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    %% Calculte the grid location of each solution
    fmax = max(PopObj,[],1);
    fmin = min(PopObj,[],1);
    d    = (fmax-fmin)/div;
    GLoc = floor((PopObj-repmat(fmin,size(PopObj,1),1))./repmat(d,size(PopObj,1),1));
    GLoc(GLoc>=div)   = div - 1;
    GLoc(isnan(GLoc)) = 0;

    %% Calculate the crowding degree of each grid
    [UniqueGLoc,~,Site] = unique(GLoc,'rows');
    CrowdG              = hist(Site,1:max(Site));

    %% Binary tournament selection
    MatingPool = zeros(1,N);
    for i = 1 : length(MatingPool)
        grid          = randi(size(UniqueGLoc,1),1,2);
        [~,best]      = min(CrowdG(grid));
        current       = find(Site==grid(best));
        MatingPool(i) = current(randi(length(current)));
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO