text:
function [stack, stack_exposure ] = SortStack( stack, stack_exposure, sorting_type )
%
%       [lin_fun, pp] = MitsunagaNayarCRF(stack, stack_exposure, N, nSamples, sampling_strategy)
%
%       This function computes camera response function using Mitsunaga and
%       Nayar method.
%
%        Input:
%           -stack: a stack of LDR images. If the stack is a single or
%           double values are assumed to be in [0,1]
%           -stack_exposure: an array containg the exposure time of each
%           image. Time is expressed in second (s)
%           -sorting_type: 'ascend' or 'descend'
%
%        Output:
%           -stack: sorted stack
%           -stack_exposure: sorted stack_exposure
%
%     Copyright (C) 2016  Francesco Banterle
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
%

if(~exist(sorting_type, 'var'))
    sorting_type =  'ascend';
end

[stack_exposure_sorted, ind] = sort(stack_exposure, sorting_type);

if(sum(abs(stack_exposure_sorted - stack_exposure)) > 0.0)
    stack_sorted = zeros(size(stack));
    for i=1:length(stack_exposure)
        stack_sorted(:,:,:,i) = stack(:,:,:,ind(i));
    end

    stack = stack_sorted;
    stack_exposure = stack_exposure_sorted;

    clear('stack_sorted');
    clear('stack_exposure_sorted');
end

end

author: banterle


text:
function h=ref_pfilt_1(f,g,a)
%REF_PFILT_1  Reference PFILT implementation by FFT
%
%   This is the old reference pfilt from before the struct filters where
%   introduced.

[L W]=size(f);

g=fir2long(g,L);

% Force FFT along dimension 1, since we have permuted the dimensions
% manually
if isreal(f) && isreal(g)
  h=ifftreal(fftreal(f,L,1).*repmat(fftreal(g,L,1),1,W),L,1);
else
  h=ifft(fft(f,L,1).*repmat(fft(g,L,1),1,W),L,1);
end;

h=h(1:a:end,:);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function determ = rutis5_determinant ( )

%*****************************************************************************80
%
%% RUTIS5_DETERMINANT returns the determinant of the RUTIS5 matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real DETERM, the determinant.
%
  determ = 1.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% POP_PROP - plot the properties of a channel or of an independent
%              component.
% Usage:
%   >> pop_prop( EEG);           % pops up a query window
%   >> pop_prop( EEG, typecomp); % pops up a query window
%   >> pop_prop( EEG, typecomp, chanorcomp, winhandle,spectopo_options);
%
% Inputs:
%   EEG        - EEGLAB dataset structure (see EEGGLOBAL)
%
% Optional inputs:
%   typecomp   - [0|1] 1 -> display channel properties
%                0 -> component properties {default: 1 = channel}
%   chanorcomp - channel or component number[s] to display {default: 1}
%
%   winhandle  - if this parameter is present or non-NaN, buttons
%                allowing the rejection of the component are drawn.
%                If non-zero, this parameter is used to back-propagate
%                the color of the rejection button.
%   spectopo_options - [cell array] optional cell array of options for
%                the SPECTOPO function.
%                For example { 'freqrange' [2 50] }
%
% Note: Type "set(gcf, ''renderer'', ''painter'')" before saving the figure
% in postscript (epsc) or jpg format
%
% Author: Arnaud Delorme, CNL / Salk Institute, 2001
%
% See also: POP_RUNICA, EEGLAB

% Copyright (C) 2001 Arnaud Delorme, Salk Institute, arno@salk.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

% hidden parameter winhandle

% 01-25-02 reformated help & license -ad
% 02-17-02 removed event index option -ad
% 03-17-02 debugging -ad & sm
% 03-18-02 text settings -ad & sm
% 03-18-02 added title -ad & sm

function com = pop_prop(EEG, typecomp, chanorcomp, winhandle, spec_opt)

com = '';
if nargin < 1
	help pop_prop;
	return;
end
if nargin < 5
	spec_opt = {};
end
if nargin == 1
	typecomp = 1;    % defaults
        chanorcomp = 1;
end
if typecomp == 0 && isempty(EEG.icaweights)
   error('No ICA weights recorded for this dataset -- first run ICA on it');
end
chanoristr = '';
if nargin == 2
    commandchan = 'tmpchanlocs = EEG(1).chanlocs; [tmp tmpval] = pop_chansel({tmpchanlocs.labels}, ''withindex'', ''on'', ''selectionmode'', ''single''); set(findobj(gcbf, ''tag'', ''chan''), ''string'',tmpval); clear tmp tmpchanlocs tmpval';

    uitext = { { 'style' 'text' 'string' fastif(typecomp,'Channel index(ices) to plot:','Component index(ices) to plot:') } ...
               { 'style' 'edit' 'string' '1', 'tag', 'chan' } ...
               { 'style' 'pushbutton' 'string'  '...', 'enable' fastif(~isempty(EEG(1).chanlocs) && typecomp, 'on', 'off') 'callback' commandchan } ...
               { 'style' 'text' 'string' 'Spectral options (see spectopo() help):' } ...
               { 'style' 'edit' 'string' '''freqrange'', [2 50]' } {} };
    uigeom = { [2 1 0.5 ] [2 1 0.5] };
    result = inputgui('geometry', uigeom, 'uilist', uitext, 'helpcom', 'pophelp(''pop_prop'');', ...
							 'title', fastif( typecomp, 'Channel properties - pop_prop()', 'Component properties - pop_prop()'));
	if size( result, 1 ) == 0 return; end

    chanoristr = result{1};
    chanorcomp   = eeg_decodechan(EEG.chanlocs, result{1} );
    spec_opt     = eval( [ '{' result{2} '}' ] );
end

if ischar(chanorcomp)
    chanoristr = chanorcomp;
    chanorcomp = eeg_decodechan(EEG.chanlocs, chanorcomp );
elseif isempty(chanoristr)
    chanoristr = int2str(chanorcomp);
end

% plotting several component properties
% -------------------------------------
if length(chanorcomp) > 1
    for index = chanorcomp
        pop_prop(EEG, typecomp, index, 0, spec_opt);  % call recursively for each chanorcomp
    end
	com = sprintf('pop_prop( EEG, %d, [%s], NaN, %s);',...
                  typecomp, int2str(chanorcomp), vararg2str( { spec_opt } ));
    return;
end

% should test for > number of components ??? -sm.
% yes (max num components is not necessarily same as nbchan). -jri
if typecomp == 1
  maxChanorcomp = EEG.nbchan;
else
  maxChanorcomp = size(EEG.icaweights, 1);
end

if chanorcomp < 1 || chanorcomp > maxChanorcomp
   error('Component index out of range');
end

% assumed input is chanorcomp
% -------------------------
try
    icadefs;
catch
	BACKCOLOR = [0.8 0.8 0.8];
	GUIBUTTONCOLOR   = [0.8 0.8 0.8];
end
basename = [fastif(typecomp,'Channel ', 'Component ') chanoristr ];

fhandle = figure('name', ['pop_prop() - ' basename ' properties'], 'color', BACKCOLOR, 'numbertitle', 'off', 'visible', 'off');
pos     = get(fhandle,'Position');
set(fhandle,'Position', [pos(1) pos(2)-500+pos(4) 500 500], 'visible', 'on');
hh = axes('parent',fhandle);
pos      = get(hh,'position'); % plot relative to current axes
q = [pos(1) pos(2) 0 0];
s = [pos(3) pos(4) pos(3) pos(4)]./100;
axis(hh,'off');

% plotting topoplot
% -----------------
h = axes('parent',fhandle,'Units','Normalized', 'Position',[-10 60 40 42].*s+q);

%topoplot( EEG.icawinv(:,chanorcomp), EEG.chanlocs); axis square;

if isfield(EEG.chanlocs, 'theta')
    if typecomp == 1 % plot single channel locations
        topoplot( chanorcomp, EEG.chanlocs, 'chaninfo', EEG.chaninfo, ...
                 'electrodes','off', 'style', 'blank', 'emarkersize1chan', 12); axis square;
    else             % plot component map
        topoplot( EEG.icawinv(:,chanorcomp), EEG.chanlocs, 'chaninfo', EEG.chaninfo, ...
                 'shading', 'interp', 'numcontour', 3); axis square;
    end
else
    axis(h,'off');
end
basename = [fastif(typecomp,'Channel ', 'IC') chanoristr ];
% title([ basename fastif(typecomp, ' location', ' map')], 'fontsize', 14);
title(basename, 'fontsize', 14);

% plotting erpimage
% -----------------
hhh = axes('Parent', fhandle,'Units','Normalized', 'Position',[45 62 48 38].*s+q);
eeglab_options;
if EEG.trials > 1
    % put title at top of erpimage
    axis(hhh,'off');
    hh = axes('Parent', fhandle,'Units','Normalized', 'Position',[45 62 48 38].*s+q);
    EEG.times = linspace(EEG.xmin, EEG.xmax, EEG.pnts);
    if EEG.trials < 6
      ei_smooth = 1;
    else
      ei_smooth = 3;
    end
    if typecomp == 1 % plot channel
         offset = nan_mean(EEG.data(chanorcomp,:));
         erp=nan_mean(squeeze(EEG.data(chanorcomp,:,:))')-offset;
         erp_limits=get_era_limits(erp);
         erpimage( EEG.data(chanorcomp,:)-offset, ones(1,EEG.trials)*10000, EEG.times*1000, ...
                       '', ei_smooth, 1, 'caxis', 2/3, 'cbar','erp','erp_vltg_ticks',erp_limits);
    else % plot component
         icaacttmp  = eeg_getdatact(EEG, 'component', chanorcomp);
         offset     = nan_mean(icaacttmp(:));
         era        = nan_mean(squeeze(icaacttmp)')-offset;
         era_limits = get_era_limits(era);
         erpimage( icaacttmp-offset, ones(1,EEG.trials)*10000, EEG.times*1000, ...
                       '', ei_smooth, 1, 'caxis', 2/3, 'cbar','erp', 'yerplabel', '','erp_vltg_ticks',era_limits);
    end
    axes(hhh);
    title(sprintf('%s activity \\fontsize{10}(global offset %3.3f)', basename, offset), 'fontsize', 14);
else
    % put title at top of erpimage
    EI_TITLE = 'Continous data';
    axis(hhh,'off');
    hh = axes('Parent', fhandle,'Units','Normalized', 'Position',[45 62 48 38].*s+q);
    ERPIMAGELINES = 200; % show 200-line erpimage
    while size(EEG.data,2) < ERPIMAGELINES*EEG.srate
       ERPIMAGELINES = 0.9 * ERPIMAGELINES;
    end
    ERPIMAGELINES = round(ERPIMAGELINES);
    if ERPIMAGELINES > 2   % give up if data too small
        if ERPIMAGELINES < 10
            ei_smooth = 1;
        else
            ei_smooth = 3;
        end
      erpimageframes = floor(size(EEG.data,2)/ERPIMAGELINES);
      erpimageframestot = erpimageframes*ERPIMAGELINES;
      eegtimes = linspace(0, erpimageframes-1, length(erpimageframes)); % 05/27/2014 Ramon: length(erpimageframes) by EEG.srate/1000  in eegtimes = linspace(0, erpimageframes-1, EEG.srate/1000);
      if typecomp == 1 % plot channel
           offset = nan_mean(EEG.data(chanorcomp,:));
           % Note: we don't need to worry about ERP limits, since ERPs
           % aren't visualized for continuous data
           erpimage( reshape(EEG.data(chanorcomp,1:erpimageframestot),erpimageframes,ERPIMAGELINES)-offset, ones(1,ERPIMAGELINES)*10000, eegtimes , ...
                         EI_TITLE, ei_smooth, 1, 'caxis', 2/3, 'cbar');
      else % plot component
         icaacttmp = eeg_getdatact(EEG, 'component', chanorcomp);
         offset = nan_mean(icaacttmp(:));
         erpimage(reshape(icaacttmp(:,1:erpimageframestot),erpimageframes,ERPIMAGELINES)-offset,ones(1,ERPIMAGELINES)*10000, eegtimes , ...
                    EI_TITLE, ei_smooth, 1, 'caxis', 2/3, 'cbar','yerplabel', '');
      end
    else
            axis(hh,'off');
            text(0.1, 0.3, [ 'No erpimage plotted' 10 'for small continuous data']);
    end
    axes(hhh);
end

% plotting spectrum
% -----------------
if ~exist('winhandle')
    winhandle = NaN;
end
if ishandle(winhandle)
	h = axes('Parent', fhandle,'units','normalized', 'position',[5 10 95 35].*s+q);
else
	h = axes('Parent', fhandle,'units','normalized', 'position',[5 0 95 40].*s+q);
end
%h = axes('units','normalized', 'position',[45 5 60 40].*s+q);
try
	eeglab_options;
	if typecomp == 1
		[spectra, freqs] = spectopo( EEG.data(chanorcomp,:), EEG.pnts, EEG.srate, spec_opt{:} );
	else
		if option_computeica
			[spectra, freqs] = spectopo( EEG.icaact(chanorcomp,:), EEG.pnts, EEG.srate, 'mapnorm', EEG.icawinv(:,chanorcomp), spec_opt{:} );
        else
    		icaacttmp = (EEG.icaweights(chanorcomp,:)*EEG.icasphere)*reshape(EEG.data(EEG.icachansind,:,:), length(EEG.icachansind), EEG.trials*EEG.pnts);
			[spectra, freqs] = spectopo( icaacttmp, EEG.pnts, EEG.srate, 'mapnorm', EEG.icawinv(:,chanorcomp), spec_opt{:} );
		end
	end
    % set up new limits
    % -----------------
    %freqslim = 50;
	%set(gca, 'xlim', [0 min(freqslim, EEG.srate/2)]);
    %spectra = spectra(find(freqs <= freqslim));
	%set(gca, 'ylim', [min(spectra) max(spectra)]);

	%tmpy = get(gca, 'ylim');
    %set(gca, 'ylim', [max(tmpy(1),-1) tmpy(2)]);
	set( get(h, 'ylabel'), 'string', 'Power 10*log_{10}(\muV^{2}/Hz)', 'fontsize', 14);
	set( get(h, 'xlabel'), 'string', 'Frequency (Hz)', 'fontsize', 14);
	title('Activity power spectrum', 'fontsize', 14);
catch err
	axis off;
    text(0.1, 0.3, [ 'Error: no spectrum plotted' 10 err.message 10]);
%     lasterror
% 	text(0.1, 0.3, [ 'Error: no spectrum plotted' 10 ' make sure you have the ' 10 'signal processing toolbox']);
end

% display buttons
% ---------------
if ishandle(winhandle)
	COLREJ = '[1 0.6 0.6]';
	COLACC = '[0.75 1 0.75]';
	% CANCEL button
	% -------------
	h  = uicontrol(fhandle, 'Style', 'pushbutton', 'backgroundcolor', GUIBUTTONCOLOR, 'string', 'Cancel', 'Units','Normalized','Position',[-10 -10 15 6].*s+q, 'callback', 'close(gcf);');

	% VALUE button
	% -------------
	hval  = uicontrol(fhandle, 'Style', 'pushbutton', 'backgroundcolor', GUIBUTTONCOLOR, 'string', 'Values', 'Units','Normalized', 'Position', [15 -10 15 6].*s+q);

	% REJECT button
	% -------------
    if ~isempty(EEG.reject.gcompreject)
    	status = EEG.reject.gcompreject(chanorcomp);
    else
        status = 0;
    end
	hr = uicontrol(fhandle, 'Style', 'pushbutton', 'backgroundcolor', eval(fastif(status,COLREJ,COLACC)), ...
				'string', fastif(status, 'REJECT', 'ACCEPT'), 'Units','Normalized', 'Position', [40 -10 15 6].*s+q, 'userdata', status, 'tag', 'rejstatus');
	command = [ 'set(gcbo, ''userdata'', ~get(gcbo, ''userdata''));' ...
				'if get(gcbo, ''userdata''),' ...
				'     set( gcbo, ''backgroundcolor'',' COLREJ ', ''string'', ''REJECT'');' ...
				'else ' ...
				'     set( gcbo, ''backgroundcolor'',' COLACC ', ''string'', ''ACCEPT'');' ...
				'end;' ];
	set(hr, 'callback', command);

	% HELP button
	% -------------
	h  = uicontrol(fhandle, 'Style', 'pushbutton', 'backgroundcolor', GUIBUTTONCOLOR, 'string', 'HELP', 'Units','Normalized', 'Position', [65 -10 15 6].*s+q, 'callback', 'pophelp(''pop_prop'');');

	% OK button
	% ---------
 	command = [ 'global EEG;' ...
 				'tmpstatus = get( findobj(''parent'', gcbf, ''tag'', ''rejstatus''), ''userdata'');' ...
 				'EEG.reject.gcompreject(' num2str(chanorcomp) ') = tmpstatus;' ];
	if winhandle ~= 0
        if VERS < 8.04
            command = [ command ...
                sprintf('if tmpstatus set(%3.15f, ''backgroundcolor'', %s); else set(%3.15f, ''backgroundcolor'', %s); end;', ...
                winhandle, COLREJ, winhandle, COLACC)];
        elseif VERS >= 8.04
            command = [ command ...
                sprintf('if tmpstatus set(findobj(''Tag'', ''%s''), ''backgroundcolor'', %s); else set(findobj(''Tag'',''%s''), ''backgroundcolor'', %s); end;', ...
                winhandle.Tag, COLREJ, winhandle.Tag, COLACC)];
        end
    end
	command = [ command 'close(gcf); clear tmpstatus' ];
	h  = uicontrol(fhandle, 'Style', 'pushbutton', 'string', 'OK', 'backgroundcolor', GUIBUTTONCOLOR, 'Units','Normalized', 'Position',[90 -10 15 6].*s+q, 'callback', command);

	% draw the figure for statistical values
	% --------------------------------------
	index = num2str( chanorcomp );
	command = [ ...
		'figure(''MenuBar'', ''none'', ''name'', ''Statistics of the component'', ''numbertitle'', ''off'');' ...
		'' ...
		'pos = get(gcf,''Position'');' ...
		'set(gcf,''Position'', [pos(1) pos(2) 340 340]);' ...
		'pos = get(gca,''position'');' ...
		'q = [pos(1) pos(2) 0 0];' ...
		's = [pos(3) pos(4) pos(3) pos(4)]./100;' ...
		'axis off;' ...
		''  ...
		'txt1 = sprintf(''(\n' ...
						'Entropy of component activity\t\t%2.2f\n' ...
					    '> Rejection threshold \t\t%2.2f\n\n' ...
					    ' AND                 \t\t\t----\n\n' ...
					    'Kurtosis of component activity\t\t%2.2f\n' ...
					    '> Rejection threshold \t\t%2.2f\n\n' ...
					    ') OR                  \t\t\t----\n\n' ...
					    'Kurtosis distribution \t\t\t%2.2f\n' ...
					    '> Rejection threshold\t\t\t%2.2f\n\n' ...
					    '\n' ...
					    'Current thresholds sujest to %s the component\n\n' ...
					    '(after manually accepting/rejecting the component, you may recalibrate thresholds for future automatic rejection on other datasets)'',' ...
						'EEG.stats.compenta(' index '), EEG.reject.threshentropy, EEG.stats.compkurta(' index '), ' ...
						'EEG.reject.threshkurtact, EEG.stats.compkurtdist(' index '), EEG.reject.threshkurtdist, fastif(EEG.reject.gcompreject(' index '), ''REJECT'', ''ACCEPT''));' ...
		'' ...
		'uicontrol(gcf, ''Units'',''Normalized'', ''Position'',[-11 4 117 100].*s+q, ''Style'', ''frame'' );' ...
		'uicontrol(gcf, ''Units'',''Normalized'', ''Position'',[-5 5 100 95].*s+q, ''String'', txt1, ''Style'',''text'', ''HorizontalAlignment'', ''left'' );' ...
		'h = uicontrol(gcf, ''Style'', ''pushbutton'', ''string'', ''Close'', ''Units'',''Normalized'', ''Position'', [35 -10 25 10].*s+q, ''callback'', ''close(gcf);'');' ...
		'clear txt1 q s h pos;' ];
	set( hval, 'callback', command);
	if isempty( EEG.stats.compenta )
		set(hval, 'enable', 'off');
	end

	com = sprintf('pop_prop( EEG, %d, %d, 0, %s);', typecomp, chanorcomp, vararg2str( { spec_opt } ) );
else
	com = sprintf('pop_prop( EEG, %d, %d, NaN, %s);', typecomp, chanorcomp, vararg2str( { spec_opt } ) );
end
disp('Note: Type "set(gcf, ''renderer'', ''painter'')" before saving the figure in postscript (epsc) or jpg format');

return;

function out = nan_mean(in)

    nans = find(isnan(in));
    in(nans) = 0;
    sums = sum(in);
    nonnans = ones(size(in));
    nonnans(nans) = 0;
    nonnans = sum(nonnans);
    nononnans = find(nonnans==0);
    nonnans(nononnans) = 1;
    out = sum(in)./nonnans;
    out(nononnans) = NaN;

function era_limits=get_era_limits(era)
%function era_limits=get_era_limits(era)
%
% Returns the minimum and maximum value of an event-related
% activation/potential waveform (after rounding according to the order of
% magnitude of the ERA/ERP)
%
% Inputs:
% era - [vector] Event related activation or potential
%
% Output:
% era_limits - [min max] minimum and maximum value of an event-related
% activation/potential waveform (after rounding according to the order of
% magnitude of the ERA/ERP)

mn=min(era);
mx=max(era);
mn=orderofmag(mn)*round(mn/orderofmag(mn));
mx=orderofmag(mx)*round(mx/orderofmag(mx));
era_limits=[mn mx];

function ord=orderofmag(val)
%function ord=orderofmag(val)
%
% Returns the order of magnitude of the value of 'val' in multiples of 10
% (e.g., 10^-1, 10^0, 10^1, 10^2, etc ...)
% used for computing erpimage trial axis tick labels as an alternative for
% plotting sorting variable

val=abs(val);
if val>=1
    ord=1;
    val=floor(val/10);
    while val>=1
        ord=ord*10;
        val=floor(val/10);
    end
    return;
else
    ord=1/10;
    val=val*10;
    while val<1
        ord=ord/10;
        val=val*10;
    end
    return;
end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function mv = mv_glmPanel(mv)
%
% Attach a uipanel to a MultiVoxelUI with controls for visualizing
% GLMs. This initializes the panel to 'off', not visible. The panel is
% made visible when the visualize GLM option is selected.
%
%
% ras, 10/06
if notDefined('mv'), mv = get(gcf, 'UserData'); end
if ishandle(mv), mv = get(mv, 'UserData'); end

mv.ui.glmPanel = mrvPanel('right',  .2);
set(mv.ui.glmPanel,  'BackgroundColor',  'w');

% voxel slider
nVoxels = size(mv.tSeries,  2);
mv.ui.glmVoxel = mrvSlider([.1 .8 .8 .1], 'Voxel', 'Parent', mv.ui.glmPanel,  ...
                           'Range',  [1 nVoxels], 'IntFlag', 1, ...
                           'MaxLabelFlag', 1, ...
                           'Callback', 'mv_visualizeGlm(gcf,  round(val));');

% buttons to grab the current voxel time course, multi voxel
cb = ['tmp = get(gcf, ''UserData''); ' ...
      'iTmp = get(tmp.ui.glmVoxel.sliderHandle, ''Value''); ' ...
      'mv_selectSubset(tmp, iTmp, ''voxels'', 2); ' ...
      'clear tmp iTmp'];
uicontrol('Parent', mv.ui.glmPanel, 'Style', 'pushbutton', ...
          'Units', 'norm', 'Position', [.1 .6 .8 .06], ...
          'String', 'Cur Voxel Timecourse', 'FontSize', 9, ...
          'BackgroundColor', 'w', 'HorizontalAlignment', 'center', ...
          'Callback', cb);

cb = ['tmp = get(gcf, ''UserData''); ' ...
      'iTmp = get(tmp.ui.glmVoxel.sliderHandle, ''Value''); ' ...
      'mv_selectSubset(tmp, iTmp, ''voxels'', 1); ' ...
      'clear tmp iTmp'];
uicontrol('Parent', mv.ui.glmPanel, 'Style', 'pushbutton', ...
          'Units', 'norm', 'Position', [.1 .54 .8 .06], ...
          'String', 'Cur Voxel MVUI', 'FontSize', 9, ...
          'BackgroundColor', 'w', 'HorizontalAlignment', 'center', ...
          'Callback', cb);

% text describing voxel location
mv.ui.glmVoxCoords = uicontrol('Parent', mv.ui.glmPanel, 'Style', 'text', ...
          'Units', 'Normalized', 'Position', [.1 .4 .8 .06], ...
          'String', '', 'BackgroundColor', 'w');

% toggle buttons for time course traces
cb = ['tmp = {''off'' ''on''}; val = get(gcbo, ''Value'')+1; ' ...
      'set(get(gcbo,''UserData''), ''Visible'', tmp{val}); ' ...
      'clear tmp val '];

h1 = uicontrol('Parent', mv.ui.glmPanel, 'Style', 'checkbox', ...
               'Units', 'normalized', ...
               'Position', [.1 .2 .8 .05], 'String', 'Time Course', ...
               'BackgroundColor', 'w', 'Value', 1, 'Callback', cb);
h2 = uicontrol('Parent', mv.ui.glmPanel, 'Style', 'checkbox', ...
               'Units', 'normalized', ...
               'Position', [.1 .15 .8 .05], 'String', 'Predictors', ...
               'BackgroundColor', 'w', 'Value', 0, 'Callback', cb);
h3 = uicontrol('Parent', mv.ui.glmPanel, 'Style', 'checkbox', ...
               'Units', 'normalized', ...
               'Position', [.1 .1 .8 .05], 'String', 'Best Fit', ...
               'BackgroundColor', 'w', 'Value', 1, 'Callback', cb);
h4 = uicontrol('Parent', mv.ui.glmPanel, 'Style', 'checkbox', ...
               'Units', 'normalized', ...
               'Position', [.1 .05 .8 .05], 'String', 'Residual', ...
               'BackgroundColor', 'w', 'Value', 0, 'Callback', cb);

mv.ui.glmToggles = [h1 h2 h3 h4];

set( get(mv.ui.glmPanel, 'Children'), 'Units', 'normalized' );

mrvPanelToggle(mv.ui.glmPanel,  'off');

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function hdr = empty_hdr
% Create an empty NIFTI-1 header
% FORMAT hdr = empty_hdr
% _______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

%
% $Id: empty_hdr.m 1143 2008-02-07 19:33:33Z spm $

org = niftistruc;
for i=1:length(org),
    hdr.(org(i).label) = org(i).def;
end;

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function clean()

clear all;

setenv('ACADOS_INSTALL_DIR', '');

ocl_path  = fileparts(which('ocl'));

% clean export
if exist(fullfile(ocl_path,'Workspace','export'), 'dir')
  addpath(fullfile(ocl_path,'Workspace','export'));
  rmpath(fullfile(ocl_path,'Workspace','export'));

  success = rmdir(fullfile(ocl_path,'Workspace','export'), 's');
  if ~success
    error('Could not remove export folder. Restart Matlab and run clean again.');
  end
end

% clean acados lib
if exist(fullfile(ocl_path,'Lib','acados'), 'dir')
  addpath(fullfile(ocl_path,'Lib','acados'));
  rmpath(fullfile(ocl_path,'Lib','acados'));

  success = rmdir(fullfile(ocl_path,'Lib','acados'), 's');
  if ~success
    error('Could not remove acados binaries folder. Restart Matlab and run clean again.');
  end
end

author: OpenOCL

repo: OpenOCL

save_path: github-repos/MATLAB/OpenOCL-OpenOCL


text:
% Project the dictionary D to satisfy the constraint that the norm of each
% column should be smaller or euqal to 1

function D = projectionDic(D)

% for i = 1: size(D,2)
%     norm_col_squared = sum(D(:,i).^2);
%     if norm_col_squared > 1
%         D(:,i) = D(:,i)/sqrt(norm_col_squared);
%     end
% end

% vectorize
norm_col_squared = sum(D.^2,1);
if norm_col_squared > 1
    D(:,norm_col_squared > 1) = D(:,norm_col_squared > 1)./repmat(sqrt(norm_col_squared(1,norm_col_squared > 1)),size(D,1),1);
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
m = 100;
n = 1000;
k = 12;
A = spx.dict.simple.gaussian_dict(m, n);
gen = spx.data.synthetic.SparseSignalGenerator(n, k);
% create a sparse vector
x =  gen.biGaussian();
b = A*x;
A  = double(A);
result = spx.fast.omp(A, b, k, 1e-12);
cmpare = spx.commons.SparseSignalsComparison(x, result, k);
cmpare.summarize();

author: indigits

repo: sparse-plex


text:
function [Estimation_X0] = notREKF_update(Estimation_X0, CameraMeasurementThis, Sigma_OB)

NumberOfLandmarksObInThisStep = size(CameraMeasurementThis,1)/3;

% initialise the IndexOfFeature if possible
if size(Estimation_X0.landmarks,2) > 0
    IndexOfFeature = Estimation_X0.landmarks(4,:)';
else
    IndexOfFeature = [];
end

IndexObservedAlreadyThis = [];
IndexObservedNew = [];
for i = 1:NumberOfLandmarksObInThisStep
    % check whether the feature is observed before or not
    M = find( IndexOfFeature== CameraMeasurementThis(3*i,2) );
    if isempty(M)
        IndexObservedNew = [IndexObservedNew;CameraMeasurementThis(3*i,2) ];
    else
        IndexObservedAlreadyThis = [IndexObservedAlreadyThis;CameraMeasurementThis(3*i,2)];
    end
end

Estimation_X0.IndexObservedNew=IndexObservedNew;
Estimation_X0.IndexObservedAlreadyThis=IndexObservedAlreadyThis;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% IndexObservedNew=[78; 97; 18] indicates that the
% robot firstly observes landmarks 78 97 18 in this step
% IndexObservedAlreadyThis=[19; 20; 53] indicates
% that the robot observes again landmarks 19 20 53 in this step
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

orientation = Estimation_X0.orientation;
position    = Estimation_X0.position;
cov         = Estimation_X0.cov;

NumberOfFeature = size( IndexOfFeature,1);
NumberOfOldFeatureInThisStep = size(IndexObservedAlreadyThis,1);
NumberOfNewFeatureInThisStep = size(IndexObservedNew,1);

% update state and covariance
if ~isempty(IndexObservedAlreadyThis)
    Z = zeros( NumberOfOldFeatureInThisStep*3 , 1);
    Y = zeros( NumberOfOldFeatureInThisStep*3 , 1);
    H = zeros(3*NumberOfOldFeatureInThisStep, 6+3*NumberOfFeature);

    temp = repmat({eye(3)}, NumberOfOldFeatureInThisStep,1 );
    R = blkdiag(temp{:});

    % update old features
    for i = 1:NumberOfOldFeatureInThisStep
        ind = find(IndexOfFeature == IndexObservedAlreadyThis(i));
        fi  = Estimation_X0.landmarks(1:3,ind);
        Y(3*i-2:3*i,1) = ObservationModel( orientation, position, fi );

        ind2 = find(CameraMeasurementThis(:,2) == IndexObservedAlreadyThis(i));
        Z(3*i-2:3*i,1 ) = CameraMeasurementThis(ind2,1);

        H(3*i-2:3*i, 4:6) = orientation';
        H(3*i-2:3*i, 6+3*ind-2:6+3*ind) = -orientation';
        R(3*i-2:3*i,3*i-2:3*i) = diag(CameraMeasurementThis(ind2,1).^2)*Sigma_OB^2;
    end

    % question @RomaTeng, different computaton scheme
    z = Z-Y;
    S = H*cov*H'+R;
    K = cov*H'*inv(S);
    s = K*z;

    Estimation_X0 = special_add_right_not(Estimation_X0,-s);
    cov = ( eye(6+3*NumberOfFeature) -K*H )*cov;
    Estimation_X0.cov = cov;

    % @todo @RomaTeng, right Jacobian % No need for consistency
    % Estimation_X0.cov=JJJr(s)*cov*(JJJr(s))';
end

% update state vector and covariance by considering
% new feature into state and covariance
if ~isempty(IndexObservedNew)
    % copy previous covariance
    temp    = repmat({eye(3)}, NumberOfNewFeatureInThisStep, 1 );
    tempKK  = blkdiag(temp{:});
    Sigma   = blkdiag(Estimation_X0.cov,tempKK);
    KK      = eye(6+3*(NumberOfFeature+NumberOfNewFeatureInThisStep));

    % add new features
    for i = 1:NumberOfNewFeatureInThisStep
        indNewf = IndexObservedNew(i);
        Estimation_X0.landmarks(4,NumberOfFeature+i) = indNewf;
        m2 = find( CameraMeasurementThis(:,2) == indNewf );
        nf = CameraMeasurementThis( m2, 1 );

        Estimation_X0.landmarks(1:3,NumberOfFeature+i) = Estimation_X0.orientation*nf+Estimation_X0.position;
        KK( 6+3*NumberOfFeature+3*i-2:6+3*NumberOfFeature+3*i,4:6 ) = eye(3);
        KK( 6+3*NumberOfFeature+3*i-2:6+3*NumberOfFeature+3*i,6+3*NumberOfFeature+3*i-2:6+3*NumberOfFeature+3*i ) = Estimation_X0.orientation;

        tempKK(3*i-2:3*i,3*i-2:3*i)=diag(nf.^2)*Sigma_OB^2;
    end
        Sigma   = blkdiag(Estimation_X0.cov,tempKK);
    Estimation_X0.cov = KK*Sigma*KK';
end

end

author: RomaTeng

repo: EKF-SLAM-on-Manifold


text:
%
% Written by M. Harper Langston - 5/10/00
% harper@cims.nyu.edu
%
% For this sparse LU solver of Ax = b, A has band s.
% For example, the following matrix has band s = 1
%
%   x x 0 0 0 0
%   x x x 0 0 0
%   0 x x x 0 0
%   0 0 x x x 0
%   0 0 0 x x x
%   0 0 0 0 x x
%
%
function [x] = Band_Solve(A,b)
% Here, s is the band number
[m,n] = size(A);
if m~=n
   error('Matrix not sqaure') % Want square matrix for simplicity
end
% Find the band of A
for l = 1:n
   if A(m,l)~=0
      s = n-l;
      break;
   end
end
% Do the sparse LU decomposition for a matrix of band s
for j = 1:m
   if j >= m-s
      L(j:m,j) = A(j:m,j)./A(j,j);
		U(j,j:m) = A(j,j:m);
   	A(j:m,j:m) = A(j:m,j:m) - L(j:m,j)*U(j,j:m);
   else
   	L(j:j+s,j) = A(j:j+s,j)./A(j,j);
   	U(j,j:j+s) = A(j,j:j+s);
   	A(j:j+s,j:j+s) = A(j:j+s,j:j+s) - L(j:j+s,j)*U(j,j:j+s);
   end
end
L = sparse(L);
U = sparse(U);
% Call backsolve routine, which I wrote to pay heed to sparsity.
[y] = Back_Solve(L,b);
[x] = Back_Solve(U,y);
%
% Written by M. Harper Langston - 5/10/00
%

author: Sable

repo: mcbench-benchmarks


text:
function vis_wp(img, opt, heatmap, center, scale, cad, dict)

    img_crop = cropImage(img,center,scale);
    S = bsxfun(@plus,opt.R*opt.S,[opt.T;0]);
    model = fullShape(S,cad,dict.kpt_id);
    mesh2d = model.vertices(:,1:2)'*200/size(heatmap,2);

    % visualization
    nplot = 4;
    h = figure('position',[100,100,nplot*300,300]);
    % cropped image
    subplot('position',[0 0 1/nplot 1]);
    imshow(img_crop); hold on;
    % heatmap
    subplot('position',[1/nplot 0 1/nplot 1]);
    response = sum(heatmap,3);
    max_value = max(max(response));
    mapIm = imresize(mat2im(response, jet(100), [0 max_value]),[200,200],'nearest');
    imToShow = mapIm*0.5 + (single(img_crop)/255)*0.5;
    imagesc(imToShow); axis equal off
    % project cad model on image
    subplot('position',[2/nplot 0 1/nplot 1]);
    imshow(img_crop); hold on;
    patch('vertices',mesh2d','faces',model.faces,'FaceColor','blue','FaceAlpha',0.3,'EdgeColor','none');
    % object viewpoint
    subplot('position',[3/nplot 0 1/nplot 1]);
    h2 = subplot('position',[3/nplot 0 1/nplot 1]);
    trisurf(model.faces,model.vertices(:,1),model.vertices(:,2),model.vertices(:,3),'EdgeColor','none');axis equal;
    view(0,-90);
    set(h2,'XTick',[],'YTick',[]);
    set(h2,'visible','off')

end

author: geopavlakos


text:
function showPrepDefaults(EEG, FID)
% Displays information about the defaults used for PREP pipeline steps
%
% Parameters:
%     EEG         a structure compatible with EEGLAB EEG structure
%     FID         (optional) an open file descriptor indication which
%                 file to output to. By default, this is the console
%                 or command window

%%
types = {'boundary', 'resample', 'globaltrend', 'detrend', ...
                'linenoise', 'reference', 'report', 'postprocess'};
if nargin < 1 || ~isstruct(EEG)
    error('showPrepDefaults:NotEnoughArguments', ...
         'first argument must be a structure');
elseif nargin < 2 || ~exist('FID', 'var')
    FID = 1;
end

indent = '  ';
for k = 1:length(types)
    theseDefaults = getPrepDefaults(EEG, types{k});
    theseFields = fieldnames(theseDefaults);
    fprintf(FID, '\nDefaults for %s:\n', types{k});
    for j = 1:length(theseFields)
        printNext(theseFields{j}, theseDefaults.(theseFields{j}));
    end
end

    function printNext(theField, theStruct)
       fprintf(FID, '%s%s:\n', indent, theField);
       fprintf(FID, '%s%svalue: %s\n', indent, indent, ...
           getString(theStruct.value));
       fprintf(FID, '%s%sclasses: %s\n', indent, indent, ...
           getString(theStruct.classes));
       fprintf(FID, '%s%sattributes: %s\n', indent, indent, ...
           getString(theStruct.attributes));
       fprintf(FID, '%s%sdescription:\n', indent, indent);
       printLines(FID, theStruct.description, 80, [indent indent indent]);
    end
end

author: VisLab

repo: EEG-Clean-Tools


text:
function elemType = updatetype(elemType,brother)
%% UPDATETYPE update the type of triangles.
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

if (size(brother,1)==0), return; end
oldNT = length(elemType);
newNT = max(brother(:,2)); % brother(t,1) < brother(t,2)
if oldNT < newNT % coarse grid to fine grid
    elemType(newNT) = 0;
    firstBrother = zeros(oldNT,1,'uint32');
    idx = (brother(:,1) <= oldNT);
    t = brother(idx,1);
    tt = brother(idx,2);
    firstBrother(t(end:-1:1)) = tt(end:-1:1);
    % update first bisected elements from 1:oldNT first
    elemType(firstBrother(t)) = ~elemType(t);
    elemType(t) = ~elemType(t);
    % updated second bisected elements
    idx = (brother(:,2) > (oldNT + length(find(firstBrother))));
    elemType(brother(idx,2)) = ~elemType(brother(idx,1));
    elemType(brother(idx,1)) = elemType(brother(idx,2));
else % fine grid to coarse grid
    elemType(brother(:,1)) = ~elemType(brother(:,1));
    elemType(brother(:,2)) = [];
end

repo: ifem


text:
function mat_time = jd2MatTime (jdate)

% convert Julian date to Gregorian (calendar) date

% input

%  jdate = julian day

% output

%  month = calendar month [1 - 12]
%  day   = calendar day [1 - 31]
%  year  = calendar year [yyyy]

%  note: day may include fractional part

% Orbital Mechanics with Matlab

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

jd = jdate;

z = fix(jd + .5);
fday = jd + .5 - z;

if (fday < 0)
   fday = fday + 1;
   z = z - 1;
end

if (z < 2299161)
   a = z;
else
   alpha = floor((z - 1867216.25) / 36524.25);
   a = z + 1 + alpha - floor(alpha / 4);
end

b = a + 1524;
c = fix((b - 122.1) / 365.25);
d = fix(365.25 * c);
e = fix((b - d) / 30.6001);
day = b - d - fix(30.6001 * e) + fday;

if (e < 14)
   month = e - 1;
else
   month = e - 13;
end

if (month > 2)
   year = c - 4716;
else
   year = c - 4715;
end

mat_time = datenum(year, month, day, 0, 0, 0);

author: goGPS-Project


text:
% PQ_SIZE returns the number of elements in the priority queue
%
% SYNTAX
% sz = pq_size(pq)
%
% INPUT PARAMETERS
%   pq: a pointer to the priority queue
%
% OUTPUT PARAMETERS
%   sz:  the number of elements in the priority queue
%
% DESCRIPTION
% Queries the priority queue for the number of elements that it contains.
% This number is not the "capacity" or the maximum number of elements which
% is possible to insert but rather the number of elements CURRENTLY in the
% priority queue
%
% See also:
% PQ_DEMO, PQ_CREATE, PQ_PUSH, PQ_POP, PQ_SIZE, PQ_TOP, PQ_DELETE
%
% References:
% Gormen, T.H. and Leiserson, C.E. and Rivest, R.L., "introduction to
% algorithms", 1990, MIT Press/McGraw-Hill, Chapter 6.

% Copyright (c) 2008 Andrea Tagliasacchi
% All Rights Reserved
% email: andrea.tagliasacchi@gmail.com
% $Revision: 1.0$  Created on: May 22, 2009

author: Cloud-CV

repo: object-proposals


text:
function sparse_unique_index = sparse_grid_mixed_unique_index ( dim_num, ...
  level_max, rule, alpha, beta, tol, point_num, point_total_num )

%*****************************************************************************80
%
%% SPARSE_GRID_MIXED_UNIQUE_INDEX maps nonunique points to unique points.
%
%  Discussion:
%
%    The sparse grid usually contains many points that occur in more
%    than one product grid.
%
%    When generating the point locations, it is easy to realize that a point
%    has already been generated.
%
%    But when it's time to compute the weights of the sparse grids, it is
%    necessary to handle situations in which weights corresponding to
%    the same point generated in multiple grids must be collected together.
%
%    This routine generates ALL the points, including their multiplicities,
%    and figures out a mapping from them to the collapsed set of unique points.
%
%    This mapping can then be used during the weight calculation so that
%    a contribution to the weight gets to the right place.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 March 2011
%
%  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.
%
%    Input, integer RULE(DIM_NUM), the rule in each dimension.
%     1, "CC",  Clenshaw Curtis, Closed Fully Nested rule.
%     2, "F2",  Fejer Type 2, Open Fully Nested rule.
%     3, "GP",  Gauss Patterson, Open Fully Nested rule.
%     4, "GL",  Gauss Legendre, Open Weakly Nested rule.
%     5, "GH",  Gauss Hermite, Open Weakly Nested rule.
%     6, "GGH", Generalized Gauss Hermite, Open Weakly Nested rule.
%     7, "LG",  Gauss Laguerre, Open Non Nested rule.
%     8, "GLG", Generalized Gauss Laguerre, Open Non Nested rule.
%     9, "GJ",  Gauss Jacobi, Open Non Nested rule.
%    10, "GW",  Golub Welsch, (presumed) Open Non Nested rule.
%    11, "CC_SE", Clenshaw Curtis Slow Exponential, Closed Fully Nested rule.
%    12, "F2_SE", Fejer Type 2 Slow Exponential, Closed Fully Nested rule.
%    13, "GP_SE", Gauss Patterson Slow Exponential, Closed Fully Nested rule.
%    14, "CC_ME", Clenshaw Curtis Moderate Exponential, Closed Fully Nested rule.
%    15, "F2_ME", Fejer Type 2 Moderate Exponential, Closed Fully Nested rule.
%    16, "GP_ME", Gauss Patterson Moderate Exponential, Closed Fully Nested rule.
%    17, "CCN", Clenshaw Curtis Nested, Linear, Closed Fully Nested rule.
%
%    Input, real ALPHA(DIM_NUM), BETA(DIM_NUM), parameters used for
%    Generalized Gauss Hermite, Generalized Gauss Laguerre, and Gauss Jacobi rules.
%
%    Input, real TOL, the tolerance for point equality.
%
%    Input, integer POINT_NUM, the number of unique points in the grid.
%
%    Input, integer POINT_TOTAL_NUM, the total number of points in the grid.
%
%    Output, integer SPARSE_UNIQUE_INDEX(POINT_TOTAL_NUM), lists, for each
%    (nonunique) point, the corresponding index of the same point in the
%    unique listing.
%

%
%  Special cases.
%
  if ( level_max < 0 )
    sparse_unique_index = [];
    return
  end

  if ( level_max == 0 )
    sparse_unique_index(1) = 1;
    return
  end
%
%  Generate SPARSE_TOTAL_ORDER and SPARSE_TOTAL_INDEX arrays for the TOTAL set of points.
%
  sparse_total_order = zeros ( dim_num, point_total_num );
  sparse_total_index = zeros ( dim_num, point_total_num );

  point_total_num2 = 0;
%
%  The outer loop generates values of LEVEL.
%
  level_min = max ( 0, level_max + 1 - dim_num );

  for level = level_min : level_max
%
%  The middle loop generates a GRID,
%  based on the next partition that adds up to LEVEL.
%
    level_1d = [];
    more_grids = 0;
    h = 0;
    t = 0;

    while ( 1 )

      [ level_1d, more_grids, h, t ] = comp_next ( level, dim_num, level_1d, ...
        more_grids, h, t );

      order_1d = level_to_order_default ( dim_num, level_1d, rule );
%
%  The inner loop generates a POINT of the GRID of the LEVEL.
%
      point_index = [];
      more_points = 0;

      while ( 1 )

        [ point_index, more_points ] = vec_colex_next3 ( dim_num, order_1d, ...
          point_index, more_points );

        if ( ~more_points )
          break
        end

        point_total_num2 = point_total_num2 + 1;
        sparse_total_order(1:dim_num,point_total_num2) = order_1d(1:dim_num);
        sparse_total_index(1:dim_num,point_total_num2) = point_index(1:dim_num);

      end

      if ( ~more_grids )
        break
      end

    end

  end
%
%  Now compute the coordinates of the TOTAL set of points.
%
  sparse_total_point = zeros ( dim_num, point_total_num );
  sparse_total_point(1:dim_num,1:point_total_num) = r8_huge ( );

  for dim = 1 : dim_num

    for level = 0 : level_max

      order = level_to_order_default ( 1, level, rule(dim) );

      if ( rule(dim) == 1 )
        points = clenshaw_curtis_compute_points ( order );
      elseif ( rule(dim) == 2 )
        points = fejer2_compute_points ( order );
      elseif ( rule(dim) == 3 )
        points = patterson_lookup_points ( order );
      elseif ( rule(dim) == 4 )
        points = legendre_compute_points ( order );
      elseif ( rule(dim) == 5 )
        points = hermite_compute_points ( order );
      elseif ( rule(dim) == 6 )
        points = gen_hermite_compute_points ( order, alpha(dim) );
      elseif ( rule(dim) == 7 )
        points = laguerre_compute_points ( order );
      elseif ( rule(dim) == 8 )
        points = gen_laguerre_compute_points ( order, alpha(dim) );
      elseif ( rule(dim) == 9 )
        points = jacobi_compute_points ( order, alpha(dim), beta(dim) );
      elseif ( rule(dim) == 10 )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'SPARSE_GRID_MIXED_SIZE - Fatal error!\n' );
        fprintf ( 1, '  Do not know how to assign points for rule 10.\n' );
        error ( 'SPARSE_GRID_MIXED_SIZE - Fatal error!' );
      elseif ( rule(dim) == 11 )
        points = clenshaw_curtis_compute_points ( order );
      elseif ( rule(dim) == 12 )
        points = fejer2_compute_points ( order );
      elseif ( rule(dim) == 13 )
        points = patterson_lookup_points ( order );
      elseif ( rule(dim) == 14 )
        points = clenshaw_curtis_compute_points ( order );
      elseif ( rule(dim) == 15 )
        points = fejer2_compute_points ( order );
      elseif ( rule(dim) == 16 )
        points = patterson_lookup_points ( order );
      elseif ( rule(dim) == 17 )
        points = ccn_compute_points ( order );
      else
        fprintf ( 1, '\n' );
        fprintf ( 1, 'SPARSE_GRID_MIXED_SIZE - Fatal error!\n' );
        fprintf ( 1, '  Unexpected value of RULE = %d\n', rule(dim) );
        error ( 'SPARSE_GRID_MIXED_SIZE - Fatal error!' );
      end

      index = find ( sparse_total_order(dim,1:point_total_num) == order );

      sparse_total_point(dim,index) = points ( sparse_total_index(dim,index) );

    end

  end
%
%  Now determine the mapping from nonunique points to unique points.
%  We can not really use the UNDX output right now.
%
  [ undx, sparse_unique_index ] = r8col_undex ( dim_num, point_total_num, ...
    sparse_total_point, point_num, tol );

  return
end


text:
function s = obj2segs(obj)

% OBJ2SEGS Get segments matrix from 3D object.
%   OBJ2SEGS(OBJ) returns a 3-by-N matrix of segments from the 3D
%   information in structure OBJ, stored in the N-by-3 matrix OBJ.vert,
%   with the object vertices, and the M-by-2 matrix OBJ.seg, with the
%   couples of indices in OBJ.vert defining each segment.
%
%   For example, the square:
%
%     1 +-----+ 2
%       |     |
%       |     |
%     3 +-----+ 4
%
%   is defined in OBJ as
%       OBJ.vert = [0 0 0; 1 0 0; 0 1 0; 1 1 0]
%       OBJ.seg  = [1 2; 2 4; 4 3; 3 1]
%   then the produced segments matrix SEG = OBJ2SEGS(OBJ) is
%       SEG = [ 0 1 1 0
%               0 0 1 1
%               0 0 0 0
%               1 1 0 0
%               0 1 1 0
%               0 0 0 0 ]

%   Copyright 2009 Joan Sola @ LAAS-CNRS.

p = obj2pnts(obj);    % endpoints in 3-by-N format
M = size(obj.seg,1);  % number of segments M
s = zeros(6,M);       % initialize output segments matrix

for i = 1:M  % for each segment ...

    j = obj.seg(i,1); % first endpoint's index
    k = obj.seg(i,2); % second endpoint's index
    s(:,i) = [p(:,j);p(:,k)]; % the segment!

end

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB 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.
%
%   SLAMTB 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 SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function mustBeZmapCatalog(x)
    % error unless provided value is a valid shape
    if ~isa(x,'ZmapCatalog') || ~isvalid(x)
        error('value must be a ZmapCatalog')
    end
end

author: CelsoReyes


text:
function B=rotmat(A,n,dim)
%ROTMAT  Rotate matrix along specified dimension.
%   ROTMAT(X), returns X.
%   ROTMAT(X,N), moves elements in matrix X with N steps in
%      a wrapping/rotating fashion. Rotation will be done row-wise.
%      Positive N will mean that elements are moved towards higher
%      row indicies. Last element(s) will be moved to the first row,
%      aso. Negative N will make the elements move upwards towards
%      lower row indicies. First element(s) will be moved to the last
%      row, aso. If N is zero or empty, no rotation will occur
%      and X will be returned.
%   ROTMAT(X,N,DIM), rotates matrix elements along dimension DIM.
%      If DIM is a singleton dimension, X will be returned.
%      Likewise if DIM < 1 or DIM > ndims(X), nothing will happen.
%      N must be scalar valued.
%
%   Examples:
%      X=rand(3,3,3)          %X is 3x3x3.
%      rotmat(X,1)            %rotates one step along rows.
%      rotmat(X,-2)           %rotates two steps along rows, opposite direction.
%      rotmat(X,1,2)          %rotates one step along columns.
%      rotmat(X,1,3)          %rotates one step along third dimension.
%
%   See also DELMAT, INSMAT, SHIFTMAT, REPMAT, RESHAPE, FLIPDIM.

% Copyright (c) 2003-10-28, B. Rasmus Anthin.
% Revision 2003-10-29.
% GPL license, freeware.

error(nargchk(1,3,nargin))
if nargin<2, n=0;end                    %do nothing
if isempty(n), n=0;end
if prod(size(n))~=1, error('N must be a scalar.'),end
n=round(n);
if nargin<3, dim=1;end                  %rotate row vectors
if dim<1, dim=ndims(A)+1;end            %if less than one, do nothing
dims=[dim:max(ndims(A),dim) 1:dim-1];
if dims(1)<=ndims(A) & n
   A=permute(A,dims);                             %move dimension to front
   sizA=size(A);                                  %get dimension sizes
   A=reshape(A,[sizA(1) prod(sizA(2:end))]);      %reshape to a 2-D matrix
   n=rem(n-1,sizA(1))+1;                          %wrapping
   idx=1:sizA(1);                                 %the row indices
   if n<0
      idx=[idx(1-n:end) idx(1:-n)];               %rotate row indices (negative rotation)
   else
      idx=[idx(end-n+1:end) idx(1:end-n)];        %rotate row indices (positive rotation)
   end
   B=A(idx,:);
   B=reshape(B,[size(B,1) sizA(2:end)]);          %back to previous shape minus removed vectors
   B=ipermute(B,dims);                            %replace dimensions to initial location
else
   B=A;
end


text:
function [out1] = getAppDir(inp1)
out1 = [];
ns = getNamespace();
if isempty(ns)
    return;
end
if strcmp(ns, 'AtlasViewerGUI')
    if nargin == 0
        [out1] = getAppDir_AtlasViewerGUI();
    elseif nargin == 1
        [out1] = getAppDir_AtlasViewerGUI(inp1);
    end
elseif strcmp(ns, 'Homer3')
    if nargin == 0
        [out1] = getAppDir_Homer3();
    elseif nargin == 1
        [out1] = getAppDir_Homer3(inp1);
    end
elseif strcmp(ns, 'DataTreeClass')
    if nargin == 0
        [out1] = getAppDir_DataTreeClass();
    elseif nargin == 1
        [out1] = getAppDir_DataTreeClass(inp1);
    end
end

% ---------------------------------------------------------
function [dirname] = getAppDir_AtlasViewerGUI(isdeployed_override)

if ~exist('isdeployed_override','var')
    isdeployed_override = 'notdeployed';
end

if isdeployed() || strcmp(isdeployed_override, 'isdeployed')
    if ispc()
        dirname = 'c:/users/public/atlasviewer/';
    else
        currdir = pwd;
        cd ~/;
        dirnameHome = pwd;
        dirname = [dirnameHome, '/atlasviewer/'];
        cd(currdir);
    end
else
    dirname = fileparts(which('AtlasViewerGUI.m'));
end

dirname(dirname=='\') = '/';
if dirname(end) ~= '/'
    dirname(end+1) = '/';
end

% ---------------------------------------------------------
function [dirname] = getAppDir_Homer3(isdeployed_override)

if ~exist('isdeployed_override','var')
    isdeployed_override = 'notdeployed';
end

if isdeployed() || strcmp(isdeployed_override, 'isdeployed')
    if ispc()
        dirname = 'c:/users/public/homer3/';
    else
        currdir = pwd;
        cd ~/;
        dirnameHome = pwd;
        dirname = [dirnameHome, '/homer3/'];
        cd(currdir);
    end
else
    dirname = fileparts(which('Homer3.m'));
end

dirname(dirname=='\') = '/';
if dirname(end) ~= '/'
    dirname(end+1) = '/';
end

% ---------------------------------------------------------
function [dirname] = getAppDir_DataTreeClass(isdeployed_override)

if ~exist('isdeployed_override','var')
    isdeployed_override = 'notdeployed';
end

if isdeployed() || strcmp(isdeployed_override, 'isdeployed')
    if ispc()
        dirname = 'c:/users/public/datatree/';
    else
        currdir = pwd;
        cd ~/;
        dirnameHome = pwd;
        dirname = [dirnameHome, '/datatree/'];
        cd(currdir);
    end
else
    dirname = fileparts(which('DataTreeClass.m'));
end

dirname(dirname=='\') = '/';
if dirname(end) ~= '/'
    dirname(end+1) = '/';
end

author: BUNPC


text:
function result = fit_maxwell_pdf( x,y,W,hAx )
% fit_maxwell_pdf - Non Linear Least Squares fit of the maxwellian distribution.
%                   given the samples of the histogram of the samples, finds the
%                   distribution parameter that fits the histogram samples.
%
%    fits data to the probability of the form:
%        p(r) = sqrt(2/pi)*(a^(-3/2))*(r^2)*exp(-(r^2)/(2*a))
%    with parameter: a
%
% format:   result = fit_maxwell_pdf( x,y,W,hAx )
%
% input:    y   - vector, samples of the histogram to be fitted
%           x   - vector, position of the samples of the histogram (i.e. y = f(x,a))
%           W   - matrix or scalar, a square weighting matrix of the size NxN where
%                 N = length(y), or 0 to indicate no weighting is needed.
%           hAx - handle of an axis, on which the fitted distribution is plotted
%                 if h is given empty, a figure is created.
%
% output:   result  - structure with the fields
%                      a   - fitted parameter
%                      VAR - variance of the estimation
%                      type- weighted LS or not weighted LS
%                      iter- number of iteration for the solution
%

%
% Algorithm
% ===========
%
% We use the WLS algorithm to estimate the PDF from the samples.%
% The maxwell distribution is given by:
%
%    p(x,a) = sqrt(2/pi)*(a^(-3/2))*(x.^2).*exp(-(x.^2)/(2*a))
%           = Const * (a^(-3/2)) .* exp(-(x.^2)/(2*a))
%
%    note that X is known and therefore, considered a constant vector
%
% The non liner WLS estimator is given by:
%
%    a(n+1) = a(n) + inv(H'*W*H)*(H') * (y-h) = a(n) + G * err
%
%    where:   h = p(x,a)
%             H = diff( p(x,a) ) with respect to "a"
%             W = weighting matrix of size NxN  (N = length(y))
%             a = a single parameter to be estimated
%
% The error estimation is given by:
%
%    VAR( a ) = G * VAR( err ) * (G')
%
%       or when W=I and the noise is a gaussian noise
%
%    VAR( a ) = inv( H' * H )
%

if (nargin<3)
    error( 'fit_maxwell_pdf - insufficient input arguments' );
end

a       = x(find(y==max(y)))^2;         % initial guess
y       = y(:);                         % both should be column vectors !
x       = x(:);
x2      = x.^2;                         % save computation time
C       = sqrt(2/pi)*x2;                % a constant vector
thresh  = 0.995;                        % convergence threshold for the loop
last_cnt= inf;
iter    = 0;

% check weight matrix input
if (size(W,1)==length(y)) & (size(W,2)==length(y))
    weights_flag    = 1;
    type            = 'WLS';
else
    weights_flag    = 0;
    type            = 'LS';
end

% Estimation
% =============
if (weights_flag)
    % loop for convergence (with weighting matrix)
    % =============================================
    while (1)
        iter    = iter + 1;
        h       = C*(a^(-1.5)).*exp(-x2/(2*a));
        H       = h.*( x2/(2*a^2) - 3/(2*a) );
        HTW     = H'*W;
        e       = inv( HTW * H ) * HTW * (y-h);
        a       = a + e;
        control = e*e;
        if ( control > (last_cnt * thresh) )
            break;
        else
            last_cnt = control;
        end
    end

    % summarize results
    h           = C*(a^(-1.5)).*exp(-x2/(2*a));
    H           = h.*( x2/(2*a^2) - 3/(2*a) );
    HTW         = H'*W;
    G           = inv( HTW * H ) * HTW;
    err         = ( y - h );
    result.a    = a;
    result.VAR  = G * var( err ) * (G');
    result.RMS  = sqrt( (err')*err/ (x(2)-x(1))^2 / (length(err)-1) );
    result.iter = iter;
    result.type = type;
else

    % loop for convergence (without a weighting matrix) - assume white noise
    % ========================================================================
    while (1)
        iter    = iter + 1;
        h       = C*(a^(-1.5)).*exp(-x2/(2*a));
        H       = h.*( x2/(2*a^2) - 3/(2*a) );
        HT      = H';
        control = inv( HT * H );
        a       = a + control * HT * (y-h);
        if ( control>(last_cnt * thresh) )
            break;
        else
            last_cnt = control;
        end
    end

	% summarize results
    h           = C*(a^(-1.5)).*exp(-x2/(2*a));
	H           = h.*( x2/(2*a^2) - 3/(2*a) );
    err         = ( y - h );
    result.a    = a;
    result.VAR  = inv( (H') * H );
    result.RMS  = sqrt( (err')*err/ (x(2)-x(1))^2 / (length(err)-1) );
    result.iter = iter;
    result.type = type;
end

% plot distribution if asked for
% ===============================
if (nargin>3)
    if ishandle( hAx )
        plot_maxwell( x,result,hAx,2 );
    else
        figure;
        plot_maxwell( x,result,gca,2 );
    end
end

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function [fig]=mvg_drawWindows(img,windows,fig)
% function [fig]=mvg_drawWindows(img,windows,numBest) draws the given windows
% on the given input image.
%
% Inputs:
% img, imgRow*imgCol*1 or 3, double, is the input image.
% windows, numWindows*4, double, is a matrix containing the windows.
%                                   Each row corresponds to one window in format:
%                                   windows(i,:)=[xmin,ymin,xmax,ymax];
% fig, 1*1, double, is a figure handle pointing to figure where image and
%                   windows are drawn. If not given, new figure will be
%                   created.
%
% Outputs:
% fig, 1*1, double, is a figure handle pointing to the figure where image
%                   and windows were drawn.
%

% 2011 MVG, Oulu, Finland, Esa Rahtu and Juho Kannala
% 2011 VGG, Oxford, UK, Matthew Blaschko

%% User defined parameters
linewidth = 3;
base_color = [1 0 0];

%% Open figure if needed
if nargin<3
    fig=figure;
end

%% Draw image
figure(fig); clf;
imshow(img);

%% Draw windows
hold on;
for idx = 1:size(windows,1)
    plot(windows(idx,[1 3 3 1 1]),windows(idx,[2 2 4 4 2]),'Color',base_color,'linewidth',linewidth);
end


text:
%   AUTHORSHIP
%   Primary Developer: Stephen Meehan <swmeehan@stanford.edu>
%   Math Lead & Secondary Developer:  Connor Meehan <connor.gw.meehan@gmail.com>
%   Bioinformatics Lead:  Wayne Moore <wmoore@stanford.edu>
%   Provided by the Herzenberg Lab at Stanford University
%   License: BSD 3 clause
%

classdef StringArray < handle
    properties(SetAccess=private)
        strings;
        N;
    end

    methods
        function this=StringArray(strings)
            this.strings=strings;
            this.N=length(strings);
        end
        function idx=lastIndexOf(this, search)
            idx=0;
            for i=this.N:-1:1
                if strcmp(this.strings{i}, search)
                    idx=i;
                    break;
                end
            end
        end

        function ok=contains(this,search)
            ok=this.indexOf(search)>0;
        end

        function idx=indexOf(this, search)
            idx=0;
            for i=1:this.N
                if strcmp(this.strings{i}, search)
                    idx=i;
                    break;
                end
            end
        end

        function idx=endsWith(this, with)
            idx=0;
            for i=1:this.N
                if endsWith(this.strings{i}, with)
                    idx=i;
                    break;
                end
            end
        end

        function indices=indicesOf(this, search)
            indices=[];
            for i=1:this.N
                if strcmp(this.strings{i}, search)
                    indices(end+1)=i;
                end
            end
        end

        function idx=indexOfI(this, search)
            idx=0;
            for i=1:this.N
                if strcmpi(this.strings{i}, search)
                    idx=i;
                    break;
                end
            end
        end

    end

    methods(Static)
        function ok=IsOk(x)
            ok=true;
            if iscell(x)
                N=length(x);
                for i=1:N
                    if ~ischar(x{i})
                        ok=false;
                        break;
                    end
                end
            elseif ischar(x)
            else
                ok=false;
                warning('Expecting string or string array');
            end
        end

        function strs=Trim(strs, limit)
            if length(strs)>limit
                strs{limit}=[num2str(length(strs)-limit) ' more...'];
                strs(limit+1:end)=[];
            end
        end
        function c=Num2Str(nums, prefix)
            if nargin<2
                prefix='#';
            end
            N_=length(nums);
            c=cell(1, N_);
            for i_=1:N_
                c{i_}=[prefix num2str(nums(i_))];
            end
        end

        function out=Intersection(left, right)
            if isempty(left)
                out=right;
            elseif isempty(right)
                out=left;
            else
                out={};
                r=StringArray.Set(right);
                N=length(left);
                for i=1:N
                    try
                        if r.contains(java.lang.String(left{i}))
                            out{end+1}=left{i};
                        end
                    catch ex
                    end
                end
            end
        end

        function out=Sort(in, canon)
            if nargin==1
                out=StringArray.Cell(StringArray.Set(in));
            else
                N=length(in);
                this=StringArray(canon);
                map=java.util.TreeMap;
                for i=1:N
                    str=in{i};
                    idx=this.indexOf(str);
                    if idx>0
                        map.put(idx, str);
                    else
                        map.put(N+i, str);
                    end
                end
                out=StringArray.Cell(map.values);
            end
        end

        function c=Cell(javaCollection)
            N=javaCollection.size;
            c=cell(1,N);
            it=javaCollection.iterator;
            for i=1:N
                c{i}=it.next;
            end
        end

        function out=JavaString2Chars(in)
            N=length(in);
            out=cell(1,N);
            for i=1:N
                out{i}=char(in{i});
            end
        end

        function c=JavaArrayToCell(javaArr)
            if iscell(javaArr)
                c=javaArr;
            else
                N=length(javaArr);
                c=cell(1,N);
                for i=1:N
                    c{i}=char(javaArr(i));
                end
            end
        end

        function lst=List(strs)
            lst=java.util.ArrayList;
            N=length(strs);
            for i=1:N
                lst.add(strs{i});
            end
        end

        function lst=Vector(strs)
            lst=java.util.Vector;
            N=length(strs);
            for i=1:N
                lst.add(strs{i});
            end
        end

        function lst=ListFromJavaArray(strs)
            lst=java.util.ArrayList;
            N=length(strs);
            for i=1:N
                lst.add(strs(i));
            end
        end
        function set=Set(strs, set)
            if nargin<2
                set=java.util.TreeSet;
            end
            N=length(strs);
            for i=1:N
                if ~isempty(strs{i})
                    set.add(java.lang.String(strs{i}));
                end
            end
        end

        function out=RemoveDuplicates(strs, set)
            if nargin<2
                set=StringArray.Set(strs);
            else
                set=StringArray.Set(strs, set);
            end
            N_=set.size;
            out=cell(1,N_);
            it=set.iterator;
            for i=1:N_
                out{i}=it.next;
            end
        end

        function [out, dupCnt]=ForceUnique(strs, prefix, suffix)
            dupCnt=0;
            map=java.util.HashMap;
            N_=length(strs);
            out=cell(1,N_);
            for i=1:N_
                if map.containsKey(strs{i})
                    cnt=map.get(strs{i})+1;
                else
                    cnt=1;
                end
                map.put(strs{i}, cnt);
                if cnt>1
                    out{i}=[strs{i} prefix num2str(cnt) suffix];
                    dupCnt=dupCnt+1;
                else
                    out{i}=strs{i};
                end
            end
        end

        function ok=Contains(searchedStrings, search)
            sa=StringArray(searchedStrings);
            ok=sa.indexOf(search)>0;
        end
        function idx=IndexOf(searchedStrings, search)
            sa=StringArray(searchedStrings);
            idx=sa.indexOf(search);
        end

        function idxs=IndexOfEmpties(searchedStrings)
            N=length(searchedStrings);
            idxs=[];
            for i=1:N
                if isempty(searchedStrings{i})
                    idxs(end+1)=i;
                end
            end
        end

        function idx=IndexOfIgnoreCase(searchedStrings, search)
            sa=StringArray(searchedStrings);
            idx=sa.indexOfI(search);
        end

        function idxs=FirstListIndexes(searched, searches)
            N=length(searches);
            idxs=zeros(1,N);
            for i=1:N
                idxs(i)=searched.indexOf(searches{i});
            end
            idxs=idxs+1;
        end

        function [idxs, foundAll]=EndsWith(searchedStrings, searches)
            idxs=[];
            sa=StringArray(searchedStrings);
            N=length(searches);
            for i=1:N
                idx=sa.endsWith(searches{i});
                if idx>0
                    idxs=[idxs idx+1];
                else
                    disp('huh');
                end
            end
            foundAll=length(idxs)>=length(searchedStrings);
        end

        %Not sure if a faster search than this is available at
        %https://www.mathworks.com/matlabcentral/answers/335930-most-efficient-way-to-search-in-text-arrays
        %I think that the the problem with find(strcmp(a,s),1) is that
        %strcmp does not STOP after its first successful strcmp
        %the Containers.map approach WILL ... but has a high startup
        %cost for under 1000 strings
        function idxs=IndexesOf(searchedStrings, searches)
            idxs=[];
            sa=StringArray.List(searchedStrings);
            N=length(searches);
            for i=1:N
                idx=sa.indexOf(searches{i});
                if idx>=0
                    idxs=[idxs idx+1];
                else
                    disp('huh');
                end
            end
        end

        function idxs=IndexesOf2(searchedStrings, searches)
            sa=StringArray.List(searchedStrings);
            N=length(searches);
            idxs=zeros(1,N);
            for i=1:N
                idxs(i)=sa.indexOf(searches{i})+1;
            end
        end

        function idxs=NotFound(searchedStrings, searches)
            idxs=[];
            sa=StringArray.List(searchedStrings);
            N=length(searches);
            for i=1:N
                idx=sa.indexOf(searches{i});
                if idx<0
                    idxs=[idxs i];
                end
            end
        end

        function [ok, sameButReordered]=AreSameOrEmpty(d1, d2)
            ok=isequal(d1, d2);
            sameButReordered=false;
            if ~ok
                if length(d1)==length(d2)
                    nf=StringArray.NotFound(d2, d1);
                    if isempty(nf)
                        sameButReordered=true;
                    elseif length(nf)==1
                        % is FMO?
                        if isempty(d1{nf}) || isempty(d2{nf})
                            d1_=d1;
                            d2_=d2;
                            d1_{nf}='dummy';
                            d2_{nf}='dummy';
                            ok=isequal(d1_,d2_);
                        end
                    end
                end
            end
        end

        function [allFound, idxs]=Find(d1, d2, allowOneEmptyLikeFMO)
            if nargin<3
                allowOneEmptyLikeFMO=false;
            end
            allFound=true;
            N=length(d1);
            idxs=zeros(1,N);

            sa2=StringArray(d2);
            for i=1:N
                ii=sa2.indexOf(d1{i});
                if ii<1
                    if allowOneEmptyLikeFMO
                        allowOneEmptyLikeFMO=false;
                        [emptyCount, ii]=StringArray.GetEmpties(d2);
                        if emptyCount ~= 1
                            [emptyCount, i2]=StringArray.GetEmpties(d1);
                            if emptyCount ~= 1
                                allFound=false;
                                ii=0;
                            else
                                ii=0;
                                if i2==i
                                    sa1=StringArray(d1);
                                    N2=length(d2);
                                    for j=1:N2
                                        jj=sa1.indexOf(d2{j});
                                        if jj<1
                                            ii=j;
                                            break;
                                        end
                                    end
                                end
                                if ii<1
                                    allFound=false;
                                end
                            end
                        end
                    else
                        allFound=false;
                    end
                end
                if ~isempty(ii)
                    idxs(i)=ii(1);
                end
            end
        end

        function str=toString(stringArray, delimiter, useQuotes)
            if nargin<3
                useQuotes=false;
                if nargin<2
                    delimiter=', ';
                end
            end
            N=length(stringArray);
            str='';
            for i=1:N
                v=stringArray{i};
                if ~ischar(v)
                    v=String.toString(v);
                end
                if useQuotes
                    v=['"' v '"'];
                end
                if (i>1)
                    str=[str delimiter v];
                else
                    str=v;
                end
            end
        end

        function ok=isEqual(stringArray1, stringArray2)
            N1=length(stringArray1);
            N2=length(stringArray2);
            ok=false;
            if N1==N2
                ok=true;
                for i=1:N1
                    if ~strcmp(stringArray1{i}, stringArray2{i})
                        ok=false;
                        break;
                    end
                end
            end
        end

        function a=Join(a, a2)
            N=length(a2);
            for i=1:N
                a{end+1}=a2{i};
            end
        end

        function stringArray=Subset(strings, idxs)
            stringArray=StringArray( strings(idxs) );
        end

        function [searchResultIdxs, unSearchedStrings, unFoundStrings]=...
                SearchSubset(searchedStrings, subsettingIndices, searchStrings)
            N=length(searchStrings);
            unFoundStrings={};
            searchResultIdxs=1:N;
            subset=StringArray.Subset(searchedStrings, subsettingIndices);
            unSearchedStrings=subset.strings;
            for i=N:-1:1
                idx=subset.indexOf(searchStrings{i});
                searchResultIdxs(i)=idx;
                if idx>0
                    %remove found string
                    unSearchedStrings(idx)=[];
                else
                    unFoundStrings{end+1}=searchStrings{i};
                end
            end
        end

        function [searchResultIdxs, unsearchedStrings, unFoundStrings]=...
                Search2Subsets(searchedStrings1, searchedStrings2,...
                subsettingIndices, searchStrings)
            N=length(searchStrings);

            unFoundStrings={};
            searchResultIdxs=1:N;
            subset1=StringArray.Subset(searchedStrings1, subsettingIndices);
            subset2=StringArray.Subset(searchedStrings2, subsettingIndices);
            temp=subset1.strings;
            for i=N:-1:1
                searchString=searchStrings{i};
                idx=subset1.indexOf(searchString);
                if idx==0
                    idx=subset2.indexOf(searchString);

                end
                searchResultIdxs(i)=idx;
                if idx>0
                    %remove found string
                    temp{idx}='';
                else
                    unFoundStrings{end+1}=searchString;
                end
            end
            unsearchedStrings={};
            for i=1:length(temp)
                if ~isempty(temp{i})
                    unsearchedStrings{end+1}=temp{i};
                end
            end
        end

         function[files]=getField(handles, fileField)
             files={};
             if isfield(handles, fileField)
                 temp=getfield(handles,fileField);
                 if ~strcmp('cell', class(temp))
                     files={temp};
                 else
                     files=temp;
                 end
             end
         end

         function array=toArray(something)
             if ~strcmp('cell', class(something))
                 array={something};
             else
                 array=something;
             end
         end

         function s=encodeCell(c, quotes)
             if nargin==1
                 quotes=false;
             end
             N=length(c);
             s='';
             for i=1:N
                 if i>1
                     s=[s ', '];
                 end
                 if quotes
                     s=[s '''' c{i} ''''];
                 else
                     s=[s c{i}];
                 end
             end
             if quotes
                 s=['{' s '}'];
             end
         end

         function c=decodeCell(s)
             if String.StartsWith(s, '{''')
                 c=eval(s);
             else
                 c=regexp(s, '[^,\s]*', 'match');
             end
         end

         function idxs=StartsWith(sa, str)
             idxs=[];
             N=length(sa);
             for i=1:N
                 if startsWith(str, sa{i})
                     idxs(end+1)=i;
                 end
             end
         end

         function ok=EqualsIndexed(these, those, idxs)
             try
                 N_=length(idxs);
                 N2=length(these);
                 N3=length(those);
                 ok=true;
                 for i=1:N_
                     idx=idxs(i);
                     if idx>N2 || idx>N3
                         continue;
                     end
                     if ~strcmp(these{idx}, those{idx})
                         ok=false;
                         break;
                     end
                 end
             catch exception
                 warning(exception.message);
                 ok=false;
             end
         end

         function ok=Equals(these, those)
             N1=length(these);
             N2=length(those);
             if N1 ~= N2
                 ok=false;
             else
                 ok=StringArray.EqualsIndexed(these, those, [1:N1]);
             end
         end

         function found=HasAny(superSet, subSet, superSetValues)
             found=false;
             for i=1:length(subSet)
                 idx=StringArray.IndexOf(superSet, subSet{i});
                 if idx>0
                     if nargin<2 || ~isempty(superSetValues{idx})
                         found=true;
                         break;
                     end
                 end
             end
         end

         function found=HasAll(superSet, subSet, superSetValues)
             found=true;
             for i=1:length(subSet)
                 idx=StringArray.IndexOf(superSet, subSet{i});
                 if idx>0
                     if nargin>1 && isempty(superSetValues{idx})
                         found=false;
                         break;
                     end
                 else
                     found=false;
                     break;
                 end
             end
         end

         function [arrayOfArrayIdx, arrayIdx]=FindAofA(arrayOfArrays, searchValue)
             arrayOfArrayIdx=0;
             arrayIdx=0;
             N=length(arrayOfArrays);
             for i=1:N
                 j=StringArray.IndexOf(arrayOfArrays{i}, searchValue);
                 if j>0
                     arrayIdx=j;
                     arrayOfArrayIdx=i;
                     return;
                 end
             end
         end

         function ok=ContainsInAofA(arrayOfArrays, arrayIdx)
             ok=false;
             N=length(arrayOfArrays);
             if N>=arrayIdx
                 j=StringArray.IndexOf(arrayOfArrays{arrayIdx}, searchValue);
                 if j>0
                   ok=true;
                 end
             end
         end

         function [lbl, sz]=GetLongest(strs)
             N=length(strs);
             sz=0;
             lbl=[];
             for i=1:N
                 ln=length(strs{i});
                 if ln>sz
                     sz=ln;
                     lbl=strs{i};
                 end
             end
         end

         function out=ToLower(strs)
             N=length(strs);
             out=cell(1,N);
             for i=1:N
                 out{i}=lower(strs{i});
             end
         end

         function [cnt, idxs]=GetEmpties(strs)
             cnt=0;
             idxs=[];
             N=length(strs);
             for i=1:N
                 if isempty(strs{i})
                     cnt=cnt+1;
                     idxs(end+1)=i;
                 end
             end
         end

         function ok=IsValid(in)
             if iscell(in)
                 ok=true;
                 N=length(in);
                 for i=1:N
                     if ~ischar(in{i})
                         ok=false;
                         break;
                     end
                 end
             else
                 ok=false;
             end

         end
    end
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function sphere_cubed_grid_lines_display ( point_num, xyz, line_num, ...
  line_data, filename )

%*****************************************************************************80
%
%% SPHERE_CUBED_GRID_LINES_DISPLAY displays the lines on a cubed sphere grid.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 May 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer LINE_NUM, the number of lines.
%
%    Input, real LINE_DATA(LINE_NUM,3,2), for each line I, the X/Y/Z
%    coordinates of the start and end of a line segment on the grid.
%
%    Input, string FILENAME, the name of the file into which a copy of the
%    graphics information should be stored.

%
%  Get the scale.
%
  for i = 1 : 3
    xyz_min(i) = min ( min ( line_data(:,i,:) ) );
    xyz_max(i) = max ( max ( line_data(:,i,:) ) );
  end

  xyz_range(1:3) = xyz_max(1:3) - xyz_min(1:3);

  margin = 0.025 * max ( xyz_range(1), ...
                   max ( xyz_range(2), xyz_range(3) ) );

  x_min = xyz_min(1) - margin;
  x_max = xyz_max(1) + margin;
  y_min = xyz_min(2) - margin;
  y_max = xyz_max(2) + margin;
  z_min = xyz_min(3) - margin;
  z_max = xyz_max(3) + margin;
%
%  Draw the picture.
%  As usual, trying to get Matlab's graphics routines to cooperate is
%  a nightmare.
%
  figure ( )
  clf
  hold on
  point_size = 50;
  point_color = [ 0.0, 1.0, 0.0 ];
  scatter3 ( xyz(:,1), xyz(:,2), xyz(:,3), point_size, 'k', 'filled' );

  for i = 1 : line_num
    xx = [ line_data(i,1,1); line_data(i,1,2) ];
    yy = [ line_data(i,2,1); line_data(i,2,2) ];
    zz = [ line_data(i,3,1); line_data(i,3,2) ];
    line ( xx, yy, zz, 'LineWidth', 2 );
  end
%
%  We want to include a sphere in the picture.  But if the sphere is close to
%  the correct radius, lines that should lie on the surface because they are
%  curved will actually tunnel through the surface because we draw them straight.
%
  [ x, y, z ] = sphere ( 40 );
  x = 0.95 * x;
  y = 0.95 * y;
  z = 0.95 * z;

  c = z ./ z;

  surf ( x, y, z, c, 'EdgeColor', 'None' );

  axis equal
  grid on
  xlabel ( '--X axis--' )
  ylabel ( '--Y axis--' )
  zlabel ( '--Z axis--' )
  title ( 'Sphere cubed grid' );
  view ( 3 )
  hold off

  print ( '-dpng', filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Saving plot as "%s".\n', filename );

  return
end

author: johannesgerer


text:
function pass = test_points(pref)

% Get preferences:
if ( nargin < 1 )
    pref = chebfunpref();
end

n = 10;
tol = 10*eps;
r = rand(n,1);
dom = [0, 1];

for k = 1:2
    f = chebfun(r, dom, 'chebkind', k);
    pass(k) = norm(f.points - chebpts(n, dom, k)) < tol;
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


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

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

% Last Modified by GUIDE v2.5 28-Jun-2012 15:16:11

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

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

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

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

axes(handles.axLogo);
path = which('splashlogo.jpg');
im = imread(path);
image(im);
axis image
axis off

% SLASH = fastif(isunix,'/','\');

% fid = fopen([fileparts(which('StartSIFT.m')) filesep 'resources' filesep 'version.txt']);
% version = fscanf(fid,'%s');
% fclose(fid);

splashTextFile = [fileparts(which('StartSIFT.m')) filesep 'Readme.txt'];

set(handles.txtSplashText,'String',...
    importdata(splashTextFile,'\n'));

% set(handles.txtSplashText,'String',...
%     {sprintf('Welcome to the Source Information Flow Toolbox version %s',version), ...
%      'Author: Tim Mullen (tim@sccn.ucsd.edu)' ...
%      'With important contributions from Arnaud Delorme and Christian Kothe' ...
%      '    ' ...
%      'DISCLAIMER:  THIS IS AN EXPERIMENTAL *ALPHA* VERSION OF SIFT. THIS VERSION IS NOT A SUPPORTED RELEASE AND IS INTENDED FOR EDUCATIONAL PURPOSES. USE AT YOUR OWN RISK! CHECK http://sccn.ucsd.edu/wiki/SIFT FOR UPDATES ON THE OFFICIAL BETA RELEASE.' ...
%      '    ' ...
%      'If you find this toolbox useful for your research, PLEASE include the following citations with any publications and/or presentations which make use of SIFT:' ...
%      '    ' ...
%      '(1) Mullen, T, Delorme, A. Kothe, C, Makeig, S (2010) "An Electrophysiological Information Flow Toolbox for EEGLAB" Society for Neuroscience 2010, San Diego, CA' ...
%      '    ' ...
%      '(2) Delorme, A., Mullen, T., Kothe C., Akalin Acar, Z., Bigdely Shamlo, N., Vankov, A., Makeig, S. (2011) "EEGLAB, SIFT, NFT, BCILAB, and ERICA: New tools for advanced EEG/MEG processing." Computational Intelligence and Neuroscience vol. 2011, Article ID 130714, 12 pages.' ...
%      });

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes gui_splashscreen wait for user response (see UIRESUME)
uiwait(handles.gui_splash);

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

% Get default command line output from handles structure
try,
varargout{1} = handles.output;
close(hObject);
catch
varargout{1} = [];
end

function gui_splash_CreateFcn(hObject, eventdata, handles)

function cmdOK_Callback(hObject, eventdata, handles)

uiresume(handles.gui_splash);

% --- Executes during object creation, after setting all properties.
function txtSplashText_CreateFcn(hObject, eventdata, handles)
% hObject    handle to txtSplashText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

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

% Hints: get(hObject,'String') returns contents of txtSplashText as text
%        str2double(get(hObject,'String')) returns contents of txtSplashText as a double

% --- Executes during object creation, after setting all properties.
function txtSplashScreen_CreateFcn(hObject, eventdata, handles)
% hObject    handle to txtSplashText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

repo: REST


text:
function [ fm ] = gsp_modulate( G,f,k )
%GSP_MODULATE Generalized modulation of the signal f to the frequency k
%   Usage: fm = gsp_modulate( G,f,k );
%
%   Input parameters
%       G   : Graph
%       f   : Signal (column)
%       k   : Indices of frequencies (int)
%   Output parameters
%       fm  : Modulated signal
%
%   This function modulate the column vector *f* onto the node i. If f is a
%   matrix, the modulation will be applicated to each column.
%

% Author: Nathanael Perraudin
% Date  : 09.12.2013

nt = size(f,2);

fm = sqrt(G.N)*repmat(f,1,nt).*repmat(G.U(:,k+1),1,nt);

end

repo: gspbox


text:
function gqu_sparse_test ( )

%*****************************************************************************80
%
%% GQU_SPARSE_TEST uses the GQU function to build a sparse grid.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 April 2014
%
%  Author:
%
%    Original MATLAB version by Florian Heiss, Viktor Winschel.
%
%  Local parameters:
%
%    Local, integer D, the spatial dimension.
%
%    Local, integer MAXK, the maximum level to check.
%
  d = 6;
  maxk = 10;
  func = 'prod ( exp(-(x/2).^2/2)/2/sqrt(2*pi), 2 ) ';
  trueval = fu_integral ( d );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'GQU_SPARSE_TEST:\n' );
  fprintf ( 1, '  GQU sparse grid:\n' );
  fprintf ( 1, '  Sparse Gaussian unweighted quadrature over [0,1].\n' );
  fprintf ( 1, '  Exact integral is %g\n', trueval );
  fprintf ( 1, '\n' );
  fprintf ( 1, '   D  Level   Nodes    SG error    MC error\n' );
  fprintf ( 1, '\n' );

  for k = 1 : maxk
%
%  Compute sparse grid estimate.
%
    [ x, w ] = nwspgr ( 'gqu', d, k );
    fx = eval ( func );
    SGappr = w' * fx;
    SGerror = sqrt ( ( SGappr - trueval ) .^ 2 ) / trueval;
%
%  Average 1000 Monte Carlo estimates.
%
    numnodes = length ( w );
    sim = zeros ( 1000, 1 );
    for r = 1 : 1000
      x = rand ( numnodes, d );
      fx = eval ( func );
      sim(r) = mean ( fx );
    end
    simerror = sqrt ( mean ( ( sim - trueval ) .^ 2 ) ) / trueval;

    fprintf( '  %2d     %2d  %6d  %10.5g  %10.5g\n', ...
      d, k, numnodes, SGerror, simerror )

  end

  return
end

repo: jburkardt-m


text:
% SDP relaxation for structured total least squares (STLS)
%
% Let k,m,n integers (m<=n)
% Let SS: R^k -> R^{m x n} be an affine map
% Let u1 be a vector in R^k
% Consider the optimization problem
%
% min_u     |u - u1|^2
% s.t.      SS(u) is rank deficient
%
% More generally, the cost can be a quadratic function
%           (u - u1)' * W * (u - u1)
% for some symmetric weight matrix W of size k x k.
% If W is diagonal and the entry W_ii is equal to zero then
% the respective u1(i) is ignored (as if it was missing).
%
% The function also accepts problems with complex data.
% They are internally converted to real valued problems.
%
% Input:
% S - matrix of size (k+1)m x n describing the affine map SS
% u1 - a vector of length k, possibly containing 'nan' entries
%
% Optional inputs:
% W - weight matrix of size k x k (default W_ii=1 except for 'nan' entries)
% solver - SDP solver (default 'sdpt3')
% quiet - set to false to display information (default true)
%
% Output:
% opt - optimal value of SDP relaxation
% u - the minimizer of the problem
% U - the matrix SS(u)
% z - a vector in the left kernel of SS(u)
% X - the PSD matrix (relaxation is exact if rank(X)=1)
%
% Usage:
% [opt,u,U,z,X] = sdp_stls(S,u1,W,solver,quiet)
% The optional arguments can be omitted or set to empty values.
% For instance: sdp_stls(S,u1,[],solver)

function [opt,u,U,z,X] = sdp_stls(S,u1,W,solver,quiet)

narginchk(2,5);
if nargin<3||isempty(W); W = diag(~isnan(u1)); end
if nargin<4||isempty(solver); solver = 'sdpt3'; end
if nargin<5||isempty(quiet); quiet = false; end

k = length(u1);
n = size(S,2);
m = size(S,1)/(k+1);
if floor(m)~=m; error('mismatch in dimensions'); end
u1(isnan(u1)) = 0;
u1 = reshape(u1,[k,1]);

iscomplex = ~(isreal(S) && isreal(u1) && isreal(W));
if iscomplex; [k,m,n,S,u1,W] = data2real(k,m,n,S,u1,W); end
W = .5*(W+W');

% fprintf('PSD matrix size: %d\n',(k+1)*m)

g = [eye(k) -u1];
G0 = g'*W*g;
E0 = sparse(k+1,k+1,1);

Im = eye(m);
G = kron(G0,Im);
E = kron(E0,Im);

[opt, X] = primal_cvx(k,m,n,S,G,E,solver,quiet);

[~,E] = eig(full(X));
r=nnz(diag(E)>1e-4);
if (~iscomplex && r>1) || (iscomplex && r>2)
    warning('sdp might not be exact');
end

J0 = k*m+1:(k+1)*m;
z = recoverSol(X(J0,J0));
u = zeros(k,1);
for i = 1:k
    Ji = (i-1)*m+1:i*m;
    u(i) = trace(X(J0,Ji));
end
U = applyAffineMap(S,u);

if iscomplex
    [u,U,z] = data2complex(k,m,n,u,U,z);
end

function [k,m,n,S,u1,W] = data2real(k,m,n,S,u1,W)
toReal = @(a) [real(a) -imag(a); imag(a) real(a)];
u1 = [real(u1); imag(u1)];
if norm(W-W')>1e-4; warning('W forced to be hermitian'); end
W = toReal(W);
SS = mat2cell(S,m*ones(k+1,1),n);
B = SS(1:k); A = SS(k+1);
iB = cellfun(@(x) 1i*x, B,'Unif',0);
SS = cellfun(toReal, [B; iB; A],'Unif',0);
S = cell2mat(SS);
k = 2*k; m = 2*m; n = 2*n;

function [u,U,z] = data2complex(k,m,n,u,U,z)
k = k/2; m = m/2; n = n/2;
u = u(1:k) + 1i * u(k+1:2*k);
U = U(1:m,1:n) + 1i * U(m+1:2*m,1:n);
z = z(1:m) + 1i * z(m+1:2*m);

% dual sdp
function [opt, X] = dual_cvx(k,m,n,S,G,E,solver,quiet)
N = (k+1)*m;
SC = cell(n,1);

if quiet
cvx_begin sdp quiet
else
cvx_begin sdp
end
    cvx_solver(solver)
    variable t(1,1);
    variable Cvec(N,n)
    dual variable X
    for i=1:n
        Ci = Cvec(:,i);
        si = S(:,i);
        SC{i} = si*Ci';
    end
    sSC = sum(cat(3,SC{:}),3);
    sSC = blksym(k+1,m,sSC);
    maximize(t);
    Q = G - t*E + sSC;
    X: Q >= 0;
cvx_end

opt = cvx_optval;

% primal sdp
function [opt, X] = primal_cvx(k,m,n,S,G,E,solver,quiet)
N = (k+1)*m;

e = speye(N);
SC = zeros(N,N,n*N);
for i=1:n
    for j=1:N
        SC(:,:,(i-1)*N+j) = S(:,i)*e(j,:);
    end
end
SC = blksym(k+1,m,SC);
A = zeros(n*N,N*(N+1)/2);
for l=1:n*N
    A(l,:) = smat2vec(SC(:,:,l));
end
A = sparse(A);

if quiet
cvx_begin sdp quiet
else
cvx_begin sdp
end
    cvx_solver(solver)
    variable X(N,N) symmetric
    dual variable Q
    y = smat2vec(X);
    minimize(smat2vec(G)'*y);
    smat2vec(E)'*y == 1;
    A*y == 0;
    Q: X >= 0;
cvx_end

opt = cvx_optval;

% recover minimizer from moment matrix
function [x,e] = recoverSol(X)
N = size(X,1);
if any(isnan(X))
    x = nan(N,1);
    e = inf;
else
    [V,E] = eig(full(X));
    e=diag(E);
    x = sqrt(e(N))*V(:,N);
    e = e(N-1);
end

% vector to symmetric matrix
function M = vec2smat(v)

N = length(v);
k = (-1+sqrt(1+8*N))/2;

M = repmat(0*v(1:k),[1,k]);
I = triu(true(k,k),0);
I2 = triu(true(k,k),1);
M(I) = v/2;
M(I2) = M(I2)*sqrt(2);
M = M + M.';

% symmetric matrix to vector
function v = smat2vec(M)

k = size(M,1);
I = triu(true(k,k),0);
I2 = triu(true(k,k),1);
M(I2) = M(I2)*sqrt(2);
v = M(I);

% Block symmetrization of a matrix
% Input: mk x mk matrix A
% Output: symmetric matrix As such that
%         each m x m block is also symmetric (there are n^2 blocks)
function As = blksym(n,m,A)

As = mysym(A);
if isempty(A); return; end

for i=1:n
    for j=1:n
        I = m*(i-1) + (1:m);
        J = m*(j-1) + (1:m);
        As(I,J,:) = mysym(As(I,J,:));
    end
end

% symmetrize matrix
function As = mysym(A)
At = permute(A,[[2,1],3:ndims(A)]);
As = .5*(A+At);

author: MIT-SPARK

repo: CertifiablyRobustPerception


%% BVAR tutorial: Minnesota prior and heteroskedasticity weights
% Author:   Filippo Ferroni and  Fabio Canova
% Date:     09/14/2020, revised 16/12/2020

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1) baseline  forecasting
% 2) forecasting  with  heteroskedastic  weights
% 3) forecasting  with  optimal  heteroskedastic  weights
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

close all; clc; clear;

addpath ../../cmintools/
addpath ../../bvartools/

%% read the data
if exist('DataCovid.mat','file') == 2
    load DataCovid
else
    url = 'https://fred.stlouisfed.org/';
    c = fred(url);
    series    = {'PAYEMS','UNRATE','PCE','INDPRO','CPIAUCSL','PCEPILFE'};
    startdate = datenum('12/01/1988', 'mm/dd/yyyy');
    enddate   = datenum('07/01/2020', 'mm/dd/yyyy');
    time      = 1988+11/12 : 1/12 : 2020+6/12;

    y=zeros(380,size(series,2)); y0=y; y1=y;
    for kk=1:size(series,2)
        tmp        = fetch(c,series{kk},startdate,enddate);
        if strcmp('UNRATE',series{kk})==1
            y(:,kk)    = tmp.Data(:,2);
            y0(:,kk)    = tmp.Data(:,2);
            y1(:,kk)  =(y0(:,kk));
        else
            y0(:,kk)    = tmp.Data(:,2);
            y1(:,kk)    = log(y0(:,kk));
            % rebase the log variable, transform to a index (2020m1=100)
            y(:,kk)= y1(:,kk) * 100 / y1(time==2019,kk);
        end
    end
    save DataCovid y y1 y0 time
end

% y, y0, y1 are  in DataCovid.

%% Forecast post COVID-19 (July 2020 - last insample data). Minnesota
lags                 = 13;
options.priors.name  = 'Minnesota';
options.fhor         = 24;
options.K            = 1000;
bvar0                = bvar_(y,lags,options);

options.nplots       = [2 3];
% start of the forecast plot - default first date in-sample data
options.time_start   = 2019;
% Titles for subplot
options.varnames = {'PAYMS','UNRATE','PCE','INDPRO','CPIAUCSL','PCEPILFE'};
% multiple credible set - default .68
options.conf_sig_2   = 0.9;
plot_frcst_(bvar0.forecasts.with_shocks,y,time,options)
pause;

% Forecast post COVID-19 (July 2020 - last insample data).
% Minnesota + Heteroskedasticity weights

lags = 13;
tstar  = find(time==2020) + 2; % pick march 2020
% scale the variables by factor <1 in the periods that
% characterize the COVID-19 induced recession
st                   = ones(size(y,1),1);
st(tstar:tstar+2 ,:) = [10 10 10]; % March, April, May
st(1:lags)           = [];

options.heterosked_weights = st;
options.priors.name  = 'Minnesota';
options.fhor         = 24;
bvar1                = bvar_(y,lags,options);

% adding  forecast  without  weights
options.add_frcst = [y; mean(bvar0.forecasts.no_shocks,3)];
plot_frcst_(bvar1.forecasts.with_shocks,y,time,options)
pause;

% Forecast  with  optimal  weights

hyperpara(1)    = 3;		  % tau
hyperpara(2)    = 0.5;		  % decay
hyperpara(3)    = 1;		  % lambda
hyperpara(4)    = 1;		  % mu
hyperpara(5)    = 2;		  % omega
hyperpara(6)    = 2; % s0: scale march 2020
hyperpara(7)    = 2; % s1: scale april 2020
hyperpara(8)    = 2; % s2: scale may   2020
% setting the options
options.index_est          = [1 6:8]; % hyper-parameter over which maximize
options.objective_function = 'bvar_opt_heterosked';
options.tstar              = find(time==2020) + 2; %march 2020
[postmode,logmlike,HH] = bvar_max_hyper(hyperpara,y,lags,options);

disp('weights for 2020:3, 2020:4, 2020:5')
disp(postmode(2:4))
heterosked_esse             = postmode(2:end); % s0, s1, s2
esse                        = ones(size(y,1),1);
esse(options.tstar : options.tstar+2 ,:)   = heterosked_esse;
esse(1:lags)                = [];

options.heterosked_weights = esse;
options.minn_prior_tau = postmode(1);
bvar2                 = bvar_(y,lags,options);
plot_frcst_(bvar2.forecasts.with_shocks,y,time,options)