text:
% pop_chanevent() - import event latencies from the rising and/or falling 'edge'
%                   latencies of a specified event-marker channel in EEG.data
% Usage:
%   >> OUTEEG = pop_chanevent( INEEG ); % select parameters via a pop-up window
%   >> OUTEEG = pop_chanevent( INEEG, chanindices, 'key', 'val' ... ); % no pop-up
%
% Graphic interface:
%   "Event channel(s)" - [edit box] indices of event channel(s) to import.
%                  Command line equivalent: chanindices.
%   "Preprocessing transform" - [edit box] apply this preprocessing
%                  formula or function to the selected data channel(s) X,
%                  transforming X into the command output before edge
%                  extraction. Command line equivalent 'oper'.
%   "Transition to extract" - [list box] extract events when the event
%                  channel values go up ('leading'), down ('trailing')
%                  or both ('both'). Command line equivalent: 'edge'.
%   "Transition length" - [edit box] Increase this number to avoid having
%                  events very close to each other due to a not perfectly
%                  straight edge. Command line equivalent: 'edgelen'.
%   "Assign duration to events?" - [checkbox] . Assign duration to each
%                  extracted event.  This option can only be used when
%                  extracting events on leading edges. Event will last
%                  until next trailing edge (down) event. Command line
%                  equivalent: 'duration'.
%   "Delete event channel(s)" - [checkbox] check to delete the event channel
%                  after events have been extracted from it.
%                  Command line equivalent: 'delchan'.
%   "Delete old events if any" - [checkbox] check this checkbox to
%                  remove any prior events in the dataset. Otherwise
%                  imported events are appended to old events. Command
%                  line equivalent: 'delevent'.
%   "Only one event type" - [checkbox] check this checkbox to assign
%                  all transitions in the event channel to one event
%                  type. Else, one type is assigned for each non-zero
%                  channel value. Command line equivalent: 'nbtype'.
% Inputs:
%   INEEG          - input dataset structure
%   chanindices    - index|(indices) of the event channel(s)
%
% Optional inputs:
%   'edge'         - ['leading'|'trailing'|'both'] extract events when values
%                    in the event channel go up ('leading'), down ('trailing')
%                    or both ('both'). {Default is 'both'}.
%   'edgelen'      - [integer] maximum edge length (for some data edge do not
%                    take whole value and it takes a few sample points for
%                    signal to rise. Default is 1 (perfect edges).
%   'oper'         - [string] prior to extracting edges, preprocess data
%                    channel(s) using the string command argument.
%                    In this command, the data channel(s) are designated by
%                    (capital) X. For example, 'X>3' will test the value of X
%                    at each time point (returning 1 if the data channel value
%                    is larger than 3, and 0 otherwise). You may also use
%                    any function (Ex: 'myfunction(X)').
%   'duration'     - ['on'|'off'] extract event duration. This option can only be
%                    used when extracting events on leading edges. Event will last
%                    until next trailing-edge (down) event { 'off' }.
%   'delchan'      - ['on'|'off'] delete channel from data { 'on' }.
%   'delevent'     - ['on'|'off'] delete old events if any { 'on' }.
%   'nbtype'       - [1|NaN] setting this to 1 will force the program to
%                    consider all events to have the same type. {Default is NaN}.
%                    If set (1), all transitions are considered the same event type
%			         If unset (NaN), each (transformed) event channel value following a
%                    transition determines an event type (Ex: Detecting leading-edge
%                    transitions of 0 0 1 0 2 0 ...  produces event types 1 and 2).
%   'typename'     - [string] event type name. Only relevant if 'nbtype' is 1
%                    or if there is only one event type in the event channel.
%                    {Default is 'chanX', X being the index of
%                    the selected event channel}.
% Outputs:
%   OUTEEG         - EEGLAB output data structure
%
% Author: Arnaud Delorme, CNL / Salk Institute, 29 July 2002
%
% See also: eeglab()

% Copyright (C) 2002 Arnaud Delorme, Salk Institute, arno@salk.edu
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

function [EEG, command] = pop_chanevent(EEG, chan, varargin);
command = '';

if nargin < 1
    help pop_chanevent;
    return;
end;

if nargin < 2
	geometry = { [1.5 1 1] [1] [1.5 1 1] [1.5 1 1] [1.5 1 1] [1.5 0.2 0.36 0.84] ...
                 [1] [1.5 0.21 1] [1.5 0.21 1] [1.5 0.21 1] };

    % callback from listbox to disable duration checkbox (if leading event is not selected)
    % --------------------------------------------------
    cb_list = [ 'if get(gcbo, ''value'') == 1,' ...
                '  set(findobj(gcbf, ''tag'', ''dur''), ''enable'', ''on'');' ...
                'else,' ...
                '  set(findobj(gcbf, ''tag'', ''dur''), ''enable'', ''off'', ''value'', 0);' ...
                'end;' ];

	strgui = { { 'style' 'text' 'string' 'Event channel(s)' 'tooltipstring' 'indexes of event channels' } ...
			   { 'style' 'edit' 'string' '' } { } ...
               {} ...
			   { 'style' 'text' 'string' 'Preprocessing transform (data=''X'')' 'tooltipstring'  ...
                 [ 'For example, ''X>3'' will test the value of X' 10 ...
                  'at each time point (returning 1 if the data channel value' 10 ...
                  'is larger than 3, and 0 otherwise).' ] } ...
			   { 'style' 'edit' 'string' '' } { 'style' 'text' 'string' 'Optional. Ex: X>3' } ...
			   { 'style' 'text' 'string' 'Transitions to extract? (up|down)' 'tooltipstring' ...
				 [ 'Extract events whenever values in the (transformed) event channel(s) shift up' 10 ...
				   '(''leading''), down (''trailing'') or either (''both'').' 10 ...
				   'AFTER SCROLLING CLICK TO SELECT' ] } ...
			   { 'style' 'listbox' 'string' 'up (leading)|both|down (trailing)' 'value' 1 'callback' cb_list } ...
               { 'style' 'text' 'string' '(click to select)'} ...
			   { 'style' 'text' 'string' 'Transition length (1=perfect edges)' 'tooltipstring'  ...
                 [ 'Increase this number to avoid having events very close to each other due.' 10 ...
                   'to a not perfectly straight edge' ] } ...
			   { 'style' 'edit' 'string' '0' } { } ...
			   { 'style' 'text' 'string' 'Assign duration to each events?' 'tag' 'dur' 'tooltipstring' ...
				 [ 'You may assign an event duration to each event if you select to detect' 10 ...
				   'event on the leading edge above. Event will last as long as the signal is non-0.' ] } ...
			   { 'style' 'checkbox' 'string' '' 'value' 0 'tag' 'dur'} { } ...
               { 'style' 'text' 'string' '(set=yes)' } ...
               {} ...
               { 'style' 'text' 'string' 'Delete event channel(s)? ' } ...
			   { 'style' 'checkbox' 'value' 1 } { 'style' 'text' 'string' '        (set = yes)'} ...
			   { 'style' 'text' 'string' 'Delete old events if any? ' } ...
			   { 'style' 'checkbox' 'value' 1 } { } ...
			   { 'style' 'text' 'string' 'All events of same type? ' 'tooltipstring' ...
			    ['If set, all transitions are considered the same event type,' 10 ...
				 'If unset, each (transformed) event channel value following a transition' 10 ...
                 'determines an event type (Ex: Detecting leading-edge transitions of' 10 ...
                 '0 0 1 0 2 0 ...  produces event types 1 and 2).' ] } ...
			   { 'style' 'checkbox' 'value' 0 } { } };
	result       = inputgui( geometry, strgui, 'pophelp(''pop_chanevent'');', 'Extract event from channel(s) - pop_chanevent()');

	if length(result) == 0 return; end;
	chan   = eval( [ '[' result{1} ']' ] );
    options = {};
	if ~isempty(result{2}), options = { options{:} 'oper' result{2} }; end;
	switch result{3},
		case 1, options = { options{:} 'edge' 'leading' };
		case 2, options = { options{:} 'edge' 'both' };
		case 3, options = { options{:} 'edge' 'trailing' };
	end;
    options = { options{:} 'edgelen' eval( [ '[' result{4} ']' ] ) };
    if result{5},  options = { options{:} 'duration' 'on' }; end;
	if ~result{6}, options = { options{:} 'delchan'  'off'}; end;
	if ~result{7}, options = { options{:} 'delevent' 'off'}; end;
	if result{8},  options = { options{:} 'nbtype'  1}; end;
else
    options = varargin;
end;
listcheck = { 'edge'      'string'     { 'both' 'leading' 'trailing'}     'both';
              'edgelen'   'integer'    [1 Inf]                            1;
              'delchan'   'string'     { 'on' 'off' }                     'on';
              'oper'      'string'     []                                 '';
              'delevent'  'string'     { 'on' 'off' }                     'on';
              'duration'  'string'     { 'on' 'off' }                     'off';
              'typename'  'string'     []                                 [ 'chan' int2str(chan) ];
              'nbtype'    'integer'    [1 NaN]                             NaN };
g = finputcheck( options, listcheck, 'pop_chanedit');
if isstr(g), error(g); end;

% check inut consistency
% ----------------------
if strcmpi(g.duration, 'on') & ~strcmpi(g.edge, 'leading')
    error('Must detect leading edge to extract event duration');
end;

% process events
% --------------
fprintf('pop_chanevent: importing events from data channel %d ...\n', chan);
counte = 1; % event counter
events(10000).latency = 0;
if isnan(g.nbtype)
    if length(unique(EEG.data(chan, :))) == 2, g.nbtype = 1; end;
end;

for ci = chan
    X = EEG.data(ci, :);

    % apply preprocessing
    % -------------------
    if ~isempty(g.oper)
        try, eval( [ 'X = ' g.oper ';' ]);
        catch, error('pop_chanevent: error executing preprocessing string');
        end;
    end;

    % extract edges
    % -------------
    tmpdiff =  diff(abs([ X X(end) ]));
    switch g.edge
     case 'both'    , tmpevent1 = find( tmpdiff > 0)-1; tmpevent2 = find( tmpdiff < 0);
     case 'trailing', tmpevent2 = find( tmpdiff < 0);
     case 'leading' , tmpevent1 = find( tmpdiff > 0)-1; tmpdur = find( tmpdiff < 0);
    end;

    % fuse close events if necessary
    % ------------------------------
    if exist('tmpevent1')
        tmpclose = find( tmpevent1(2:end)-tmpevent1(1:end-1) < g.edgelen)+1;
        tmpevent1(tmpclose) = [];
        tmpevent    = tmpevent1+1;
        tmpeventval = tmpevent1+2;
    end;
    if exist('tmpevent2')
        tmpclose = find( tmpevent2(2:end)-tmpevent2(1:end-1) < g.edgelen); % not +1
        tmpevent2(tmpclose) = [];
        tmpevent    = tmpevent2+1;
        tmpeventval = tmpevent2;
    end;
    if exist('tmpevent1') & exist('tmpevent2')
        tmpevent    = sort([ tmpevent1+1 tmpevent2+1]);
        tmpeventval = sort([ tmpevent1+2 tmpevent2]);
    end;

    % adjust edges for duration  if necessary
    % ---------------------------------------
    if strcmpi(g.duration, 'on')
        tmpclose = find( tmpdur(2:end)-tmpdur(1:end-1) < g.edgelen); % not +1 (take out the first)
        tmpdur(tmpclose) = [];
        if tmpdur(1) < tmpevent(1), tmpdur(1) = []; end;
        if length(tmpevent) > length(tmpdur), tmpdur(end+1) = EEG.pnts; end;
        if length(tmpevent) ~= length(tmpdur)
            error([ 'Error while attempting to extract event durations' 10 ...
                    'Maybe edges are not perfectly defined, try increasing edge length' ]);
        end;
    end;

    if isempty(tmpevent),
        fprintf('No event found for channel %d\n', ci);
    else
        for tmpi = 1:length(tmpevent)
            if ~isnan(g.nbtype)
                events(counte).type    = g.typename;
            else
                events(counte).type    = X(tmpeventval(tmpi));
            end;
            events(counte).latency = tmpevent(tmpi);
            if strcmpi(g.duration, 'on')
                events(counte).duration = tmpdur(tmpi) - tmpevent(tmpi);
            end;
            counte = counte+1;
        end;
    end;
    events = events(1:counte-1);
end;

% resort events
% --------------
if strcmp(g.delevent, 'on')
	EEG.event = events;
    if EEG.trials > 1
        for index = 1:length(events)
            EEG.event(index).epoch = 1+floor((EEG.event(index).latency-1) / EEG.pnts);
        end;
    end;
else
	for index = 1:length(events)
		EEG.event(end+1).type  = events(index).type;
		EEG.event(end).latency = events(index).latency;
        if EEG.trials > 1 | isfield(EEG.event, 'epoch');
            EEG.event(end).epoch = 1+floor((EEG.event(end).latency-1) / EEG.pnts);
        end;
	end;
    if EEG.trials > 1
        EEG = pop_editeventvals( EEG, 'sort', {  'epoch' 0 'latency', [0] } );
    else
        EEG = pop_editeventvals( EEG, 'sort', { 'latency', [0] } );
    end;
end;
if isfield(EEG.event, 'urevent'), EEG.event = rmfield(EEG.event, 'urevent'); end;
EEG = eeg_checkset(EEG, 'eventconsistency');
EEG = eeg_checkset(EEG, 'makeur');

% delete channels
% ---------------
if strcmp(g.delchan, 'on')
	EEG = pop_select(EEG, 'nochannel', chan);
end;

if nargin < 2
    command = sprintf('%s = pop_chanevent(%s, %s);', inputname(1), inputname(1), ...
                      vararg2str({ chan options{:} }));
end;
return;

author: PatternRecognition

repo: OpenBMI


text:
function ll = mlpLogLikelihood(model)

% MLPLOGLIKELIHOOD Multi-layer perceptron log likelihood.
% FORMAT
% DESC computes the log likelihood of a multi-layer perceptron
% model. For single hidden layer models this is done by wrapping
% the mlperr command.
% ARG model : the model structure for computing the log likelihood.
% RETURN ll : the model log likelihood.
%
% SEEALSO : modelLogLikeihood, mlperr
%
% COPYRIGHT : Neil D. Lawrence, 2006, 2007

% MLTOOLS

if length(model.hiddenDim) == 1
  ll = -mlperr(model, model.X, model.y);
else
  Y = mlpOut(model, model.X);
  ll = -0.5*sum(sum((model.Y - Y).^2));
end

ll = ll - size(model.X, 1)/2*log(2*pi);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function D = diffmat(n, p)
%DIFFMAT  Trigonometric differentiation matrix.
%   D = DIFFMAT(N) is the matrix that maps function values at N equally-spaced
%   points in [0 2*pi) to values of the derivative of the Fourier interpolant
%   at those points.
%
%   D = DIFFMAT(N, K) is the same, but for the Kth derivative.
%
% See also DIFF.

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% References for the rectangular differentiation matrices:
%
% J.A.C. Weideman and S.C. Reddy, A MATLAB differentiation matrix suite, ACM
% Transcations on Mathematical Software, Vol. 26, No. 4, Page 465--519, 2000.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% grid spacing:
h = 2*pi/n;

% sign of entries in a column:
sgn = (-1).^(1:n-1).';

% indices:
n1 = floor((n-1)/2);
n2 = ceil((n-1)/2);

% grid points on [-1 0):
v = (1:n2)'*h/2;

% For different p:
if ( p == 0 )           % Zeroth-derivative
    D = eye(n);
else
    switch p
        case 1          % 1st-derivative

            % Forming the first column by 'flipping trick':
            if ( rem(n, 2) )           % n is odd
                tmp = csc(v);
                col = [0; (pi/2)*sgn.*[tmp; flipud(tmp(1:n1))]];
            else                      % n is even
                tmp = cot(v);
                col = [0; (pi/2)*sgn.*[tmp; -flipud(tmp(1:n1))]];
            end

            % Form the first row:
            row = -col;

        case 2         % 2nd-derivative

            % Form the first column by 'flipping trick':
            if ( rem(n, 2) )           % n is odd
                tmp = csc(v).*cot(v);
                col = pi^2*[-pi^2/(3*h^2)+1/12; ...
                    -0.5*sgn.*[tmp; -flipud(tmp(1:n1))]];
            else                      % n is even
                tmp = csc(v).^2;
                col = pi^2*[-pi^2/(3*h^2)-1/6; ...
                    -0.5*sgn.*[tmp; flipud(tmp(1:n1))]];
            end

            % Form the first row:
            row = col;

        otherwise      % pth-derivative for p >= 3

            % Form the first column using FFT:
            n3 = (-n/2)*rem(p+1,2)*ones(rem(n+1,2));
            waveNumber = 1i*[(0:n1) n3 (-n1:-1)];
            col = pi^p*real(ifft(waveNumber.^p.*fft(eye(1,n))));

            % Form the first row:
            if ( rem(p,2) )      % p is odd
                col = [0 col(2:n)]';
                row = -col;
            else                 % p is even
                row = col;
            end
    end

    % Form the differentiation matrix which is toeplitz:
    D = toeplitz(col, row);

end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
  function sh = de_ftab_s_iter(fit, fh, varargin)
%|function sh = de_ftab_s_iter(fit, fh, [options])
%| estimate s from fh by iterative LS
%| sh = argmin{s} | fh - fm(s) |^2
%| in
%|	fit			from de_ftab_fit() / de_ftab_curv()
%|	fh	[(Nd) M]	estimates of f (nonlinear BH function)
%| option
%|	'niter'			# of iterations (default: 1 - not enough!)
%|				(if 0, then just initial estimate returned)
%|	'init'	[(Nd) L]	initial estimates (default: linear inverse)
%|	'ctype'	char		curvature type (see de_ftab_curv.m)
%|				default '': inherit from de_ftab_curv.m
%| out
%|	sh	[(Nd) L]	estimates of s (component density integrals)
%|
%| Copyright 2006-05-21, Jeff Fessler, University of Michigan

if nargin == 1 && streq(fit, 'test'), de_ftab_s_iter_test, return, end
if nargin < 2, ir_usage, end

arg.ctype = '';
arg.niter = 1;
arg.init = [];
arg.strue = []; % optional for testing
arg.nprint = 100; % every 100th iteration show rms
arg = vararg_pair(arg, varargin);

LL = fit.LL;
Nd = size(fh); MM = Nd(end); Nd = Nd(1:end-1);
fh = reshapee(fh, [], MM); % [*Nd M]
if MM ~= fit.MM, error 'M', end

if ~isempty(arg.strue)
	arg.strue = reshapee(arg.strue, prod(Nd), LL); % [*Nd L]
end

% initialize
if isempty(arg.init)
	sh = (pinv(fit.mac_eff) * fh')'; % linear inverse (perfect if mono)

	if ~isempty(arg.strue)
		pr 'rms(sh - arg.strue)' % report initial error, for testing!
	end
else
	sh = max(arg.init, 0);
	sh = reshape(sh, prod(Nd), LL); % [*Nd L]
end

if arg.niter < 1, return, end

% curvature
switch arg.ctype
case ''
%	fstep = @(sh, fm, fh) 1 ./ sum(fit.curv2(sh, fm, fh), 2); % [*Nd 1]
	fstep = @(sh, fm, fh) 1 ./ sum(fit.ls_curv(sh, fh, fm), 2); % [*Nd 1]

case 'old'
	[curv1 curv2] = de_ftab_curv_ub(fit); % each is [MM 1]
	fcurv2 = @(fm, fh) max(fh - fm, 0) * curv2; % [*Nd M] * [M 1] -> [*Nd 1]
%	fcurv2 = @(fm, fh) max(fh, 0) * curv2; % [*Nd M] * [M 1] -> [*Nd 1] TOO SLOW!?
	fstep = @(sh, fm, fh) 1 ./ (sum(curv1) + fcurv2(fm, fh)); % [*Nd 1]
	step = fstep(0, 0, fh); % [*Nd 1] precomputed!
	minmax(step)
	step = repmat(step, [1 LL]);

otherwise
	fstep = @(sh, fm, fh) 1 ./ ...
		sum(fit.ls_curv(sh, fh, fm, 'ctype', arg.ctype), 2); % [*Nd 1]
%	fail('bad ctype: %s', arg.ctype)
end

ticker reset
for ii=1:arg.niter % iterate
	ticker(mfilename, ii, arg.niter)

	fm = fit.fmfun(sh); % [*Nd M]
%	cost = mean(col(fh - fm).^2);

	if 1
		fgrad = fit.fgrad(sh); % [*Nd L M]
		tmp = repmat(fh - fm, [1 1 LL]); % [*Nd M L]
		tmp = permute(tmp, [1 3 2]); % [*Nd L M]
		fgrad = fgrad .* tmp;
	else
%		fgrad = fit.ls_grad(fh, sh);
	end

	step = fstep(sh, fm, fh); % [*Nd 1]
	step = repmat(step, [1 LL]); % [*Nd L]
%	minmax(step)
	sh = sh + step .* sum(fgrad, 3);

	if ~isempty(arg.strue) && (~rem(ii, arg.nprint) || ii==arg.niter)
		pr 'rms(sh - arg.strue)' % report error, for testing!
	end
end

sh = reshapee(sh, Nd, LL); % [(Nd) L]

% de_ftab_curv_ub()
% curv* are [M 1] upper bounds
%
function [curv1, curv2] = de_ftab_curv_ub(fit)

MM = fit.MM;
new = fit.fgrad(zeros([1 fit.LL])); % [1 L M]
curv1 = squeeze(sum(new.^2,2)); % [1 1 M] -> [M 1]

neghess = -fit.fhess(zeros([1 fit.LL])); % [1 L L M]
curv2 = zeros(MM,1);
for mm=1:MM
	h0 = squeeze(neghess(1,:,:,mm)); % [1 L L 1] -> [L L]
	curv2(mm) = norm(h0);
end

return

% old way:

oldcurv1 = zeros(MM,1);
oldcurv2 = zeros(MM,1);

for mm=1:MM
	alf = fit.coef{mm}; % [ne 1]
	mac = fit.mac{mm}; % [ne L]
	g0 = mac' * alf; % gradient of fit at s=0 (largest point)
	h0 = mac' * diag(alf) * mac - g0 * g0'; % hessian at s=0 (largest?)
	oldcurv1(mm) = norm(g0)^2;
	oldcurv2(mm) = norm(h0);
end

equivs(curv1, oldcurv1)
equivs(curv2, oldcurv2)

% de_ftab_s_iter_test1
% L=1 M=1 case for sanity checking
%
function de_ftab_s_iter_test1
stype = 'ps1';
stype = 'poly1,100';
xrs = xray_read_spectra(stype);
mtype = 'water';
mas = xray_read_mac(mtype);
s1 = linspace(0, 50, 26)';
fm = de_ftab_fm(s1, mas.mac(xrs.en), xrs.Ide);
fit = de_ftab_fit({s1}, fm, 'type', 'exp', 'mtype', mtype);
ctype = 'newt';
%ctype = 'pre10'; % slow!
fit = de_ftab_curv(fit, 'ctype', ctype);

%fm = fit.fmfun(s1); % test with model fm, not original

% iterate
sh = de_ftab_s_iter(fit, fm, 'strue', s1, 'ctype', '', ...
	'niter', 5, 'nprint', 1);

if im
	clf, plot(s1, s1, 'c-', s1, sh, 'y.-')
	grid, axis equal, axis square
	xlabel 'true', ylabel 'noiseless estimate'
	title(fit.ctype)
prompt
end

% de_ftab_s_iter_test2
%
function de_ftab_s_iter_test2
sl{1} = linspace(0, 50, 26);
sl{2} = linspace(0, 30, 31);
xrs = xray_read_spectra('ps1');
mtype = {'water', 'bone'};
mtype = {'water', 'aluminum'};
mas = xray_read_mac(mtype);
sll = ndgrid_jf('mat', sl{:});
fm = de_ftab_fm(sll, mas.mac(xrs.en), xrs.Ide);
fit = de_ftab_fit(sl, fm, 'type', 'exp', 'mtype', mtype);
fit = de_ftab_curv(fit, 'ctype', 'newt');
%fit = de_ftab_curv(fit, 'ctype', 'pre2');

fm = fit.fmfun(sll); % test with model fm, not original

if 1 % initialize with polynomial inverse
	mac = xray_make_mac(xrs, mas);
	T = pinv(mac.bar);
	inv2 = de_ftab_inv2(fit, sl, 'T', T);
	s0 = inv2.fun(fm);
	s0 = reshapee(s0, [], 2);
else
	% initializing with linear inv.
	s0 = de_ftab_s_iter(fit, fm, 'niter', 0);
end

strue = reshapee(sll, [], 2); % true
pr rms(s0 - strue) % initial error - large!

if 1 % picture of error of initial guess vs (s_1,s_2)
	pr fit.mac_eff
	pr cond(fit.mac_eff)
	tmp = reshape(sqrt(mean((s0 - strue).^2, 2)), size(sll(:,:,1)));
	if im
		clf, im(sl{:}, tmp), cbar % error is smallest at (0,0)
		xlabel(mtype{1}), ylabel(mtype{2})
	end
prompt
end

% todo: run profiler, or pre-tabulate the inverse...
% iterate; initializing with linear inv.
sh = de_ftab_s_iter(fit, fm, 'init', s0, 'niter', 90, ...
		'strue', sll, 'nprint', 10);
pr 'rms(reshapee(sh, [], 2) - strue)'

%sh = reshape(sh, size(sll));
if im
	clf, plot(sl{1}, sh(:,:,1), 'c', sl{2}, sh(:,:,2), 'y')
	grid, axis equal, axis square
	xlabel 'true', ylabel 'noiseless estimate'
prompt
%	plot(fm(:,:,1), fm(:,:,2), '.')
end

% de_ftab_s_iter_test
%
function de_ftab_s_iter_test
de_ftab_s_iter_test2
de_ftab_s_iter_test1

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [varargout]=quiverCurve(varargin)

%% Parse input

switch nargin
    case 1
        Vg=varargin{1};
        dirOpt=1;
        colorOpt='k';
        alphaLevel=1;
    case 2
        Vg=varargin{1};
        dirOpt=varargin{2};
        colorOpt='k';
        alphaLevel=1;
    case 3
        Vg=varargin{1};
        dirOpt=varargin{2};
        colorOpt=varargin{3};
        alphaLevel=1;
    case 4
        Vg=varargin{1};
        dirOpt=varargin{2};
        colorOpt=varargin{3};
        alphaLevel=varargin{4};
end

%%

if size(Vg,2)==2
    Vg(:,3)=0;
end

switch dirOpt
    case 1 %Forward
        V=Vg(1:end-1,:);
        U=diff(Vg,1,1);
    case 2 %Backward
        V=Vg(2:end,:);
        U=flipud(diff(flipud(Vg),1,1));
    case 3 %Central
        V=Vg;
        Uf=[diff(Vg,1,1); nan(1,size(Vg,2))];
        Ub=-[nan(1,size(Vg,2)); flipud(diff(flipud(Vg),1,1))];
        U=Uf;
        U(:,:,2)=Ub;
        U=gnanmean(U,3);
end
D=sqrt(sum(U.^2,2));
a=[min(D(:)) max(D(:))];

[F,V,C]=quiver3Dpatch(V(:,1),V(:,2),V(:,3),U(:,1),U(:,2),U(:,3),[],a);

if isempty(colorOpt) %If empty use colormapping
    h=gpatch(F,V,C,'none',alphaLevel);
else %else use specified which could be 'k'
    h=gpatch(F,V,colorOpt,'none',alphaLevel);
end

varargout{1}=h;

%%
% _*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 [mask3MU, zValues] = getMask3D(structNum,planC)
%function [mask3MU, zValues] = getMask3D(structNum,planC)
%Assemble the 3-D CT-registered mask (type UINT8).
%JOD.
%Latest modifications:  JOD, 19 Feb 03, added zValues output.
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

indexS = planC{end};

scanNum = getStructureAssociatedScan(structNum, planC);

ROIImageSize   = [planC{indexS.scan}(scanNum).scanInfo(1).sizeOfDimension1  planC{indexS.scan}(scanNum).scanInfo(1).sizeOfDimension2];

numSlices = length(planC{indexS.scan}(scanNum).scanInfo);

mask3MU = false(ROIImageSize(1),ROIImageSize(2),numSlices);

zValues = [];

for sliceNum = 1 : numSlices

  z    = planC{indexS.scan}(scanNum).scanInfo(sliceNum).zValue;

  zValues = [zValues, z];

  [segmentsM, planC, isError] = getRasterSegments(structNum, planC);
%  segmentsM = planC{indexS.structures}(structNum).rasterSegments;

  indV = find(segmentsM(:,1) == z);  %mask values on this slice

  segmentsM = segmentsM(indV(:),7:9);     %segments

  %reconstruct the mask:

  for j = 1 : size(segmentsM,1)
    mask3MU(segmentsM(j,1),segmentsM(j,2):segmentsM(j,3),sliceNum) = 1;
  end

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function K = covPERiso(cov, hyp, x, z, i)

% Stationary periodic covariance function for an isotropic stationary covariance
% function k0 such as covMaterniso, covPPiso, covRQiso and covSEiso.
% Isotropic stationary means that the covariance function k0(x,z) depends on the
% data points x,z only through the squared distance
% dxz = (x-z)'*inv(P)*(x-z) where the P matrix is ell^2 times the unit matrix.
% The covariance function is parameterized as:
%
% k(x,z) = k0(u(x),u(z)), u(x) = [sin(pi*x/p); cos(pi*x/p)]
%
% where the period p belongs to covPERiso and hyp0 belong to k0:
%
% hyp = [ log(p)
%         hyp0 ]
%
% The first hyperparameter of k0 is the log lengthscale hyp0(1) = log(ell).
% Note that for k0 = covSEiso and D = 1, a faster alternative is covPeriodic.
%
% Copyright (c) by Hannes Nickisch, 2013-10-15.
%
% See also COVFUNCTIONS.M.

nocov = false;                   % default case when no cov argument is provided
if nargin==0, cov = {@covSEiso}; nocov = true; end                % default case
if isnumeric(cov)       % detect old version where the cov parameter was missing
  % i <- z, z <- x, x <- hyp, hyp <- cov
  if nargin>3, i = z; end
  if nargin>2, z = x; end
  if nargin>1, x = hyp; end
  hyp = cov; cov = {@covSEiso}; nocov = true;
end

if nocov && nargin<2 || ~nocov && nargin<3         % report number of parameters
  K = ['(1+',feval(cov{:}),')']; return
end
if nocov && nargin<3 || ~nocov && nargin<4, z = []; end    % make sure, z exists
xeqz = isempty(z); dg = strcmp(z,'diag');                       % determine mode

[n,D] = size(x);
p = exp(hyp(1));

if nocov && nargin<4 || ~nocov && nargin<5
  [x,z] = u(x,z,p,dg);                      % apply the embedding u:IR^D->IR^2*D
  K = feval(cov{:},hyp(2:end),x,z);
else
  if i==1
    if dg                                                   % compute distance d
      d = zeros([n,1,D]);
    else
      if xeqz                                             % symmetric matrix Kxx
        d = repmat(reshape(x,n,1,D),[1,n, 1])-repmat(reshape(x,1,n, D),[n,1,1]);
      else                                               % cross covariances Kxz
        nz = size(z,1);
        d = repmat(reshape(x,n,1,D),[1,nz,1])-repmat(reshape(z,1,nz,D),[n,1,1]);
      end
    end
    d = 2*pi*d/p; dD2_dlp = -2*sum(sin(d).*d,3);        % derivative dD2/dlog(p)
    [x,z] = u(x,z,p,dg);                    % apply the embedding u:IR^D->IR^2*D
    if dg                                            % compute squared distances
      D2 = zeros(n,1);
    else
      if xeqz, D2 = sq_dist(x'); else D2 = sq_dist(x',z'); end
    end
    % reconstruct derivative w.r.t. D2 from derivative w.r.t. log(ell)
    dK_dD2 = feval(cov{:},hyp(2:end),x,z,1)./(-2*D2); dK_dD2(D2<1e-12) = 0;
    K = dK_dD2.*dD2_dlp;                                      % apply chain rule
  else
    [x,z] = u(x,z,p,dg);                    % apply the embedding u:IR^D->IR^2*D
    K = feval(cov{:},hyp(2:end),x,z,i-1);
  end
end

function [x,z] = u(x,z,p,dg)                % apply the embedding u:IR^D->IR^2*D
  x = 2*pi*x/p; x = [sin(x), cos(x)];
  if numel(z)>0 && ~dg, z = 2*pi*z/p; z = [sin(z), cos(z)]; end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function UpdateCondPopupmenu(handles)
global maingui

if isempty(handles)
    return;
end
if isempty(maingui)
    return;
end
if maingui.dataTree.IsEmpty()
    return;
end

iG = maingui.dataTree.GetCurrElemIndexID();
CondNames = maingui.dataTree.groups(iG).GetConditions();
CondNamesCurrElem = maingui.dataTree.currElem.GetConditionsActive();
for jj=1:length(CondNames)
    k = find(strcmp(['-- ', CondNames{jj}], CondNamesCurrElem));
    if ~isempty(k)
        CondNames{jj} = ['-- ', CondNames{jj}];
    end
end
if isempty(CondNames)
   CondNames = {' '};
   set(handles.popupmenuConditions, 'enable', 'off');
else
    set(handles.popupmenuConditions, 'enable', 'on');
end
set(handles.popupmenuConditions, 'string', CondNames);
maingui.condition = get(handles.popupmenuConditions, 'value');

author: BUNPC


text:
%-------------------------------------------------------------------------------------------------------------------%
%
% IB2d is an Immersed Boundary Code (IB) for solving fully coupled non-linear
% 	fluid-structure eraction models. This version of the code is based off of
%	Peskin's Immersed Boundary Method Paper in Acta Numerica, 2002.
%
% Author: Nicholas A. Battista
% Email:  nick.battista@unc.edu
% Date Created: May 27th, 2015
% Institution: UNC-CH
%
% This code is capable of creating Lagrangian Structures using:
% 	1. Springs
% 	2. Beams (*torsional springs)
% 	3. Target Pos
%	4. Muscle-Model (combined Force-Length-Velocity model, "Hill+(Length-Tension)")
%
% One is able to update those Lagrangian Structure parameters, e.g., spring ants, resting lengths, etc
%
% There are a number of built in Examples, mostly used for teaching purposes.
%
% If you would like us to add a specific muscle model, please let Nick (nick.battista@unc.edu) know.
%
%--------------------------------------------------------------------------------------------------------------------%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: updates the beam attributes!
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function beams_info = update_nonInv_Beams(dt,current_time,beams_info)

% beams_info:   col 1: 1ST PT.
%               col 2: MIDDLE PT. (where force is exerted)
%               col 3: 3RD PT.
%               col 4: beam stiffness
%               col 5: x-curvature
%               col 6: y-curvature

% GEOMETRIC PARAMETERS
pi = 4*atan(1);
L1 = 8;                              % length of computational domain (m)
N1 = 256;                            % number of Cartesian grid meshwidths at the finest level of the AMR grid
bell_length = 2;                     % bell length (m)
npts_bell = ceil(2.0*(bell_length/L1)*N1);  % number of pos along the length of the bell
npts_circ = 1;                       %number of pos along the circumference (if in 3D)
npts = npts_bell*npts_circ;	         % total number pos
ds1 = bell_length/(npts_bell-1);     % mesh spacing(m) along length of bell

% Values from Alben, Peng, and Miller
betao = 0.5;
betam = 0.3;
to = 0.4;
Zs = L1/8;
xRef = L1/4;

 %These are used to keep track of cycle number and time into the cycle
pulse_time = current_time-floor(current_time); % determine time since beginning of first pulse

    % GIVE BELL STATES
    if (pulse_time<to)  %contract bell
        beta = betao+(betam-betao)*(pulse_time/to);
    else                % expand bell
        beta = betam+(betao-betam)*((pulse_time-to)/(1-to));
    end

    %---------------------
    %Xb and Yb are calculated here and will be used to determine new curvatures
    s1=0;
    zl = Zs;
    ro = xRef;

    % Pre-allocate memory for speed
    Xb_lam = zeros(npts,1);
    Yb_lam = zeros(npts,1);

    %top of bell
    Xb_lam(1)=ro;
    Yb_lam(1)=zl;

    %right side of bell
    for s1 = 2:ceil(npts_bell/2)
        thetaj = -1.55*(1-exp(-(s1)*ds1/beta));
        zl = zl + ds1*sin(thetaj);
        ro = ro + ds1*cos(thetaj);
        Xb_lam(s1)=ro;
        Yb_lam(s1)=zl;
    end

    zl = Zs;
    ro = xRef;

    %left side of bell
    for s1 = (ceil(npts_bell/2))+1:npts_bell
        s2=s1-(ceil(npts_bell/2));
        thetaj = -1.55*(1-exp(-(s2)*ds1/beta));
        zl = zl + ds1*sin(thetaj);
        ro = ro - ds1*cos(thetaj);
        Xb_lam(s1)=ro;
        Yb_lam(s1)=zl;
    end

% CHANGE CURVATURE BTWN SIDES OF JELLYFISH BELL
for s = 1:length( beams_info(:,1) )

    s1 = beams_info(s,2); % gives lag_idx of the middle node

    %top of bell
    if (s1==1)
         beams_info(s,5) = Xb_lam(ceil(npts/2)+1)+Xb_lam(2)-2*Xb_lam(s1);
         beams_info(s,6) = Yb_lam(ceil(npts/2)+1)+Yb_lam(2)-2*Yb_lam(s1);
    elseif ( s1<ceil(npts/2) )
    % right side of bell
        beams_info(s,5) = Xb_lam(s1-1)+Xb_lam(s1+1)-2*Xb_lam(s1);
        beams_info(s,6) = Yb_lam(s1-1)+Yb_lam(s1+1)-2*Yb_lam(s1);
    elseif ( s1==ceil(npts/2)+1 )
    % top-left side of bell
         beams_info(s,5) = Xb_lam(s1+1)+Xb_lam(1)-2*Xb_lam(s1);
         beams_info(s,6) = Yb_lam(s1+1)+Yb_lam(1)-2*Yb_lam(s1);
    elseif (s1 > ceil(npts/2)+1)
    % left side of bell
         beams_info(s,5) = Xb_lam(s1+1)+Xb_lam(s1-1)-2*Xb_lam(s1);
         beams_info(s,6) = Yb_lam(s1+1)+Yb_lam(s1-1)-2*Yb_lam(s1);
    end

end

author: nickabattista


text:
function q = fromTransfMatrixToPosQuat(H)

    % FROMTRANSFMATRIXTOQUATERNIONS computes a link pose using position
    %                               + quaternion rapresentation. The input is
    %                               the link pose represented by a transformation matrix.
    %
    % FORMAT: q = fromTransfMatrixToPosQuat(H)
    %
    % INPUT:  - H = [4 * 4] transf matrix representing link pose
    %
    % OUTPUT: - q = [7 * 1] position + quaternions representing link pose
    %
    % Authors: Daniele Pucci, Marie Charbonneau, Gabriele Nava
    %
    %          all authors are with the Italian Istitute of Technology (IIT)
    %          email: name.surname@iit.it
    %
    % Genoa, Dec 2017
    %

    %% --- Initialization ---

    % Separate the rotation matrix
    R       = H(1:3,1:3);
    qt_b    = zeros(4,1);

    % min. value to treat a number as zero
    epsilon = 1e-12;

    % Compute the corresponding (unit) quaternion from a given rotation matrix R:
    %
    % The transformation uses the computational efficient algorithm of Stanley.
    % To be numerically robust, the code determines the set with the maximum
    % divisor for the calculation.
    %
    % For further details about the Stanley Algorithm, see:
    %   [1] Optimal Spacecraft Rotational Maneuvers, John L. Junkins & James D. Turner, Elsevier, 1986, pp. 28-29, eq. (2.57)-(2.59).
    %   [2] Theory of Applied Robotics: Kinematics, Dynamics, and Control, Reza N. Jazar, 2nd Edition, Springer, 2010, p. 110, eq. (3.149)-(3.152).
    % Note: There exist also an optimized version of the Stanley method and is the fastest
    %       possible computation method for Matlab, but it does not cover all special cases.
    % Further details about the fast calculation can be found at:
    %   [3] Modelling and Control of Robot Manipulators, L. Sciavicco & B. Siciliano, 2nd Edition, Springer, 2008,
    %       p. 36, formula (2.30).
    %
    tr = R(1,1) + R(2,2) + R(3,3);

    if (tr > epsilon)

            % scalar part:
            qt_b(1,1) = 0.5*sqrt(tr + 1);
            s_inv     = 1/(qt_b(1,1)*4);

            % vector part:
            qt_b(2,1) = (R(3,2) - R(2,3))*s_inv;
            qt_b(3,1) = (R(1,3) - R(3,1))*s_inv;
            qt_b(4,1) = (R(2,1) - R(1,2))*s_inv;
    else
        % if tr <= 0, find the greatest diagonal element for calculating
        % the scale factor s and the vector part of the quaternion
        if ((R(1,1) > R(2,2)) && (R(1,1) > R(3,3)))

            qt_b(2,1) = 0.5*sqrt(R(1,1) - R(2,2) - R(3,3) + 1);
            s_inv     = 1/(qt_b(2,1)*4);

            qt_b(1,1) = (R(3,2) + R(2,3))*s_inv;
            qt_b(3,1) = (R(2,1) + R(1,2))*s_inv;
            qt_b(4,1) = (R(3,1) + R(1,3))*s_inv;

        elseif (R(2,2) > R(3,3))

            qt_b(3,1) = 0.5*sqrt(R(2,2) - R(3,3) - R(1,1) + 1);
            s_inv     = 1/(qt_b(3,1)*4);

            qt_b(1,1) = (R(1,3) - R(3,1))*s_inv;
            qt_b(2,1) = (R(2,1) + R(1,2))*s_inv;
            qt_b(4,1) = (R(3,2) + R(2,3))*s_inv;
        else
            qt_b(4,1) = 0.5*sqrt(R(3,3) - R(1,1) - R(2,2) + 1);
            s_inv     = 1/(qt_b(4,1)*4);

            qt_b(1,1) = (R(2,1) - R(1,2))*s_inv;
            qt_b(2,1) = (R(3,1) + R(1,3))*s_inv;
            qt_b(3,1) = (R(3,2) + R(2,3))*s_inv;
        end
    end

    % Final state converted into quaternion rapresentation
    q = [H(1:3,4); qt_b];
end

author: robotology

repo: whole-body-controllers


text:
function [ step_init ] = bb_init(problem, w)
%
% Barzilai-Borwein step-size initialization:
%
%
% This file is part of GDLibrary.
%
% This file originally comes from https://github.com/bodono/apg.
% Modifeid by H.Kasai on Apr. 17, 2017

    grad = problem.full_grad(w);
    step = 1 / norm(grad);
    w_hat = w - step*grad;
    grad_hat = problem.full_grad(w_hat);

    s = w - w_hat;
    y = grad - grad_hat;
    step_init = abs(s'*y / (norm(y)^2));
end

author: hiroyuki-kasai

repo: SGDLibrary


text:
function [phi] = evalPhi( r );
%
% [phi] = evalPhi( r );
%
%  Calculates phi( r ), phi the discrete delta function (without scaling)
%
%  Returns:
%     phi = value of unscaled discrete delta function at each point, r
%
%  Input:
%     r   = the points to evaluate the discrete delta function at
%
%  Note:
%     This code doesn't handle wrapping across periodic boundaries.
%
%
%  License: This code is free to use for any purposes, provided
%           any publications resulting from the use of this code
%           reference the original code/author.
%
%  Author:  Samuel Isaacson (isaacson@math.utah.edu)
%  Date:    11/2007
%
%  Please notify the author of any bugs, and contribute any
%  modifications or bug fixes back to the original author.
%
%  Disclaimer:
%   This code is provided as is. The author takes no responsibility
%   for its results or effects.

phi  = zeros(length(r), 1);
aR   = abs( r );
idx1 = find( aR < 1 );
idx2 = find( aR >= 1 );

rT        = aR(idx1);
phi(idx1) = 3 - 2*rT + sqrt( 1 + 4*rT - 4*rT.*rT );

rT        = aR(idx2);
phi(idx2) = 5 - 2*rT - sqrt( -7 + 12*rT - 4*rT.*rT );

phi = phi / 8;

repo: Aero-matlab


text:
% ismember_bc2 - ismember backward compatible with Matlab versions prior to 2013a

function [C,IA] = ismember_bc2(A,B,varargin)
% errorFlag = false;
v = version;
indp = find(v == '.');
v = str2num(v(1:indp(2)-1));
if v > 7.19
        v = floor(v) + rem(v,1)/10;
end
if nargin > 2
        ind = find(ismember(varargin, 'legacy'));
        if ~isempty(ind)
                varargin(ind) = [];
        end
end
if v >= 7.14
        %disp('legacy')
        [C,IA] = ismember(A,B,varargin{:},'legacy');
        %if true
        %        [C2,IA2] = ismember(A,B,varargin{:});
        %        if (~isequal(C, C2) || ~isequal(IA, IA2))
        %                warning('Warn:BckwrdcompTest','ismember_bc2: Backward compatibility issue with call to ismember function');
        %        end
        %end
else
        %disp('old')
        [C,IA] = ismember(A,B,varargin{:});
end

author: ucdavis

repo: erplab

save_path: github-repos/MATLAB/ucdavis-erplab


text:
function rank = subset_colex_rank ( n, t )

%*****************************************************************************80
%
%% SUBSET_COLEX_RANK computes the colexicographic rank of a subset.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 August 2011
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Donald Kreher, Douglas Simpson,
%    Combinatorial Algorithms,
%    CRC Press, 1998,
%    ISBN: 0-8493-3988-X,
%    LC: QA164.K73.
%
%  Parameters:
%
%    Input, integer N, the number of items in the master set.
%    N must be positive.
%
%    Input, integer T(N), the subset.  If T(I) = 0, item I is
%    not in the subset; if T(I) = 1, item I is in the subset.
%
%    Output, integer RANK, the rank of the subset.
%

%
%  Check.
%
  subset_check ( n, t );

  rank = 0;

  for i = 1 : n

    if ( t(i) == 1 )
      rank = rank + 2^( i - 1 );
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function varargout = ft_transform_vol(varargin)

% This function is a backward compatibility wrapper for existing MATLAB scripts
% that call a function that is not part of the FieldTrip toolbox any more.
%
% Please update your code to make it future-proof.

oldname = mfilename;
newname = 'ft_transform_geometry';

ft_warning('%s is only a backward compatibility wrapper, which will soon be removed. Please call %s instead.', upper(oldname), upper(newname));

[varargout{1:nargout}] = feval(newname, varargin{:});

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function mia_Ypixval(arg1, arg2)
% mia_Ypixval display information about image pixels in coronal image
% in the case of 3Dcursor mode
%
% Matlab library function for mia_gui utility.
% University of Debrecen, PET Center/LB 2003

%   PIXVAL ON turns on interactive display of information about image pixels
%   in the current figure. PIXVAL installs a black bar at the bottom of the
%   figure which displays the (x,y) coordinates for whatever pixel the
%   cursor is currently over, and the color information for that pixel. If
%   the image is binary or intensity, the color information is a single
%   intensity value. If the image is indexed or RGB, the color information
%   is an RGB triplet. The values displayed are the actual data values,
%   regardless of the class of the image array, or whether the data is in
%   normal image range.
%
%   If you click on the image and hold down the mouse button while you move
%   the cursor, PIXVAL also displays the Euclidean distance between the
%   point you clicked on and current cursor location. PIXVAL draws a line
%   between these points to indicate the distance being measured. When you
%   release the mouse button, the line and the distance display disappear.
%
%   You can move the display bar by clicking on it and dragging it to
%   another place in the figure.
%
%   PIXVAL OFF turns interactive display off in the current figure. You can
%   also turn off the display by clicking the button on the right side of
%   the display bar.
%
%   PIXVAL toggles interactive display on or off in the current figure.
%
%   PIXVAL(FIG,OPTION) applies the PIXVAL command to the figure specified
%   by FIG. OPTION is string containing 'on' or 'off'.
%
%   See also IMPIXEL, IMPROFILE.

%   Copyright 1993-2001 The MathWorks, Inc.
%   $Revision: 1.19 $  $Date: 2001/01/18 15:30:19 $

error(nargchk(0,2,nargin));

if nargin==0
   arg1 = gcf;
end

if ischar(arg1) % execute a callback
   switch lower(arg1)
   case 'pixvalmotionfcn'
      PixvalMotionFcn;
   case 'buttondownonimage'
      ButtonDownOnImage;
   case 'backtonormalpixvaldisplay'
      BackToNormalPixvalDisplay;
   case 'movedisplaybar'
      MoveDisplayBar;
   case 'on'
      mia_Ypixval(gcf, 'on');
   case 'off'
      mia_Ypixval(gcf, 'off');
   case 'deletedisplaybar'
      DeleteDisplayBar;
   otherwise
      error(['Unknown PIXVAL command or callback: ' arg1 '.']);
   end

elseif ishandle(arg1)
   images = findobj(arg1, 'type', 'image');
   if isempty(images)
      return
   end

   figureHandle = get(get(images(1), 'parent'), 'parent');

   if nargin>1   % mia_Ypixval(handle, 'on|off')
      if ischar(arg2)
         if strcmp(lower(arg2), 'off')
            action = 'destroy';
         elseif strcmp(lower(arg2), 'on')
            action = 'create';
         else
            error(['Unknown PIXVAL input string: ' arg2]);
         end
      else
         error('Invalid input arguments');
      end
   else  % mia_Ypixval(handle) or mia_Ypixval
      % Toggle the state
      displayBar = findobj(figureHandle, 'tag', 'pixel value display bar');
      if isempty(displayBar)
         action = 'create';
      else
         action = 'destroy';
      end
   end

   % Put things in the Images' UserData so that pixel value display works
   switch action
   case 'destroy'
      displayBar = findobj(figureHandle, 'tag', 'pixel value display bar');
      if ~isempty(displayBar)
         DeleteDisplayBar(displayBar);
      end

   case 'create'
      otherDisplayBar = findobj(figureHandle, 'tag', 'pixel value display bar');
      if ~isempty(otherDisplayBar)
         % mia_Ypixval is already ON for this figure.
         return
      end

      % Save the interactive state of the figure. UICLEARMODE now
      % does everything that UISUSPEND used to do. It calls
      % SCRIBECLEARMODE to take care of SCRIBE.  We WILL be
      % restoring the state - this is a change.  UICLEARMODE also
      % registers a way to disable mia_Ypixval
      dbud.oldFigureUIState = uiclearmode(figureHandle,'mia_Ypixval',figureHandle,'off');
      AxisHandle = findobj(figureHandle, 'type', 'axes');
      dbud.xLim = get(AxisHandle,'XLim');% for zoom mode
      dbud.yLim = get(AxisHandle,'YLim');% for zoom mode
      dbud.old_pt = 0;% for zoom mode

      % Make sure the SCRIBE toolbar doesn't go away
      plotedit(figureHandle,'locktoolbarvisibility');

      % Make sure that the Motion function doesn't interrupt itself
      set(figureHandle, 'Interruptible', 'off', 'busyaction', 'queue');

      % Make sure that the Image's Button Down & Up functions will queue
      set(images, 'Interruptible', 'off', 'BusyAction', 'Queue');

      % dbud is the UserData structure which is going into the display bar
      dbud.figureHandle = figureHandle;
      dbud.displayMode = 'normal';

      figPos = get(dbud.figureHandle, 'position');
      tempStr = sprintf(' %g, %g = %6.4f,%6.4f,%6.4f  dist = %3.3f', ...
        -0.9999,-0.9999,pi/10,pi/10,pi/10,pi*100);

      % Calculate the position vector of the display bar
      tempTxt = uicontrol('Parent', dbud.figureHandle, 'horiz', 'left', ...
         'Style', 'text', 'Units', 'pixels', 'visible', 'off', 'fontname', ...
         'FixedWidth','String', tempStr, 'position', [0 0 figPos(3) 20]);
      extent= get(tempTxt, 'extent');
      txtHt = extent(4);
      txtWid = extent(3)+txtHt; % Leave room for the close button
      delete(tempTxt);

      % Create position vectors for Text, Frame, and Close Button
      pos = [0 0 (min(txtWid,figPos(3))-txtHt)*1/2 txtHt];
      btnSide = pos(4)-4;  % In pixels
      btnPos = [pos(1)+pos(3)+2 pos(2)+2 btnSide btnSide];
      frmPos = [pos(1)+pos(3) 0 txtHt txtHt];

      % Create a frame for the display bar. I do this because of a bug in
      % MATLAB 5.2 where the text field pops to the front of the button each
      % time it is updated.  It should stay behind in which case, all we would
      % need is the text uicontrol and the button.  We can take the frame out
      % if the bug gets fixed.  Make sure the frame doesn't fall behind the
      % text control, otherwise you can't execute the text control's buttondown
      % function (another MATLAB 5.2 bug).

      dbud.DisplayFrame = uicontrol('Parent', figureHandle, ...
         'style', 'frame', ...
         'Units', 'pixels', ...
         'Position', frmPos, ...
         'Background', [0 0 0], ...
         'BusyAction', 'queue', ...
         'Interruptible', 'off');

      % Create the display bar
      dbud.buttonDownImage = 0;  % Image 0 never exists
      DisplayBar = uicontrol('Parent', figureHandle, ...
         'Style','text', ...
         'Units','pixels', ...
         'Position',pos, ...
         'Foreground', [1 1 .5], ...
         'Background', [0 0 0], ...
         'Horiz','left', ...
         'Tag', 'pixel value display bar', ...
         'String','', ...
         'fontname', 'FixedWidth', ...
         'BusyAction', 'queue', ...
         'enable', 'inactive', ...
         'Interruptible', 'off');
%          'ButtonDownFcn', 'mia_Ypixval(''MoveDisplayBar'')', ...
%          'DeleteFcn', 'mia_Ypixval(''DeleteDisplayBar'')', ...

      % Create the close button
      dbud.CloseButton = uicontrol('Parent', figureHandle, ...
         'style', 'pushbutton', ...
         'units', 'pixels', ...
         'position', btnPos, ...
         'BusyAction', 'queue', ...
         'Interruptible', 'off');
%          'String', 'X', ...
%          'Callback', 'mia_Ypixval(''DeleteDisplayBar'')', ...
%          'DeleteFcn', 'mia_Ypixval(''DeleteDisplayBar'')', ...

      % Register the motion function and button up function
      set(figureHandle, 'WindowButtonMotionFcn',...
         'mia_Ypixval(''PixvalMotionFcn'')')

      set(DisplayBar, 'UserData', dbud);
      PixvalMotionFcn(DisplayBar);
   end

else
   error('Invalid input argument.')
end

%%%
%%%  Sub-function - PixvalMotionFcn
%%%

function PixvalMotionFcn(displayBar)

if nargin < 1
   displayBar = findobj(gcbf, 'Tag', 'pixel value display bar');
end

if isempty(displayBar)
   % Pixval is in a half-broken state.  Another function (like zoom) has
   % restored pixval callbacks and PIXVAL has already been uninstalled.
   % Call uisuspend to gracefully get rid of the callbacks.

   % Note 7/21/98 - Since we are now using UICLEARMODE instead of
   % UISUSPEND, I think that we should never get into this
   % state.  It will only happen if a user writes a function
   % which calls UIRESTORE without ever calling UICLEARMDOE or
   % SCRIBECLEARMODE.  We'll leave this code here just to be safe.
   uisuspend(gcbf);
   return
end

dbud = get(displayBar, 'UserData');
pt = get(dbud.figureHandle, 'CurrentPoint');
x = pt(1,1); y = pt(1,2);

if strcmp(dbud.displayMode, 'normal')
   % See if we're above the displayBar
   dbpos = get(displayBar, 'Position');
   left   = dbpos(1);
   bottom = dbpos(2);
   width  = dbpos(3);
   height = dbpos(4);
   % On the line below, we look to see if x<=left+width+height because
   % there is a frame next to the text object, see note above about the
   % MATLAB 5.2 bug we're working around
   if x>=left & x<=left+width+height & y>=bottom & y<=bottom+height
      % We're hovering above the display bar

      % Look to see if we are above the Close button
      cp = get(dbud.CloseButton, 'Position');
      if  x>=cp(1) & x<=cp(1)+cp(3) & y>=cp(2) & y<=cp(2)+cp(4)
         if isfield(dbud, 'oldPointer') % Switch back to the default pointer
            set(dbud.figureHandle, 'Pointer', dbud.oldPointer, ...
               'PointerShapeCData', dbud.oldPointerShapeCData, ...
               'PointerShapeHotSpot', dbud.oldPointerShapeHotSpot);
            dbud = rmfield(dbud, {'oldPointer', ...
                  'oldPointerShapeCData', 'oldPointerShapeHotSpot'});
         end
      else
         % Save the default pointer if there is one
         if ~isfield(dbud, 'oldPointer') % We still have the default pointer
            dbud.oldPointer = get(dbud.figureHandle, 'Pointer');
            dbud.oldPointerShapeCData = get(dbud.figureHandle, 'PointerShapeCData');
            dbud.oldPointerShapeHotSpot = get(dbud.figureHandle, 'PointerShapeHotSpot');
         end
         setptr(dbud.figureHandle, 'hand');
      end
   else
      % We're not hovering over the display bar
      [imageHandle,imageType,img,x,y] = OverImage(dbud.figureHandle);
      if imageHandle ~= 0
         % Save old pointer and get a new one (use Custom pointer for non-Windows machine)
         if ~isfield(dbud, 'oldPointer') % We still have the default pointer
            dbud.oldPointer = get(dbud.figureHandle, 'Pointer');
            dbud.oldPointerShapeCData = get(dbud.figureHandle, 'PointerShapeCData');
            dbud.oldPointerShapeHotSpot = get(dbud.figureHandle, 'PointerShapeHotSpot');
         end
         if strcmp(computer, 'PCWIN')
            set(dbud.figureHandle, 'Pointer', 'cross');
         else
            [pointerShape, pointerHotSpot] = CreatePointer;
            set(dbud.figureHandle, 'Pointer', 'custom', ...
               'PointerShapeCData', pointerShape, ...
               'PointerShapeHotSpot', pointerHotSpot);
         end

         % Update the Pixel Value display
         UpdatePixelValues(imageHandle, imageType, displayBar,img,x,y );
      else
         if isfield(dbud, 'oldPointer') % Switch back to the default pointer
            set(dbud.figureHandle, 'Pointer', dbud.oldPointer, ...
               'PointerShapeCData', dbud.oldPointerShapeCData, ...
               'PointerShapeHotSpot', dbud.oldPointerShapeHotSpot);
            dbud = rmfield(dbud, {'oldPointer', ...
                  'oldPointerShapeCData', 'oldPointerShapeHotSpot'});
            set(displayBar, 'String', '');
         end
      end
   end
elseif strcmp(dbud.displayMode, 'distance')
   % If we're in distance mode and in another image, clean up a bit.
   [imageHandle,imageType,img,x,y]= OverImage(dbud.figureHandle);
   if imageHandle~=0
      if imageHandle==dbud.buttonDownImage
         UpdatePixelValues(imageHandle, imageType, displayBar,img,x,y);
      end
   end
end
set(displayBar, 'UserData', dbud);

%%%
%%%  Sub-function - OverImage
%%%

function [imageHandle,imageType,img,x,y] = OverImage(figHandle)
% Return the index of which image we are over, and return a 0 if we
% aren't above an image.

images = findobj(figHandle, 'type', 'image');
if isempty(images)
   imageHandle=0; imageType=''; img=[]; x=0; y=0;
   return
end
% Make sure that the Image's Button Down & Up functions will queue
set(images, 'ButtonDownFcn', 'mia_Ypixval(''ButtonDownOnImage'')', ...
   'Interruptible', 'off', 'BusyAction', 'Queue');
axHandles = get(images, {'Parent'});
axPositions = get([axHandles{:}], {'Position'});
axCurPt = get([axHandles{:}], {'CurrentPoint'});

% Loop over the axes, see if we are above any of them
imageHandle = 0;
for k=1:length(axHandles)
   XLim = get(axHandles{k}, 'XLim');
   YLim = get(axHandles{k}, 'YLim');
   pt = axCurPt{k};
   x = pt(1,1); y = pt(1,2);
   if x>=XLim(1) & x<=XLim(2) & y>=YLim(1) & y<=YLim(2)
      imageHandle = images(k);
      break;
   end
end

% Figure out image type
if imageHandle ~= 0
   [img,flag] = getimage(imageHandle);

   switch flag
   case 1
      imageType = 'indexed';
   case {2,3}  % Grayscale or binary
      imageType = 'intensity';
   case 4
      imageType = 'rgb';
   otherwise
      error(['Invalid image, GETIMAGE returned flag = ' flag '.']);
   end
else
   imageHandle=0; imageType=''; img=[]; x=0; y=0;
end

%%%
%%%  Sub-function - UpdatePixelValues
%%%

function UpdatePixelValues(imageHandle, imageType, displayBar, img, x, y)
%   This is the motion callback for when we are displaying
%   pixels.  Either we are in automatic display mode and the
%   mouse pointer is moving around or we are in normal mode
%   and there has been a button-down but not yet a button up.
%   I get the current point and update the string.

global old_pt_zoom_y x0_pan_y;%temp for zoom and pan

dbud = get(displayBar, 'UserData');
[rows, cols, colors] = size(img);
y = min(y,rows);% y > rows happens in some case /LB
rp = axes2pix(rows, get(imageHandle,'YData'),y);
cp = axes2pix(cols, get(imageHandle,'XData'),x);
r = min(rows, max(1, round(rp)));
c = min(cols, max(1, round(cp)));

if strcmp(imageType,'indexed')
   map=get(dbud.figureHandle, 'Colormap');
   idx = img(r,c);
   if ~isa(idx, 'double'),
      idx = double(idx)+1;
   end
   idx = round(idx);
   if idx<=size(map,1)
      pixel = map(idx,:);
   else
      pixel = map(end,:);
   end
else
   pixel = double(img(r,c,:));
end

% figure out the new string
figureHandle = get(get(imageHandle, 'parent'), 'parent');
switch dbud.displayMode
case 'normal'   % Just display intensity information
	FigDataIn = get(figureHandle,'Userdata');
    ImVolume = get(FigDataIn.ImaHandlerX,'UserData');
    %pixel(1) = ImVolume.PixInt;

	if isa(img, 'uint8')
     pixval_str = sprintf(' %g,%g,%g = %3d',round(y),round(ImVolume.Y), round(x),pixel(1));
	elseif isa(img, 'int16')
     pixval_str = sprintf(' %g,%g,%g = %5d',round(y),round(ImVolume.Y), round(x), pixel(1));
	else
     pixval_str = sprintf(' %g,%g,%g = %5.1f',round(y),round(ImVolume.Y), round(x),pixel(1));
	end
    set(displayBar, 'String', pixval_str, 'UserData', dbud);
case 'distance'
    mouseclick = get(gcf,'SelectionType');
    if strcmp(mouseclick,'normal')
		FigDataIn = get(figureHandle,'Userdata');
		fighandlerZ = FigDataIn.ImaHandlerZ;
		fighandlerX = FigDataIn.ImaHandlerX;
		FigData = FigDataIn.CData;

        set(fighandlerZ,'CData',squeeze(FigData(:,:,size(FigData,3) - round(rp) + 1)) );
        set(fighandlerX,'CData',rot90(squeeze(FigData(:,round(cp),:))));

        ImVolume = get(FigDataIn.ImaHandlerX,'UserData');
        ImVolume.Z = y;
        ImVolume.X = x;
        ImVolume.PixInt = pixel(1);
        set(FigDataIn.ImaHandlerX,'UserData',ImVolume);

		if isa(img, 'uint8')
         pixval_str = sprintf(' %g,%g,%g = %3d',round(y),round(ImVolume.Y), round(x),pixel(1));
		elseif isa(img, 'int16')
         pixval_str = sprintf(' %g,%g,%g = %5d', round(y),round(ImVolume.Y), round(x), pixel(1));
		else
         pixval_str = sprintf(' %g,%g,%g = %5.1f',round(y),round(ImVolume.Y), round(x),pixel(1));
		end
        set(displayBar, 'String', pixval_str, 'UserData', dbud);
    elseif strcmp(mouseclick,'extend') % zoom mode
        CurrentAxes = findobj(figureHandle, 'type', 'axes');
        %set(CurrentAxes,'Units','pixel');
        new_pt = get(0,'PointerLocation');
        dy = (new_pt(2) - old_pt_zoom_y(2))/abs(old_pt_zoom_y(2))*5;
        zoomFactor = (1-dy);
        xLim=get(CurrentAxes,'XLim');
        yLim=get(CurrentAxes,'YLim');
     	xLimNew = [0 zoomFactor*diff(xLim)] + xLim(1) + (1-zoomFactor)*diff(xLim)/2;
        yLimNew = [0 zoomFactor*diff(yLim)] + yLim(1) + (1-zoomFactor)*diff(yLim)/2;
		set(CurrentAxes,'XLim',xLimNew,'YLim',yLimNew);
        old_pt_zoom_y  = new_pt;
    elseif strcmp(mouseclick,'alt') % pan mode
        currentPoint = get(figureHandle,'CurrentPoint');
        CurrentAxes = findobj(figureHandle, 'type', 'axes');
        set(CurrentAxes,'Units','pixel');
        dx = currentPoint - x0_pan_y;
        x0_pan_y = currentPoint;
        ap = get(CurrentAxes,'Position');
        xLim = get(CurrentAxes,'XLim');
        yLim = get(CurrentAxes,'YLim');
        set(CurrentAxes,'XLim',xLim-(diff(xLim)*dx(1)/ap(3)), ...
           'YLim',yLim+(diff(yLim)*dx(2)/ap(4)));
       set(CurrentAxes,'Units','normalized');
    elseif strcmp(mouseclick,'open')% restore the original image
        xLimNew = dbud.xLim;
        yLimNew = dbud.yLim;
        CurrentAxes = findobj(gcf, 'type', 'axes');
        set(CurrentAxes,'XLim',xLimNew,'YLim',yLimNew);
    end
end

%%%
%%%  Sub-function - ButtonDownOnImage
%%%

function ButtonDownOnImage

global old_pt_zoom_y x0_pan_y;%temp for zoom

%imageHandle = gcbo;
currentHandle = gcbo;
if strcmp(get(currentHandle,'type'),'image')
    imageHandle = currentHandle;
else
    currentParentHandle = get(currentHandle,'Parent');
    currentChildrensHandles = get(currentParentHandle,'children');
    imageHandle = currentChildrensHandles(3);
end
figureHandle = get(get(imageHandle,'Parent'),'Parent');
displayBar = findobj(figureHandle, 'Tag', 'pixel value display bar');
dbud = get(displayBar, 'UserData');
axesHandle = get(imageHandle, 'Parent');

mouseclick = get(gcf,'SelectionType');
if strcmp(mouseclick,'normal')% case of distance measuring
	%
elseif strcmp(mouseclick,'alt') % pan mode
    setptr(figureHandle,'closedhand');
elseif strcmp(mouseclick,'extend') % zoom mode
    setptr(figureHandle,'glass');
end

%
% Turn on the PlotEdit options on X,Y slice images
%
FigDataIn = get(figureHandle,'Userdata');

%delete the positioning line
ImVolume = get(FigDataIn.ImaHandlerX,'UserData');
if ishandle(ImVolume.Xxline)
    delete(ImVolume.Xxline);
    delete(ImVolume.Xyline);
    delete(ImVolume.Yxline);
    delete(ImVolume.Yyline);
    delete(ImVolume.Zxline);
    delete(ImVolume.Zyline);
end
set(FigDataIn.ImaHandlerX,'UserData',ImVolume);

% Set the initial point (x0,y0)

%pt = get(axesHandle, 'CurrentPoint');
%dbud.x0 = pt(1,1);
%dbud.y0 = pt(1,2);
%dbud.line = line('Parent', axesHandle, ...
%  'erasemode', 'xor', ...
%   'color', [1 0 0], ...
%   'Xdata', [dbud.x0 dbud.x0], ...
%   'Ydata', [dbud.y0 dbud.y0]);
dbud.displayMode = 'distance';
old_pt_zoom_y = get(0,'PointerLocation');% for zoom mode
x0_pan_y = get(figureHandle,'CurrentPoint');% for pan mood
dbud.buttonDownImage = imageHandle;
set(displayBar, 'UserData', dbud);
set(dbud.figureHandle, 'WindowButtonUpFcn', 'mia_Ypixval(''BackToNormalPixvalDisplay'')');
PixvalMotionFcn(displayBar);

%%%
%%%  Sub-function - BackToNormalPixvalDisplay
%%%

function BackToNormalPixvalDisplay

displayBar = findobj(gcbo, 'Tag', 'pixel value display bar');
dbud = get(displayBar, 'UserData');
imageHandle = dbud.buttonDownImage;
%
% Turn off the PlotEdit options on X,Y slice images and
% write the last X,Y,Z coordinates to the FigDataIn struct
%
[imageHandleTmp,imageType,img,x,y]= OverImage(dbud.figureHandle);
figureHandle = get(get(imageHandle,'Parent'),'Parent');
FigDataIn = get(figureHandle,'Userdata');
axesHandleX = get(FigDataIn.ImaHandlerX, 'Parent');
axesHandleY = get(FigDataIn.ImaHandlerY, 'Parent');
axesHandleZ = get(FigDataIn.ImaHandlerZ, 'Parent');

%draw the positioning line
Xyrange = get(axesHandleX,'Ylim');
Xxrange = get(axesHandleX,'Xlim');
Yyrange = get(axesHandleY,'Ylim');
Yxrange = get(axesHandleY,'Xlim');
Zyrange = get(axesHandleZ,'Ylim');
Zxrange = get(axesHandleZ,'Xlim');

ImVolume = get(FigDataIn.ImaHandlerX,'UserData');
if ishandle(ImVolume.Xxline)
    delete(ImVolume.Xxline);
    delete(ImVolume.Xyline);
    delete(ImVolume.Yxline);
    delete(ImVolume.Yyline);
    delete(ImVolume.Zxline);
    delete(ImVolume.Zyline);
end

LineWidthCur = 2;
%lines in saggital slice
ImVolume.Xxline=line('Parent', axesHandleX,'color', [0 1 1],'EraseMode','xor', ....
    'LineWidth',LineWidthCur,'Xdata',[ImVolume.Y ImVolume.Y],'Ydata',[0 Xyrange(2)], ...
    'ButtonDownFcn','mia_Xpixval(''ButtonDownOnImage'')');
ImVolume.Xyline=line('Parent', axesHandleX,'color', [0 1 1],'EraseMode','xor', ...
'LineWidth',LineWidthCur,'Xdata',[0 Xxrange(2)],'Ydata', [ImVolume.Z ImVolume.Z], ...
    'ButtonDownFcn','mia_Xpixval(''ButtonDownOnImage'')');

%lines in axial slice
ImVolume.Zxline=line('Parent', axesHandleZ,'color', [0 1 1],'EraseMode','xor', ....
    'LineWidth',LineWidthCur,'Xdata',[ImVolume.X ImVolume.X],'Ydata',[0 Zyrange(2)], ...
    'ButtonDownFcn','mia_Zpixval(''ButtonDownOnImage'')');
ImVolume.Zyline=line('Parent', axesHandleZ,'color', [0 1 1],'EraseMode','xor', ...
    'LineWidth',LineWidthCur,'Xdata',[0 Zxrange(2)],'Ydata', [ImVolume.Y ImVolume.Y], ...
    'ButtonDownFcn','mia_Zpixval(''ButtonDownOnImage'')');

% set(FigDataIn.ImaHandlerX,'UserData',ImVolume);
%lines in coronal slice
ImVolume.Yxline=line('Parent', axesHandleY,'color', [0 1 1],'EraseMode','xor', ....
    'LineWidth',LineWidthCur,'Xdata',[ImVolume.X ImVolume.X],'Ydata',[0 Xyrange(2)], ...
    'ButtonDownFcn','mia_Ypixval(''ButtonDownOnImage'')');
ImVolume.Yyline=line('Parent', axesHandleY,'color', [0 1 1],'EraseMode','xor', ...
    'LineWidth',LineWidthCur,'Xdata',[0 Yxrange(2)],'Ydata', [ImVolume.Z ImVolume.Z], ...
    'ButtonDownFcn','mia_Ypixval(''ButtonDownOnImage'')');

set(FigDataIn.ImaHandlerX,'UserData',ImVolume);

%delete(dbud.line);
%dbud.line = [];
%dbud.x0 = []; dbud.y0 = [];
set(dbud.figureHandle, 'WindowButtonUpFcn', '');
dbud.displayMode = 'normal';
dbud.buttonDownImage = 0;
set(displayBar, 'UserData', dbud);
PixvalMotionFcn( displayBar);

%%%
%%%  Sub-function - MoveDisplayBar
%%%

function MoveDisplayBar

displayBar = gcbo;
dbud = get(displayBar, 'UserData');
oldWindowMotionFcn = get(dbud.figureHandle, 'WindowButtonMotionFcn');
set(dbud.figureHandle, 'WindowButtonMotionFcn', '');
oldPointer = get(dbud.figureHandle, 'Pointer');
oldPointerShapeCData = get(dbud.figureHandle, 'PointerShapeCData');
oldPointerShapeHotSpot = get(dbud.figureHandle, 'PointerShapeHotSpot');
setptr(dbud.figureHandle,'closedhand');

txtpos = get(displayBar, 'Position');
frmPos = get(dbud.DisplayFrame, 'Position');
position = txtpos + [0 0 frmPos(3) 0];
position2 = dragrect(position);
dx = position2(1)-position(1);
dy = position2(2)-position(2);

txtpos = txtpos + [dx dy 0 0];
set(displayBar, 'Position', txtpos);

frmPos = frmPos + [dx dy 0 0];
set(dbud.DisplayFrame, 'Position', frmPos);

btnPos = get(dbud.CloseButton, 'Position');
btnPos = btnPos + [dx dy 0 0];
set(dbud.CloseButton, 'Position', btnPos);

set(dbud.figureHandle, 'Pointer', oldPointer, ...
   'PointerShapeCData', oldPointerShapeCData, ...
   'PointerShapeHotSpot', oldPointerShapeHotSpot)
set(dbud.figureHandle, 'WindowButtonMotionFcn', oldWindowMotionFcn);

%%%
%%%  Sub-function - DeleteDisplayBar
%%%

function DeleteDisplayBar(displayBar)
%  Take apart the pixel value display - get rid of the text
%  uicontrol and clean up the image objects (remove their
%  UserData's and get rid of their ButtonDown & Delete Fcn's)

if nargin<1
   displayBar = findobj(gcbf, 'Tag', 'pixel value display bar');
end

if ~isempty(displayBar) % if it's empty, there's no work to do.
   dbud = get(displayBar, 'UserData');
   set(displayBar, 'DeleteFcn', '');  % Make sure there's no recursion
   if ishandle(displayBar)
      delete(displayBar);
   end
   if ishandle(dbud.CloseButton)
      delete(dbud.CloseButton);
   end
   if ishandle(dbud.DisplayFrame)
      delete(dbud.DisplayFrame);
   end

   % We are once again restoring the figure state, since
   % UIRESTORE also will reactivate the SCRIBE toolbar.
   uirestore(dbud.oldFigureUIState);

   % This is what we did for IPT 2.1/ML 5.2 :
   %   uisuspend(dbud.figureHandle);
end

%%%
%%%  Sub-function - CreatePointer
%%%

function [pointerShape, pointerHotSpot] = CreatePointer

pointerHotSpot = [8 8];
pointerShape = [ ...
      NaN NaN NaN NaN NaN   1   2 NaN   2   1 NaN NaN NaN NaN NaN NaN
   NaN NaN NaN NaN NaN   1   2 NaN   2   1 NaN NaN NaN NaN NaN NaN
   NaN NaN NaN NaN NaN   1   2 NaN   2   1 NaN NaN NaN NaN NaN NaN
   NaN NaN NaN NaN NaN   1   2 NaN   2   1 NaN NaN NaN NaN NaN NaN
   NaN NaN NaN NaN NaN   1   2 NaN   2   1 NaN NaN NaN NaN NaN NaN
   1   1   1   1   1   1   2 NaN   2   1   1   1   1   1   1   1
   2   2   2   2   2   2   2 NaN   2   2   2   2   2   2   2   2
   NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
   2   2   2   2   2   2   2 NaN   2   2   2   2   2   2   2   2
   1   1   1   1   1   1   2 NaN   2   1   1   1   1   1   1   1
   NaN NaN NaN NaN NaN   1   2 NaN   2   1 NaN NaN NaN NaN NaN NaN
   NaN NaN NaN NaN NaN   1   2 NaN   2   1 NaN NaN NaN NaN NaN NaN
   NaN NaN NaN NaN NaN   1   2 NaN   2   1 NaN NaN NaN NaN NaN NaN
   NaN NaN NaN NaN NaN   1   2 NaN   2   1 NaN NaN NaN NaN NaN NaN
   NaN NaN NaN NaN NaN   1   2 NaN   2   1 NaN NaN NaN NaN NaN NaN
   NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN];

author: Sable

repo: mcbench-benchmarks


text:
function [] = plotHammDist( jobIDs, taskIDs, jobNames, objIDs, doHourly )
% View trace plots of hamming distance for sampler state at each
%    recorded iteration of the sampler chain.  Useful for diagnosing
%    convergence and mixing rates.
% ________________________________________________________________________

if ~exist( 'objIDs','var')
    objIDs = 1:6;
end

figure;
set( gcf, 'Units', 'normalized', 'Position', [0.5 0.5 0.5 0.5] );
if exist( 'jobNames', 'var' ) && ~isempty( jobNames )
    plotColors  = get(0,'defaultAxesColorOrder');
    hold on;
else
    plotColors = get(0,'defaultAxesColorOrder');
    hold all;
end

for jobID = jobIDs
    taskNames = {};

    doFirstTask = 1;
    for taskID = taskIDs
        DATA = loadSamplerOutput( jobID, taskID , {'times', 'A'} );
        if isnumeric(DATA) && DATA == -1
            continue;
        end

        times = DATA.times.Psi;

        Ts = DATA.A(end).HDist.Ts( objIDs );
        Hdist = zeros( 1, length(DATA.A ) );
        for ss = 1:length( DATA.A )
            objDists = DATA.A(ss).HDist.obj( objIDs );
            Hdist(ss) = sum( objDists .* Ts )/sum(Ts);
        end

        if exist( 'jobNames', 'var' )  && ~isempty( jobNames )
            jj = 1 + mod( find( jobID == jobIDs )-1, size(plotColors,1) );
            curColor = plotColors(jj,:);
            if doFirstTask;
                taskVis = 'on';
                doFirstTask = 0;
            else
                taskVis = 'off';
            end
        else
            taskNames{end+1} = num2str( taskID );
            jj = 1 + mod( find( taskID == taskIDs )-1, size(plotColors,1) );
            curColor = plotColors(jj,:);
            taskVis = 'on';
        end

        styleStr = '.-';

        while length( times ) > 200
           times = times(1:2:end);
           Hdist = Hdist(1:2:end);
        end

        if doHourly
            times = times/3600;
        end

        plot( times, Hdist, styleStr, 'MarkerSize', 15, 'LineWidth', 2, 'HandleVisibility', taskVis, 'Color', curColor );

    end
end

if exist( 'jobNames', 'var' ) && ~isempty( jobNames )
    legend( jobNames , 'Location', 'SouthEast' );
else
    legend( taskNames, 'Location', 'SouthEast' );
end

ylabel( 'Hamming Dist', 'FontSize', 18 );

if doHourly
xlabel ('CPU time (hours)', 'FontSize', 18);
else
xlabel ('CPU time (sec)', 'FontSize', 18);
end
grid on;
set( gca, 'FontSize', 16 );
end % main function

author: michaelchughes

repo: NPBayesHMM


text:
classdef ProgressBar < handle

properties
    width
    back, msg
    i, N
end

methods

    % ===========================
    % Initialize the progress bar
    % ===========================
	function obj = ProgressBar( N )
        obj.width = 25;
        obj.i     = 1;
        obj.N     = N;
        % display empty bar
        obj.msg = [ '   [' repmat(' ',1,obj.width) '] ' ];
        fprintf([obj.back, obj.msg]);
        obj.back = repmat('\b',1,length(obj.msg));
    end

    % ==========================================================
    % Update the counter by 1 and print the updated progress bar
    % ==========================================================
    function obj = update( obj )
        if ( obj.i < 1 || obj.i > obj.N )
            return
        end

        if ( mod( obj.i, obj.width ) == 0 || obj.i==obj.N )
            p = floor(obj.i/obj.N*obj.width);
            obj.msg = [ '   [' repmat('=',1,p) repmat(' ',1,obj.width-p) '] ' ];
            fprintf([obj.back, obj.msg]);
            obj.back = repmat('\b',1,length(obj.msg));
        end
        obj.i = obj.i + 1;
    end

    % ======================
    % Close the progress bar
    % ======================
    function close( obj )
        fprintf('\n');
        obj.i    = 0;
        obj.N    = 0;
        obj.msg  = '';
        obj.back = '';
    end

end

end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function output = calllindo_nlp(interfacedata)

global MY_LICENSE_FILE

persistent iEnv

% Instead of calling lindo, we define the parameters we need. This is
% needed to speed up repeated calls
%lindo
LSERR_NO_ERROR                                               = 0000;
LS_IPARAM_NLP_PRINTLEVEL                                     = 203;
LS_IPARAM_NLP_SOLVER                                         = 201;
LS_IPARAM_NLP_MAXLOCALSEARCH                                 = 221;
LS_STATUS_OPTIMAL                                            = 1;
LS_STATUS_BASIC_OPTIMAL                                      = 2;
LS_STATUS_INFEASIBLE                                         = 3;
LS_STATUS_UNBOUNDED                                          = 4;
LS_STATUS_FEASIBLE                                           = 5;
LS_STATUS_INFORUNB                                           = 6;
LS_STATUS_NEAR_OPTIMAL                                       = 7;
LS_STATUS_LOCAL_OPTIMAL                                      = 8;
LS_STATUS_LOCAL_INFEASIBLE                                   = 9;
LS_STATUS_CUTOFF                                             = 10;
LS_STATUS_NUMERICAL_ERROR                                    = 11;
LS_STATUS_UNKNOWN                                            = 12;
LS_STATUS_UNLOADED                                           = 13;
LS_STATUS_LOADED                                             = 14;
LS_METHOD_FREE                                               = 0;
LS_METHOD_PSIMPLEX                                           = 1;
LS_METHOD_DSIMPLEX                                           = 2;
LS_METHOD_BARRIER                                            = 3;
LS_METHOD_NLP                                                = 4;
LS_NMETHOD_FREE                                              = 4;
LS_NMETHOD_CONOPT                                            = 7;
LS_NMETHOD_MSW_GRG                                           = 9;

if isempty(iEnv)
    % This call is mighty slow, so we do it only once, unless uses clears
    % everything
    [MY_LICENSE_KEY,Err] = mxlindo('LSloadLicenseString',MY_LICENSE_FILE);
    [iEnv,nErr]=mxlindo('LScreateEnv',MY_LICENSE_KEY);
    if nErr ~= LSERR_NO_ERROR;output = returnempty(-5); return; end;
end

% Retrieve needed data
options = interfacedata.options;
F_struc = interfacedata.F_struc;
c       = interfacedata.c;
K       = interfacedata.K;
x0      = interfacedata.x0;
Q       = interfacedata.Q;
lb      = interfacedata.lb;
ub      = interfacedata.ub;
monomtable = interfacedata.monomtable;

% Do some pre-calc to be used in callbacks
nonlinearindicies = find(interfacedata.variabletype~=0);
nonlinearindicies = unionstripped(nonlinearindicies,interfacedata.evalVariables);
linearindicies    = find(interfacedata.variabletype==0);
linearindicies    = setdiff1D(linearindicies,nonlinearindicies);
interfacedata.nonlinearindicies = nonlinearindicies;
interfacedata.linearindicies    = linearindicies;

% % Move nonlinear bounds to constraints
if ~isempty(lb)
    finite = find(~isinf(lb(nonlinearindicies)));
    if ~isempty(finite)
        temp = F_struc(1:K.f,:);
        F_struc(1:K.f,:) = [];
        n = length(c);

        for i = 1:length(finite)
            j = nonlinearindicies(i);
            F_struc = [-lb(j) sparse(1,j,1,1,n);F_struc];
        end
        K.l = K.l + length(finite);
        F_struc = [temp;F_struc];
        interfacedata.K = K;
        interfacedata.F_struc = F_struc;
    end
end
if ~isempty(ub)
    finite = find(~isinf(ub(nonlinearindicies)));
    if ~isempty(finite)
        temp = F_struc(1:K.f,:);
        F_struc(1:K.f,:) = [];
        n = length(c);

        for i = 1:length(finite)
            j = nonlinearindicies(i);
            F_struc = [ub(j) -sparse(1,j,1,1,n);F_struc];
        end
        K.l = K.l + length(finite);
        F_struc = [temp;F_struc];
        interfacedata.K = K;
        interfacedata.F_struc = F_struc;
    end
end
%
% Init model size
m  = K.l + K.f;
n  = length(c);

% Specifying variable types...
vtype = repmat('C',1,length(c(linearindicies)));
vtype(interfacedata.integer_variables) = 'I';

oshift = interfacedata.f;

if m == 0
    interfacedata.F_struc = [1e6 -ones(1,length(c))];
    K.l = 1;
    F_struc = [1e6 -ones(1,length(c))];
    m = 1;
    csense = [repmat('E',1,K.f) repmat('L',1,K.l)];
end

[Nbegcol,Nlencol,Nrowndx,Nobjcnt,Nobjndx,Apatt] = jacSparsity(interfacedata);
A = -F_struc(:,1+linearindicies);
b = full(F_struc(:,1));
csense = [repmat('E',1,K.f) repmat('L',1,K.l)];
A = A.*(~Apatt);
b(any(Apatt,2)) = 0;

%[MY_LICENSE_KEY,nErr] = mxlindo('LSloadLicenseString',MY_LICENSE_FILE);

%[iEnv,nErr]=mxlindo('LScreateEnv',MY_LICENSE_KEY);
%if nErr ~= LSERR_NO_ERROR;output = returnempty(-5); return; end;
[iModel,nErr]=mxlindo('LScreateModel',iEnv);
if nErr ~= LSERR_NO_ERROR;output = returnempty(11); return; end;
constant_data = setup_fmincon_params(interfacedata);
constant_data.F_struc = F_struc;
lindo_fun([],[],[],[],[],[],constant_data);
[nErr] = mxlindo('LSsetFuncalc', iModel, 'lindo_fun',constant_data);
if nErr ~= LSERR_NO_ERROR;output = returnempty(11); return; end;
[nErr] = mxlindo('LSsetModelIntParameter', iModel, LS_IPARAM_NLP_PRINTLEVEL, options.verbose+1);
if nErr ~= LSERR_NO_ERROR;output = returnempty(11); return; end;

% Set NLP solver
[nErr] = mxlindo('LSsetModelIntParameter', iModel, LS_IPARAM_NLP_SOLVER, eval(options.lindo.LS_IPARAM_NLP_SOLVER));
[nErr] = mxlindo('LSsetModelIntParameter', iModel, LS_IPARAM_NLP_MAXLOCALSEARCH,options.lindo.LS_IPARAM_NLP_MAXLOCALSEARCH);

% Load the LP portion of  model
[nErr] = mxlindo('LSXloadLPData', iModel, 1, 0, c(linearindicies), b, csense,sparse(A), lb(linearindicies), ub(linearindicies));
if nErr ~= LSERR_NO_ERROR;output = createoutput(11); return; end;

nErr = mxlindo('LSloadVarType',iModel,vtype);
if nErr ~= LSERR_NO_ERROR;output = createoutput(11); return; end;

% Load the NLP portion of the model
[nErr] = mxlindo('LSloadNLPData', iModel, Nbegcol, Nlencol,[], Nrowndx, Nobjcnt,Nobjndx,[]);
if nErr ~= LSERR_NO_ERROR;output = createoutput(11); return; end;

% Optimize model
solvertime = tic;

if isempty(interfacedata.integer_variables)
    solver = 2;
else
    solver = 1;
end
solvertime = tic;
switch solver
    case 1
        [solstat,nErr] = mxlindo('LSsolveMIP', iModel);
        if ~ismember(solstat,[2009 LS_STATUS_INFEASIBLE])
            [x,nErr] = mxlindo('LSgetMIPPrimalSolution',iModel);
        else
            x = zeros(length(linearindicies),1);
        end
    case 2
        [solstat,nErr] = mxlindo('LSoptimize', iModel, eval(options.lindo.LS_METHOD));
        if ~ismember(solstat,[2009 LS_STATUS_INFEASIBLE])
            [x,nErr] = mxlindo('LSgetPrimalSolution',iModel);
        else
            x = zeros(length(linearindicies),1);
        end
    case 3
        [solStatus,nErr] = mxlindo('LSsolveGOP', iModel);
        [x,nErr] = mxlindo('LSgetPrimalSolution',iModel);
    otherwise
end
solvertime = toc(solvertime);

w = zeros(length(c),1);w(linearindicies) =x;
y = [];

%[nErr]=mxlindo('LSdeleteEnv',iEnv);
[nErr]=mxlindo('LSdeleteModel',iModel);

switch solstat
    case {LS_STATUS_OPTIMAL,LS_STATUS_BASIC_OPTIMAL,7,8}
        problem = 0;
    case {LS_STATUS_INFEASIBLE,LS_STATUS_LOCAL_INFEASIBLE}
        problem = 1;
    case {LS_STATUS_UNBOUNDED}
        problem = 2;
    otherwise
        problem = 11;
end
infostr = yalmiperror(problem,'LINDO-QP');

% Save all data sent to solver?
if options.savesolverinput
    solverinput.solstat = solstat;
    solverinput.nErr = nErr;
    solverinput.x = x;
else
    solverinput = [];
end

% Save all data from the solver?
if options.savesolveroutput
    solveroutput.x = x;
    solveroutput.fmin = fmin;
    solveroutput.flag = flag;
    solveroutput.output=output;
    solveroutput.lambda=lambda;
else
    solveroutput = [];
end

% Standard interface
output = createoutput(w,y,[],problem,'LINDO',solverinput,solveroutput,solvertime);

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function str = showdetail(item)

% function str = showdetail(item)
% Display details for a cfg_exbranch and all of its options.
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: showdetail.m 4867 2012-08-30 13:04:51Z volkmar $

rev = '$Rev: 4867 $'; %#ok

str = showdetail(item.cfg_branch);
str{end+1} = 'Class  : cfg_exbranch';
str{end+1} = 'Functions associated with this item:';
% Display details about function handles
str{end+1} = sprintf('.prog  : %s', func2str(item.prog));
if ~isempty(item.vout)
    str{end+1} = sprintf('.vout  : %s', func2str(item.vout));
elseif ~isempty(item.vfiles)
    str{end+1} = sprintf('.vfiles: %s', func2str(item.vfiles));
end

author: spm


text:
function files = mrtrix_build_files(fname_trunk,lmax)
% Builds a structure with the names of the files that the MRtrix commands
% will generate and need.
%
% files = mrtrix_build_files(fname_trunk,lmax)
%
% Franco Pestilli, Ariel Rokem, Bob Dougherty Stanford University

% Convert the raw dwi data to the mrtrix format:
files.dwi = strcat(fname_trunk,'_dwi.mif');

% This file contains both bvecs and bvals, as per convention of mrtrix
files.b     = strcat(fname_trunk, '.b');

% Convert the brain mask from mrDiffusion into a .mif file:
files.brainmask = strcat(fname_trunk,'_brainmask.mif');

% Generate diffusion tensors:
files.dt = strcat(fname_trunk, '_dt.mif');

% Get the FA from the diffusion tensor estimates:
files.fa = strcat(fname_trunk, '_fa.mif');

% Generate the eigenvectors, weighted by FA:
files.ev = strcat(fname_trunk, '_ev.mif');

% Estimate the response function of single fibers:
files.sf = strcat(fname_trunk, '_sf.mif');
files.response = strcat(fname_trunk, '_response.txt');

% Create a white-matter mask, tracktography will act only in here.
files.wm    = strcat(fname_trunk, '_wm.mif');

% Compute the CSD estimates:
files.csd = strcat(fname_trunk, sprintf('_csd_lmax%i.mif',lmax));

end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function flatten_hmm(hmm,filename)

K = hmm.K;
train = hmm.train;
prior = hmm.prior;
Dir_alpha = hmm.Dir_alpha;
Pi = hmm.Pi;
Dir2d_alpha = hmm.Dir2d_alpha;
P = hmm.P;

for k = 1:K
    if isfield(hmm.state(k),'W') && isfield(hmm.state(k).W,'Mu_W') && ~isempty(hmm.state(k).W.Mu_W)
        eval(['state_' num2str(k-1) '_Mu_W = hmm.state(k).W.Mu_W;']);
        eval(['state_' num2str(k-1) '_S_W = hmm.state(k).W.S_W;']);
        if isfield(hmm.state(k).W,'iS_W')
            eval(['state_' num2str(k-1) '_iS_W = hmm.state(k).W.iS_W;']);
        end
    end
    if isfield(hmm.state(k),'Omega') && isfield(hmm.state(k).Omega,'Gam_rate')
        eval(['state_' num2str(k-1) '_Omega_Gam_rate = hmm.state(k).Omega.Gam_rate;']);
        eval(['state_' num2str(k-1) '_Omega_Gam_shape = hmm.state(k).Omega.Gam_shape;']);
        if isfield(hmm.state(k).Omega,'Gam_irate')
            eval(['state_' num2str(k-1) '_Omega_Gam_irate = hmm.state(k).Omega.Gam_irate;']);
        end
        eval(['state_' num2str(k-1) '_prior_Omega_Gam_rate = hmm.state(k).prior.Omega.Gam_rate;']);
        eval(['state_' num2str(k-1) '_prior_Omega_Gam_shape = hmm.state(k).prior.Omega.Gam_shape;']);
    end
end

if isfield(hmm,'Omega') && isfield(hmm.Omega,'Gam_rate')
    Omega_Gam_rate = hmm.Omega.Gam_rate;
    Omega_Gam_shape = hmm.Omega.Gam_shape;
    if isfield(hmm.Omega,'Gam_irate')
        Omega_Gam_irate = hmm.Omega.Gam_irate;
    end
    prior_Omega_Gam_rate = hmm.prior.Omega.Gam_rate;
    prior_Omega_Gam_shape = hmm.prior.Omega.Gam_shape;
end

save_string = strcat("save('",filename,"','K','train','Dir_alpha','Pi','Dir2d_alpha','P'");

for k = 1:K
    if isfield(hmm.state(k),'W') && isfield(hmm.state(k).W,'Mu_W') && ~isempty(hmm.state(k).W.Mu_W)
        save_string = strcat(save_string,",'state_",num2str(k-1),"_Mu_W'");
        save_string = strcat(save_string,",'state_",num2str(k-1),"_S_W'");
        if isfield(hmm.state(k).W,'iS_W')
            save_string = strcat(save_string,",'state_",num2str(k-1),"_iS_W'");
        end
    end
    if isfield(hmm.state(k),'Omega') && isfield(hmm.state(k).Omega,'Gam_rate')
        save_string = strcat(save_string,",'state_",num2str(k-1),"_Omega_Gam_rate'");
        save_string = strcat(save_string,",'state_",num2str(k-1),"_Omega_Gam_shape'");
        if isfield(hmm.state(k).Omega,'Gam_irate')
            save_string = strcat(save_string,",'state_",num2str(k-1),"_Omega_Gam_irate'");
        end
        save_string = strcat(save_string,",'state_",num2str(k-1),"_prior_Omega_Gam_rate'");
        save_string = strcat(save_string,",'state_",num2str(k-1),"_prior_Omega_Gam_shape'");
    end
end

if isfield(hmm,'Omega') && isfield(hmm.Omega,'Gam_rate')
    save_string = strcat(save_string,",'Omega_Gam_rate'");
    save_string = strcat(save_string,",'Omega_Gam_shape'");
    if isfield(hmm.Omega,'Gam_irate')
        save_string = strcat(save_string,",'Omega_Gam_irate'");
    end
    save_string = strcat(save_string,",'prior_Omega_Gam_rate'");
    save_string = strcat(save_string,",'prior_Omega_Gam_shape'");
end
save_string = strcat(save_string,")");

eval(save_string)

end

author: OHBA-analysis

repo: HMM-MAR


text:
% GET_BOUNDARY_INDICES: given a set of indices (elements or dofs) in the
%  whole domain, collect those that belong to a given boundary.
%
%   bnd_indices = get_boundary_indices (iside, size_dir, indices)
%
% INPUT:
%
%   iside:    number of the side where to compute the indices
%   size_dir: the size in each direction (usually given by ndof_dir or nel_dir)
%   indices:  set of indices in the whole domain
%
% OUTPUT:
%
%   bnd_indices: set of indices that belong to the boundary, with the
%     local numbering of the boundary.
%
% Copyright (C) 2015 Rafael Vazquez
%
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.

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

function bnd_indices = get_boundary_indices (iside, size_dir, indices)

  ndim = numel (size_dir);
  ind2 = ceil (iside/2);
  ind = setdiff (1:ndim, ind2);

  if (mod(iside,2) == 1)
    boundary_ind = 1;
  else
    boundary_ind = size_dir (ind2);
  end

  indsub = cell (1, ndim);
  [indsub{:}] = ind2sub (size_dir, indices);
  aux = find (indsub{ind2} == boundary_ind);
  ppp = cellfun (@(x) x(aux), indsub(ind), 'UniformOutput', false);
  bnd_indices = sub2ind ([size_dir(ind), 1], ppp{:});

end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function [v] = tapas_mh_mc3g_arc(ollh, olpp, nllh, nlpp, ...
    ratio, T)
%% Acceptance rejection criterion for metropolis hastings in the context of
% population mcmc generalized for hierarchical models.
%
% aponteeduardo@gmail.com
% copyright (C) 2016
%

v = nllh .* T + nlpp - (ollh .* T + olpp) + ratio;

nansv = isnan(v);
v(nansv) = -inf;

v = rand(size(v)) < exp(min(v, 0));

assert(all(-inf < nllh(v) + nlpp(v)), 'tapas:mh', ...
    '-inf value in the new samples');

end

author: translationalneuromodeling

repo: tapas


text:
function test_triangulation_test02 ( )

%*****************************************************************************80
%
%% TEST02 tests P00_TEST_NUM and P00_SAMPLE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 April 2009
%
%  Author:
%
%    John Burkardt
%
  m = 2;
  seed_start = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST02\n' );
  fprintf ( 1, '  P00_TEST_NUM reports the number of problems.\n' );
  fprintf ( 1, '  P00_SAMPLE returns sample points from the region.\n' );

  test_num = p00_test_num ( );

  for test = 1 : test_num

    seed = seed_start;

    title = p00_title ( test );

    n = 20;

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Test number         =  %d\n', test );
    fprintf ( 1, '  Title:              = "%s"\n', title );
    fprintf ( 1, '  Spatial dimension M =  %d\n', m );
    fprintf ( 1, '  Number of samples N =  %d\n', n );
    fprintf ( 1, '  Initial SEED:       =  %d\n', seed );
    fprintf ( 1, '\n' );

    [ point, seed ] = p00_sample ( test, m, n, seed );

    r8mat_transpose_print ( m, n, point, '  The sample points:' );

  end

  return
end


text:
function [mbar] = mmHg2mbar(mmHg)
% Convert pressure from millimeters of mercury at 0 degrees C to millibars
% Chad Greene 2012
mbar = mmHg*1.33322;


text:
function im = ifftc(d)
% Function performs a centered ifft
im = ifftshift(ifft(ifftshift(d)));


text:
function replaceMainCatalog(otherCatalog)
    % replaceMainCatalog(otherCatalog)
    % use this instead of a=something
    % protects the main catalog from corruption
    report_this_filefun();
    if isempty(otherCatalog)
        otherCatalog=ZmapCatalog();
    end
    if isa(otherCatalog,'ZmapCatalogView')
        otherCatalog=otherCatalog.Catalog();
    end
    assert(isa(otherCatalog,'ZmapCatalog'));

    ZG = ZmapGlobal.Data;
    ZG.primeCatalog = copy(otherCatalog);
end

author: CelsoReyes


text:
ops.chanMap             = 'D:\GitHub\KiloSort2\configFiles\neuropixPhase3A_kilosortChanMap.mat';
% ops.chanMap = 1:ops.Nchan; % treated as linear probe if no chanMap file

% sample rate
ops.fs = 30000;

% frequency for high pass filtering (150)
ops.fshigh = 150;

% minimum firing rate on a "good" channel (0 to skip)
ops.minfr_goodchannels = 0.1;

% threshold on projections (like in Kilosort1, can be different for last pass like [10 4])
ops.Th = [10 4];

% how important is the amplitude penalty (like in Kilosort1, 0 means not used, 10 is average, 50 is a lot)
ops.lam = 10;

% splitting a cluster at the end requires at least this much isolation for each sub-cluster (max = 1)
ops.AUCsplit = 0.9;

% minimum spike rate (Hz), if a cluster falls below this for too long it gets removed
ops.minFR = 1/50;

% number of samples to average over (annealed from first to second value)
ops.momentum = [20 400];

% spatial constant in um for computing residual variance of spike
ops.sigmaMask = 30;

% threshold crossings for pre-clustering (in PCA projection space)
ops.ThPre = 8;
%% danger, changing these settings can lead to fatal errors
% options for determining PCs
ops.spkTh           = -6;      % spike threshold in standard deviations (-6)
ops.reorder         = 1;       % whether to reorder batches for drift correction.
ops.nskip           = 25;  % how many batches to skip for determining spike PCs

ops.GPU                 = 1; % has to be 1, no CPU version yet, sorry
% ops.Nfilt               = 1024; % max number of clusters
ops.nfilt_factor        = 4; % max number of clusters per good channel (even temporary ones)
ops.ntbuff              = 64;    % samples of symmetrical buffer for whitening and spike detection
ops.NT                  = 64*1024+ ops.ntbuff; % must be multiple of 32 + ntbuff. This is the batch size (try decreasing if out of memory).
ops.whiteningRange      = 32; % number of channels to use for whitening each channel
ops.nSkipCov            = 25; % compute whitening matrix from every N-th batch
ops.scaleproc           = 200;   % int16 scaling of whitened data
ops.nPCs                = 3; % how many PCs to project the spikes into
ops.useRAM              = 0; % not yet available

%%

author: MouseLand

repo: Kilosort

save_path: github-repos/MATLAB/MouseLand-Kilosort


text:
classdef MUSIC < handle
    % Implements MUSIC-MMV algorithm for sparse approximation

    properties
        % Indicates if log messages should be printed at each iteration
        Verbose = false
    end

    properties(SetAccess=private)
        % The dictionary
        Dict
        % Ambient signal dimensions
        N
        % Number of atoms in dictionary
        D
        % Sparsity level of representations (may be negative)
        K
        % Result of a solver
        result
    end

    methods
        function self  = MUSIC(Dict, K)
            % We assume that all the columns in dictionary are normalized.
            if isa(Dict, 'spx.dict.Operator')
                self.Dict = Dict;
            elseif ismatrix(Dict)
                self.Dict = spx.dict.MatrixOperator(Dict);
            else
                error('Unsupported operator.');
            end
            [self.N, self.D] = size(Dict);
            if nargin < 2
                error('Sparsity level must be specified.');
            end
            self.K = K;
        end

        function result  = solve(self, Y)
            % Initialization
            % Solves approximation problem using OMP
            d = self.D;
            n = self.N;
            % The number of signals being approximated.
            s = size(Y, 2);
            dict = self.Dict;
            % Active indices
            omega = [];
            % Estimate
            Z = zeros(d, s);
            % Orthonormalize Y
            U = orth(Y); % Y is n x s. U is n x r.
            % Verify its rank
            k = self.K;
            r = rank(U);
            if r < k
                error('MUSIC principle cannot be directly applied.');
            end
            % Prepare the projection operator
            P = eye(n) - U*U'; % P is n x n
            % Compute the projection of each atom on
            % orthogonal complement of R(U)
            % compute the norms of projections
            projection_norms = zeros(1, d);
            for i=1:d
                projection_norms(i) = norm(P * (dict.column(i)) );
            end
            % sort the norms
            [sorted_norms, indices] = sort(projection_norms);
            % disp(sorted_norms);
            % Pick first K indices as the solution support.
            omega = indices(1:self.K);
            % Solve least squares problem
            subdict = columns(dict, omega);
            tmp = linsolve(subdict, Y);
            % Updated solution
            Z(omega, :) = tmp;
            % Let us update the residual.
            R = Y - dict.apply(Z);
            resNorm = norm(R, 'fro');
            % Solution vector
            result.Z = Z;
            % Residual obtained
            result.R = R;
            % Solution support
            result.support = omega;
            % residual Frobenius norm
            result.residual_frobenius_norm = resNorm;
            self.result = result;
        end
    end
end

author: indigits

repo: sparse-plex


text:
function data=clust_normalize(data,method);

%   method     description
%   'var'      Variance is normalized to one (linear operation).
%   'range'    Values are normalized between [0,1] (linear operation).

data.Xold=data.X;
if strcmp(method,'range')
     data.min=min(data.X);
     data.max=max(data.X);
     data.X=(data.X-repmat(min(data.X),size(data.X,1),1))./(repmat(max(data.X),...
         size(data.X,1),1)-repmat(min(data.X),size(data.X,1),1));
 elseif strcmp(method,'var')
     data.X=(data.X-repmat(mean(data.X),size(data.X,1),1))./(repmat(std(data.X),size(data.X,1),1));
     data.mean=mean(data.X);
     data.std=std(data.X);
 else
     error('Unknown method given')
end


text:
function c=ref_fcdgt(f,g,a,M,m_t,m_f,w_t,w_f)
%REF_CDGT  Reference centered DGT
%   Usage:  c=ref_dgtiv(f,g,a,M,c_t,c_f,c_w);
%
%   Linear algebra version of the algorithm. Create big matrix
%   containing all the basis functions and multiply with the transpose.
%
%   For easy work, m_t,m_f,w_t,w_f are all just 0/1 indicator variables.

L=size(f,1);
N=L/a;
b=L/M;

m_t=m_t*.5;
w_t=w_t*floor(a/2);
w_f=w_f*ceil(b/2);

F=zeros(L,M*N);

l=(0:L-1).';

if m_f==0

  for n=0:N-1
    for m=0:M-1
      F(:,M*n+m+1)=exp(2*pi*i*(m*b+w_f)*(l+m_t)/L).*circshift(g,n*a+w_t);
    end;
  end;

else

  for n=0:N-1
    for m=0:M-1
      F(:,M*n+m+1)=exp(2*pi*i*(m*b+.5+w_f)*(l+m_t-n*a)/L).*circshift(g,n*a+w_t);
    end;
  end;

end;

c=F'*f;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat