text:
function params = mrInitGUI_analParams(params);
% Assign analysis parameters for a new mrVista session (mrInit2).
%
% params = mrInitGUI_analParams(params);
%
% This GUI allows the user to assign GLM ("event-related") and
% coherence ("block") analysis parameters to scans, assign parfiles
% to scans, and group scans. It returns params for mrInit2 which include
% the updated parameters.
%
%
% ras, 07/2007.
if notDefined('params'),    params = mrInitDefaultParams;       end

%% GUI update function switch
% the callbacks will call this function with special calls to sub-functions
% in lieu of a params struct. Parse this:
if ischar(params)
	switch lower(params)
		case 'assignglmparams', mrInitGUI_assignGlmParams;
		case 'assigncoparams', mrInitGUI_assignCoParams;
		case 'assignparfiles', mrInitGUI_assignParfiles;
		case 'assignparfilepattern', mrInitGUI_assignParfilePattern;
		case 'groupscans', mrInitGUI_groupScans;
		otherwise, error(sprintf('Invalid argument for %s', mfilename));
	end
	return
end

% --- the code below makes the GUI --- %
%%%%% params
fsz = 10;  % fontsize for controls
bgColor = [.9 .9 .9];  % background color for figure

%%%%% Create the main GUI figure
hFig = figure('Units', 'norm', 'Position', [.15 .3 .7 .35], ...
	'UserData', params, ...
	'Name', 'Assign Analysis Parameters', 'NumberTitle', 'off', ...
	'Color', bgColor, 'CloseRequestFcn', 'OK = 0; uiresume;');

%%%%% listboxes
%% scans listbox
annotations = {};
for i = 1:length(params.annotations)
	annotations{i} = sprintf('%i: %s', i, params.annotations{i});
end
hScans = uicontrol('Style', 'listbox', 'String', annotations, ...
	'Min', 0, 'Max', 2, 'Tag', 'ScanList', ... % multiple selection
	'Units', 'norm', 'Position', [.05 .2 .25 .7], ...
	'BackgroundColor', 'w', 'ForegroundColor', 'k', 'FontSize', fsz);

% text label
uicontrol('Style', 'text', 'String', 'Scans', 'FontSize', fsz+2, ...
	'Units', 'norm', 'Position', [.05 .93 .25 .05], ...
	'HorizontalAlignment', 'center', 'BackgroundColor', bgColor);

%% parfiles listbox
for i = 1:length(annotations)  % nScans
	if (length(params.parfile) >= i) & (~isempty(params.parfile{i}))
		parfiles{i} = sprintf('%i: %s', i, params.parfile{i});
	else
		parfiles{i} = sprintf('%i: (none)', i);
	end
end
hParfiles = uicontrol('Style', 'listbox', 'String', parfiles, ...
	'Min', 0, 'Max', 2, 'Tag', 'ParfilesList', ... % multiple selection
	'Units', 'norm', 'Position', [.35 .2 .15 .7], 'Value', [], ...
	'BackgroundColor', bgColor, 'ForegroundColor', 'k', 'FontSize', fsz);

% text label
uicontrol('Style', 'text', 'String', 'Parfiles', 'FontSize', fsz+2, ...
	'Units', 'norm', 'Position', [.35 .93 .15 .05],  ...
	'HorizontalAlignment', 'center', 'BackgroundColor', bgColor);

%% scan groups listbox
for i = 1:length(annotations)  % nScans
	scanGroups{i} = sprintf('%i: (none)', i);
end

for i = 1:length(params.scanGroups)
	for j = params.scanGroups{i}
		scanGroups{j} = sprintf('%i: %s', i, num2str(params.scanGroups{i}));
	end
end

hScanGroups = uicontrol('Style', 'listbox', 'String', scanGroups, ...
	'Min', 0, 'Max', 2, 'Tag', 'ScanGroupsList', ... % multiple selection
	'Units', 'norm', 'Position', [.55 .2 .15 .7], 'Value', [], ...
	'BackgroundColor', bgColor, 'ForegroundColor', 'k', 'FontSize', fsz);

% text label
uicontrol('Style', 'text', 'String', 'Scan Groups', 'FontSize', fsz+2, ...
	'Units', 'norm', 'Position', [.55 .93 .12 .05], ...
	'HorizontalAlignment', 'center', 'BackgroundColor', bgColor);

%%%%% buttons
%% assign co params button
uicontrol('Style', 'pushbutton', 'String', 'Assign Coherence Params', ...
	'Units', 'norm', 'Position', [.75 .8 .2 .12], ...
	'BackgroundColor', [.8 .8 .9], 'ForegroundColor', 'k', 'FontSize', fsz, ...
	'Callback', 'mrInitGUI_analParams(''assignCoParams''); ');

%% assign GLM params button
uicontrol('Style', 'pushbutton', 'String', 'Assign GLM Params', ...
	'Units', 'norm', 'Position', [.75 .65 .2 .12], ...
	'BackgroundColor', [.8 .8 .9], 'ForegroundColor', 'k', 'FontSize', fsz, ...
	'Callback', 'mrInitGUI_analParams(''assignGlmParams''); ');

%% assign parfiles button
uicontrol('Style', 'pushbutton', 'String', 'Assign Parfiles', ...
	'Units', 'norm', 'Position', [.75 .5 .2 .12], ...
	'UserData', hParfiles, ...
	'BackgroundColor', [.8 .8 .9], 'ForegroundColor', 'k', 'FontSize', fsz, ...
	'Callback', 'mrInitGUI_analParams(''assignParfiles''); ');

%% assign parfiles pattern button
uicontrol('Style', 'pushbutton', 'String', 'Assign Parfile Pattern', ...
	'Units', 'norm', 'Position', [.75 .35 .2 .12], ...
	'UserData', hParfiles, ...
	'BackgroundColor', [.8 .8 .9], 'ForegroundColor', 'k', 'FontSize', fsz, ...
	'Callback', 'mrInitGUI_analParams(''assignParfilePattern''); ');

%% group scans button
uicontrol('Style', 'pushbutton', 'String', 'Group Scans', ...
	'Units', 'norm', 'Position', [.75 .2 .2 .12], ...
	'UserData', hScanGroups, ...
	'BackgroundColor', [.8 .8 .9], 'ForegroundColor', 'k', 'FontSize', fsz, ...
	'Callback', 'mrInitGUI_analParams(''groupScans''); ');

%%%%%  GO / Cancel Buttons
uicontrol('Style', 'pushbutton', 'Units', 'norm', 'Position', [.1 .02 .2 .08], ...
		  'BackgroundColor', [.9 .8 .8], 'String', 'Cancel', ...
		  'Callback', 'OK = 0; uiresume;');

uicontrol('Style', 'pushbutton', 'Units', 'norm', 'Position', [.5 .02 .2 .08], ...
		  'BackgroundColor', [.4 .8 .5], 'String', 'OK', ...
		  'Callback', 'OK = 1; uiresume;');

%%%%% get the params: do a uiwait/uiresume
% (uiresume is in the callbacks for this GUI)
uiwait;

%% we uiresumed: get the updated values from the figure
OK = evalin('base', 'OK');
evalin('base', 'clear OK'); % clear from base workspace
if OK==1
	% get modified params
	params = get(gcf, 'UserData');
end

delete(hFig);

return
% /------------------------------------------------------------------/ %

% /------------------------------------------------------------------/ %
function mrInitGUI_assignCoParams;
% put up a dialog to set the coherence params for the selected scans,
% then modify the stored coParams field to match the edited params.
params = get(gcf, 'UserData');

hScans = findobj('Parent', gcf, 'Tag', 'ScanList');
selScans = get(hScans, 'Value');

% get initial values: take from the first selected scan
% if it's assigned in params, use that; else use defaults
s = selScans(1);
if (length(params.coParams) >= s)  &  (~isempty(params.coParams{s}))
	coParams = params.coParams{s};
else
	coParams = coParamsDefault;
end

% have user edit the values
[coParams ok] = coParamsEdit(coParams);
if ~ok, return; end

% update the params field
for s = selScans
	params.coParams{s} = coParams;
end

set(gcf, 'UserData', params);
return
% /------------------------------------------------------------------/ %

% /------------------------------------------------------------------/ %
function mrInitGUI_assignGlmParams;
% put up a dialog to set the GLM analysis params for the selected scans,
% then modify the stored glmParams field to match the edited params.
params = get(gcf, 'UserData');

hScans = findobj('Parent', gcf, 'Tag', 'ScanList');
selScans = get(hScans, 'Value');

% get initial values: take from the first selected scan
% if it's assigned in params, use that; else use defaults
s = selScans(1);
if (length(params.coParams) >= s)  &  (~isempty(params.coParams{s}))
	glmParams = params.glmParams{s};
else
	glmParams = er_defaultParams;
end

% have user edit the values
[glmParams ok] = er_editParams(glmParams);
if ~ok, return; end

% update the params field
for s = selScans
	params.glmParams{s} = glmParams;
end

set(gcf, 'UserData', params);
return
% /------------------------------------------------------------------/ %

% /------------------------------------------------------------------/ %
function mrInitGUI_assignParfiles;
% prompt the user to select parfiles to assign to the selected scans,
% and update the params and GUI accordingly.
params = get(gcf, 'UserData');

hScans = findobj('Parent', gcf, 'Tag', 'ScanList');
selScans = get(hScans, 'Value');

% first, we need to make sure we have a parfiles directory:
% this is either stim/parfiles or Stimuli/parfiles in the session dir:
parDir = fullfile(params.sessionDir, 'Stimuli', 'parfiles');
if ~exist(parDir, 'dir')
	if exist( fullfile(params.sessionDir, 'stim', 'parfiles'), 'dir' )
		parDir = fullfile(params.sessionDir, 'stim', 'parfiles');
	end
end

if ~exist(parDir, 'dir')
	myWarnDlg(['You don''t have a Stimuli/parfiles directory ' ...
			   'in your session directory. Parfiles should be kept there.']);
   return
end

parList = dir( fullfile(parDir, '*.par') );
parList = {parList.name};

% also add an option to remove the assignment of a parfile:
parList = [{'(none)'} parList];

% put up the dialog
ttl = sprintf('Select parfiles to assign to scans: %s', num2str(selScans));

[sel, ok] = listdlg('PromptString', ttl, 'ListSize', [400 600], ...
					'ListString', parList, 'InitialValue', 1, ...
					'OKString', 'OK');
if ~ok  return;  end

% ensure the selected parfiles match the number of selected scans.
% The strategy here is simple: if sel is shorter than selScans,
% all scans in between get assigned the last selected parfile.
% This way, you can assign one parfile to many scans easily.
if length(sel) < length(selScans)
	nExtra = length(selScans) - length(sel);
	sel = [sel repmat(sel(end), [1 nExtra])];
end

% update the parfiles listbox
hParfiles = findobj('Parent', gcf, 'Tag', 'ParfilesList');
str = get(hParfiles, 'String');
for s = 1:length(selScans)
	str{selScans(s)} = sprintf('%i: %s', s, parList{sel(s)});
end
set(hParfiles, 'String', str);

% update the params field
for s = 1:length(selScans)
	if sel(s)==1	% selected 'none', no parfile to this scan
		params.parfile{selScans(s)} = '';
	else
		params.parfile{selScans(s)} = parList{sel(s)};
	end
end
set(gcf, 'UserData', params);

return
% /------------------------------------------------------------------/ %

% /------------------------------------------------------------------/ %
function mrInitGUI_assignParfilePattern;
% a separate dialog providing a more flexible way to assign parfiles
% to scans: if the parfiles for several scans have a common pattern
% in their names -- e.g., 'run1.par' through 'run5.par' -- this dialog
% lets the user specify the pattern, and the sequence of numbers which
% correspond to the selected scans.
%
% Note that the sequence of numbers need not run from 1 to N, nor
% be linear increasing. For instance, suppose I ran my scripts in a
% counterbalanced order, going forwards (1-10) in one session, and
% backwards (10-1) in another session. And suppose a scanning glitch
% caused run 1 to precede run 2. I might assign it thusly:
%
% pattern: run%i.par
% runs:  [10:-1:3 1 2]
%
% Since this code uses sprintf formatting, the user needs to put the
% symbol '%i' where the number should go.
%
% The code then updates the params and the GUI with the new parfiles.
params = get(gcf, 'UserData');

hScans = findobj('Parent', gcf, 'Tag', 'ScanList');
selScans = get(hScans, 'Value');

% build the dialog
dlg(1).fieldName = 'pattern';
dlg(1).style	 = 'edit';
dlg(1).string	 = ['Parfile pattern: insert a %i where the run number ' ...
					'should go. (E.g. ''run%i.par'') '];
dlg(1).value	 = '%i.par';

dlg(2).fieldName = 'runVals';
dlg(2).style	 = 'edit';
dlg(2).string	 = ['Run values (will be pasted in %i position); ' ...
					'E.g. ''1:8'', ''8:-1:1'', ''[1 3 5:9]'' '];
dlg(2).value	 = '';

% put up the dialog and get a response
resp = generalDialog(dlg, 'Assign Parfile Pattern', 'center');

% generate the new parfile list
runVals = str2num( resp.runVals );
for s = runVals
	parList{s} = sprintf( resp.pattern, s );
end

% update the parfiles listbox
hParfiles = findobj('Parent', gcf, 'Tag', 'ParfilesList');
str = get(hParfiles, 'String');
for s = selScans
	str{s} = sprintf('%i: %s', s, parList{s});
end
set(hParfiles, 'String', str);

% update the params field
for s = 1:length(selScans)
	params.parfile{selScans(s)} = parList{s};
end
set(gcf, 'UserData', params);

return
% /------------------------------------------------------------------/ %

% /------------------------------------------------------------------/ %
function mrInitGUI_groupScans;
% group the selected scans in the scans list. This modified the
% params.scanGroupField, and updates the scan group listbox to
% reflect that change.
params = get(gcf, 'UserData');

hScans = findobj('Parent', gcf, 'Tag', 'ScanList');
selScans = get(hScans, 'Value');

params.scanGroups{end+1} = selScans;

% update the scan groups listbox
hScanGroup = findobj('Parent', gcf, 'Tag', 'ScanGroupsList');
str = get(hScanGroup, 'String');
for s = selScans
	str{s} = sprintf('%i: [%s]', s, num2str(selScans));
end
set(hScanGroup, 'String', str);

% set updated params
set(gcf, 'UserData', params);
return

author: vistalab

repo: vistasoft


text:
function ADEM_lorenz_entropy
% This demo shows how structure can be instilled from an environment. It
% uses an agent that optimises its recognition density over successive
% epochs of exposure to an environment. This environment causes the agent
% to flow on a Lorenz attractor with random perturbations. As the agent
% learns the causal regularities in its environment, it is better able to
% predict them and act to oppose the random effect. The result is that it
% is more robust to random forces and therefore exhibits states with lower
% entropy. This routine takes several minutes to run.
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: ADEM_lorenz_entropy.m 4804 2012-07-26 13:14:18Z karl $

% generative process (environment)
%==========================================================================

% set up
%--------------------------------------------------------------------------
G(1).E.s = 1/4;                             % smoothness
G(1).E.n = 4;                               % embedding
G(1).E.d = 2;                               % embedding

% dynamics
%--------------------------------------------------------------------------
fG      = '[v + a; 0; 0] + [-P(1) P(1) 0; P(3) -1 -x(1); x(2) 0 P(2)]*x/64';
fM      = '[v    ; 0; 0] + [-P(1) P(1) 0; P(3) -1 -x(1); x(2) 0 P(2)]*x/64';

% parameters
%--------------------------------------------------------------------------
PG      = [10; -8/3; 32];                   % Target parameters
PM      = [10;    0;  0];                   % initial parameters

% level 1
%--------------------------------------------------------------------------
G(1).x  = [1; 1; 24];
G(1).f  = inline(fG ,'x','v','a','P');
G(1).g  = inline('x','x','v','a','P');
G(1).pE = PG;
G(1).V  = exp(8);                           % error precision
G(1).W  = exp(16);                          % error precision

% level 2
%--------------------------------------------------------------------------
G(2).a  = 0;                                % action variables
G(2).v  = 0;                                % inputs
G(2).V  = exp(16);
G       = spm_ADEM_M_set(G);

% % plot flow fields and equilibrium densities
% %==========================================================================
% spm_figure('GetWin','Figure 1');
%
% x{1}    = linspace(-20,20,32);
% x{2}    = linspace(-32,32,32);
% x{3}    = linspace(  10,40,8);
%
%
% % controlled flow (P0)
% %--------------------------------------------------------------------------
% subplot(3,2,1)
% spm_fp_display_density(G,x);
% xlabel('position','Fontsize',12)
% ylabel('velocity','Fontsize',12)
% title('controlled','Fontsize',16)

% recognition model: learn the controlled environmental dynamics
%==========================================================================

% make a naive model (M)
%--------------------------------------------------------------------------
M       = G;
M(1).f  = inline(fM ,'x','v','P');
M(1).g  = inline('x','x','v','P');
M(1).pE = PM;
M(1).pC = eye(3)*4;
M(1).V  = exp(8);
M(1).W  = exp(8);
M       = spm_DEM_M_set(M);

% teach naive model by exposing it to a controlled environment (G)
%==========================================================================

% length of realization
%--------------------------------------------------------------------------
n     = 256;
U     = sparse(n,1);

DEM.M = M;
DEM.G = G;
DEM.C = U;
DEM.U = U;

% optimise recognition model in epochs
%--------------------------------------------------------------------------
DEM.M(1).E.nE = 1;
H     = sparse(128,128);
for i = 1:32

        % random perturbations
        %------------------------------------------------------------------
        DEM.C       = spm_conv(randn(n,1)*8,8);

        % integrate and update priors
        %------------------------------------------------------------------
        DEM         = spm_ADEM(DEM);
        DEM.M(1).pE = DEM.qP.P{1};
        DEM.M(1).x  = DEM.qU.x{1}(:,end);
        DEM.M(2).a  = DEM.qU.a{2}(:,end);

        % save free-energy and conditional expectations
        %------------------------------------------------------------------
        qP(:,i) = DEM.qP.P{1};
        F(i)    = DEM.F;

        % display states at the beginning and end
        %------------------------------------------------------------------
        spm_figure('GetWin','Figure 1');

        if i < 8
            subplot(3,2,3)
            x = DEM.pU.v{1}(2,:);
            y = DEM.pU.v{1}(3,:);
            plot(x,y); hold on
            title('states before','FontSize',16)
            a = axis;

        elseif i > 24
            subplot(3,2,4)
            x = DEM.pU.v{1}(2,:);
            y = DEM.pU.v{1}(3,:);
            H = H + sparse(64 + fix(x)',64 + fix(y)',1,128,128);
            plot(x,y); hold on
            title('and after','FontSize',16)
            axis(a)

        end
    end

% graphics
%==========================================================================
spm_figure('GetWin','Figure 1');

% plot free-energy
%--------------------------------------------------------------------------
t   = 1:length(F);
subplot(3,1,1)
plot(t,-F)
xlabel('time (epochs)','FontSize',12)
title('Free-energy','FontSize',16)
axis square

% and underlying conditional expectations
%--------------------------------------------------------------------------
subplot(3,2,5)
plot(t,qP,t,kron(ones(1,length(t)),PG),'-.')
xlabel('time (epochs)','FontSize',12)
title('conditional parameters','FontSize',16)
axis square

subplot(3,2,6)
imagesc(H')
xlabel('first state','FontSize',12)
ylabel('second state','FontSize',12)
title('ensemble density','FontSize',16)
axis square xy

author: spm


text:
function n=v_voicebox_update(fn,m)
%V_VOICEBOX_UPDATE update v_voicebox calls by prefixing with 'v_'
%  Inputs: fn   file or folder name to update [default: current folder]
%          m    'r' to search recursively through the folder tree
%               's' silent mode: do not print out statistics
%               'v' verbose: list files as they are updated
%               'V' extra verbose: list files even if not updated
%               'b' do not backup old file
%
% Outputs: n    number of files updated
%
% This routine is not all that smart. Words that match the function names
% included in the 'vbn' list below will be updated even if they are inside
% strings (except for the last string on a line). Local variables that
% share the name of a v_voicebox function will also be modified.
%
% vbp lists function names that may have an argument list
% vbn lists function names that may be called without arguments
vbp='activlev|activlevg|atan2sc|axisenlarge|bark2frq|berk2prob|besselratio|besselratioi|bitsprec|cblabel|ccwarpf|cent2frq|cep2pow|choosenk|choosrnk|convfft|correlogram|distchar|distchpf|disteusq|distisar|distispf|distitar|distitpf|ditherq|dlyapsq|dualdiag|dypsa|enframe|entropy|erb2frq|estnoiseg|estnoisem|ewgrpdel|fig2emf|fig2pdf|figbolden|filtbankm|filterbank|finishat|fopenmkd|frac2bin|fram2wav|frq2bark|frq2cent|frq2erb|frq2mel|frq2midi|fxpefac|fxrapt|gammabank|gammalns|gausprod|gaussmix|gaussmixd|gaussmixg|gaussmixk|gaussmixm|gaussmixp|gaussmixt|glotlf|glotros|gmmlpdf|histndim|hostipinfo|huffman|hypergeom1f1|imagehomog|importsii|irdct|irfft|kmeanhar|kmeanlbg|lambda2rgb|ldatrace|lin2pcma|lin2pcmu|lognmpdf|logsum|lpcaa2ao|lpcaa2dl|lpcaa2rf|lpcao2rf|lpcar2am|lpcar2cc|lpcar2db|lpcar2ff|lpcar2fm|lpcar2im|lpcar2ls|lpcar2pf|lpcar2pp|lpcar2ra|lpcar2rf|lpcar2rr|lpcar2zz|lpcauto|lpcbwexp|lpccc2ar|lpccc2cc|lpccc2db|lpccc2ff|lpccc2pf|lpcconv|lpccovar|lpccw2zz|lpcdb2pf|lpcdl2aa|lpcff2pf|lpcfq2zz|lpcifilt|lpcim2ar|lpcis2rf|lpcla2rf|lpclo2rf|lpcls2ar|lpcpf2cc|lpcpf2ff|lpcpf2rr|lpcpp2cw|lpcpp2pz|lpcpz2zz|lpcra2ar|lpcra2pf|lpcra2pp|lpcrand|lpcrf2aa|lpcrf2ao|lpcrf2ar|lpcrf2is|lpcrf2la|lpcrf2lo|lpcrf2rr|lpcrr2am|lpcrr2ar|lpcss2zz|lpcstable|lpczz2ar|lpczz2cc|lpczz2ss|m2htmlpwd|maxfilt|maxgauss|meansqtf|mel2frq|melbankm|melcepst|midi2frq|minspane|mintrace|modspect|momfilt|mos2pesq|nearnonz|normcdflog|overlapadd|paramsetch|pcma2lin|pcmu2lin|pdfmoments|peak2dquad|permutes|pesq2mos|phon2sone|polygonarea|polygonwind|polygonxline|potsband|pow2cep|prob2berk|psycdigit|psycest|psycestu|psychofunc|qrabs|qrdivide|qrdotdiv|qrdotmult|qrmult|qrpermute|quadpeak|randfilt|randiscr|randvec|rdct|readaif|readau|readcnx|readflac|readhtk|readsfs|readsph|readwav|rectifyhomog|regexfiles|rfft|rhartley|rnsubset|rotation|rotax2qr|roteu2qr|roteu2ro|rotmc2qc|rotmr2qr|rotpl2ro|rotqc2mc|rotqc2qr|rotqr2ax|rotqr2eu|rotqr2mr|rotqr2qc|rotqr2ro|rotqrmean|rotqrvec|rotro2eu|rotro2pl|rotro2qr|rsfft|sapisynth|schmitt|sigalign|skew3d|snrseg|sone2phon|soundspeed|specsub|specsubm|spendred|spgrambw|sphrharm|sprintcpx|sprintsi|ssubmmse|ssubmmsev|stdspectrum|stoi2prob|teager|texthvc|tilefigs|txalign|unixwhich|upolyhedron|usasi|vadsohn|voicebox|vonmisespdf|wavread|wavwrite|winenvar|writehtk|writewav|xticksi|xyzticksi|yticksi|zerocros|zerotrim|zoomfft';
if nargin<2 || isempty(m)
    m='';
end
lm=lower(m);
if nargin<1 || isempty(fn)
    fn='.'; % do files in current folder
end
if exist(fn,'dir')
    ff=v_regexfiles('\.m$',fn,char('n'+('r'-'n')*any(m=='r')));
elseif exist(fn,'file')
    ff={fn};
else
    error('first argument must be a file or folder name');
end
n=0;
for i=1:length(ff)
    fid=fopen(ff{i},'rt');
    s=fread(fid,'*char')';
    fclose(fid);
    t=regexprep(s,['\<(' vbp ')\>(?![^''\n]*''[^''\n]*\n)'],'v_$1'); % keyword unless followed by a single quote somewhere on the same line
    if numel(t)>numel(s)
        if ~any(m=='b')
            if ~movefile(ff{i},[ff{i} '.old'])
                error(sprintf('Cannot rename %s\n',ff{i}));
            end
        end
        fid=fopen(ff{i},'wt');
        if fid<0
            error(sprintf('Cannot write to %s\n',ff{i}));
        end
        fwrite(fid,t,'*char');
        fclose(fid);
        n=n+1;
        if any(lm=='v')
            nch=floor((length(t)-length(s))/2);
            if nch>1
                fprintf(1,'%d changes in %s\n',nch,ff{i});
            else
                fprintf(1,'%d change in %s\n',nch,ff{i});
            end
        end
    elseif any(m=='V')
        fprintf(1,'No changes in %s\n',ff{i});
    end
end
if ~any(m=='s')
    if any(lm=='v') && n>0
        fprintf(1,'===========\n');
    end
    if length(ff)>1
        fprintf(1,'%d of %d files updated\n',n,length(ff));
    else
        fprintf(1,'%d of %d file updated\n',n,length(ff));
    end
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
  function [mat nx ny nb na] = wtf_read(file, varargin)
%|function [mat nx ny nb na] = wtf_read(file, [options])
%| read aspire sparse matrix file (usually file.wtf)
%| and return matlab sparse matrix
%| option
%|	'parse'	0|1	if 1, do it by old matlab parsing way (old wtfread.m)
%| Copyright 2008-9-26, Jeff Fessler, University of Michigan

if ~nargin, ir_usage, end
if nargin == 1 && streq(file, 'test'), wtf_read_test, return, end

arg.chat = 0;
arg.parse = 0;
arg = vararg_pair(arg, varargin);

if ~has_mex_jf
	fail('cannot read .wtf due to mex problem')
end

if arg.parse
	[mat nx ny nb na mask] = wtf_read_parse(file, arg.chat);
return
end

[buff nx ny nb na] = wtfmex('asp:read', file, int32(arg.chat));
[mat nx ny nb na is_tranpose] = wtfmex('asp:mat', buff, int32(arg.chat));
if is_tranpose
	mat = mat';
end

%
% wtf_read_parse()
% reads in a sparse matrix from Aspire format
% A is the returned sparse matrix
% n is a structure containing dimensions
% mask is the nx by ny support mask
%

function [A, nx, ny, nb, na, mask, n] = wtf_read_parse(file, chat)

machine = 'native';		% this should work
%machine = 'ieee-be';		% but try this if it doesn't...
fp = fopen(file, 'r', machine);
if (fp == -1), error fopen, end

%
% read until two form feeds
%
while (1)
	c = fread(fp, 1, 'char');
	if chat
		fprintf(1, '%c', char(c))	% echo header
	end
	if isempty(c)
		error eof
	end
	if (c == 12)	% form feed
		c == fread(fp, 1, 'char');
		if (c ~= 12)
			error only one form feed?
		else
			break;
		end
	end
end

%
% read binary header, extract some dimensions from it
%
tmp = fread(fp, 128/4, 'int');
type.group = tmp(1);
if (type.group ~= 1), error(sprintf('type.group = %d', type.group)), end
type.index = tmp(2);
if (type.index ~= 0), error 'type.index', end
type.value = tmp(3);
if (type.value ~= 0), error 'type.value', end
n.wt = tmp(6);
nx = tmp(12);
ny = tmp(13);
nb = tmp(14);
na = tmp(15);
if chat
	printf('nxy=%d,%d nba=%d,%d nwt=%d', ...
		nx, ny, nb, na, n.wt)
end

%
% read support
%
mask = fread(fp, nx*ny, 'uchar');
mask = reshape(mask, nx, ny);
if chat
	imagesc(mask');
end

%
% read length and offset
%
length = fread(fp, nx*ny, 'uint32');
offset = fread(fp, nx*ny, 'uint32');
if any(diff(offset) ~= length(1:end-1)), error bug, end

%
% read index and value arrays
%
index = fread(fp, n.wt, 'uint16');
value = fread(fp, n.wt, 'float32');

c = fread(fp, 1, 'char');
if ~isempty(c), error 'extra stuff in file?', end

%
% close file
%
if (fclose(fp)), error fclose, end

%
% generate sparse matrix
%
i = 1 + index;
j = zeros(n.wt,1);
n.col = nx * ny;
h = waitbar(0, 'Sparse matrix formation');
for ii=1:n.col
	t = [1:length(ii)] + offset(ii);
	j(t) = ii;
	waitbar(ii/n.col)
end
close(h)
if any(i <= 0) || any(j <= 0) || any(i > nb*na) || any(j > n.col)
	error 'bad indices'
end
if chat
	printf('wt value range [%g,%g]', min(value), max(value))
end
A = sparse(i, j, value, nb*na, n.col);

%
% wtf_read_test()
%
function wtf_read_test
nx = 6;
ny = 4;
nb = 8;
na = 5;
rng(0)
A = rand(nb*na, nx*ny);
A = double(single(A));

file = [test_dir 't.wtf'];
for row_grouped = 0:1
	pr row_grouped
	delete(file)
	wtf_write(file, A, nx, ny, nb, na, ...
		'chat', 0, 'row_grouped', row_grouped)
%	eval(['!wt test ' file])

	if 0 % test old "load" use
		[B mx my mb ma is_transpose] = wtfmex('load', file);
		pr is_transpose
		if is_transpose, B = B'; end
		jf_equal(A, B)
		jf_equal([mx my mb ma], [nx ny nb na])
	end

	if 1 % test new "asp:load" use
		[B mx my mb ma is_transpose] = wtfmex('asp:load', file);
		pr is_transpose
		if is_transpose, B = B'; end
		jf_equal(A, B)
		jf_equal([mx my mb ma], [nx ny nb na])
	end

	if 1 % test wtf_read
		[B mx my mb ma] = wtf_read(file);
		jf_equal(A, B)
		jf_equal([mx my mb ma], [nx ny nb na])
	end
end

ig = image_geom('nx', 22, 'ny', 20, 'dx', 2);
sg = sino_geom('par', 'nb', 24, 'na', 18, 'dr', 1.8);
ig.mask = ig.circ(ig.fov/2) > 0;

list = {'col', 'row'};
for ii=1:2
	gtype = list{ii};
	tmp = aspire_pair(sg, ig, 'support', 'array');
	buff = wtfmex('asp:gensys', tmp', gtype, uint8(ig.mask), int32(0));
	[mat mx my mb ma is_transpose] = wtfmex('asp:mat', buff, int32(0));
	jf_equal([mx my mb ma], [ig.nx ig.ny sg.nb sg.na])
	pr is_transpose
	if is_transpose, mat = mat'; end
	jf_equal(size(mat), [sg.nb*sg.na ig.nx*ig.ny])
end

% todo: test parse version?

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function cS = titanite
% sphene

%
cs_Ttn = crystalSymmetry('2/m',[0.753 1 0.854], [90, 119.72, 90]*degree);
N = Miller({0,1,0},{0,0,1},{1,1,1},{0,1,1},cs_Ttn);
dist = [0.89, 1.39, 1, 1.4];
cS = crystalShape(N./dist);

author: mtex-toolbox

repo: mtex

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


text:
function p = getNMSPenalty(B,b)

p = -0.5*(getMaxIncFloat(B',b)+getIOUFloat(B',b));

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
classdef SPX_MatchedFilter < handle

    properties(SetAccess=private)
        % The signals to be used by the filter
        Signals
    end

    methods
        function self  = SPX_MatchedFilter(signals)
            self.Signals  = signals;
        end

        function [ matchedFilterStatistic ] = apply(self, receivedSequence)
            %APPLY Computes matched filter based statistic on sequence
            signals = self.Signals;

            % Length of each signal to compare (N) and number of signals to compare (L)
            % Each signal is stored in one column
            % There are S such columns
            [N, S] = size(signals); % NxS
            % Number of received samples
            numReceivedSamples = length(receivedSequence);
            % Number of received bits
            B = round(numReceivedSamples/N);
            % We initialize output data
            % Number of rows = number of received bits
            % Each row contains results for each of the signals
            % Number of columns = number of signals
            % We reshape the received sequence
            % each column contains received sequence for one sample
            receivedSequence = reshape(receivedSequence, N, B); % NxB
            matchedFilterStatistic = receivedSequence' * signals;
        end
    end

    methods(Static)
        function [ matchedFilterStatistic ] = filter( receivedSequence,...
            signals )
            mf = SPX_MatchedFilter(signals);
            matchedFilterStatistic = mf.apply(receivedSequence);
        end
    end
end

author: indigits

repo: sparse-plex


text:
%INTBDR Integation of expression over boundaries.
%
%   [ VAL ] = INTBDR( S_EXPR, PROB, IND_B, I_CUB, SOLNUM, IND_S )
%   Integrates the expression S_EXPR over the boundaries indicated in
%   IND_B. PROB is a valid finite element problem struct, and I_CUB
%   specifies the numerical integration rule. IND_S optionally
%   specifies which subdomain to use as reference for internal/
%   interior boundaries (normals point out from these subdomains).
%
%       Input       Value/[Size]           Description
%       -----------------------------------------------------------------------------------
%       s_expr      string                 Expression to integrate
%       prob        struct                 Finite element problem struct
%       ind_b       [1,n_bdr]              Boundary numbers (default all)
%       i_cub       scalar                 Numerical integration rule (default 2)
%       solnum      scalar {n_sols}        Solution number/time to evaluate
%       ind_s       integer array          Integration subdomains for internal boundaries
%                                                                                .
%       Output      Value/[Size]           Description
%       -----------------------------------------------------------------------------------
%       val         scalar                 Result of integration
%
%   See also INTSUBD, MINMAXSUBD, MINMAXBDR

% Copyright 2013-2022 Precise Simulation, Ltd.

author: precise-simulation

repo: featool-multiphysics


text:
function [] = assertEqual(a, b)
  if (a != b)
    testFailed;
  end

function [] = testFailed()
  [ST, I] = dbstack(2);
  disp(strcat("FAILED: ",  ST(1).name));

author: ricedsp


text:
function [L,S,obj,err,iter] = igc(A,C,lambda,opts)

% Reference: Chen, Yudong, Sujay Sanghavi, and Huan Xu. Improved graph clustering.
% IEEE Transactions on Information Theory 60.10 (2014): 6440-6455.
%
% min_{L,S} ||L||_*+lambda*||C \cdot S||_1, s.t. A=L+S, 0<=L<=1.
%
% ---------------------------------------------
% Input:
%       A       -    d*n matrix
%       C       -    d*n matrix
%       lambda  -    >0, parameter
%       opts    -    Structure value in Matlab. The fields are
%           opts.tol        -   termination tolerance
%           opts.max_iter   -   maximum number of iterations
%           opts.mu         -   stepsize for dual variable updating in ADMM
%           opts.max_mu     -   maximum stepsize
%           opts.rho        -   rho>=1, ratio used to increase mu
%           opts.DEBUG      -   0 or 1
%
% Output:
%       L       -    d*n matrix
%       S       -    d*n matrix
%       obj     -    objective function value
%       err     -    residual
%       iter    -    number of iterations
%
% version 1.0 - 19/06/2016
%
% Written by Canyi Lu (canyilu@gmail.com)
%

tol = 1e-8;
max_iter = 500;
rho = 1.1;
mu = 1e-4;
max_mu = 1e10;
DEBUG = 0;

if ~exist('opts', 'var')
    opts = [];
end
if isfield(opts, 'tol');         tol = opts.tol;              end
if isfield(opts, 'max_iter');    max_iter = opts.max_iter;    end
if isfield(opts, 'rho');         rho = opts.rho;              end
if isfield(opts, 'mu');          mu = opts.mu;                end
if isfield(opts, 'max_mu');      max_mu = opts.max_mu;        end
if isfield(opts, 'DEBUG');       DEBUG = opts.DEBUG;          end

C = abs(C);
[d,n] = size(A);

L = zeros(d,n);
S = L;
Z = L;
Y1 = L;
Y2 = L;

iter = 0;
for iter = 1 : max_iter
    Lk = L;
    Sk = S;
    Zk = Z;
    % first super block {L,S}
    [L,nuclearnormL] = prox_nuclear(Z-Y2/mu,1/mu);
    S = prox_l1(-Z+A-Y1/mu,C*(lambda/mu));

    % second super block {Z}
    Z = project_box((-S+A+L+(Y2-Y1)/mu)/2,0,1);

    dY1 = Z+S-A;
    dY2 = L-Z;
    chgL = max(max(abs(Lk-L)));
    chgS = max(max(abs(Sk-S)));
    chgZ = max(max(abs(Zk-Z)));
    chg = max([chgL chgS chgZ max(abs(dY1(:))) max(abs(dY2(:)))]);
    if DEBUG
        if iter == 1 || mod(iter, 10) == 0
            obj = nuclearnormL+lambda*sum(sum(C.*abs(S)));
            err = sqrt(norm(dY1,'fro')^2+norm(dY2,'fro')^2);
            disp(['iter ' num2str(iter) ', mu=' num2str(mu) ...
                    ', obj=' num2str(obj) ', err=' num2str(err)]);
        end
    end

    if chg < tol
        break;
    end
    Y1 = Y1 + mu*dY1;
    Y2 = Y2 + mu*dY2;
    mu = min(rho*mu,max_mu);
end
obj = nuclearnormL+lambda*sum(sum(C.*abs(S)));
err = sqrt(norm(dY1,'fro')^2+norm(dY2,'fro')^2);

author: canyilu

repo: LibADMM-toolbox


text:
function [ids,fullfilenames] = GenerateListOfFilesTBA(ext,datadir,recursive)
%
%[ids,fullfilenames] = GenerateListOfFilesTBA(ext, datadir, recursive)
%
%   OBJECTIVE:
%       This function takes the extension ext and finds all relevant
%       files within the working directory and the subfolders of the
%       working directory. Works with records that are stored with a
%       filename in the format '####.ext'
%
%   INPUT:  ext       - string variable of the extension of filetype
%           	        example: 'mat' 'qrs'
%           datadir   - string variable representing the name of the
%           	        directory of the data
%           recursive -
%
%   OUTPUT: ids          - vector of identifiers (doubles) listing all records
%           fulfilenames - full file name for each record(strings)
%
%	REPO:
%       https://github.com/cliffordlab/PhysioNet-Cardiovascular-Signal-Toolbox
%   ORIGINAL SOURCE AND AUTHORS:
%       Main script written by Adriana N. Vest
%       Dependent scripts written by various authors
%       (see functions for details)
%	COPYRIGHT (C) 2016
%   LICENSE:
%       This software is offered freely and without warranty under
%       the GNU (v3 or later) public license. See license file for
%       more information
%

% Check for folders inside the data directory so that a decision can be
% made about whether to include subfolders or not

if strcmp(computer,'GLNXA64')
    if recursive
    	[~, sysout] = system(['ls ' datadir filesep '*.',ext]);
    else
        [~, sysout] = system(['ls ' datadir filesep '*.',ext]);
        % ANV modified this line to work in linux for MIT NSR files on box
        % May need to revisit
    end
    allfiles = strsplit(sysout);
elseif strcmp(computer,'MACI64') || strcmp(computer,'PCWIN64')
    allfiles1 = [];
    if recursive
        allfiles1 = dir(fullfile(datadir, '*', filesep, ['*.' ext]));
    end
	allfiles2 = dir(fullfile(datadir, ['*.' ext]));
    allfiles = [allfiles1;allfiles2];
else
    error('no files or records found');
end

% Initialize arrays to store data
ids = cell(length(allfiles), 1); % Create ids vector to store subject identifiers
fullfilenames = cell(length(ids),1); % Create cell array to store filenames

for i = 1:length(fullfilenames)
    try
        if strcmp(computer,'GLNXA64')
            filename = dir(allfiles{i});
            fullfilenames{i} = [allfiles{i}];

            filename = filename.name;
            filename = strrep(filename, ['.' ext], '');
            ids{i} = filename;
        elseif strcmp(computer,'MACI64')
            fullfilenames{i} = [allfiles(i).folder filesep allfiles(i).name];
            filename = allfiles(i).name;
            filename = strrep(filename, ['.' ext], '');
            ids{i} = filename;
        elseif strcmp(computer,'PCWIN64')
            fullfilenames{i} = [datadir filesep allfiles(i).name];
            filename = allfiles(i).name;
            filename = strrep(filename, ['.' ext], '');
            ids{i} = filename;
        end
    catch

    end
end

if strcmp(computer,'GLNXA64')
	%idx_rem = find((isempty(ids))); %% This line isn't working on Linux
    emptyCells = cellfun(@isempty,ids);
    idx_rem = find(emptyCells);
    ids(idx_rem) = [];
	fullfilenames(idx_rem) = [];
end

if isnumeric(ids{1})
    for k = 1:length(ids)
        ids{k} = num2str(ids{k});
    end
end

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox


text:
%% Mean-Shift Algorithm
% by Sylvain Bernhardt
% July 2008
%% Description
% Computes the mask of a Parzen window
% and its gradient in respect of the x-axis
% and the y-axis.
% The different types of kernel are:
% {Uniform,Triangular,Epanechnikov,Gaussian}
%
% [k,gx,gy] = Parzen_window(H,W,R,type,graph)
% with:
% k - the mask
% gx,gy - its gradients
% (H,W) - size of the mask
% type - its type
% graph - plot the masks if graph=1

function [k,gx,gy] = Parzen_window(H,W,R,type,graph)
k = zeros(H,W);

%% ----Uniform----
if strcmp(type,'Uniform')==1
    for i=1:H
        for j=1:W
            if (((2*i)/H-1)/R)^2+(((2*j)/W-1)/R)^2 <= 1
                k(i,j) = 1;
            end
        end
    end
end

%% ----Triangular----
if strcmp(type,'Triangular')==1
    Max = max(H,W);
    for z=1:round(R*Max/2)-1
        h = zeros(H,W);
        for i=1:H
            for j=1:W
                if ((i-(H/2))/(R*H/2-z*H/Max))^2+...
                        ((j-(W/2))/(R*W/2-z*W/Max))^2 <= 1
                    h(i,j) = 2/(Max*R);
                end
            end
        end
        k = k+h;
    end
end

%% ----Epanechnikov
if strcmp(type,'Epanechnikov')==1
    for i=1:H
        for j=1:W
            k(i,j) = (1-(2*i/(R*H)-1/R)^2-...
                (2*j/(R*W)-1/R)^2);
            if k(i,j) < 0
                k(i,j) = 0;
            end
        end
    end
end

%% ----Gaussian----
if strcmp(type,'Gaussian')==1
    sigmaH = (R*H/2)/3;
    sigmaW = (R*W/2)/3;
    % sigma = x/3 as a gaussian is almost equal to 0
    % from 3*sigma.
    for i=1:H
        for j=1:W
            k(i,j) = exp(-.5*((i-.5*H)^2/sigmaH^2+...
                (j-.5*W)^2/sigmaW^2));
        end
    end
end

%% Gradient of kernel
[gx,gy] = gradient(-k);

%% Plotting the window
if graph==1
    figure (4)
    scrsz = get(0,'ScreenSize');
    set(4,'Position',[scrsz(3)/4 scrsz(4)/4 ...
        scrsz(3)/1.5 scrsz(4)/1.5])
    subplot(2,2,1)
    mesh(k);
    surf(k);
    shading interp
    subplot(2,2,2)
    mesh(gx);
    surf(gx);
    shading interp
    subplot(2,2,3)
    mesh(gy);
    surf(gy);
    shading interp
end

author: Sable

repo: mcbench-benchmarks


text:
function [ bic ] = BIC_f(loglik,k,nbParamK,N)
%BIC

bic = -2*loglik + k*nbParamK*log(N);

end

author: epfl-lasa


text:
function response = Cautious(n,p)

if n==1 || n==20
    response = 'defect';
elseif p==0 || p==3
    response = 'cooperate';
else
    response = 'defect';
end


text:
function [vect] = getPointingVectFromPitchHeading(rVect, pitchAngle, headingAngle)
%getPointingVectFromPitchHeading Summary of this function goes here
%   Detailed explanation goes here

    zHat = [0,0,1];

    up = rVect/norm(rVect);
    north = cross(cross(up,zHat),up);
    north = north/norm(north);
    west = cross(up,north);
    west = west/norm(west);

    hP = 2*pi - headingAngle;
    vect = cos(pitchAngle)*(cos(hP)*north + sin(hP)*west) + sin(pitchAngle)*up;
    vect = vect/norm(vect);
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
%% AMG TEST V: DIFFERENT INTERPOLATION OPERATORS
%
% We consider the effect of using different interpolation operators to
% interpolate fine grid values by coarse grid values. It is tested through
% the following choices in option.interpolation
%
% * 's' standard interpolation. Use the matrix A_fc as a weighted average
% of all connected coarse nodes.
% * 't' two-points interpolation. Use at most two connected coarse
% nodes.
% * 'a' aggegration (one-point) interpolation. Use the strongest connected
% coarse node.

%%
clear all; close all;
%% Unstructured mesh in 2-D
load lakemesh
showmesh(node,elem);
%%
% Standard interpolation
option.interpolation = 's';
[N,itStep,time,err] = amgtest(node,elem,[],option);
%%
colHeaders = {'Unknowns','Iterations','Time (sec)','Error'};
makeHtmlTable([N itStep time err],[],[],colHeaders,[],6);
%%
clf;
r = showrate(N,time,1);
xlabel('N'); ylabel('Time');
title(['Complexity is N^{' num2str(r) '}'],'Fontsize', 14);
%%
% Two-points interpolation
load lakemesh
option.interpolation = 't';
[N,itStep,time,err] = amgtest(node,elem,[],option);
%%
colHeaders = {'Unknowns','Iterations','Time (sec)','Error'};
makeHtmlTable([N itStep time err],[],[],colHeaders,[],6);
%%
clf;
r = showrate(N,time,1);
xlabel('N'); ylabel('Time');
title(['Complexity is N^{' num2str(r) '}'],'Fontsize', 14);
%%
% One point interpolation
load lakemesh
option.interpolation = 'a';
[N,itStep,time,err] = amgtest(node,elem,[],option);
%%
colHeaders = {'Unknowns','Iterations','Time (sec)','Error'};
makeHtmlTable([N itStep time err],[],[],colHeaders,[],6);
%%
clf;
r = showrate(N,time,1);
xlabel('N'); ylabel('Time');
title(['Complexity is N^{' num2str(r) '}'],'Fontsize', 14);

%% Unstructured mesh in 3-D
clear all; close all
load bunny;
showboundary3(node,elem);
view([-179 74]);
%%
% Standard interpolation
option.interpolation = 's';
[N,itStep,time,err] = amgtest3(node,elem,1,option);
%%
colHeaders = {'Unknowns','Iterations','Time (sec)','Error'};
makeHtmlTable([N itStep time err],[],[],colHeaders,[],6);
%%
clf;
r = showrate(N,time,1);
xlabel('N'); ylabel('Time');
title(['Complexity is N^{' num2str(r) '}'],'Fontsize', 14);
%%
% Two-points interpolation
load bunny
option.interpolation = 't';
[N,itStep,time,err] = amgtest3(node,elem,1,option);
%%
colHeaders = {'Unknowns','Iterations','Time (sec)','Error'};
makeHtmlTable([N itStep time err],[],[],colHeaders,[],6);
%%
clf;
r = showrate(N,time,1);
xlabel('N'); ylabel('Time');
title(['Complexity is N^{' num2str(r) '}'],'Fontsize', 14);
%%
% One point interpolation
load bunny
option.interpolation = 'a';
[N,itStep,time,err] = amgtest3(node,elem,1,option);
%%
colHeaders = {'Unknowns','Iterations','Time (sec)','Error'};
makeHtmlTable([N itStep time err],[],[],colHeaders,[],6);
%%
clf;
r = showrate(N,time,1);
xlabel('N'); ylabel('Time');
title(['Complexity is N^{' num2str(r) '}'],'Fontsize', 14);

%% Conclusion
% In general, if we use more coarse grids, we get more accurate
% interpolation and require less iteration steps. But the matrix on coarse
% level becomes denser, see nnz/Nc, which requires, indeed, more
% computational time in both smoother and coarse grid solver. In the
% extremly case, the one-point interpolation almost keeps the sparsity but
% iteration steps increase as level increases in the speed of J. The
% average using all neighboring points keeps the iteration steps but the
% sparsity increase a lot, especially in 3-D. The two-points interpolation
% seems a good balance.

repo: Aero-matlab


text:
function p=gaussian(x,m,C);
% p=gaussian(x,m,C);
%
% Evaluate the multi-variate density with mean vector m and covariance
% matrix C for the input vector x.
%
% p=gaussian(X,m,C);
%
% Vectorized version: Here X is a matrix of column vectors, and p is
% a vector of probabilities for each vector.
% Jianbo Shi, 1997
d=length(m);

if size(x,1)~=d
   x=x';
end
N=size(x,2);

detC = det(C);
if rcond(C)<eps
%   fprintf(1,'Covariance matrix close to singular. (gaussian.m)\n');
   p = zeros(N,1);
else
   m=m(:);
   M=m*ones(1,N);
   denom=(2*pi)^(d/2)*sqrt(abs(detC));
   mahal=sum(((x-M)'*inv(C)).*(x-M)',2);   % Chris Bregler's trick
   numer=exp(-0.5*mahal);
   p=numer/denom;
end

author: jwyang


text:
function [ns_RESULT, LastError] = ns_GetLastErrorMsg();

%ns_GetLastErrorMsg   Retrieves the extended error message for the last
%error
%
%   Usage:
%      [ns_RESULT, LastError] = ns_GetLastErrorMsg
%
%   Description:
%       Returns the last error message in text form to LastError. This
%       function should be called immediately following a function whose
%       return value indicates that an error occurred. The maximum size of
%       the error message text is 256 characters.
%
%   Return Values:
%       LastError	Variable to receive the extended last error message.
%       ns_RESULT   This function returns ns_OK if the file is successfully
%                   opened. Otherwise one of the following error codes is
%                   generated:
%
%                       ns_LIBERROR	Library error
%
%   Copyright (C) 2003 Neuroshare Project
%   Author: Almut Branner
%   Last modification: 2/21/2003

[ns_RESULT, LastError] = mexprog(17);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%% Orientation Dependent Functions
%
%%
% An orientation dependent function is a function that assigns to each
% rotation or orientation a numerical value. An import example of a
% rotational function is the <ODFTheory.html orientation density function
% (ODF)> that assignes to each crystal orientation the probability of its
% occurence within a specimen. Other examples are the Schmidt or the Taylor
% factor as a function of the crystal orientation.
%
%% Definition of a orientation dependent function
%
% Within MTEX a rotational function is represented by a variable of type
% <SO3Fun.SO3Fun.html |SO3Fun|>. Let us consider as an example the function
% that takes an orientation and returns it rotational angle modulo cubic
% crystal symmetry. In MTEX the rotational angle is computed by the command
% <orientation.angle.html |angle(ori)|>. In order to turn this
% correspondence into a |SO3Fun| we use the command @SO3FunHandle and pass
% the angle command as an
% <https://de.mathworks.com/help/matlab/matlab_prog/anonymous-functions.html
% anonymous function>.

% define the crystal symmetry
cs = crystalSymmetry('cubic');

% construct the SO3Fun
SO3F = SO3FunHandle(@(ori) angle(ori) ./ degree, cs)

%%
% Many more methods for defining orientation dependent functions are
% discussed <SO3FunDefinition.html here>.
%
%%
% The entire information abot the orientation dependent function is now
% stored in the variable |SO3F|. In order to determine its value for a
% specific orientation |ori| the function <SO3Fun.eval.html |eval(ori)|> is
% used

ori = orientation.rand(cs)
SO3F.eval(ori)

%% Plotting an orientation Dependent Function
%
% Orientation dependent functions are most of visualized by sections
% according to the third Euler angle $\varphi_2$
%

plotSection(SO3F)

%%
% The plot tells us for which Euler angles the the resulting rotational
% angle is large and for which Euler angles it is low. The plot of this
% "angle function" |SO3F| becomes trivial if represented in an axis angle
% sections

plotSection(SO3F,'axisAngle','upper')
mtexColorbar
mtexColorMap parula

%%
% as obviously, the function value is constant in each section.
% Many more methods for visualizing orientation dependent functions are
% discussed <SO3FunVisualization.html here>.
%
%% Computing with orientation dependent functions
%
% The power of representing an orientation dependent functions as a
% variables of type @SO3Fun is that we may apply to it a
% <SO3FunOperations.html large number of analysis tools>. In particular,
% one can add, subtract and mutiply orientation dependent functions, plot
% them in various projections or detect the local minima or maxima. In the
% case of our example function the local maxima refers to the orientations
% with maximum rotational angle in cubic symmetry. We may compute them by
% the command <SO3Fun.max.html |max|>

[value,ori] = max(SO3F,'numLocal',10)

%%
% We observe that there are exactly six symetrically not equivalent
% orientations that realize an orientation angle of about 62.77 degree and
% form the vertices of the fundamental region in orienation space

plot(ori.symmetrise,'axisAngle','filled','markerSize',15,'restrict2FundamentalRegion')

%% Representations of Rotational Functions
%
% Internally MTEX represents rotational functions in different ways:
%
% || by a harmonic series expansion || <SO3FunHarmonicRepresentation.html SO3FunHarmonic> ||
% || as superposition of radial function || <RadialODFs.html SO3FunRBF> ||
% || as superposition of fibre elements || <FibreODFs.html SO3FunCBF> ||
% || as Bingham distribution || <BinghamODFs.html SO3FunBingham> ||
% || as sum of different components || @SO3FunComposition ||
% || explicitely given by a formula || @SO3FunHandle ||
%
% All representations allow the same operations which are specified for
% the abstact class |@SO3Fun|. In particular it is possible
% to calculate with $SO(3)$ functions as with ordinary numbers, i.e., you
% can add, multiply arbitrary functions, take the mean, integrate them or
% compute gradients, see <SO3FunOperations.html Operations>.
%
%% Generalizations of Rotational Functions
%
% || rotational vector fields || <SO3FunVectorField.html SO3VectorField> ||
% || radial rotational functions || <SO3Kernels.html SO3Kernel> ||
%


text:
% collect the training data annotations (from the menpo challenge)
clear
train_data_loc = 'C:\Users\tbaltrus\Documents\menpo_data_orig/';

dataset_locs = {[train_data_loc, '/train/'];};

left_to_frontal_map = [17,28; 18,29; 19,30; 20,31;
                       21,34; 22,32; 23,39; 24,38; 25,37; 26,42; 27,41;
                       28,52; 29,51; 30,50; 31,49; 32,60; 33,59; 34,58;
                       35,63; 36,62; 37,61; 38,68; 39,67];

right_to_frontal_map = [17,28; 18,29; 19,30; 20,31;
                       21,34; 22,36; 23,44; 24,45; 25,46; 26,47; 27,48;
                       28,52; 29,53; 30,54; 31,55; 32,56; 33,57; 34,58;
                       35,63; 36,64; 37,65; 38,66; 39,67];

all_pts_left = [];
all_pts_right = [];
all_pts_orig_left = [];
all_pts_orig_right = [];

for i=1:numel(dataset_locs)
    landmarkLabels = dir([dataset_locs{i} '\*.pts']);
    landmarkImgs = dir([dataset_locs{i} '\*.jpg']);

    for p=1:numel(landmarkLabels)
        landmarks = importdata([dataset_locs{i}, landmarkLabels(p).name], ' ', 3);
        img = imread([dataset_locs{i}, landmarkImgs(p).name]);
        landmarks = landmarks.data;
        landmark_labels = -ones(68,2);

        if(size(landmarks,1) == 39)
            % Determine if the points are clock-wise or counter clock-wise
            % Clock-wise points are facing left, counter-clock-wise right
            sum = 0;
            for k=1:11
                step = (landmarks(k+1,1) - landmarks(k,1)) * (landmarks(k+1,2) + landmarks(k,2));
                sum = sum + step;
            end

            if(sum > 0)
                % First need to resample the face outline as there are 9
                % points in the near-frontal and 10 points in profile for
                % the outline of the face

                outline = iterate_piece_wise(landmarks(1:10,:), 9);
                brow = iterate_piece_wise(landmarks(13:16,:), 5);
                landmark_labels(1:9,:) = outline;
                landmark_labels(18:22,:) = brow;
                landmark_labels(left_to_frontal_map(:,2),:) = landmarks(left_to_frontal_map(:,1),:);

                all_pts_orig_right = cat(3, all_pts_orig_right, landmarks);
                all_pts_right = cat(3, all_pts_right, landmark_labels);
            else
                outline = iterate_piece_wise(landmarks(10:-1:1,:), 9);
                brow = iterate_piece_wise(landmarks(16:-1:13,:), 5);

                landmark_labels(9:17,:) = outline;
                landmark_labels(23:27,:) = brow;

                landmark_labels(right_to_frontal_map(:,2),:) = landmarks(right_to_frontal_map(:,1),:);

                all_pts_orig_left = cat(3, all_pts_orig_left, landmarks);
                all_pts_left = cat(3, all_pts_left, landmark_labels);
            end
        end

        if(mod(p,50) == 0)
            fprintf('%d/%d\n', p, numel(landmarkLabels))
        end
    end
end

%%

save('menpo_68_pts_train_profile', 'all_pts_orig_right', 'all_pts_right', 'all_pts_orig_left', 'all_pts_left');

author: TadasBaltrusaitis

repo: OpenFace


text:
function [rad] = deg2rad(deg)
%deg2rad Summary of this function goes here
%   Detailed explanation goes here
    rad=deg*(pi/180);
end


text:
% example_mri_rf_spsp
% This script is an example of spectral-spatial pulse design for
% through-plane phase precompensatory slice selection for T2*-weighted
% functional MR, based on our publication:
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Spectral-spatial RF pulse design script, based on "Spectral-spatial RF
% pulse design for through-plane phase precompensatory slice selection for
% T2*-weighted functional MRI", Chun-yu Yip et al, Magnetic Resonance in
% Medicine, May 2009.  doi 10.1002/mrm.21938
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Written by Chun-yu Yip, University of Michigan, Ann Arbor, 4/1/2009
%
% Running this script requires that the entire image reconstruction toolbox
% by Professor Jeffrey A. Fessler be installed first:
%
% http://web.eecs.umich.edu/~fessler/code/index.html
%
% You can conveniently download the whole package at
%
% http://www.eecs.umich.edu/~fessler/irt/fessler.tgz
%
% and use the unix command 'tar' to retrieve the files and folders.
% Use the "setup.m" command in the toolbox to set up paths first.

% Load pulse design parameters into matlab structures kp, rfp, iop.
if ~isvar('kp'), printm 'Loading parameters...'
	[kp, iop, rfp] = example_spsp_param1();
	iop.writetofile_scanner = false; % comment out to create scanner files
%	iop.writetofile_sim = true;
end

% Design z-gradient waveform, based on parameters in kp.
if ~isvar('gz'), printm 'Designing z gradient waveform...'
	[kp gz kz kf t] = compute_gz_spsp(kp);

	if 1 % display gradient waveform
		clf, pl = @(i) subplot(340+i);
		pl(1)
		plot(t*1000, gz), axis tight
		xlabel('time (msec)')
		ylabel('g/cm')
		title('Gz waveform')
		grid
	drawnow
	end
end

% Design complex-valued RF waveform iteratively using conjugate gradient
if ~isvar('b'), printm 'Computing RF pulse waveform...'
	[b d f z mm] = compute_rf_spsp_mgh(kp, rfp, gz, kz, kf);

	ir_mri_rf_spsp_plot(b, gz, d, f, z, mm, t) % display results
	drawnow

	% Write computed waveforms to files for simulation and/or scanner.
	if iop.writetofile_scanner || iop.writetofile_sim
		iop.waveformfilespath = './o/'; % put output in a subdirectory
		if ~exist(iop.waveformfilespath, 'dir')
			fail('must create output directory "%s"', iop.waveformfilespath)
		end
		write2files_spsp(kp, rfp, iop, gz, b);
	end
end

% Perform Bloch simulation in SPSP space
if ~isvar('mresult'), printm 'Performing Bloch simulation...'
	mresult = dosim7_spsp(kp, rfp, iop);
end


text:
function soln = directCollocation(problem)
% soln = directCollocation(problem)
%
% OptimTraj utility function
%
% This function is designed to be called by either "trapezoid" or
% "hermiteSimpson". It actually calls FMINCON to solve the trajectory
% optimization problem.
%
% Analytic gradients are supported.
%
% NOTES:
%
%   If analytic gradients are used, then the sparsity pattern is returned
%   in the struct: soln.info.sparsityPattern. View it using spy().
%

%To make code more readable
G = problem.guess;
B = problem.bounds;
F = problem.func;
Opt = problem.options;

nGrid = length(F.weights);

flagGradObj = strcmp(Opt.nlpOpt.GradObj,'on');
flagGradCst = strcmp(Opt.nlpOpt.GradConstr,'on');

% Print out notice about analytic gradients
if Opt.verbose > 0
    if flagGradObj
        fprintf('      - using analytic gradients of objective function\n');
    end
    if flagGradCst
        fprintf('      - using analytic gradients of constraint function\n');
    end
    fprintf('\n');
end

% Interpolate the guess at the grid-points for transcription:
guess.tSpan = G.time([1,end]);
guess.time = linspace(guess.tSpan(1), guess.tSpan(2), nGrid);
guess.state = interp1(G.time', G.state', guess.time')';
guess.control = interp1(G.time', G.control', guess.time')';

[zGuess, pack] = packDecVar(guess.time, guess.state, guess.control);

if flagGradCst || flagGradObj
    gradInfo = grad_computeInfo(pack);
end

% Unpack all bounds:
tLow = linspace(B.initialTime.low, B.finalTime.low, nGrid);
xLow = [B.initialState.low, B.state.low*ones(1,nGrid-2), B.finalState.low];
uLow = B.control.low*ones(1,nGrid);
zLow = packDecVar(tLow,xLow,uLow);

tUpp = linspace(B.initialTime.upp, B.finalTime.upp, nGrid);
xUpp = [B.initialState.upp, B.state.upp*ones(1,nGrid-2), B.finalState.upp];
uUpp = B.control.upp*ones(1,nGrid);
zUpp = packDecVar(tUpp,xUpp,uUpp);

%%%% Set up problem for fmincon:
if flagGradObj
    P.objective = @(z)( ...
        myObjGrad(z, pack, F.pathObj, F.bndObj, F.weights, gradInfo) );   %Analytic gradients
    [~, objGradInit] = P.objective(zGuess);
    sparsityPattern.objective = (objGradInit~=0)';  % Only used for visualization!
else
    P.objective = @(z)( ...
        myObjective(z, pack, F.pathObj, F.bndObj, F.weights) );   %Numerical gradients
end
if flagGradCst
    P.nonlcon = @(z)( ...
        myCstGrad(z, pack, F.dynamics, F.pathCst, F.bndCst, F.defectCst, gradInfo) ); %Analytic gradients
    [~,~,cstIneqInit,cstEqInit] = P.nonlcon(zGuess);
    sparsityPattern.equalityConstraint = (cstEqInit~=0)';  % Only used for visualization!
    sparsityPattern.inequalityConstraint = (cstIneqInit~=0)';  % Only used for visualization!
else
    P.nonlcon = @(z)( ...
        myConstraint(z, pack, F.dynamics, F.pathCst, F.bndCst, F.defectCst) ); %Numerical gradients
end

P.x0 = zGuess;
P.lb = zLow;
P.ub = zUpp;
P.Aineq = []; P.bineq = [];
P.Aeq = []; P.beq = [];
P.options = Opt.nlpOpt;
P.solver = 'fmincon';

%%%% Call fmincon to solve the non-linear program (NLP)
tic;
[zSoln, objVal,exitFlag,output] = fmincon(P);
[tSoln,xSoln,uSoln] = unPackDecVar(zSoln,pack);
nlpTime = toc;

%%%% Store the results:

soln.grid.time = tSoln;
soln.grid.state = xSoln;
soln.grid.control = uSoln;

soln.interp.state = @(t)( interp1(tSoln',xSoln',t','linear',nan)' );
soln.interp.control = @(t)( interp1(tSoln',uSoln',t','linear',nan)' );

soln.info = output;
soln.info.nlpTime = nlpTime;
soln.info.exitFlag = exitFlag;
soln.info.objVal = objVal;

if flagGradCst || flagGradObj  % Then return sparsity pattern for visualization
    if flagGradObj
        [~, objGradInit] = P.objective(zSoln);
        sparsityPattern.objective = (objGradInit~=0)';
    end
    if flagGradCst
        [~,~,cstIneqInit,cstEqInit] = P.nonlcon(zSoln);
        sparsityPattern.equalityConstraint = (cstEqInit~=0)';
        sparsityPattern.inequalityConstraint = (cstIneqInit~=0)';
    end
    soln.info.sparsityPattern = sparsityPattern;
end

soln.problem = problem;  % Return the fully detailed problem struct

end

%%%%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%%%%
%%%%                          SUB FUNCTIONS                            %%%%
%%%%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%%%%

function [z,pack] = packDecVar(t,x,u)
%
% This function collapses the time (t), state (x)
% and control (u) matricies into a single vector
%
% INPUTS:
%   t = [1, nTime] = time vector (grid points)
%   x = [nState, nTime] = state vector at each grid point
%   u = [nControl, nTime] = control vector at each grid point
%
% OUTPUTS:
%   z = column vector of 2 + nTime*(nState+nControl) decision variables
%   pack = details about how to convert z back into t,x, and u
%       .nTime
%       .nState
%       .nControl
%

nTime = length(t);
nState = size(x,1);
nControl = size(u,1);

tSpan = [t(1); t(end)];
xCol = reshape(x, nState*nTime, 1);
uCol = reshape(u, nControl*nTime, 1);

indz = reshape(2+(1:numel(u)+numel(x)),nState+nControl,nTime);

% index of time, state, control variables in the decVar vector
tIdx = 1:2;
xIdx = indz(1:nState,:);
uIdx = indz(nState+(1:nControl),:);

% decision variables
% variables are indexed so that the defects gradients appear as a banded
% matrix
z = zeros(2+numel(indz),1);
z(tIdx(:),1) = tSpan;
z(xIdx(:),1) = xCol;
z(uIdx(:),1) = uCol;

pack.nTime = nTime;
pack.nState = nState;
pack.nControl = nControl;
pack.tIdx = tIdx;
pack.xIdx = xIdx;
pack.uIdx = uIdx;

end

%%%%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%%%%

function [t,x,u] = unPackDecVar(z,pack)
%
% This function unpacks the decision variables for
% trajectory optimization into the time (t),
% state (x), and control (u) matricies
%
% INPUTS:
%   z = column vector of 2 + nTime*(nState+nControl) decision variables
%   pack = details about how to convert z back into t,x, and u
%       .nTime
%       .nState
%       .nControl
%
% OUTPUTS:
%   t = [1, nTime] = time vector (grid points)
%   x = [nState, nTime] = state vector at each grid point
%   u = [nControl, nTime] = control vector at each grid point
%

nTime = pack.nTime;
nState = pack.nState;
nControl = pack.nControl;

t = linspace(z(1),z(2),nTime);

x = z(pack.xIdx);
u = z(pack.uIdx);

% make sure x and u are returned as vectors, [nState,nTime] and
% [nControl,nTime]
x = reshape(x,nState,nTime);
u = reshape(u,nControl,nTime);

end

%%%%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%%%%

function cost = myObjective(z,pack,pathObj,bndObj,weights)
%
% This function unpacks the decision variables, sends them to the
% user-defined objective functions, and then returns the final cost
%
% INPUTS:
%   z = column vector of decision variables
%   pack = details about how to convert decision variables into t,x, and u
%   pathObj = user-defined integral objective function
%   endObj = user-defined end-point objective function
%
% OUTPUTS:
%   cost = scale cost for this set of decision variables
%

[t,x,u] = unPackDecVar(z,pack);

% Compute the cost integral along trajectory
if isempty(pathObj)
    integralCost = 0;
else
    dt = (t(end)-t(1))/(pack.nTime-1);
    integrand = pathObj(t,x,u);  %Calculate the integrand of the cost function
    integralCost = dt*integrand*weights;  %Trapazoidal integration
end

% Compute the cost at the boundaries of the trajectory
if isempty(bndObj)
    bndCost = 0;
else
    t0 = t(1);
    tF = t(end);
    x0 = x(:,1);
    xF = x(:,end);
    bndCost = bndObj(t0,x0,tF,xF);
end

cost = bndCost + integralCost;

end

%%%%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%%%%

function [c, ceq] = myConstraint(z,pack,dynFun, pathCst, bndCst, defectCst)
%
% This function unpacks the decision variables, computes the defects along
% the trajectory, and then evaluates the user-defined constraint functions.
%
% INPUTS:
%   z = column vector of decision variables
%   pack = details about how to convert decision variables into t,x, and u
%   dynFun = user-defined dynamics function
%   pathCst = user-defined constraints along the path
%   endCst = user-defined constraints at the boundaries
%
% OUTPUTS:
%   c = inequality constraints to be passed to fmincon
%   ceq = equality constraints to be passed to fmincon
%

[t,x,u] = unPackDecVar(z,pack);

%%%% Compute defects along the trajectory:
dt = (t(end)-t(1))/(length(t)-1);
f = dynFun(t,x,u);
defects = defectCst(dt,x,f);

%%%% Call user-defined constraints and pack up:
[c, ceq] = collectConstraints(t,x,u,defects, pathCst, bndCst);

end

%%%%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%%%%
%%%%            Additional Sub-Functions for Gradients                 %%%%
%%%%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%%%%

%%%% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ %%%%

function gradInfo = grad_computeInfo(pack)
%
% This function computes the matrix dimensions and indicies that are used
% to map the gradients from the user functions to the gradients needed by
% fmincon. The key difference is that the gradients in the user functions
% are with respect to their input (t,x,u) or (t0,x0,tF,xF), while the
% gradients for fmincon are with respect to all decision variables.
%
% INPUTS:
%   nDeVar = number of decision variables
%   pack = details about packing and unpacking the decision variables
%       .nTime
%       .nState
%       .nControl
%
% OUTPUTS:
%   gradInfo = details about how to transform gradients
%

nTime = pack.nTime;
nState = pack.nState;
nControl = pack.nControl;
nDecVar = 2 + nState*nTime + nControl*nTime;

zIdx = 1:nDecVar;
gradInfo.nDecVar = nDecVar;
[tIdx, xIdx, uIdx] = unPackDecVar(zIdx,pack);
gradInfo.tIdx = tIdx([1,end]);
gradInfo.xuIdx = [xIdx;uIdx];

%%%% Compute gradients of time:
% alpha = (0..N-1)/(N-1)
% t = alpha*tUpp + (1-alpha)*tLow
alpha = (0:(nTime-1))/(nTime-1);
gradInfo.alpha = [1-alpha; alpha];

if (gradInfo.tIdx(1)~=1 || gradInfo.tIdx(end)~=2)
    error('The first two decision variables must be the initial and final time')
end
gradInfo.dtGrad = [-1; 1]/(nTime-1);

%%%% Compute gradients of state
gradInfo.xGrad = zeros(nState,nTime,nDecVar);
for iTime=1:nTime
    for iState=1:nState
        gradInfo.xGrad(iState,iTime,xIdx(iState,iTime)) = 1;
    end
end

%%%% For unpacking the boundary constraints and objective:
gradInfo.bndIdxMap = [tIdx(1); xIdx(:,1); tIdx(end); xIdx(:,end)];

end

%%%% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ %%%%

function [c, ceq, cGrad, ceqGrad] = grad_collectConstraints(t,x,u,defects, defectsGrad, pathCst, bndCst, gradInfo)
% [c, ceq, cGrad, ceqGrad] = grad_collectConstraints(t,x,u,defects, defectsGrad, pathCst, bndCst, gradInfo)
%
% OptimTraj utility function.
%
% Collects the defects, calls user-defined constraints, and then packs
% everything up into a form that is good for fmincon. Additionally, it
% reshapes and packs up the gradients of these constraints.
%
% INPUTS:
%   t = time vector
%   x = state matrix
%   u = control matrix
%   defects = defects matrix
%   pathCst = user-defined path constraint function
%   bndCst = user-defined boundary constraint function
%
% OUTPUTS:
%   c = inequality constraint for fmincon
%   ceq = equality constraint for fmincon
%

ceq_dyn = reshape(defects,numel(defects),1);
ceq_dynGrad = grad_flattenPathCst(defectsGrad);

%%%% Compute the user-defined constraints:
if isempty(pathCst)
    c_path = [];
    ceq_path = [];
    c_pathGrad = [];
    ceq_pathGrad = [];
else
    [c_pathRaw, ceq_pathRaw, c_pathGradRaw, ceq_pathGradRaw] = pathCst(t,x,u);
    c_path = reshape(c_pathRaw,numel(c_pathRaw),1);
    ceq_path = reshape(ceq_pathRaw,numel(ceq_pathRaw),1);
    c_pathGrad = grad_flattenPathCst(grad_reshapeContinuous(c_pathGradRaw,gradInfo));
    ceq_pathGrad = grad_flattenPathCst(grad_reshapeContinuous(ceq_pathGradRaw,gradInfo));
end
if isempty(bndCst)
    c_bnd = [];
    ceq_bnd = [];
    c_bndGrad = [];
    ceq_bndGrad = [];
else
    t0 = t(1);
    tF = t(end);
    x0 = x(:,1);
    xF = x(:,end);
    [c_bnd, ceq_bnd, c_bndGradRaw, ceq_bndGradRaw] = bndCst(t0,x0,tF,xF);
    c_bndGrad = grad_reshapeBoundary(c_bndGradRaw,gradInfo);
    ceq_bndGrad = grad_reshapeBoundary(ceq_bndGradRaw,gradInfo);
end

%%%% Pack everything up:
c = [c_path;c_bnd];
ceq = [ceq_dyn; ceq_path; ceq_bnd];

cGrad = [c_pathGrad;c_bndGrad]';
ceqGrad = [ceq_dynGrad; ceq_pathGrad; ceq_bndGrad]';

end

%%%% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ %%%%

function C = grad_flattenPathCst(CC)
%
% This function takes a path constraint and reshapes the first two
% dimensions so that it can be passed to fmincon
%
if isempty(CC)
    C = [];
else
    [n1,n2,n3] = size(CC);
    C = reshape(CC,n1*n2,n3);
end

end

%%%% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ %%%%

function CC = grad_reshapeBoundary(C,gradInfo)
%
% This function takes a boundary constraint or objective from the user
% and expands it to match the full set of decision variables
%

CC = zeros(size(C,1),gradInfo.nDecVar);
CC(:,gradInfo.bndIdxMap) = C;

end

%%%% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ %%%%

function grad = grad_reshapeContinuous(gradRaw,gradInfo)
% grad = grad_reshapeContinuous(gradRaw,gradInfo)
%
% OptimTraj utility function.
%
% This function converts the raw gradients from the user function into
% gradients with respect to the decision variables.
%
% INPUTS:
%   stateRaw = [nOutput,nInput,nTime]
%
% OUTPUTS:
%   grad = [nOutput,nTime,nDecVar]
%

if isempty(gradRaw)
    grad = [];
else
    [nOutput, ~, nTime] = size(gradRaw);

    grad = zeros(nOutput,nTime,gradInfo.nDecVar);

    % First, loop through and deal with time.
    timeGrad = gradRaw(:,1,:); timeGrad = permute(timeGrad,[1,3,2]);
    for iOutput=1:nOutput
        A = ([1;1]*timeGrad(iOutput,:)).*gradInfo.alpha;
        grad(iOutput,:,gradInfo.tIdx) = permute(A,[3,2,1]);
    end

    % Now deal with state and control:
    for iOutput=1:nOutput
        for iTime=1:nTime
            B = gradRaw(iOutput,2:end,iTime);
            grad(iOutput,iTime,gradInfo.xuIdx(:,iTime)) = permute(B,[3,1,2]);
        end
    end
end

end

%%%% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ %%%%

function [cost, costGrad] = myObjGrad(z,pack,pathObj,bndObj,weights,gradInfo)
%
% This function unpacks the decision variables, sends them to the
% user-defined objective functions, and then returns the final cost
%
% INPUTS:
%   z = column vector of decision variables
%   pack = details about how to convert decision variables into t,x, and u
%   pathObj = user-defined integral objective function
%   endObj = user-defined end-point objective function
%
% OUTPUTS:
%   cost = scale cost for this set of decision variables
%

%Unpack the decision variables:
[t,x,u] = unPackDecVar(z,pack);

% Time step for integration:
dt = (t(end)-t(1))/(length(t)-1);
dtGrad = gradInfo.dtGrad;
nTime = length(t);
nState = size(x,1);
nControl = size(u,1);
nDecVar = length(z);

% Compute the cost integral along the trajectory
if isempty(pathObj)
    integralCost = 0;
    integralCostGrad = zeros(nState+nControl,1);
else

    % Objective function integrand and gradients:
    [obj, objGradRaw] = pathObj(t,x,u);
    nInput = size(objGradRaw,1);
    objGradRaw = reshape(objGradRaw,1,nInput,nTime);
    objGrad = grad_reshapeContinuous(objGradRaw,gradInfo);

    % integral objective function
    unScaledIntegral = obj*weights;
    integralCost = dt*unScaledIntegral;

    % Gradient of integral objective function
    dtGradTerm = zeros(1,nDecVar);
    dtGradTerm(1) = dtGrad(1)*unScaledIntegral;
    dtGradTerm(2) = dtGrad(2)*unScaledIntegral;
    objGrad = reshape(objGrad,nTime,nDecVar);
    integralCostGrad = ...
        dtGradTerm + ...
        dt*sum(objGrad.*(weights*ones(1,nDecVar)),1);
end

% Compute the cost at the boundaries of the trajectory
if isempty(bndObj)
    bndCost = 0;
    bndCostGrad = zeros(1,nDecVar);
else
    t0 = t(1);
    tF = t(end);
    x0 = x(:,1);
    xF = x(:,end);
    [bndCost, bndCostGradRaw] = bndObj(t0,x0,tF,xF);
    bndCostGrad = grad_reshapeBoundary(bndCostGradRaw,gradInfo);
end

% Cost function
cost = bndCost + integralCost;

% Gradients
costGrad = bndCostGrad + integralCostGrad;

end

%%%%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%%%%

function [c, ceq, cGrad, ceqGrad] = myCstGrad(z,pack,dynFun, pathCst, bndCst, defectCst, gradInfo)
%
% This function unpacks the decision variables, computes the defects along
% the trajectory, and then evaluates the user-defined constraint functions.
%
% INPUTS:
%   z = column vector of decision variables
%   pack = details about how to convert decision variables into t,x, and u
%   dynFun = user-defined dynamics function
%   pathCst = user-defined constraints along the path
%   endCst = user-defined constraints at the boundaries
%
% OUTPUTS:
%   c = inequality constraints to be passed to fmincon
%   ceq = equality constraints to be passed to fmincon
%

%Unpack the decision variables:
[t,x,u] = unPackDecVar(z,pack);

% Time step for integration:
dt = (t(end)-t(1))/(length(t)-1);
dtGrad = gradInfo.dtGrad;

% Gradient of the state with respect to decision variables
xGrad = gradInfo.xGrad;

%%%% Compute defects along the trajectory:
[f, fGradRaw] = dynFun(t,x,u);
fGrad = grad_reshapeContinuous(fGradRaw,gradInfo);

[defects, defectsGrad] = defectCst(dt,x,f,...
    dtGrad, xGrad, fGrad);

% Compute gradients of the user-defined constraints and then pack up:
[c, ceq, cGrad, ceqGrad] = grad_collectConstraints(t,x,u,...
    defects, defectsGrad, pathCst, bndCst, gradInfo);

end

author: MatthewPeterKelly

repo: OptimTraj


text:
function box = intersectBoxes3d(box1, box2)
%INTERSECTBOXES3D Intersection of two 3D bounding boxes
%
%   RES = intersectBoxes3d(BOX1, BOX2)
%
%   Example
%   box1 = [5 20 5 30 10 50];
%   box2 = [0 15 0 15 0 20];
%   intersectBoxes3d(box1, box2)
%   ans =
%       5 15 5 15 10 20
%
%   See also
%   boxes3d, drawBox3d, mergeBoxes3d
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2010-07-26,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2010 INRA - Cepia Software Platform.

% unify sizes of data
if size(box1,1) == 1
    box1 = repmat(box1, size(box2,1), 1);
elseif size(box2, 1) == 1
    box2 = repmat(box2, size(box1,1), 1);
elseif size(box1,1) ~= size(box2,1)
    error('Bad size for inputs');
end

% compute extreme coords
mini = min(box1(:,2:2:end), box2(:,2:2:end));
maxi = max(box1(:,1:2:end), box2(:,1:2:end));

% concatenate result into a new box structure
box = [maxi(:,1) mini(:,1) maxi(:,2) mini(:,2) maxi(:,3) mini(:,3)];

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function cvx_expert_check( fname, varargin )

global cvx___
if ~cvx___.expert && ( isempty(varargin) || any( cellfun( 'isclass', varargin, 'cvx' ) ) ),
    fprintf( 1, [ 'CVX Warning:\n', ...
        '   Models involving "%s" or other functions in the log, exp, and entropy\n', ...
        '   family are solved using an experimental successive approximation method.\n', ...
        '   The method requires multiple calls to the solver, so it can be slow; and\n', ...
        '   in certain cases it fails to converge. See Appendix D of the the user''s\n', ...
        '   guide for more information about this method, and for instructions on how\n', ...
        '   to suppress this warning in the future.\n' ], fname );
    cvx___.expert = true;
end

repo: REST


text:
function [ obj ] = bold_grad_cd( obj, n )
% Calculate jacobian matrix of predicted fMRI BOLD signal with respect to
% DCM parameters using central difference method.
%
% This is a protected method of the tapas_Huge class. It cannot be called
% from outside the class.
%

% Author: Yu Yao (yao@biomed.ee.ethz.ch)
% Copyright (C) 2019 Translational Neuromodeling Unit
%                    Institute for Biomedical Engineering,
%                    University of Zurich and ETH Zurich.
%
% This file is part of TAPAS, which is released under the terms of the GNU
% General Public Licence (GPL), version 3. For further details, see
% <https://www.gnu.org/licenses/>.
%
% This software is provided "as is", without warranty of any kind, express
% or implied, including, but not limited to the warranties of
% merchantability, fitness for a particular purpose and non-infringement.
%
% This software is intended for research only. Do not use for clinical
% purpose. Please note that this toolbox is under active development.
% Considerable changes may occur in future releases. For support please
% refer to:
% https://github.com/translationalneuromodeling/tapas/issues
%

% prediction error for current mean DCM parameter
mu_n = obj.posterior.mu_n(n,:);
data = obj.data(n);
obj.aux.epsilon{n} = obj.bold_gen(mu_n, data, obj.inputs(n), ...
    obj.options.hemo, obj.R, obj.L, obj.idx );

% gradient
obj.aux.G{n} = zeros(numel(obj.data(n).bold), obj.idx.P_c + obj.idx.P_h);
data.bold(:) = 0;
% central difference
for p = 1:obj.idx.P_c + obj.idx.P_h
    for s = -1:2:1
        % perturb parameters
        mu_n = obj.posterior.mu_n(n,:);
        mu_n(p) = mu_n(p) + s*obj.options.delta;
        % generate bold signal
        pred = obj.bold_gen(mu_n, data, obj.inputs(n), obj.options.hemo,...
            obj.R, obj.L, obj.idx );
        % calculate difference
        obj.aux.G{n}(:,p) = obj.aux.G{n}(:,p) - s*pred(:);
    end
end
% normalize by step size
obj.aux.G{n} = obj.aux.G{n}/obj.options.delta/2;

end

author: translationalneuromodeling

repo: tapas


text:
function point_num = sparse_grid_gl_size ( dim_num, level_max )

%*****************************************************************************80
%
%% SPARSE_GRID_GL_SIZE sizes a sparse grid of Gauss-Legendre points.
%
%  Discussion:
%
%    The grid is defined as the sum of the product rules whose LEVEL
%    satisfies:
%
%      LEVEL_MIN <= LEVEL <= LEVEL_MAX.
%
%    where LEVEL_MAX is user specified, and
%
%      LEVEL_MIN = max ( 0, LEVEL_MAX + 1 - DIM_NUM ).
%
%    The grids are only very weakly nested, since Gauss-Legendre rules
%    only have the origin in common.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 July 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Fabio Nobile, Raul Tempone, Clayton Webster,
%    A Sparse Grid Stochastic Collocation Method for Partial Differential
%    Equations with Random Input Data,
%    SIAM Journal on Numerical Analysis,
%    Volume 46, Number 5, 2008, pages 2309-2345.
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer LEVEL_MAX, the maximum value of LEVEL.
%
%    Output, integer POINT_NUM, the number of points in the grid.
%

%
%  Special case.
%
  if ( level_max == 0 )
    point_num = 1;
    return;
  end
%
%  The outer loop generates LEVELs from LEVEL_MIN to LEVEL_MAX.
%
  point_num = 0;

  level_min = max ( 0, level_max + 1 - dim_num );

  for level = level_min : level_max
%
%  The middle loop generates the next partition that adds up to LEVEL.
%
    level_1d = [];
    more = false;
    h = 0;
    t = 0;

    while ( 1 )

      [ level_1d, more, h, t ] = comp_next ( level, dim_num, level_1d, more, h, t );
%
%  Transform each 1D level to a corresponding 1D order.
%
      order_1d = level_to_order_open ( dim_num, level_1d );

      for dim = 1 : dim_num
%
%  If we can reduce the level in this dimension by 1 and
%  still not go below LEVEL_MIN.
%
        if ( level_min < level & 1 < order_1d(dim) )
          order_1d(dim) = order_1d(dim) - 1;
        end

      end

      point_num = point_num + prod ( order_1d(1:dim_num) );

      if ( ~more )
        break
      end

    end
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function y = tt_blockwise(arr)
%% TT-Toolbox 2.2.1, 2009-2015
%This is TT Toolbox, written by Ivan Oseledets et al.
%Institute of Numerical Mathematics, Moscow, Russia
%webpage: http://spring.inm.ras.ru/osel

%for questions on this particular code email Alexey Boyko, Skolkovo Institute of Science and Technology,
% a.boyko@skoltech.ru or alexey.boyko@skolkovotech.ru
%Moscow, Russia, 2016

%%
%INPUTS
% arr {2d cell array} of <tt_matrix>  variables
%OUTPUTS
% y <tt_matrix>
%
%this function assembles a block matrix in tt-format from a given 2d cell array
%of blocks in tt-format

    [n,m]=size(arr);

    block=zeros(n,m); block(1,1)=1;
    %since current version of TT-Toolbox doesn't support
    %init like M=0; M=M+<tt_matrix>, one should construct
    %a valid tt_matrix of consistent dimensions as an initial variable.

    y=tkron(arr{1,1},tt_matrix(block));

    for i=1:n
        for j=1:m
%             full(a{i,j})
            if (i~=1)||(j~=1)
                block=zeros(n,m);
                block(i,j)=1;
                y=y+tkron(arr{i,j},tt_matrix(block));
            end
        end
    end
return
end

author: oseledets

repo: TT-Toolbox


text:
function [cluster, total] = findcluster(onoff, spatdimneighbstructmat, varargin)

% FINDCLUSTER returns all connected clusters in a 3 dimensional matrix
% with a connectivity of 6.
%
% Use as
%   [cluster, num] = findcluster(onoff, spatdimneighbstructmat, minnbchan)
% or ar
%   [cluster, num] = findcluster(onoff, spatdimneighbstructmat, spatdimneighbselmat, minnbchan)
% where
%   onoff                   is a 3D boolean matrix with size N1xN2xN3
%   spatdimneighbstructmat  defines the neighbouring channels/combinations, see below
%   minnbchan               the minimum number of neighbouring channels/combinations
%   spatdimneighbselmat     is a special neighbourhood matrix that is used for selecting
%                           channels/combinations on the basis of the minnbchan criterium
%
% The neighbourhood structure for the first dimension is specified using
% spatdimneighbstructmat, which is a 2D (N1xN1) matrix. Each row and each column corresponds
% to a channel (combination) along the first dimension and along that row/column, elements
% with "1" define the neighbouring channel(s) (combinations). The first dimension of
% onoff should correspond to the channel(s) (combinations).
%
% See also SPM_BWLABEL (spm toolbox)

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

spatdimlength = size(onoff, 1);
nfreq = size(onoff, 2);
ntime = size(onoff, 3);

if length(size(spatdimneighbstructmat))~=2 || ~all(size(spatdimneighbstructmat)==spatdimlength)
  ft_error('invalid dimension of spatdimneighbstructmat');
end

minnbchan=0;
if length(varargin)==1
    minnbchan=varargin{1};
end
if length(varargin)==2
    spatdimneighbselmat=varargin{1};
    minnbchan=varargin{2};
end

if minnbchan>0
    % For every (time,frequency)-element, it is calculated how many significant
    % neighbours this channel has. If a significant channel has less than minnbchan
    % significant neighbours, then this channel is removed from onoff.

    if length(varargin)==1
        selectmat = single(spatdimneighbstructmat | spatdimneighbstructmat');
    end
    if length(varargin)==2
        selectmat = single(spatdimneighbselmat | spatdimneighbselmat');
    end
    nremoved=1;
    while nremoved>0
        nsigneighb=reshape(selectmat*reshape(single(onoff),[spatdimlength (nfreq*ntime)]),[spatdimlength nfreq ntime]);
        remove=(onoff.*nsigneighb)<minnbchan;
        nremoved=length(find(remove.*onoff));
        onoff(remove)=0;
    end
end

% for each channel (combination), find the connected time-frequency clusters
labelmat = zeros(size(onoff));
total = 0;
if nfreq*ntime>1
  for spatdimlev=1:spatdimlength
    [labelmat(spatdimlev, :, :), num] = spm_bwlabel(double(reshape(onoff(spatdimlev, :, :), nfreq, ntime)), 6); % the previous code contained a '4' for input
    labelmat(spatdimlev, :, :) = labelmat(spatdimlev, :, :) + (labelmat(spatdimlev, :, :)~=0)*total;
    total = total + num;
  end
else
  labelmat(onoff>0) = 1:sum(onoff(:));
  total = sum(onoff(:));
end
% combine the time and frequency dimension for simplicity
labelmat = reshape(labelmat, spatdimlength, nfreq*ntime);

% combine clusters that are connected in neighbouring channel(s)
% (combinations). Convert inputs to uint32 as that is required by the mex
% file (and the values will be positive integers anyway).
cluster = combineClusters(uint32(labelmat), logical(spatdimneighbstructmat), uint32(total));

% reshape the output to the original format of the data
cluster = reshape(cluster, spatdimlength, nfreq, ntime);

% update the total number
total = numel(unique(cluster(:)))-1; % the value of 0 does not count

author: spm


text:
function test_suite = test_createBasisTransform
%TEST_CREATEBASISTRANSFORM  Test case for the file createBasisTransform
%
%   Test case for the file createBasisTransform

%   Example
%   test_createBasisTransform
%
%   See also
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2011-10-13,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2011 INRA - Cepia Software Platform.

test_suite = functiontests(localfunctions);

function test_Translate(testCase) %#ok<*DEFNU>
% Basic test to check the function runs

p1 = [3 4];
p2 = [10 20];
basis1 = [p1 1 0 0 1];
basis2 = [p2 1 0 0 1];

dp = p1-p2;
exp = [eye(2) dp' ; 0 0 1];

trans = createBasisTransform(basis1, basis2);
testCase.assertEqual(exp, trans, 'AbsTol', .01);

function test_NonOrthogonalBasis(testCase)
% Check consistency of transfomating to Non Ortho bases, and back to global

src = [0 0  1 0  0 1];
tgt = [1 2  3 4  -5 6];

% transform a polygon to TGT basis
poly = [10 10;30 10; 30 20; 20 20;20 40; 10 40];
trans = createBasisTransform(src, tgt);
poly2 = transformPoint(poly, trans);

% transform back to original basis
trans2 = createBasisTransform(tgt, src);
poly3 = transformPoint(poly2, trans2);

testCase.assertEqual(poly, poly3, 'AbsTol', .01);

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function y = cvx_vexity( x )
sx = size( x );
if cvx_use_sparse( sx, 0, true )
    y = sparse( sx(1), sx(2) );
else
    y = zeros( sx );
end

% Copyright 2005-2016 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function ijk = points_to_idx(map, xyz)
% convert 3d points xyz (Nx3) into indices in occupancy grid
% @author: Yiren Lu
% @date: 02/13/2017
% @input:   map         map generated by load_map
%           xyz         Nx3 3d points in actual space
% @output:  ijk         Nx3 indices in map's occupancy grid
ijk = ones(size(xyz,1),3);
if size(xyz,1) > 0
ijk(:,1) = min(max(floor((xyz(:,1) - map.boundary(1))/map.xy_res + 1),1), map.nx);
ijk(:,2) = min(max(floor((xyz(:,2) - map.boundary(2))/map.xy_res + 1),1), map.ny);
ijk(:,3) = min(max(floor((xyz(:,3) - map.boundary(3))/map.z_res + 1),1), map.nz);
end
end

author: yrlu

repo: quadrotor

save_path: github-repos/MATLAB/yrlu-quadrotor


text:
function r = cos(a)
%COS          Hessian (elementwise) cosine
%

% written  04/04/04     S.M. Rump
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 09/28/08     S.M. Rump  check for rounding to nearest improved
% modified 08/26/12     S.M. Rump  global variables removed
%

  e = 1e-30;
  if 1+e==1-e                           % fast check for rounding to nearest
    rndold = 0;
  else
    rndold = getround;
    setround(0)
  end

  K = prod(size(a.x));
  if K==1                   % scalar hessian

    r.x = cos(a.x);
    msinax = -sin(a.x);
    r.dx = msinax * a.dx;
    r.hx = msinax * a.hx - reshape( ((0.5*r.x)*a.dx) * a.dx.',size(a.hx));

  else                      % matrix hessian

    N = getappdata(0,'INTLAB_HESSIAN_NUMVAR');
    N2 = N^2;

    r.x = cos(full(a.x));
    if issparse(a.hx)               % input sparse

      ax = -sin(a.x(:));
      sizeax = length(ax);
      [ia,ja,sa] = find(a.dx);
      % check for emptyness: cures Matlab bug
      % a=sparse([],[],[],2,1), [i,j,s]=find(a), s(i).*s(:)  yields error
      if isempty(ia)
        r.dx = sparse([],[],[],N,sizeax);
        r.hx = sparse([],[],[],N2,sizeax);
      else
        adx1 = (-0.5)*r.x(:);
        if isa(a.x,'intval')          % sparse intval
          rdx = times(ax(ja),sa(:),0);
          adx1 = times(adx1(ja),sa(:),0);
          if rdx.complex
            r.dx = intval( sparse(ia,ja,rdx.mid,N,sizeax) , sparse(ia,ja,rdx.rad,N,sizeax) , 'midrad' );
          else
            r.dx = intval( sparse(ia,ja,rdx.inf,N,sizeax) , sparse(ia,ja,rdx.sup,N,sizeax) , 'infsup' );
          end
          if adx1.complex
            adx1 = intval( sparse(ia,ja,adx1.mid,N,sizeax) , sparse(ia,ja,adx1.rad,N,sizeax) , 'midrad' );
          else
            adx1 = intval( sparse(ia,ja,adx1.inf,N,sizeax) , sparse(ia,ja,adx1.sup,N,sizeax) , 'infsup' );
          end
        else                          % sparse point
          r.dx = sparse(ia,ja,ax(ja).*sa(:),N,sizeax);
          adx1 = sparse(ia,ja,adx1(ja).*sa(:),N,sizeax);
        end
        r.hx = adx2rhx(N,sizeax,adx1,a.dx);
      end
      [ia,ja,sa] = find(a.hx);        % sparse point or intval
      % check for emptyness: cures Matlab bug
      % a=sparse([],[],[],2,1), [i,j,s]=find(a), s(i).*s(:)  yields error
      if ~isempty(ia)
        if isa(a.x,'intval')
          rhx = times(ax(ja),sa(:),0);
          if rhx.complex
            r.hx = r.hx + intval( sparse(ia,ja,rhx.mid,N2,sizeax) , sparse(ia,ja,rhx.rad,N2,sizeax) , 'midrad' );
          else
            r.hx = r.hx + intval( sparse(ia,ja,rhx.inf,N2,sizeax) , sparse(ia,ja,rhx.sup,N2,sizeax) , 'infsup' );
          end
        else
          r.hx = r.hx + sparse(ia,ja,ax(ja).*sa(:),N2,sizeax);
        end
      end

    else                            % input full

      r.x = cos(a.x);
      rx = r.x(:).';
      msinax = -sin(a.x(:).');
      msinax = msinax(ones(N*N,1),:);
      r.dx = a.dx .* msinax(1:N,:);
      adx = repmat(0.5*rx,N,1) .* a.dx;
      r.hx = a.hx .* msinax - adx(repmat(1:N,N,1),:) .* a.dx(repmat(1:N,1,N),:);

    end

  end

  r = class(r,'hessian');

  if rndold
    setround(rndold)
  end

repo: OpenMAS