text:
classdef prtUiDataSetClassExploreWidgetTabPlotOptions < prtUiDataSetClassExploreWidgetTab

    properties

        titleStr = 'Plot Options';

        possibleMarkers = {'o','s','^','d','v','x','.','<','>','p','h'};

        dataCell
        handles
    end

    methods
        function self = prtUiDataSetClassExploreWidgetTabPlotOptions(varargin)
            self = self@prtUiDataSetClassExploreWidgetTab(varargin{:});

            init(self);
        end

        function init(self)

            self.handles.table = uitable('parent',self.managedHandle,...
                'units','normalized','position',[0 0 1 1]); % Dummy position

            %self.handles.jScrollPane = findjobj(self.handles.table);
            %self.handles.jTable = self.handles.jScrollPane.getViewport.getView;

            % uY, visible, legend string, shape, color, edgeColor, markerSize, edgeSize

            initData(self);

            columnFormat = {'char', 'logical', self.possibleMarkers,'char','char','numeric','numeric'};
            columnHeadings = {'String','Visible','Shape','Color','Edge Color',' Size','Edge Size'};

            self.handles.uicontextmenu = uicontextmenu;
            self.handles.uicontextmenuReset = uimenu(self.handles.uicontextmenu, 'Label', 'reset', 'Callback', @(h,e)self.reset);

            % Set the actual data in the cell
            set(self.handles.table,...
                'columnName',columnHeadings,...
                'columnFormat',columnFormat,...
                'RearrangeableColumns','on',...
                'RowName',{},...cat(1,num2cell(self.widget.plotManager.dataSet.uniqueClasses(:)),{nan}),...
                'TooltipString','Plotting Options',...
                'ColumnEditable',true(1,size(self.dataCell,2)),...
                'CellEditCallback',@(h,e)self.tableCellEditCallback(h,e),...
                'UIContextMenu',self.handles.uicontextmenu);

            %javaaddpath(fullfile(prtRoot,'+prtExternal','+ColorCell'));

            %self.handles.jTable.setModel(javax.swing.table.DefaultTableModel(self.dataCell,columnHeadings))
            %self.handles.jTable.getColumnModel.getColumn(3).setCellRenderer(ColorCellRenderer);
            %self.handles.jTable.getColumnModel.getColumn(3).setCellEditor(ColorCellEditor);

            %self.handles.jTable.getColumnModel.getColumn(4).setCellRenderer(ColorCellRenderer);
            %self.handles.jTable.getColumnModel.getColumn(4).setCellEditor(ColorCellEditor);

            % Set the column widths to take up the whole area
            %self.handles.jTable.setAutoResizeMode(self.handles.jTable.AUTO_RESIZE_SUBSEQUENT_COLUMNS);
        end

        function reset(self)
            self.widget.plotManager.setDataSetDependentProperties();
            initData(self);
        end

        function initData(self)

            uYs = self.widget.plotManager.dataSet.uniqueClasses;

            data = cell(length(uYs)+1, 7);
            data(:,1) = cat(1,self.widget.plotManager.legendStrings,{self.widget.plotManager.legendStringUnlabeled});
            data(:,2) = num2cell(logical(cat(1, self.widget.plotManager.isVisible(:), self.widget.plotManager.isVisibleUnlabeled)));
            data(:,3) = cat(1,num2cell(self.widget.plotManager.markers(:)),{self.widget.plotManager.markerUnlabeled});

            for iClass = 1:length(uYs)

                cColor = max(min(round(self.widget.plotManager.colors(iClass,:)*255),255),0);
                cEdgeColor = max(min(round(self.widget.plotManager.edgeColors(iClass,:)*255),255),0);

                data{iClass,4} = sprintf('%d,%d,%d',cColor(1), cColor(2), cColor(3));
                data{iClass,5} = sprintf('%d,%d,%d',cEdgeColor(1), cEdgeColor(2), cEdgeColor(3));
            end
            cColor = max(min(round(self.widget.plotManager.colorUnlabeled*255),255),0);
            cEdgeColor = max(min(round(self.widget.plotManager.edgeColorUnlabeled*255),255),0);

            data{end,4} = sprintf('%d,%d,%d',cColor(1), cColor(2), cColor(3));
            data{end,5} = sprintf('%d,%d,%d',cEdgeColor(1), cEdgeColor(2), cEdgeColor(3));

            data(:,6) = num2cell(cat(1,self.widget.plotManager.markerSizes(:), self.widget.plotManager.markerSizeUnlabeled));
            data(:,7) = num2cell(cat(1,self.widget.plotManager.markerLineWidths(:), self.widget.plotManager.markerLineWidthUnlabeled));

            self.dataCell = data;

            set(self.handles.table,'data',self.dataCell);
        end
        function tableCellEditCallback(self, h, e)

            inds = e.Indices;

            switch inds(2)
                case 1
                    changeLegend(self, h, e);
                case 2
                    changeVisible(self, h, e);
                case 3
                    changeMarkers(self, h, e);
                case 4
                    changeColors(self, h, e);
                case 5
                    changeEdgeColors(self, h, e);
                case 6
                    changeMarkerSize(self, h, e);
                case 7
                    changeMarkerEdgeSize(self, h, e)
            end
        end

        function changeVisible(self, h, e)
            newValue = logical(e.NewData);

            cData = self.dataCell;
            if e.Indices(1) < size(cData,1)
                self.widget.plotManager.isVisible(e.Indices(1)) = newValue;
            else
                self.widget.plotManager.isVisibleUnlabeled = newValue;
            end

            self.widget.plotManager.legendRefresh()

            self.dataCell{e.Indices(1), e.Indices(2)} = newValue;
        end
        function changeLegend(self, h, e)
            if e.Indices(1) < size(self.dataCell,1)
                self.widget.plotManager.legendStrings{e.Indices(1)} = e.NewData;
            else
                self.widget.plotManager.legendStringUnlabeled = e.NewData;
            end

            self.dataCell{e.Indices(1), e.Indices(2)} = e.NewData;
        end
        function changeMarkers(self, h, e)
            if e.Indices(1) < size(self.dataCell,1)
                self.widget.plotManager.markers(e.Indices(1)) = e.NewData;
            else
                self.widget.plotManager.markerUnlabeled = e.NewData;
            end

            self.dataCell{e.Indices(1), e.Indices(2)} = e.NewData;
        end
        function changeColors(self, h, e)
            [err, color] = testColorValueString(self, e.NewData);
            if err
                set(self.handles.table,'data',self.dataCell); % Keep what we had
                return
            end

            if e.Indices(1) < size(self.dataCell,1)
                self.widget.plotManager.colors(e.Indices(1),:) = color;
            else
                self.widget.plotManager.colorUnlabeled = color;
            end

            self.dataCell{e.Indices(1), e.Indices(2)} = e.NewData;
        end
        function changeEdgeColors(self, h, e)
            [err, color] = testColorValueString(self, e.NewData);
            if err
                set(self.handles.table,'data',self.dataCell); % Keep what we had
                return
            end

            if e.Indices(1) < size(self.dataCell,1)
                self.widget.plotManager.edgeColors(e.Indices(1),:) = color;
            else
                self.widget.plotManager.edgeColorUnlabeled = color;
            end

            self.dataCell{e.Indices(1), e.Indices(2)} = e.NewData;
        end
        function [err, color] = testColorValueString(self, val)

            err = true;
            color = [];
            try
                commaLocs = find(val==',');
                if length(commaLocs) ~= 2
                    return
                end

                color = zeros(1,3);
                color(1) = str2double(val(1:commaLocs(1)-1));
                color(2) = str2double(val((commaLocs(1)+1):(commaLocs(2)-1)));
                color(3) = str2double(val((commaLocs(2)+1):end));

                color = color/255;

                if any(color < 0) || any(color > 1) || any(isnan(color))
                    color = [];
                    return
                end
                err = false;
            catch
                err = true;
                color = [];
            end
        end

        function changeMarkerSize(self, h, e)
            try
                assert(e.NewData > 0);
                if e.Indices(1) < size(self.dataCell,1)
                    self.widget.plotManager.markerSizes(e.Indices(1)) = e.NewData;
                else
                    self.widget.plotManager.markerSizeUnlabeled = e.NewData;
                end

                self.dataCell{e.Indices(1), e.Indices(2)} = e.NewData;

            catch %#ok<CTCH>
                set(self.handles.table,'data',self.dataCell);
            end
        end
        function changeMarkerEdgeSize(self, h, e)
            try
                assert(e.NewData > 0);
                if e.Indices(1) < size(self.dataCell,1)
                    self.widget.plotManager.markerLineWidths(e.Indices(1)) = e.NewData;
                else
                    self.widget.plotManager.markerLineWidthUnlabeled = e.NewData;
                end

                self.dataCell{e.Indices(1), e.Indices(2)} = e.NewData;

            catch %#ok<CTCH>
                set(self.handles.table,'data',self.dataCell);
            end
        end
    end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
  function groups = group2d(nx, ny, ngv, mask, chat)
%|function groups = group2d(nx, ny, ngv, mask, chat)
%| build groups for GCA for 2d case
%| ngv [2]	[ngx ngy]
%|
%| Copyright Mar. 1999, Jeff Fessler

if 1 ~= exist('chat'), chat = 1; end
if 1 ~= exist('mask') || isequal(mask, 1)
	mask = ones(nx,ny);
end

	ngx = ngv(1);
	ngy = ngv(2);

	np = nx*ny;

	ng = ngx*ngy;
	groups	= logical(zeros(np, ng));
	for igx = 1:ngx
		for igy = 1:ngy
			ig = igx + (igy-1)*ngx;
			gx = [igx:ngx:nx]';
			gy = [igy:ngy:ny]';
			ggx = gx * ones(1,length(gy));
			ggy = ones(length(gx),1) * gy';
			gg = ggx + (ggy-1) * nx;
			gg = gg(:);
			groups(gg,ig) = ones(size(gg));
		end
	end

	groups = groups(find(mask(:)),:);

	% look at groups
	if chat
		im(groups', 'groups')
		for ig=1:ng
			prompt
			im(embed(groups(:,ig), mask), sprintf('group %d', ig))
		end
		prompt
		im(embed(groups * [1:ng]',mask), 'All Groups')
	end

	if any(sum(groups')) > 1, error 'overlap groups', end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function test_bug3417

% WALLTIME 00:20:00
% MEM 3gb
% DEPENDENCY

% ... works fine on his computer, which has a 2012 version of FieldTrip and Matlab 2013a.
% But when we try to run either a later version of Fieldtrip (e.g. 2015) or a later version
% of Matlab (e.g. 2017b), we get the following bug:
%
% Error using findcluster (line 50)
% invalid dimension of spatdimneighbstructmat

% figure 4

%% SL - one sample t

load(dccnpath('/home/common/matlab/fieldtrip/data/test/bug3417.mat'))

freqs=14; % or 21
times=81;

cond1=permute(SL_Fig4,[3,1,2]);
cond1=cond1(:,1:freqs,:);

cond2=cond1;
cond2(:,:,:)=0;

data1.globspctrm=cond1;
data2.globspctrm=cond2;

[subj,~,freq] = size(data1.globspctrm);
nsubj = subj;

cfg.statistic = 'depsamplesT';
cfg.design = zeros(2,2*nsubj);
cfg.design(1,:)=[ones(1,nsubj) ones(1,nsubj)+1];
cfg.design(2,:)=[1:nsubj 1:nsubj];
cfg.ivar = 1;
cfg.uvar = 2;

data1.label = {'Positive'};
data2.label = {'Negative'};
data1.dimord = 'rpt_freq_time'; % does clustering in 2 dimensions (freq, time)
data2.dimord = 'rpt_freq_time';
data1.freq = [1:freqs];
data2.freq = [1:freqs];
data1.time =[1:times];
data2.time= [1:times];

% fix the invalid data structure, see ft_datatype_freq
data1.dimord = 'rpt_chan_freq_time';
data2.dimord = 'rpt_chan_freq_time';
data1.globspctrm = reshape(data1.globspctrm, [nsubj 1 freqs times]);
data2.globspctrm = reshape(data2.globspctrm, [nsubj 1 freqs times]);
data1.label = {'TheSame'};
data2.label = {'TheSame'};

cfg.parameter = 'globspctrm'; % 'avg', 'cohspctrm'
cfg.method = 'montecarlo'; %  'montecarlo',   'analytic' , 'stats', 'crossvalidate'

cfg.neighbours = [];
cfg.channel = 'all';
cfg.correctm = 'cluster'; % 'no', 'max', 'cluster', 'bonferroni', 'holm', 'hochberg', 'fdr' (default = 'no')

cfg.clusterstatistc = 'wcm'; % 'wcm', 'maxsum', 'maxsize'
cfg.clusterthreshold = 'parametric'; %  'parametric', 'nonparametric_individual', 'nonparametric_common'
cfg.numrandomization = 1000;
cfg.clusteralpha = 0.05;
cfg.tail = 0;

stats = ft_freqstatistics(cfg,data1,data2);

figure;

subplot(2,3,1)
tres=squeeze(mean(cond1,1));
imagesc(tres)
colorbar
title('Adult SL: Mean')
ylabel('Frequencies')
xlabel('Time')

subplot(2,3,4)
tres=squeeze(stats.stat)
imagesc(tres)
colorbar
title('SL: t matrix')
ylabel('Frequencies')
xlabel('Time')

subplot(2,3,5)
Positive_cluster=squeeze(stats.posclusterslabelmat);
Positive_cluster(Positive_cluster > 1) = 0; % keep just cluster 1
imagesc(Positive_cluster)
title(['Positive cluster: p=' num2str(stats.posclusters(1, 1).prob)])
colorbar
ylabel('Frequencies')
xlabel('Time')

subplot(2,3,6)
Negative_cluster=squeeze(stats.negclusterslabelmat);
Negative_cluster(Negative_cluster > 1) = 0; % keep just cluster 1
imagesc(Negative_cluster)
title(['Negative cluster: p=' num2str(stats.negclusters(1, 1).prob)])
colorbar
ylabel('Frequencies')
xlabel('Time')

%% JA - one sample t

load(dccnpath('/home/common/matlab/fieldtrip/data/test/bug3417.mat'))

freqs=14; % or 21
times=81;

cond1=permute(JA_Fig4,[3,1,2]);
cond1=cond1(:,1:freqs,:);

data1.globspctrm=cond1;
data2.globspctrm=cond2;

[subj,~,freq] = size(data1.globspctrm);
nsubj = subj;

cfg.statistic = 'depsamplesT';
cfg.design = zeros(2,2*nsubj);
cfg.design(1,:)=[ones(1,nsubj) ones(1,nsubj)+1];
cfg.design(2,:)=[1:nsubj 1:nsubj];
cfg.ivar = 1;
cfg.uvar = 2;

data1.label = {'Positive'};
data2.label = {'Negative'};
data1.dimord = 'rpt_freq_time'; % does clustering in 2 dimensions (freq, time)
data2.dimord = 'rpt_freq_time';
data1.freq = [1:freqs];
data2.freq = [1:freqs];
data1.time =[1:times];
data2.time= [1:times];

% fix the invalid data structure, see ft_datatype_freq
data1.dimord = 'rpt_chan_freq_time';
data2.dimord = 'rpt_chan_freq_time';
data1.globspctrm = reshape(data1.globspctrm, [nsubj 1 freqs times]);
data2.globspctrm = reshape(data2.globspctrm, [nsubj 1 freqs times]);
data1.label = {'TheSame'};
data2.label = {'TheSame'};

cfg.parameter = 'globspctrm'; % 'avg', 'cohspctrm'
cfg.method = 'montecarlo'; %  'montecarlo',   'analytic' , 'stats', 'crossvalidate'

cfg.neighbours = [];
cfg.channel = 'all';
cfg.correctm = 'cluster'; % 'no', 'max', 'cluster', 'bonferroni', 'holm', 'hochberg', 'fdr' (default = 'no')

cfg.clusterstatistc = 'wcm'; % 'wcm', 'maxsum', 'maxsize'
cfg.clusterthreshold = 'parametric'; %  'parametric', 'nonparametric_individual', 'nonparametric_common'
cfg.numrandomization = 1000;
cfg.clusteralpha = 0.05;
cfg.tail = 0;

stats = ft_freqstatistics(cfg,data1,data2)

figure;

subplot(2,3,1)
tres=squeeze(mean(cond1,1));
imagesc(tres)
colorbar
title('Adult SL: Mean')
ylabel('Frequencies')
xlabel('Time')

subplot(2,3,4)
tres=squeeze(stats.stat)
imagesc(tres)
colorbar
title('SL: t matrix')
ylabel('Frequencies')
xlabel('Time')

subplot(2,3,5)
Positive_cluster=squeeze(stats.posclusterslabelmat);
Positive_cluster(Positive_cluster > 1) = 0; % keep just cluster 1
imagesc(Positive_cluster)
title(['Positive cluster: p=' num2str(stats.posclusters(1, 1).prob)])
colorbar
ylabel('Frequencies')
xlabel('Time')

subplot(2,3,6)
Negative_cluster=squeeze(stats.negclusterslabelmat);
Negative_cluster(Negative_cluster > 1) = 0; % keep just cluster 1
imagesc(Negative_cluster)
title(['Negative cluster: p=' num2str(stats.negclusters(1, 1).prob)])
colorbar
ylabel('Frequencies')
xlabel('Time')

%% SJ vs. AJ

load(dccnpath('/home/common/matlab/fieldtrip/data/test/bug3417.mat'))

freqs=14; % or 21
times=81;

cond1=permute(SL_Fig4,[3,1,2]);
cond2=permute(JA_Fig4,[3,1,2]);

cond1=cond1(:,1:freqs,:);
cond2=cond2(:,1:freqs,:);

data1.globspctrm=cond1;
data2.globspctrm=cond2;

[subj,~,freq] = size(data1.globspctrm);
nsubj = subj;

cfg.statistic = 'depsamplesT';
cfg.design = zeros(2,2*nsubj);
cfg.design(1,:)=[ones(1,nsubj) ones(1,nsubj)+1];
cfg.design(2,:)=[1:nsubj 1:nsubj];
cfg.ivar = 1;
cfg.uvar = 2;

data1.label = {'Positive'};
data2.label = {'Negative'};
data1.dimord = 'rpt_freq_time'; % does clustering in 2 dimensions (freq, time)
data2.dimord = 'rpt_freq_time';
data1.freq = [1:freqs];
data2.freq = [1:freqs];
data1.time =[1:times];
data2.time= [1:times];

% fix the invalid data structure, see ft_datatype_freq
data1.dimord = 'rpt_chan_freq_time';
data2.dimord = 'rpt_chan_freq_time';
data1.globspctrm = reshape(data1.globspctrm, [nsubj 1 freqs times]);
data2.globspctrm = reshape(data2.globspctrm, [nsubj 1 freqs times]);
data1.label = {'TheSame'};
data2.label = {'TheSame'};

cfg.parameter = 'globspctrm'; % 'avg', 'cohspctrm'
cfg.method = 'montecarlo'; %  'montecarlo',   'analytic' , 'stats', 'crossvalidate'

cfg.neighbours = [];
cfg.channel = 'all';
cfg.correctm = 'cluster'; % 'no', 'max', 'cluster', 'bonferroni', 'holm', 'hochberg', 'fdr' (default = 'no')

cfg.clusterstatistc = 'wcm'; % 'wcm', 'maxsum', 'maxsize'
cfg.clusterthreshold = 'parametric'; %  'parametric', 'nonparametric_individual', 'nonparametric_common'
cfg.numrandomization = 1000;
cfg.clusteralpha = 0.05;
cfg.tail = 0;

stats = ft_freqstatistics(cfg,data1,data2)

figure;

subplot(2,3,1)
tres=squeeze(mean(cond1,1));
imagesc(tres)
colorbar
title('Adult SL: Mean')
ylabel('Frequencies')
xlabel('Time')

subplot(2,3,2)
tres=squeeze(mean(cond2,1));
imagesc(tres)
colorbar
title('Adult JA: Mean')
ylabel('Frequencies')
xlabel('Time')

subplot(2,3,4)
tres=squeeze(stats.stat)
imagesc(tres)
colorbar
title('SL vs. JA: t matrix')
ylabel('Frequencies')
xlabel('Time')

subplot(2,3,5)
Positive_cluster=squeeze(stats.posclusterslabelmat);
Positive_cluster(Positive_cluster > 1) = 0; % keep just cluster 1
imagesc(Positive_cluster)
title(['SL > JA cluster: p=' num2str(stats.posclusters(1, 1).prob)])
colorbar
ylabel('Frequencies')
xlabel('Time')

subplot(2,3,6)
Negative_cluster=squeeze(stats.negclusterslabelmat);
Negative_cluster(Negative_cluster > 1) = 0; % keep just cluster 1
imagesc(Negative_cluster)
title(['SL < JA cluster: p=' num2str(stats.negclusters(1, 1).prob)])
colorbar
ylabel('Frequencies')
xlabel('Time')

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [pesq_mos]= pesq(ref_wav, deg_wav)

% ----------------------------------------------------------------------
%            PESQ objective speech quality measure
%
%   This function implements the PESQ measure based on the ITU standard
%   P.862 [1].
%
%
%   Usage:  pval=pesq(cleanFile.wav, enhancedFile.wav)
%
%         cleanFile.wav - clean input file in .wav format
%         enhancedFile  - enhanced output file in .wav format
%         pval          - PESQ value
%
%    Note that the PESQ routine only supports sampling rates of 8 kHz and
%    16 kHz [1]
%
%  Example call:  pval = pesq ('sp04.wav','enhanced.wav')
%
%
%  References:
%   [1] ITU (2000). Perceptual evaluation of speech quality (PESQ), and
%       objective method for end-to-end speech quality assessment of
%       narrowband telephone networks and speech codecs. ITU-T
%       Recommendation P. 862
%
%   Authors: Yi Hu and Philipos C. Loizou
%
%
% Copyright (c) 2006 by Philipos C. Loizou
% $Revision: 0.0 $  $Date: 10/09/2006 $
% ----------------------------------------------------------------------

if nargin<2
    fprintf('Usage: [pesq_mos]=pesq(cleanfile.wav,enhanced.wav) \n');
    return;
end;

global Downsample DATAPADDING_MSECS SEARCHBUFFER Fs WHOLE_SIGNAL
global Align_Nfft Window

[ref_data,sampling_rate]= audioread( ref_wav);
% if sampling_rate~=8000 & sampling_rate~=16000
%     error('Sampling frequency needs to be either 8000 or 16000 Hz');
% end

setup_global( sampling_rate);

% Window= hann( Align_Nfft, 'periodic'); %Hanning window
% Window= Window';
TWOPI= 6.28318530717959;
%for count = 0: Align_Nfft- 1
%    Window(1+ count) = 0.5 * (1.0 - cos((TWOPI * count) / Align_Nfft));
%end

count=0:Align_Nfft- 1;
Window= 0.5 * (1.0 - cos((TWOPI * count) / Align_Nfft));

ref_data= ref_data';
ref_data= ref_data* 32768;
ref_Nsamples= length( ref_data)+ 2* SEARCHBUFFER* Downsample;
ref_data= [zeros( 1, SEARCHBUFFER* Downsample), ref_data, ...
    zeros( 1, DATAPADDING_MSECS* (Fs/ 1000)+ SEARCHBUFFER* Downsample)];

deg_data= audioread( deg_wav);
deg_data= deg_data';
deg_data= deg_data* 32768;
deg_Nsamples= length( deg_data)+ 2* SEARCHBUFFER* Downsample;
deg_data= [zeros( 1, SEARCHBUFFER* Downsample), deg_data, ...
    zeros( 1, DATAPADDING_MSECS* (Fs/ 1000)+ SEARCHBUFFER* Downsample)];

maxNsamples= max( ref_Nsamples, deg_Nsamples);

ref_data= fix_power_level( ref_data, ref_Nsamples, maxNsamples);
deg_data= fix_power_level( deg_data, deg_Nsamples, maxNsamples);

standard_IRS_filter_dB= [0, -200; 50, -40; 100, -20; 125, -12; 160, -6; 200, 0;...
    250, 4; 300, 6; 350, 8; 400, 10; 500, 11; 600, 12; 700, 12; 800, 12;...
    1000, 12; 1300, 12; 1600, 12; 2000, 12; 2500, 12; 3000, 12; 3250, 12;...
    3500, 4; 4000, -200; 5000, -200; 6300, -200; 8000, -200];

ref_data= apply_filter( ref_data, ref_Nsamples, standard_IRS_filter_dB);
deg_data= apply_filter( deg_data, deg_Nsamples, standard_IRS_filter_dB);
%

% for later use in psychoacoustical model
model_ref= ref_data;
model_deg= deg_data;

[ref_data, deg_data]= input_filter( ref_data, ref_Nsamples, deg_data, ...
    deg_Nsamples);

[ref_VAD, ref_logVAD]= apply_VAD( ref_data, ref_Nsamples);
[deg_VAD, deg_logVAD]= apply_VAD( deg_data, deg_Nsamples);

crude_align (ref_logVAD, ref_Nsamples, deg_logVAD, deg_Nsamples,...
    WHOLE_SIGNAL);

utterance_locate (ref_data, ref_Nsamples, ref_VAD, ref_logVAD,...
    deg_data, deg_Nsamples, deg_VAD, deg_logVAD);

ref_data= model_ref;
deg_data= model_deg;

% make ref_data and deg_data equal length
if (ref_Nsamples< deg_Nsamples)
    newlen= deg_Nsamples+ DATAPADDING_MSECS* (Fs/ 1000);
    ref_data( newlen)= 0;
elseif (ref_Nsamples> deg_Nsamples)
    newlen= ref_Nsamples+ DATAPADDING_MSECS* (Fs/ 1000);
    deg_data( newlen)= 0;
end

pesq_mos= pesq_psychoacoustic_model (ref_data, ref_Nsamples, deg_data, ...
    deg_Nsamples );

author: vipchengrui

repo: traditional-speech-enhancement


text:
classdef CelestialBodySunRelStateDataCache < matlab.mixin.SetGet
    %CelestialBodySunRelStateDataCache Summary of this class goes here
    %   Detailed explanation goes here

    properties
        bodyInfo KSPTOT_BodyInfo

        times(:,1) double = [];
        stateVects(:,6) double = [];
        relFrame AbstractReferenceFrame

        gi %gridded interpolant
    end

    methods
        function obj = CelestialBodySunRelStateDataCache(bodyInfo)
            obj.bodyInfo = bodyInfo;
        end

        function setStateData(obj, times, stateVects, relFrame)
            obj.times = times;
            obj.stateVects = stateVects;
            obj.relFrame = relFrame;

            obj.gi = griddedInterpolant(obj.times, obj.stateVects, 'spline');
        end

        function [rVect, vVect] = getCachedBodyStateAtTime(obj, time)
            if(numel(obj.times) <= 1 || any(time > max(obj.times)) || any(time < min(obj.times)))
                maxUT = max(obj.times);
                boolMaxUT = maxUT - time < 0;
                worstMaxUTViolation = max(time(boolMaxUT));

                minUT = min(obj.times);
                boolMinUT = time - minUT < 0;
                worstMinUTViolation = min(time(boolMinUT));

                if(not(isempty(worstMaxUTViolation)) && not(isempty(worstMinUTViolation)))
                    str = sprintf('Adjust the minimum cache UT to below %0.3f sec and the maximum cache UT to above %0.3f sec.', worstMinUTViolation, worstMaxUTViolation);

                elseif(not(isempty(worstMaxUTViolation)))
                    str = sprintf('Adjust the maximum cache UT to above %0.3f sec.', worstMaxUTViolation);

                elseif(not(isempty(worstMinUTViolation)))
                    str = sprintf('Adjust the minimum cache UT to below %0.3f sec.', worstMinUTViolation);

                end

                msg = sprintf('At least one of the points in time queried are out of the bounds of the numerical integration (%0.3f sec to %0.3f sec).  %s', min(obj.times), max(obj.times), str);

                errordlg(msg, 'Body State Cache Bounds', 'modal');
                error(msg);
            else
                vq = obj.gi(time);

                rVect = vq(:,1:3)';
                vVect = vq(:,4:6)';
            end
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
% DEMVOWELS2 Model the vowels data with a 2-D FGPLVM using RBF kernel.

% FGPLVM

% Fix seeds
randn('seed', 1e5);
rand('seed', 1e5);

dataSetName = 'vowels';
experimentNo = 2;

% load data
[Y, lbls] = lvmLoadData(dataSetName);

% Set up model
options = fgplvmOptions('fitc');
options.numActive = 200;
latentDim = 2;
d = size(Y, 2);

model = fgplvmCreate(latentDim, d, Y, options);

% Optimise the model.
iters = 1000;
display = 1;

model = fgplvmOptimise(model, display, iters);

% Save the results.
capName = dataSetName;;
capName(1) = upper(capName(1));
save(['dem' capName num2str(experimentNo) '.mat'], 'model');

if exist('printDiagram') & printDiagram
  fgplvmPrintPlot(model, lbls, capName, experimentNo);
end

% Load the results and display dynamically.
fgplvmResultsDynamic(dataSetName, experimentNo, 'vector')

errors = fgplvmNearestNeighbour(model, lbls);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [u] = spm_uc_RF(a,df,STAT,R,n)
% Corrected critical height threshold at a specified significance level
% FORMAT [u] = spm_uc_RF(a,df,STAT,R,n)
% a     - critical probability - {alpha}
% df    - [df{interest} df{residuals}]
% STAT  - Statistical field
%         'Z' - Gaussian field
%         'T' - T field
%         'X' - Chi-squared field
%         'F' - F field
% R     - RESEL Count {defining search volume}
% n     - number of conjoint SPMs
%
% u     - critical height {corrected}
%
%__________________________________________________________________________
%
% spm_uc returns the corrected critical threshold at a specified significance
% level (a). If n > 1 a conjunction probability over the n values of the
% statistic is returned.
%__________________________________________________________________________
% Copyright (C) 1999-2012 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_uc_RF.m 4634 2012-02-01 19:01:09Z guillaume $

%-Find approximate value
%--------------------------------------------------------------------------
u  = spm_u((a/max(R))^(1/n),df,STAT);
du = 1e-6;

%-Approximate estimate using E{m}
%--------------------------------------------------------------------------
d  = 1;
while abs(d) > 1e-6
    [P, P, p] = spm_P_RF(1,0,u,df,STAT,R,n);
    [P, P, q] = spm_P_RF(1,0,u + du,df,STAT,R,n);
    d         = (a - p)/((q - p)/du);
    u         = u + d;
    if isinf(u), u=+Inf; return; end
end

%-Refined estimate using 1 - exp(-E{m})
%--------------------------------------------------------------------------
d  = 1;
while abs(d) > 1e-6
    p         = spm_P_RF(1,0,u,df,STAT,R,n);
    q         = spm_P_RF(1,0,u + du,df,STAT,R,n);
    d         = (a - p)/((q - p)/du);
    u         = u + d;
    if isinf(u), u=+Inf; return; end
end

author: spm


text:
function [tdiff, ac]  = timediff(clus_idx, look_ahead_time, clus, eqtimes)
    % TIMEDIFF calculates the time difference between the ith and jth event
    % works with variable eqtime from function CLUSTIME
    % gives the indices ac of the eqs not already related to cluster k1
    % eqtimes should be sorted!
    %
    % clus_idx : ith cluster (ci)
    % look_ahead : look-ahead time (tau)
    % clus: clusters (length of catalog)
    % eqtimes: datetimes for event catalog, in days  [did not use duration because of overhead]
    %
    % tdiff: is time between jth event and eqtimes(clus_idx)
    % ac: index of each event within the cluster

    comparetime = eqtimes(clus_idx);

    first_event = clus_idx + 1; % in cluster
    last_event = numel(eqtimes);
    max_elapsed = comparetime + look_ahead_time;

    if eqtimes(end) >= max_elapsed
        last_event = find(eqtimes(first_event : last_event) < max_elapsed, 1, 'last') + clus_idx;
    end

    tdiff = eqtimes(first_event : last_event) - comparetime;

    if first_event == last_event
        % no additional events were found.
        ac = [];
    else

        this_clusternum = clus(clus_idx);

        if this_clusternum == 0
            ac = first_event:last_event;
        else
            % indices of eqs not already related to this cluster
            ac = (find(clus(first_event : last_event) ~= this_clusternum)) + clus_idx;
        end
    end
    ac = ac(:);
end

author: CelsoReyes


text:
 function ob = abs(ob)
%function ob = abs(ob)

if isempty(ob.handle_abs)
	error 'no abs() method for this object'
end

ob = ob.handle_abs(ob);


text:
function chars = getCharSizes(axis)
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

axes(axis);
oldUnits = get(axis, 'units');
set(axis, 'units', 'pixels');
h = text(1, 1, 1, '||', 'margin', .01, 'interpreter', 'none', 'units', 'pixels');
extent = get(h, 'extent');
delete(h);
baseline = extent(3);

for i=0:255
    h = text(1, 1, 1, ['|' char(i) '|'], 'margin', .01, 'interpreter', 'none', 'units', 'pixels');
    extent = get(h, 'extent');
    output(i+1) = extent(3) - baseline;
    delete(h);
end
chars = output;
set(axis, 'units', oldUnits);

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function optimal = fibonacci(f,a,b,n)

% Program : fibonacci.m
%
% Purpose : find the optimal point of a given single
%           variable function within the given interval
%
% Author :  Aamir Alaud Din
% Date :    26.09.2013
%
% Inputs :  Four input arguments are required viz., function, initial and
%           final points of the interval, and function evaluation number n
%
% Syntax :  optimal = fibonacci(f,a,b)
%           where f is an inline function
%           e.g., f = inline('sin(x)','x');
%           a and b are initial and final points of the interval
%
% Example:  optimal = fibonacci(f,-3,5)
%
%           Fourth input argument is fibonacci number decided by the user.
%           Default fibonacci number is 11 which gives 10 iterations.
%           Number of iterations is one less then fibonacci number, i.e.,
%           number of iterations = fibonacci number - 1
%
% Syntax:   optimal = fibonacci(f,a,b,n)
%           f, a, and b are described above
%           n is the number of iterations decided by the user
% Example:  optimal = fibonacci(f,0,1,20)
%

if (nargin < 3)
    error('Too few input arguments');
    optimal = 'Fibonacci search is not possible with too few inputs';
    return;

elseif (nargin == 3)
    n = 11;
    fss = [1,1];
    for ii = 1:(n-1)
        fs = fss(end) + fss(end-1);
        fss = [fss,fs];
    end

    disp('Iteration          N          x1              x2             f(x1)           f(x2)');
    disp('=========        ===       ========        ========        =========       =========');
    Iter = 0;
    for jj = n:-1:2
        Iter = Iter + 1;
        L = b - a;
        x1 = a + ((fss(jj-1)/fss(jj+1))*L);
        x2 = b - ((fss(jj-1)/fss(jj+1))*L);
        if (f(x1) < f(x2))
            b = x2;
        elseif (f(x1) > f(x2))
            a = x1;
        elseif (f(x1) == f(x2))
            a = x1;
            b = x2 + 0.0001;
        end
        fprintf('%3d', Iter);
        fprintf('\t\t\t\t');
        fprintf('%4d', jj-1);
        fprintf('\t');
        fprintf('%11.4f', x1);
        fprintf('\t\t');
        fprintf('%11.4f', x2);
        fprintf('\t\t');
        fprintf('%12.4f', f(x1));
        fprintf('\t\t');
        fprintf('%12.4f', f(x2));
        fprintf('\n');
    end
    optimal = min(x1,x2);

elseif (nargin == 4)
    fss = [1,1];
    for ii = 1:(n-1)
        fs = fss(end) + fss(end-1);
        fss = [fss,fs];
    end

    disp('Iteration          N          x1              x2             f(x1)           f(x2)');
    disp('=========        ===       ========        ========        =========       =========');
    Iter = 0;
    for jj = n:-1:2
        Iter = Iter + 1;
        L = b - a;
        x1 = a + ((fss(jj-1)/fss(jj+1))*L);
        x2 = b - ((fss(jj-1)/fss(jj+1))*L);
        if (f(x1) < f(x2))
            b = x2;
        elseif (f(x1) > f(x2))
            a = x1;
        elseif (f(x1) == f(x2))
            a = x1;
            b = x2 + 0.0001;
        end
        fprintf('%3d', Iter);
        fprintf('\t\t\t\t');
        fprintf('%4d', jj-1);
        fprintf('\t');
        fprintf('%11.4f', x1);
        fprintf('\t\t');
        fprintf('%11.4f', x2);
        fprintf('\t\t');
        fprintf('%12.4f', f(x1));
        fprintf('\t\t');
        fprintf('%12.4f', f(x2));
        fprintf('\n');
    end
    optimal = min(x1,x2);

elseif (nargin > 4)
    error('Too many input arguments');
    optimal = 'Fibonacci search is not possible with too many inputs';
    return

end

author: Sable

repo: mcbench-benchmarks


text:
function [TrainX, TestX, TrainY, TestY, numSource, numTrain, numTarget, numTest] = HIDC_loadData(inputPath)
    load(inputPath);
    traindata{1,1} = TrainData;
    trainlabel{1,1} = TrainLabel;
    testdata{1,1} = TestData;
    testlabel{1,1} = TestLabel;
    numSource = size(traindata,2);                              % number of train file
    numTrain = zeros(1,numSource);
    TrainX = [];
    TrainY = [];
    for i = 1 : length(numTrain)
        numTrain(1,i) = size(traindata{1,i},2);
        TrainX = [TrainX traindata{1,i}];
        TrainY = [TrainY trainlabel{1,i}];
    end

    numTarget = size(testdata,2);                              % number of train file
    numTest = zeros(1,numTarget);
    TestX = [];
    TestY = [];
    for i = 1 : length(numTest)
        numTest(1,i) = size(testdata{1,i},2);
        TestX = [TestX testdata{1,i}];
        TestY = [TestY testlabel{1,i}];
    end
    TrainY = TrainY+1;
    TestY = TestY+1;
end

author: FuzhenZhuang

repo: Transfer-Learning-Toolkit


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
% This is a testing environment for the files in the folder kernel/data
% 1. Based on data/contents, a list of required files is generated and it
%    is verified, that all files are present; additional files are listed.
% 2. All c-files are compiled.
% 3. All files are executed.
% 4. Check administrative modul Trafo.m
% 5. Transform certain 2D and 3D grids
%==============================================================================

FAIRcheckFiles(mfilename);

%% Test module trafo, start with syntax
trafo('reset','trafo','splineTransformation2D',...
  'p',[4,5],'omega',[1,1],'m',[6,40]);
trafo('disp');
[scheme,parameter] = trafo;
trafo('clear');
%% test 2D cases

omega  = [0,1,0,2];
m      = [6,7];
xc     = getCellCenteredGrid(omega,m);
center = (omega(2:2:end)+omega(1:2:end))'/2;

% initialize various transformations
trafos = {
  'affine2D',{},...
  'affine2Dsparse',{},...
  'rigid2D',{},...
  'rotation2D',{'c',center},...
  'translation2D',{},...
  'splineTransformation2D',{'p',[4,5],'omega',omega,'m',m},...
  'splineTransformation2Dsparse',{'p',[4,5],'omega',omega,'m',m}
  };

for k=2:length(trafos)/2,
  fprintf(2,'\n\n\ntest [ %s ] --- %d of %d\n',...
    trafos{2*k-1},k,length(trafos)/2)
  optn = trafos{2*k};
  trafo('reset','trafo',trafos{2*k-1},'debug','on',optn{:});
  trafo('disp');
  wc = trafo('w0');
  fctn = @(wc) trafo(wc,xc);
  checkDerivative(fctn,wc+randn(size(wc)),'fig',1+mod(k,2));
  title(trafo)
  builtin('pause',2);
end;

%% test 3D cases
omega = [0,1,0,2,0,3]; m = [6,7,8]; xc = getCellCenteredGrid(omega,m);
trafos = {
  'affine3D',{},...
  'affine3Dsparse',{},...
  'rigid3D',{},...
  'splineTransformation3Dsparse',{'p',[4,5,6],'omega',omega,'m',m},...
  'translation3D',{}
  };

for k=1:length(trafos)/2,
  fprintf('\n\n\ntest [ %s ] --- %d of %d\n',...
    trafos{2*k-1},k,length(trafos)/2)
  optn = trafos{2*k};
  trafo('reset','trafo',trafos{2*k-1},'debug','on',optn{:});
  trafo('disp');
  wc = trafo('w0');
  fctn = @(wc) trafo(wc,xc);
  checkDerivative(fctn,wc+randn(size(wc)),'fig',1+mod(k,2));
  title(trafo)
  builtin('pause',2);
end;

testEnd;
%==============================================================================


text:
% DEMVOWELSFGPLVM2 Model the vowels data with a 2-D FGPLVM using RBF kernel.

% FGPLVM

% Fix seeds
randn('seed', 1e5);
rand('seed', 1e5);

dataSetName = 'vowels';
experimentNo = 2;

% load data
[Y, lbls] = lvmLoadData(dataSetName);

% Set up model
options = fgplvmOptions('fitc');
options.numActive = 200;
latentDim = 2;
d = size(Y, 2);

model = fgplvmCreate(latentDim, d, Y, options);

% Optimise the model.
iters = 1000;
display = 1;

model = fgplvmOptimise(model, display, iters);

% Save the results.
modelWriteResult(model, dataSetName, experimentNo);

if exist('printDiagram') & printDiagram
  lvmPrintPlot(model, lbls, dataSetName, experimentNo);
end

% Load the results and display dynamically.
lvmResultsDynamic(dataSetName, experimentNo, 'vector')

errors = lvmNearestNeighbour(model, lbls);


text:
function change_making_test01 ( )

%*****************************************************************************80
%
%% CHANGE_MAKING_TEST01 lists the problem data.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 August 2014
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CHANGE_MAKING_TEST01:\n' );
  fprintf ( 1, '  List the problem data.\n' );

  test_num = 7;

  coin_num_list = [ ...
    3, ...
    5, ...
    6, ...
    7, ...
    3, ...
    6, ...
    3 ];
  coin_value_list = { ...
    [  5,  9, 13 ], ...
    [  1,  4,  5,  8, 11 ], ...
    [  1,  5, 10, 25, 50, 100 ], ...
    [  1,  2,  6, 12, 24,  48,  60 ], ...
    [  1,  3,  4 ], ...
    [ 16, 17, 23, 24, 39,  40 ], ...
    [  6,  9, 20 ] };
  target_list = [ ...
    19, ...
    29, ...
    96, ...
    96, ...
     6, ...
   100, ...
    43 ];

  for test = 1 : test_num

    coin_num = coin_num_list(test);
    coin_value = coin_value_list{test};
    target = target_list(test);

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Test %d:\n', test );
    fprintf ( 1, '  Number of coins = %d\n', coin_num );
    fprintf ( 1, '  Values = ' );
    for i = 1 : coin_num
      fprintf ( 1, '  %d', coin_value(i) );
    end
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Target = %d\n', target );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function ans = calc_meanDose(doseBinsMidPtsV, volsHistV, volumeType)
%Calculate the mean dose for a given DVH
%  The last parameter 'volumeType' is a wash in this function again
%
%  MODIFICATION ALERT:  THIS FUNCTION IS UTILIZED BY THE DREXLER CODEBASE
%
%  LM: 6 Oct 06, JOD, corrected slight error in not taking middle of dose
%  bin.  Added warning if relative volume not close to one (0.5%
%  tolerance).
%      4 Apr 17, APA, passed the dose as min-points of bins. Hence, no need
%  to add half the binwidth to the passed dose. Call this function after
%  obtaining DVH from loadDVHMatrix.m, since it returns mid-points of dose bins.
%
% Usage: calc_meanDose(doseBinsMidPtsV, volsHistV)
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

% doseBinsMidPtsV = (doseBinsLowerPtsV(1:end-1)+doseBinsLowerPtsV(2:end))/2;
% ans = (sum(doseBinsMidPtsV.*volsHistV(1:end-1))+doseBinsLowerPtsV(end)*volsHistV(end))/sum(volsHistV);

ans = sum(doseBinsMidPtsV.*volsHistV)/sum(volsHistV);

return;


text:
% Demo: on-line design optimization for DCM of fMRI data.
% This demo simplifies the network identification of Friston et al. (2003).
% In brief, photic input enters V1, which is reciprocally connected to V5.
% The goal of the experiment is to assess whether attention modulates the
% feedforward connection from V1 to V5. This is addressed using Bayesian
% model comparison given fMRI data (i.e. model 1: attention modulates the
% V1->V5 connection; model 2: no modulatory effect). Within the session,
% each block consists of 16 seconds of stimulation and 16 of rest. The
% on-line design optimization consists in deciding, before each block,
% whether we modulate subjects' attention. This is done by comparing the
% design efficiency  of two canonical block designs (i.e. with and without
% attentional modulation).
% NB: the procedure is adaptative, because it updates the posterior on
% unknown model parameters after each block. This is why the design
% efficiency inceases across time.

close all
clear variables

% Basic settings
n_t = 16; % number of time samples per block
nblocks = 10; % number of blocks
nreg = 2; % number of regions
TR = 3e0; % sampling period (in sec)
microDT = 5e-2; % micro-time resolution (in sec)
f_fname = @f_DCMwHRF;
g_fname = @g_HRF3;

% construct atomic designs
nu = 2;
u1 = [0,ones(1,floor(n_t/2))];
u1(end:n_t) = 0;
u2 = zeros(1,n_t);
u = {[u1;u2];[u1;u1]}; % u{1}/u{2}: w/wo attentinal modulation

nm = 2; % # candidate models (see below)
truemodel = 1; % index of true model (which generates the data)

% model 1: u_att modulates connection from V1 to V5
A = [0 1
    1 0];
B{1} = zeros(nreg,nreg);
B{2} = zeros(nreg,nreg);
B{2}(2,1) = 1;
C = [1 0
     0 0];
D{1} = [0 0
    0 0];
D{2} = zeros(nreg,nreg);
D{3} = zeros(nreg,nreg);
[o4design{1},dim{1}] = getOptions4dcm(A,B,C,D,TR,microDT,n_t,1,1,1);

% model 1: no modulatory effect of attention
A = [0 1
    1 0];
nreg = size(A,1);
B{1} = zeros(nreg,nreg);
B{2} = zeros(nreg,nreg);
C = [1 0
     0 0];
D{1} = [0 0
    0 0];
D{2} = zeros(nreg,nreg);
D{3} = zeros(nreg,nreg);
[o4design{2},dim{2}] = getOptions4dcm(A,B,C,D,TR,microDT,n_t,1,1,1);

o4design{truemodel}.verbose = 0;

% simu parameters
t_A = exp([ -0.5
            -0.5
            ]);
t_Aself = -0;
t_B{1} = [];
t_B{2} = exp([ -0.5 ]);
t_C = exp([ -0.5 ]);
t_D{1} = [];
t_D{2} = [];
theta = zeros(dim{truemodel}.n_theta,1);
phi = zeros(dim{truemodel}.n_phi,1);
theta(o4design{truemodel}.inF.indA) = t_A;
for i=1:nu
    theta(o4design{truemodel}.inF.indB{i}) = t_B{i};
end
theta(o4design{truemodel}.inF.indC) = t_C;
alpha = Inf; % state noise precision
sigma = 5e-1; % measurement noise precision
x0 = zeros(dim{truemodel}.n,1); % initial conditions

% prepare graphics & montoring variables
hf = figure('color',[1 1 1]);
ha = subplot(3,2,1,'parent',hf,'nextplot','add','xlim',[1 nblocks+1]);
xlabel(ha,'time (blocks)')
ylabel(ha,'design efficiency')
ha2 = subplot(3,2,3,'parent',hf,'nextplot','add','xlim',[1 n_t*nblocks]);
xlabel(ha2,'time (scans)')
ylabel(ha2,'BOLD signal')
ha3 = subplot(3,2,2,'parent',hf,'nextplot','add','xlim',[1 nblocks+1]);
xlabel(ha3,'time (blocks)')
ylabel(ha3,'log p(y|m_1) -log p(y|m_2)')
e = zeros(2,1);
Y = zeros(dim{truemodel}.p,n_t*nblocks);
U = zeros(2,n_t*nblocks);
F = zeros(2,nblocks);
eb = zeros(1,nblocks);
vb = zeros(1,nblocks);

% on-line experiment

for tt=1:nblocks

    % 1- find best design given current information
    fprintf(1,'\n')
    fprintf(1,['-- Optimizing design (block ',num2str(tt),')...'])
    fprintf(1,'\n')
    for i=1:length(u) %loops over experimental designs
        fprintf(1,['- Candidate design #',num2str(i),'...'])
        fprintf(1,'\n')
        [e(i,tt)] = VBA_designEfficiency(f_fname,g_fname,dim,o4design,u{i},'models');
    end
    [em,ind] = max(e(:,tt));
    U(:,(tt-1)*n_t+1:tt*n_t) = u{ind};
    fprintf(1,['Optimizing design (block ',num2str(tt),')...  OK.'])
    fprintf(1,'\n')
    plot(ha,tt,e(1,tt)','ro')
    plot(ha,tt,e(2,tt)','go')
    legend(ha,{'u_{att} = off','u_{att} = on'},'Location','southeast','Orientation','horizontal')

    % 2- simulate BOLD response to chosen design under true model
    fprintf(1,['-- Simulating data (block ',num2str(tt),')...  '])
    [y,x,x0,eta,ee] = VBA_simulate (n_t,f_fname,g_fname,theta,phi,u{ind},alpha,sigma,o4design{truemodel},x0);
    Y(:,(tt-1)*n_t+1:tt*n_t) = y;
    x0 = x(:,end);

    try
        set(pl(1), 'XData', 1 : tt*n_t);
        set(pl(1), 'YData', [get(pl(1), 'YData') y(1,:)]);
        set(pl(2), 'XData', 1 : tt*n_t);
        set(pl(2), 'YData', [get(pl(2), 'YData') y(2,:)]);
    catch
        pl(1) = plot(ha2,1:n_t,y(1,:),'m');
        pl(2) = plot(ha2,1:n_t,y(2,:),'b');
    end
    legend(ha2,{'V1','V5'},'Location','southeast','Orientation','horizontal')
    fprintf(1,[' OK.'])
    fprintf(1,'\n')
    drawnow

    % 3- invert both models given new piece of dataset
    % NB: priors for the next block are updated to current posterior
    fprintf(1,['-- VB (block ',num2str(tt),'):   inverting model       '])
    for j=1:length(o4design)
        fprintf(1,repmat('\b',1,6))
        fprintf(1,[num2str(j),'/',num2str(length(o4design)),'...'])
        o4design{j}.DisplayWin = 0;
        o4design{j}.verbose = 0;
        [posterior,out] = VBA_NLStateSpaceModel(y,u{ind},f_fname,g_fname,dim{j},o4design{j});
        o4design{j}.priors = posterior;
        o4design{j}.priors.muX0 = posterior.muX(:,end);
        o4design{j}.priors.SigmaX0 = posterior.SigmaX.current{end};
        if tt > 1
            F(j,tt) = out.F + F(j,tt-1);
        else
            F(j,tt) = out.F;
        end
        OUT(j,tt).out = out;
        OUT(j,tt).posterior = posterior;
    end
    plot(ha3,tt,F(1,tt)-F(2,tt),'k*')

    eb(tt) = OUT(1,tt).posterior.muTheta(OUT(1,tt).out.options.inF.indB{2});
    vb(tt) = OUT(1,tt).posterior.SigmaTheta(OUT(1,tt).out.options.inF.indB{2},OUT(1,tt).out.options.inF.indB{2});

    fprintf(1,repmat('\b',1,24))
    fprintf(1,[' OK.'])
    fprintf(1,'\n')

end

ha4 = subplot(3,2,4,'parent',hf,'nextplot','add','xlim',[0 nblocks+1]);
xlabel(ha4,'time (blocks)')
ylabel(ha4,'modulatory effect')
plotUncertainTimeSeries(eb',1.96^2*vb',[],ha4)
hold(ha4,'on')
plot(ha4,[0,nblocks+1],[t_B{2},t_B{2}],'g--')

ha5 = subplot(3,2,5,'parent',hf);
imagesc(U,'parent',ha5)
title(ha5,'chosen (online) design')
xlabel(ha5,'time (scans)')
hold(ha5,'on')
plot(get(ha5,'xlim'),[1.5 1.5],'k')
set(ha5,'ytick',[1,2],'yticklabel',{'u1','u2'})
colormap(flipud(bone))

% invert full datasets at once
n_t = size(Y,2);
B{2}(2,1) = 1; % add modulatory effect (true model)
[OPT,DIM] = getOptions4dcm(A,B,C,D,TR,microDT,n_t,1,1,1);
[posterior,out] = VBA_NLStateSpaceModel(Y,U,f_fname,g_fname,DIM,OPT);

author: MBB-team

repo: VBA-toolbox


text:
function An = slnormalize(A, p, d)
%SLNORMALIZE Normalize the sub-arrays
%
% $ Syntax $
%   - An = slnormalize(A)
%   - An = slnormalize(A, p)
%   - An = slnormalize(A, [], d)
%   - An = slnormalize(A, p, d)
%
% $ Arguments $
%   - A:        the input array
%   - p:        the order of the norm (default = 2)
%   - d:        the dimension of subarrays (default = 1, column vectors)
%   - An:       the normalized array
%
% $ Description $
%   - An = slnormalize(A) normalizes the column vectors of A by 2nd-order.
%
%   - An = slnormalize(A, p) normalizes the column vectors of a by
%     pth-order.
%
%   - An = slnormalize(A, [], d) normalizes the subarrays of A along
%     dimensions specified by d by 2nd-order.
%
%   - An = slnormalize(A, p, d) normalizes the subarrays of A along
%     dimensions specified by d by pth-order.
%
% $ Remarks $
%   # Normalize an array by pth order means dividing the elements of
%     the array by its p-th norm.
%   # p can be inf or -inf. If p = inf, then the Lp norm is simply the
%     maximum magnitude value; while if p = -inf, then the Lp norm is the
%     minimum magnitude value.
%
% $ History $
%   - Created by Dahua Lin on Nov 19th, 2005
%   - Modified by Dahua Lin on Sep 10th, 2006
%       - replace slmul by slmulvec to increase efficiency
%

%% parse and verify input arguments
if nargin < 2 || isempty(p)
    p = 2;
end
if nargin < 3 || isempty(d)
    d = 1;
end
if ~isscalar(p)
    error('sltoolbox:notscalar', 'p should be a scalar');
end
if p == 0
    error('sltoolbox:zeroerr', 'p should not be zero');
end

%% compute
nrms = slnorm(A, p, d);
nrms = slenforce(nrms, 'positive');
An = slmulvec(A, 1 ./ nrms);

author: lmthang


text:
function [C,m]=covmatrix(X)
[K,n]=size(X);
X=double(X);
if K==1
    C=0;
    m=X;
else
    m=sum(X,1)/K;
    X=X-m(ones(K,1),:);
    C=(X'*X)/(K-1);
    m=m';
end
end


text:
% Test file for chebtech/mat2cell.m

function pass = test_mat2cell(pref)

if ( nargin < 2 )
    pref = chebtech.techPref();
end

for n = 1:2
    if ( n == 1 )
        testclass = chebtech1();
    else
        testclass = chebtech2();
    end

    f = testclass.make(@(x) [sin(x) cos(x) exp(x)], [], pref);
    g = testclass.make(@(x) sin(x), [], pref);
    h = testclass.make(@(x) [cos(x) exp(x)], [], pref);

    F = mat2cell(f, 1, [1 2]);
    pass(n, 1) = sum(F{1} - g) < 10*vscale(g)*eps;
    pass(n, 2) = all( sum(F{2} - h) < 10*max(vscale(h)*eps) );
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function frequency = convertPhaseToFrequency(phase,fs,sps)
factor =  (fs/sps)/(2*pi);
frequency = filter(ones(200,1)/200, 1, ... % Moving average
    diff(unwrap(phase(:)))*factor);
end

author: analogdevicesinc


text:
function [pp, pa, ap, aa] = tapas_sem_prosa_gen_lh(x, theta, ptheta)
%% Generate the likelihood function.
%
% Input
%   x       Time interval to plot the likelihood function.
%   theta   Parameters to use
%   ptheta  Priors
%
% Output
%   pp      Prosaccade trials, prosaccade action
%   pa      Prosaccade trial, antisaccade action
%   ap      Antisaccade trial, prosaccade action
%   aa      Antisaccade trial, antisaccade action
%

% aponteeduardo@gmail.com
% copyright (C) 2015
%

[PROSACCADE, ANTISACCADE] = tapas_sem_action_codes();

DIM_THETA = tapas_sem_prosa_ndims();

nt = numel(x);

u = struct('tt', zeros(nt, 1));
y = struct('t', x, 'a', zeros(nt, 1));

%
method = ptheta.method;

u.tt(:) = PROSACCADE;
y.a(:) = PROSACCADE;

pp = exp(method(y.t, y.a, u.tt, theta));

%

u.tt(:) = ANTISACCADE;
y.a(:) = PROSACCADE;

ap = exp(method(y.t, y.a, u.tt, theta));

%

u.tt(:) = PROSACCADE;
y.a(:) = ANTISACCADE;

pa = exp(method(y.t, y.a, u.tt, theta));

%

u.tt(:) = ANTISACCADE;
y.a(:) = ANTISACCADE;

aa = exp(method(y.t, y.a, u.tt, theta));

end

author: translationalneuromodeling

repo: tapas


text:
function r8_sind_test ( )

%*****************************************************************************80
%
%% R8_SIND_TEST tests R8_SIND.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 January 2015
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8_SIND_TEST\n' );
  fprintf ( 1, '  R8_SIND computes the sine of an angle\n' );
  fprintf ( 1, '  given in degrees.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  ANGLE    R8_SIND(ANGLE)\n' );
  fprintf ( 1, '\n' );

  for i = 0 : 15 : 360
    angle = i;
    fprintf ( 1, '  %8.2f  %14.6g\n', angle, r8_sind ( angle ) );
  end

  return
end


text:
function obs_mask= ObservMask( band )
%MASK Create the mask of obervation using the overlap of all bands.
%
% Syntax
%
%     obs_mask= ObservMask( band )
%
% Description
%
%     Create the mask of obervation using the overlap of all bands, which
%     have been calculated in LoadData.
% History
%
%     1. Create this function. (22. September, 2017)
%
% Input arguments
%
%     band     Blue band.
%
% Output arguments
%
%     obs_mask      Observation mask(outside).
%
% See also: Loaddata
%
%
% Author:  Shi Qiu (shi.qiu@ttu.edu)
% Date: 22. September, 2017

    obs_mask=band>-9999;
%
%     obs_mask=band~=0;
%
%     %% simulate Sentinel 2 scene by 100 km x 100 km.
%     % spatial resolution is 30 m for each pixel. That means there will be
%     % 3333 x 3333 30m pixels should be selected expanding from the centre.
%     s2_pixels=[3333 3333];
%     s2_pixels_half=round(s2_pixels./2);
%     obs_dim=size(obs_mask);
%     pixel_centre=round(obs_dim./2);
%
%     s2_mask=obs_mask.*0;
%     r_start=pixel_centre(1)-s2_pixels_half(1);
%     r_end=pixel_centre(1)+s2_pixels_half(1);
%     c_start=pixel_centre(2)-s2_pixels_half(2);
%     c_end=pixel_centre(2)+s2_pixels_half(2);
%     s2_mask(r_start:r_end,c_start:c_end)=1;
%     obs_mask=s2_mask&obs_mask;
end

author: GERSL

repo: Fmask

save_path: github-repos/MATLAB/GERSL-Fmask


text:
%%
% Perform softmax prediction.
%
% Thang Luong @ 2015, <lmthang@stanford.edu>
%
%%
function [cost, probs, scores, scoreIndices] = softmaxLayerForward(W, inVec, predLabels, curMask)
  mask = curMask.mask;
  unmaskedIds = curMask.unmaskedIds;

  % softmax_h -> predictions
  [probs, scores, norms] = softmax(W*inVec, mask);

  % cost
  predLabels = predLabels(unmaskedIds);
  scoreIndices = sub2ind(size(scores), predLabels, unmaskedIds); % 1 * length(tgtPredictedWords)
  cost = - sum(scores(scoreIndices)) + sum(log(norms).*mask);
end


text:
function [d,g,rr,ss]=v_sigalign(s,r,maxd,m,fs)
%V_SIGALIGN align a clean reference with a noisy signal [d,g,rr,ss]=(s,r,maxd,m,fs)
% Inputs:
%            m  mode
%                 u = unity gain
%                 g = find optimal gain [default]
%                 a = A-weight the signals
%                 b = weight signals by BS-468
%                 s = find delay to maximize the correlation coefficient between r and s [default]
%                 S = find delay to maximize the energy of the component of r in s
%                 p = plot result
%            s  test signal
%            r  reference signal
%         maxd  [+-max] or [min max] delay allowed in samples or fractions of length(r)
%               default is maximum that ensures at least 50% of r or s in the overlap
%           fs  sample frequency (only used for filtering and plotting)
%
% Outputs:
%            d = optimum delay to apply to r
%            g = optimal gain to apply to r
%           rr = g*r(* -d)  [zero padded to match s if ss output is not given]
%           ss = s truncated if necessary to match the length of rr

%      Copyright (C) Mike Brookes 2011
%      Version: $Id: v_sigalign.m 10865 2018-09-21 17:22:45Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Bugs/Suggestions
% 1. add option to calculate a DC offset
% 2. optionally find optimal fractional time shift
% 3. split long signals into chunks to reduce memory requirements

ns=length(s);
nr=length(r);
if numel(s)~=ns || numel(r)~=nr
    error('Inputs cannot be matrices');
end
s=s(:);
r=r(:);
if nargin<3
    maxd=[];
end
switch numel(maxd)
    case 0
        if nr<ns
            lmm=[-0.25*nr ns-0.75*nr];
        else
            lmm=[-0.25*ns nr-0.75*ns];
        end
    case 1
        lmm=[-maxd maxd];
    otherwise
        lmm=maxd(1:2);
end
lmm=round(lmm.*(1+(nr-1)*(abs(lmm)<1)));  % convert fractions of nr to samples
lmin=lmm(1);
lmax=lmm(2);
lags=lmax-lmin+1;
if lags<=0
    error('Invalid lag limits');
end
if nargin<4 || ~numel(m)
    m='gs';
end
if nargin<5 || ~numel(fs)
    fs=[];
else
    if any(m=='a')
        [b,a]=v_stdspectrum(2,'z',fs);
        s=filter(b,a,s);
        r=filter(b,a,r);
    elseif any(m=='b')
        [b,a]=v_stdspectrum(8,'z',fs);
        s=filter(b,a,s);
        r=filter(b,a,r);
    end
end

% now do cross correlation

rxi=max(1,1-lmin);   % first reference sample needed
rxj=min(nr,ns-lmax); % last reference sample needed
nrx=rxj-rxi+1;          % length of reference segment
if nrx<1
    error('Reference signal too short');
end
fl=2^nextpow2(lmax-lmin+nrx);
sxi=max(1,rxi+lmin); % first signal sample needed
sxj=min(ns,rxj+lmax); % last signal sample needed
rs=v_irfft(v_rfft([s(sxi:sxj); zeros(fl-sxj+sxi-1,1)]).*conj(v_rfft([r(rxi:rxj); zeros(fl-rxj+rxi-1,1)])));
rsu=rs(1:lags);
ssq=cumsum(s(sxi:sxj).^2);
ssqd=[ssq(nrx); ssq(nrx+1:lmax-lmin+nrx)-ssq(1:lmax-lmin)];
if any (m=='S') % maximize energy of common component
    [cmx,icx]=max(abs(rsu)); % maximize cross correlation
else
    [cmx,icx]=max(rsu.^2./ssqd); % maximize correlation coefficient
end
d=icx-1+lmin;
ia=max(1,d+1); % first sample of s in common region
ja=min(ns,d+nr); % last sample of s in common region
ija=ia:ja;
ijad=ija-d;
rr=r(ijad);
ss=s(ija);
if any (m=='u')
    g=1;
else
g=sum(rr.*ss)/sum(rr.^2);   % gain to apply to r
end
rr=rr*g;
if ~nargout || any(m=='p')
    xco=sum(rr.*ss)/sqrt(sum(rr.^2)*sum(ss.^2));
    snr=sum(rr.^2)/sum((rr-ss).^2);
    if numel(fs)==1
        tun='s';
    else
        tun='samples';
        fs=1;
    end
    subplot(311);
    plot(ija/fs,rr);
    pm='+-';
    title(sprintf('Ref delay = %.2g %s, %cGain = %.2g dB, Xcorr = %.2g, SNR = %.2g dB',d/fs,tun,pm(1+(g<0)),20*log10(g),xco,10*log10(snr)));
    ylabel('Reference');
    set(gca,'XLim',ija([1 end])/fs);
    axh(2)=gca;
    subplot(312);
    plot(ija/fs,ss);
    ylabel('Signal');
    set(gca,'XLim',ija([1 end])/fs);
    axh(1)=gca;
    subplot(313);
    plot(ija/fs,ss-rr);
    ylabel('Residual');
    xlabel(sprintf('Time (%s)',tun));
    set(gca,'XLim',ija([1 end])/fs);
    axh(3)=gca;
    linkaxes(axh(1:3),'x');
end
if nargout==3
    rr=[zeros(ia-1,1); rr; zeros(ns-ja,1)]; % force to be the size of s
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function [bcx,bcy] = specific_flow(xbd,ybd)
%poiseuille_flow   Reference problem 5.1 inflow condition
%   [bcx,bcy] = specific_flow(xbd,ybd);
%   input
%          xbd          x coordinate vector
%          ybd          y coordinate vector
%
%   specifies Poiseuille flow boundary condition
%   IFISS function: DJS; 6 March 2005.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage
bcx=1-ybd.*ybd; bcy=0*xbd;
return


text:
function [trans, points] = registerPoints3dAffine(points, target, varargin)
%REGISTERPOINTS3DAFFINE Fit 3D affine transform using iterative algorithm.
%
%   TRANS = registerPoints3dAffine(POINTS, TARGET)
%   Computes the affine transform that maps the shape defines by POINTS
%   onto the shape defined by the points TARGET. Both POINTS and TARGET are
%   N-by-3 array of point coordinates, not necessarily the same size.
%   The result TRANS is a 4-by-4 affine transform.
%
%   TRANS = registerPoints3dAffine(POINTS, TARGET, NITER)
%   Specifies the number of iterations for the algorithm.
%
%   [TRANS, POINTS2] = registerPoints3dAffine(...)
%   Also returns the set of transformed points.
%
%   Example
%     registerPoints3dAffine
%
%   See also
%     transforms3d, fitAffineTransform3d
%

% ------
% Author: David Legland
% E-mail: david.legland@inra.fr
% Created: 2015-02-24, using Matlab 8.4.0.150421 (R2014b)
% Copyright 2015-2022 INRA - Cepia Software Platform

nIters = 10;
if ~isempty(varargin)
    nIters = varargin{1};
end

% keep original points to transform them at each
trans = [1 0 0 0;0 1 0 0;0 0 1 0;0 0 0 1];

for i = 1:nIters
    % identify target points for each source point
    inds = findClosestPoint(points, target);
    corrPoints = target(inds, :);

    % compute transform for current iteration
    trans_i = fitAffineTransform3d(points, corrPoints);

    % apply transform, and update cumulated transform
    points = transformPoint3d(points, trans_i);
    trans = trans_i * trans;
end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function turkRemainingHITs(input_csv_fname, result_csv_fname, output_csv_fname, uniqueColumn)
% this function take a turk input.csv file and result.csv file
% and output a file remaining_input.csv
% for the HITs that haven't been done
% based on a specific column
%{
%example:
input_csv_fname    = '/data/vision/torralba/gigaSUN/inputTurk/c/cybercafe.csv';
result_csv_fname   = '/data/vision/torralba/gigaSUN/resultTurk/c/cybercafe.csv';
output_csv_fname   = '/data/vision/torralba/gigaSUN/inputTurk/c/cybercafe_2nd.csv';
uniqueColumn = 'sunID';
%}

input_csv  = readCSVfromMTurk(input_csv_fname);
result_csv = readCSVfromMTurk(result_csv_fname);

uniqueColumnIDinput = find(ismember(input_csv(1,:),uniqueColumn));
uniqueColumnIDresult = find(ismember(result_csv(1,:),['Input.' uniqueColumn]));

good_input_row = [];
for row = 2:size(input_csv,1)
    if ~any(find(ismember(result_csv(2:end,uniqueColumnIDresult)',input_csv{row,uniqueColumnIDinput})))
        good_input_row(end+1) = row;
    end
end

fp = fopen(output_csv_fname,'w');
if ~isempty(good_input_row)

    fid = fopen(input_csv_fname,'r');
    tline = fgetl(fid);
    fprintf(fp,'%s\n',tline);
    for row = 2:size(input_csv,1)
        tline = fgetl(fid);
        if ~isempty(find(good_input_row, row))
            fprintf(fp,'%s\n',tline);
        end
    end
    fclose(fid);
end
fclose(fp);

author: jianxiongxiao

repo: ProfXkit


text:
function results = vl_test_imarray(varargin)
% VL_TEST_IMARRAY
vl_test_init ;

function test_movie_rgb(s)
A = rand(23,15,3,4) ;
B = vl_imarray(A,'movie',true) ;

function test_movie_indexed(s)
cmap = get(0,'DefaultFigureColormap') ;
A = uint8(size(cmap,1)*rand(23,15,4)) ;
A = min(A,size(cmap,1)-1) ;
B = vl_imarray(A,'movie',true) ;

function test_movie_gray_indexed(s)
A = uint8(255*rand(23,15,4)) ;
B = vl_imarray(A,'movie',true,'cmap',gray(256)) ;

for k=1:size(A,3)
  vl_assert_equal(squeeze(A(:,:,k)), ...
                  frame2im(B(k))) ;
end

function test_basic(s)
M = 3 ;
N = 4 ;
width = 32 ;
height = 15 ;
for i=1:M
  for j=1:N
    A{i,j} = rand(width,height) ;
  end
end
A1 = A';
A1 = cat(3,A1{:}) ;
A2 = cell2mat(A) ;
B = vl_imarray(A1, 'layout', [M N]) ;
vl_assert_equal(A2,B) ;

author: Cloud-CV

repo: object-proposals


text:
function cs_demo (do_pause, matrixpath)
%CS_DEMO run all CXSparse demos.
%   cs_demo(0) will run all demos without pausing.
%
% Example:
%   cs_demo
% See also: cs_demo1, cs_demo2, cs_demo3

%   Copyright 2006-2007, Timothy A. Davis.
%   http://www.cise.ufl.edu/research/sparse

help cs_demo
if (nargin < 1)
    do_pause = 1 ;
end
if (nargin < 2)
    matrixpath = [] ;
end

figure (1)
clf

fprintf ('\n\n-------------------------------------------------------\n') ;
help cs_demo1 ;
cs_demo1 (matrixpath) ;

fprintf ('\n\n-------------------------------------------------------\n') ;
help cs_demo2
cs_demo2 (do_pause, matrixpath) ;

fprintf ('\n\n-------------------------------------------------------\n') ;
help cs_demo3
cs_demo3 (do_pause, matrixpath) ;

fprintf ('\n\n-------------------------------------------------------\n') ;
help private/ex_1
ex_1

fprintf ('\n\n-------------------------------------------------------\n') ;
help private/ex2
ex2

fprintf ('\n\n-------------------------------------------------------\n') ;
help private/ex3
ex3

fprintf ('\nAll CXSparse demos finished.\n') ;

author: fieldtrip


text:
function [y,Jyl,Jyr,ilocl,ilocr,evalcnt]=greedy2_cross(n, fun, tol, varargin)
% Two-site greedy cross interpolation scheme.
%   [y,Jyl,Jyr,ilocl,ilocr]=greedy2_cross(n, fun, tol, varargin)
% Tries to interpolate the tensor with mode sizes n specified by the
% function fun up to the accuracy tol using the "basic" greedy restricted
% cross interpolation algorithm in the TT format.
% The method sequentially adds one pivot to each superblock in a forward
% DMRG-type half-sweep (k=1,...,d).
%
% The input n should be a vector of mode sizes of length d,
% fun = @(ind)fun(ind) is a sought function of index ind.
% By default, ind is an array of d indices, and the function fun should
% return one value of the corresponding tensor entry.
% To speed up the computations, set the optional parameter 'vec' to true,
% and provide the function which takes ind as an array of sizes M x d, and
% returns an array of M values.
%
% In addition to the indexwise function, one may provide the value-wise
% function of another tt_tensor (funcrs style) via optional parameters
% 'aux', 'auxfun' (see below).
% The resulting tensor is the sum    y(ind)=fun(ind) + auxfun(aux(ind)).
%
% The first output parameter is the computed tt_tensor, the rest four
% parameters return the pivot indices (left-global, right-global, left-local, right-local) for
% testing purposes.
%
% Optional arguments are provided in the form
% 'PropertyName1',PropertyValue1,'PropertyName2',PropertyValue2 and so on.
% The list of option names and default values:
%       o nswp - maximal number of sweeps [20]
%       o tol_exit - stopping difference between consecutive iterations [tol]
%       o verb - verbosity level, 0-silent, 1-sweep info, 2-block info [1]
%       o vec - whether fun can accept and return vectorized values [false]
%       o aux - set of tt_tensors for auxiliary funcrs contribution []
%       o auxfun - an auxiliary function defined pointwise at the elements
%           of aux []
%       o 'locsearch' - an algorithm for the error pivoting in superblocks:
%           try a lottery of n*r random elements ('lot'), or
%           conduct two-dimensional maxvol ALS iteration ('als').
%           The first one uses less evaluations, but may be costly due to a
%           nonvectorized MATLAB loop ['lot']
%
%
% This procedure implements Algorithm 2 from
%   D. Savostyanov, Quasioptimality of maximum-volume cross interpolation of tensors, Linear Algebra Applications 458, pp 217-244, 2014.
%    http://dx.doi.org/10.1016/j.laa.2014.06.006
% Please cite this paper if your research benefits from the use of this code.
% Development of the MATLAB version: S. Dolgov.
% Please send feedback to: {dmitry.savostyanov,sergey.v.dolgov}@gmail.com
%
%---------------------------

if (~isempty(varargin))
    v1 = varargin{1};
    if (isa(v1, 'cell'))
        varargin=v1;
    end;
end;
vars = varargin;

nswp = 20;
tol_exit = tol;
verb = 1;
vec = false;
aux = []; % Extra tt_tensors to pass into this cross
auxfun = []; % the total function equals fun(ind)+auxfun(aux(ind)), since simple aux(ind) via tt_tensor/subsref suxx.
% locsearch = 'als';
locsearch = 'lot';
y0 = [];

i = 1;
while (i<length(vars))
    switch lower(vars{i})
        case 'nswp'
            nswp=vars{i+1};
        case 'tol_exit'
            tol_exit=vars{i+1};
        case 'verb'
            verb = vars{i+1};
        case 'vec'
            vec = vars{i+1};
        case 'aux'
            aux = vars{i+1};
        case 'auxfun'
            auxfun = vars{i+1};
        case 'locsearch'
            locsearch = vars{i+1};
        case 'xtru'
            xtru = vars{i+1};
        case 'y0'
            y0 = vars{i+1};
        otherwise
            warning('Option %s was not recognized', vars{i});
    end;
    i=i+2;
end;

d = numel(n);
ry = ones(d+1,1);
y = cell(d,1);

if (~isempty(aux))
    if (isa(aux, 'tt_tensor'))
        raux = aux.r;
        aux = core2cell(aux);
        phiauxl = cell(d+1,1); phiauxl{1}=1; phiauxl{d+1}=1;
        phiauxr = cell(d+1,1); phiauxr{1}=1; phiauxr{d+1}=1;
%         for i=2:d
%             phiauxl{i} = ones(1,raux(i));
%             phiauxr{i} = ones(raux(i),1);
%         end;
        Raux = 1;
    else
        Raux = size(aux, 2);
        raux = ones(d+1,Raux);
        aux = [aux; cell(d-1, Raux)];
        phiauxl = cell(d+1,Raux);
        phiauxr = cell(d+1,Raux);
        for j=1:Raux
            raux(:,j) = aux{1,j}.r;
            aux(:,j) = core2cell(aux{1,j});
            phiauxl{1,j}=1; phiauxl{d+1,j}=1;
            phiauxr{1,j}=1; phiauxr{d+1,j}=1;
%             for i=2:d
%                 phiauxl{i,j} = ones(1,raux(i,j));
%                 phiauxr{i,j} = ones(raux(i,j),1);
%             end;
        end;
    end;
end;

% Factorized inverse interpolation matrix -- in the form U^{-1}, L^{-1}
mid_inv = cell(d+1,2); mid_inv{1,1}=1; mid_inv{1,2}=1; mid_inv{d+1,1}=1; mid_inv{d+1,2}=1;

Jyr = cell(d+1,1);
Jyl = cell(d+1,1);
ilocr = cell(d+1,1);
ilocl = cell(d+1,1);

evalcnt = 0;

% Start with some rand indices
for i=1:d-1
    if (isempty(y0))
        ilocl{i+1} = rand(1,1)*ry(i)*n(i);
        ilocl{i+1} = round(ilocl{i+1});
        if (ilocl{i+1}==0)
            ilocl{i+1} = 1;
        end;
    else
        ilocl{i+1} = y0(i);
    end
    Jyl{i+1} = indexmerge(Jyl{i}, (1:n(i))');
    Jyl{i+1} = Jyl{i+1}(ilocl{i+1},:);

    if (~isempty(aux))
        for j=1:Raux
            craux1 = reshape(aux{i,j}, raux(i,j), n(i)*raux(i+1,j));
            craux1 = phiauxl{i,j}*craux1;
            craux1 = reshape(craux1, ry(i)*n(i), raux(i+1,j));
            phiauxl{i+1,j} = craux1(ilocl{i+1}, :);
        end;
    end;
end;
% Perform one sweep back-forth to fix the initial index
for i=d:-1:1
    J = indexmerge(Jyl{i}, (1:n(i))', Jyr{i+1});
    evalcnt = evalcnt + size(J,1);
    cry1 = autovecfun(fun, J, vec);
    if (~isempty(aux))
        craux = zeros(ry(i)*n(i)*ry(i+1), Raux);
        for j=1:Raux
            craux1 = reshape(aux{i,j}, raux(i,j), n(i)*raux(i+1,j));
            craux1 = phiauxl{i,j}*craux1;
            craux1 = reshape(craux1, ry(i)*n(i), raux(i+1,j));
            craux1 = craux1*phiauxr{i+1,j};
            craux(:,j) = reshape(craux1, ry(i)*n(i)*ry(i+1), 1);
        end;
        cry1 = cry1+auxfun(craux);
    end;

    if (i>1)
        cry1 = reshape(cry1, ry(i), n(i)*ry(i+1));
        [v,ilocr{i}] = max(abs(cry1.'));
        ilocr{i} = ilocr{i}(:);
        mid_inv{i,1} = 1/cry1(:,ilocr{i});
        mid_inv{i,2} = 1;
        Jyr{i} = indexmerge((1:n(i))', Jyr{i+1});
        Jyr{i} = Jyr{i}(ilocr{i},:);

        if (~isempty(aux))
            for j=1:Raux
                phiauxr{i,j} = reshape(aux{i,j}, raux(i,j)*n(i), raux(i+1,j));
                phiauxr{i,j} = phiauxr{i,j}*phiauxr{i+1,j};
                phiauxr{i,j} = reshape(phiauxr{i,j}, raux(i,j), n(i)*ry(i+1));
                phiauxr{i,j} = phiauxr{i,j}(:, ilocr{i});
            end;
        end;
    end;

    y{i} = reshape(cry1, ry(i), n(i), ry(i+1));
end;
for i=1:d
    J = indexmerge(Jyl{i}, (1:n(i))', Jyr{i+1});
    evalcnt = evalcnt + size(J,1);
    cry1 = autovecfun(fun, J, vec);
    if (~isempty(aux))
        craux = zeros(ry(i)*n(i)*ry(i+1), Raux);
        for j=1:Raux
            craux1 = reshape(aux{i,j}, raux(i,j), n(i)*raux(i+1,j));
            craux1 = phiauxl{i,j}*craux1;
            craux1 = reshape(craux1, ry(i)*n(i), raux(i+1,j));
            craux1 = craux1*phiauxr{i+1,j};
            craux(:,j) = reshape(craux1, ry(i)*n(i)*ry(i+1), 1);
        end;
        cry1 = cry1+auxfun(craux);
    end;

    if (i<d)
        cry1 = reshape(cry1, ry(i)*n(i), ry(i+1));
        [v,ilocl{i+1}] = max(abs(cry1));
        mid_inv{i+1,1} = 1/cry1(ilocl{i+1});
        mid_inv{i+1,2} = 1;
        Jyl{i+1} = indexmerge(Jyl{i}, (1:n(i))');
        Jyl{i+1} = Jyl{i+1}(ilocl{i+1},:);
        if (~isempty(aux))
            for j=1:Raux
                phiauxl{i+1,j} = reshape(aux{i,j}, raux(i,j), n(i)*raux(i+1,j));
                phiauxl{i+1,j} = phiauxl{i,j}*phiauxl{i+1,j};
                phiauxl{i+1,j} = reshape(phiauxl{i+1,j}, ry(i)*n(i), raux(i+1,j));
                phiauxl{i+1,j} = phiauxl{i+1,j}(ilocl{i+1}, :);
            end;
        end;
    end;

    y{i} = reshape(cry1, ry(i), n(i), ry(i+1));
end;

last_sweep = false;
maxy = 0;
max_dx = 0;
swp = 1;
dir = 1;
i = 1;
while (swp<=nswp)
    % Generate a random test set for new indices
    % Prepare candidate index sets--the ones with the current crosses excluded
    cind1 = (1:ry(i)*n(i))';
    cind2 = (1:n(i+1)*ry(i+2))';
    cind1(ilocl{i+1})=[];
    cind2(ilocr{i+1})=[];

    if (strcmp(locsearch, 'als'))
        testsz = min(numel(cind1), numel(cind2));
    else
        %%% % rn Lottery
        % Now draw random entries
        tind = rand(min(numel(cind1), numel(cind2)), 1)*numel(cind1)*numel(cind2);
        tind = round(tind);
        tind(tind>numel(cind1)*numel(cind2))=numel(cind1)*numel(cind2);
        tind(tind<1)=1;
        tind = unique(tind);
        testsz = size(tind,1);
    end;

    if (~isempty(aux))
        craux1 = cell(1,Raux);
        craux2 = cell(1,Raux);
        for j=1:Raux
            craux = reshape(aux{i,j}, raux(i,j), n(i)*raux(i+1,j));
            craux = phiauxl{i,j}*craux;
            craux1{j} = reshape(craux, ry(i)*n(i), raux(i+1,j));
            craux = reshape(aux{i+1,j}, raux(i+1,j)*n(i+1), raux(i+2,j));
            craux = craux*phiauxr{i+2,j};
            craux2{j} = reshape(craux, raux(i+1,j), n(i+1)*ry(i+2));
        end;
    end;

    % Check that we are not in the full rank case
    if (testsz>0)
        if (strcmp(locsearch, 'als'))
            % 2D ALS cross
            % Evaluate y at tind
            cry1 = reshape(y{i}, ry(i)*n(i), ry(i+1));
            cry2 = reshape(y{i+1}, ry(i+1), n(i+1)*ry(i+2));
            ys1 = cry1*mid_inv{i+1,1};
            ys1 = ys1(cind1,:);
            ys2 = mid_inv{i+1,2}*cry2;
            ys2 = ys2(:,cind2);
            % Full indices
            J1 = indexmerge(Jyl{i}, (1:n(i))');
            J1c = J1(cind1,:);
            J2 = indexmerge((1:n(i+1))', Jyr{i+2});
            J2c = J2(cind2,:);
            rz = 2;
            cre2 = randn(numel(cind2), rz);
            [cre2,rv]=qr(cre2,0);
            indr = maxvol2(cre2);
            Jr = J2c(indr,:);
            Ye2 = ys2(:,indr);

            J = indexmerge(J1c,Jr);
            evalcnt = evalcnt + size(J,1);
            cre1 = autovecfun(fun, J, vec);
            if (~isempty(aux))
                craux = zeros(numel(cind1)*numel(indr), Raux);
                for j=1:Raux
                    craux(:,j) = reshape(craux1{j}(cind1,:)*craux2{j}(:,cind2(indr)), [], 1);
                end;
                craux = auxfun(craux);
                cre1 = cre1+craux;
            end;
            maxy = max(maxy, max(abs(cre1)));
            cre1 = reshape(cre1, numel(cind1), rz);
            cre1 = cre1-ys1*Ye2;
            [zmax1,imaxnew]=max(abs(cre1(:)));
            imaxnew = tt_ind2sub([numel(cind1), rz], imaxnew);
            imax1 = imaxnew(1);
            [cre1,rv]=qr(cre1,0);
            indl = maxvol2(cre1);
            Jl = J1c(indl,:);
            Ye1 = ys1(indl,:);

            J = indexmerge(Jl,J2c);
            evalcnt = evalcnt + size(J,1);
            cre2 = autovecfun(fun, J, vec);
            if (~isempty(aux))
                craux = zeros(numel(indl)*numel(cind2), Raux);
                for j=1:Raux
                    craux(:,j) = reshape(craux1{j}(cind1(indl),:)*craux2{j}(:,cind2), [], 1);
                end;
                craux = auxfun(craux);
                cre2 = cre2+craux;
            end;
            maxy = max(maxy, max(abs(cre2)));
            cre2 = reshape(cre2, rz, numel(cind2));
            cre2 = cre2-Ye1*ys2;
            [zmax2,imaxnew]=max(abs(cre2(:)));
            imaxnew = tt_ind2sub([rz, numel(cind2)], imaxnew);
            imax2 = imaxnew(2);
            imax1 = cind1(imax1);
            imax2 = cind2(imax2);
            emax = max(zmax1,zmax2);
        else
            %%% % rn Lottery
            tind = tt_ind2sub([numel(cind1), numel(cind2)], tind);
            tind = [cind1(tind(:,1)), cind2(tind(:,2))];
            % Full indices
            J1 = indexmerge(Jyl{i}, (1:n(i))');
            J1c = J1(tind(:,1),:);
            J2 = indexmerge((1:n(i+1))', Jyr{i+2});
            J2c = J2(tind(:,2),:);
            J = [J1c,J2c];
            % Evaluate the function
            evalcnt = evalcnt + testsz;
            crt = autovecfun(fun, J, vec);
            if (~isempty(aux))
                craux = zeros(testsz, Raux);
                for k=1:Raux
                    for j=1:testsz
                        craux(j,k) = craux1{k}(tind(j,1),:)*craux2{k}(:,tind(j,2));
                    end;
                end;
                crt = crt+auxfun(craux);
            end;
            maxy = max(maxy, max(abs(crt)));
            % Evaluate y at tind
            cry1 = reshape(y{i}, ry(i)*n(i), ry(i+1));
            cry2 = reshape(y{i+1}, ry(i+1), n(i+1)*ry(i+2));
            % Subtract the current approx.
            cry = zeros(testsz, 1);
            % Apply the inverse interp at the middle
            cre1 = cry1*mid_inv{i+1,1};
            cre2 = mid_inv{i+1,2}*cry2;
            for j=1:testsz
                cry(j) = cre1(tind(j,1),:)*cre2(:,tind(j,2));
            end;
            cre = crt-cry;
            % Take the max-err index to enrich
            [emax,imax2] = max(abs(cre));

            % All again: at (:,tind(imax,2)), run the error maximization
            J1c = J1(cind1,:);
            J = indexmerge(J1c, J2c(imax2,:));
            evalcnt = evalcnt + size(J,1);
            crt = autovecfun(fun, J, vec);
            if (~isempty(aux))
                craux = zeros(numel(cind1), Raux);
                for j=1:Raux
                    craux(:,j) = craux1{j}(cind1,:)*craux2{j}(:,tind(imax2,2));
                end;
                crt = crt+auxfun(craux);
            end;
            maxy = max(maxy, max(abs(crt)));
            % Subtract the current approx.
            cry = cre1(cind1,:)*cre2(:,tind(imax2,2));
            cre = crt-cry;
            % Take the max-err index to enrich
            [emax,imax1] = max(abs(cre));
            imax1 = cind1(imax1);
            imax2 = tind(imax2,2);
            % imax1 samples from ry(i)*n(i)
            % imax2 samples from n(i+1)*ry(i+2)
        end;

        dx = emax/maxy;
        max_dx = max(max_dx,dx);
        if (verb>1)
            fprintf('=greedy_cross= i=%d, swp=%d, testsz=%d, emax=%3.3e, dx=%3.3e, cond=[%3.3e,%3.3e]\n', i, swp, testsz, emax, dx, cond(mid_inv{i+1,1}), cond(mid_inv{i+1,2}));
        end;

        if (dx>tol)
            % Now evaluate new factors
            J1m = J1(imax1,:);
            J2m = J2(imax2,:);
            % Generate a full index = [Jl, (1:n)', Jr]
            % Jyl to lowest index           n(i) to the middle one,                      Jyr to the senior index
            Jl = indexmerge(J1, J2m);
            Jr = indexmerge(J1m, J2);

            evalcnt = evalcnt + size(Jl,1);
            cre1 = autovecfun(fun, Jl, vec);
            if (~isempty(aux))
                craux = zeros(ry(i)*n(i), Raux);
                for j=1:Raux
                    craux(:,j) = craux1{j}*craux2{j}(:,imax2);
                end;
                cre1 = cre1+auxfun(craux);
            end;

            evalcnt = evalcnt + size(Jr,1);
            cre2 = autovecfun(fun, Jr, vec);
            if (~isempty(aux))
                craux = zeros(n(i+1)*ry(i+2), Raux);
                for j=1:Raux
                    craux(:,j) = craux1{j}(imax1,:)*craux2{j};
                end;
                cre2 = cre2+auxfun(craux);
            end;

            cre1 = reshape(cre1, ry(i)*n(i), 1);
            cre2 = reshape(cre2, 1, n(i+1)*ry(i+2));

            % Expand blocks
            y{i} = [cry1, cre1];
            y{i+1} = [cry2; cre2];
            ry(i+1)=ry(i+1)+1;
            % ilocl{i+2} is smashed now, since ry(i+1) changed. Recover...
            if (i<d-1)
                ii = floor((ilocl{i+2}-1)/(ry(i+1)-1));
                ilocl{i+2} = ilocl{i+2}+ii;
            end;
            % Expand indices
            Jyl{i+1} = [Jyl{i+1}; J1m];
            Jyr{i+1} = [Jyr{i+1}; J2m];
            ilocl{i+1} = [ilocl{i+1}; imax1];
            ilocr{i+1} = [ilocr{i+1}; imax2];
            % We know that the inverse LU may be written analytically
            uold = mid_inv{i+1,1};
            lold = mid_inv{i+1,2};
            % Expanding vectors
            erow = cry1(imax1, :);
            ecol = cre1(ilocl{i+1}(1:ry(i+1)-1), 1);
            % ..and the scalar
            eel = cre1(imax1,1);
            ecol = lold*ecol;
            erow = erow*uold;
            eel = eel - erow*ecol; % alpha-dA^{-1}c
            ecol = uold*ecol; % A^{-1} c
            erow = erow*lold; % d A^{-1}
            % New inv(U)
            mid_inv{i+1,1} = zeros(ry(i+1), ry(i+1));
            mid_inv{i+1,1}(1:ry(i+1)-1, 1:ry(i+1)-1)=uold;
            mid_inv{i+1,1}(1:ry(i+1)-1, ry(i+1)) = -ecol/eel;
            mid_inv{i+1,1}(ry(i+1), ry(i+1))=1/eel;
            % New inv(L)
            mid_inv{i+1,2} = zeros(ry(i+1), ry(i+1));
            mid_inv{i+1,2}(1:ry(i+1)-1, 1:ry(i+1)-1)=lold;
            mid_inv{i+1,2}(ry(i+1), 1:ry(i+1)-1) = -erow;
            mid_inv{i+1,2}(ry(i+1), ry(i+1))=1;

            y{i} = reshape(y{i}, ry(i), n(i), ry(i+1));
            y{i+1} = reshape(y{i+1}, ry(i+1), n(i+1), ry(i+2));
        end;
    end;

    if (~isempty(aux))
        for j=1:Raux
            phiauxl{i+1,j} = craux1{j}(ilocl{i+1}, :);
            phiauxr{i+1,j} = craux2{j}(:, ilocr{i+1});
        end;
    end;

    i = i+dir;

    % Check the convergence, restart
    if ((i==d)||(i==0))
        if (verb>0)
            if (exist('xtru'))
                ytest=formtensor(y,mid_inv,d,ry,n);
                fprintf('=greedy_cross= swp=%d, max_dx=%3.3e, max_rank=%d, cum#evals=%d, err_tru=%3.3e\n', swp, max_dx, max(ry), evalcnt, norm(ytest-xtru)/norm(xtru));
            else
                fprintf('=greedy_cross= swp=%d, max_dx=%3.3e, max_rank=%d, cum#evals=%d\n', swp, max_dx, max(ry), evalcnt);
            end;
        end;

        if (max_dx<tol_exit)
            break;
        end;

        max_dx = 0;
        i=1;
        swp = swp+1;
    end;
end;

y=formtensor(y,mid_inv,d,ry,n);
end

function [y]=formtensor(y,mid_inv,d,ry,n)
% Merge mid_inv: y:=inv(L)*y*inv(U)
for i=1:d
    y{i} = reshape(y{i}, ry(i), n(i)*ry(i+1));
    y{i} = mid_inv{i,2}*y{i};
    y{i} = reshape(y{i}, ry(i)*n(i), ry(i+1));
    y{i} = y{i}*mid_inv{i+1,1};
    y{i} = reshape(y{i}, ry(i), n(i), ry(i+1));
end;
y = cell2core(tt_tensor,y);
end

function [y]=autovecfun(fun, J, vec)
if (vec)
    y = fun(J); % User can cast M x d -> M x 1
else
    % They can't -- run a loop manually
    sz = size(J,1);
    y = zeros(sz, 1);
    for j=1:sz
        y(j) = fun(J(j,:));
    end;
end;
end

function [J]=indexmerge(varargin)
% Merges two or three indices in the little-endian manner
sz1 = max(size(varargin{1},1),1);
sz2 = max(size(varargin{2},1),1);
sz3 = 1;
if (nargin>2) % Currently allows only 3
    sz3 = max(size(varargin{3}, 1), 1);
end;
% J1 goes to the fastest index, just copy it
J1 = repmat(varargin{1}, sz2*sz3, 1);
% J2 goes to the middle
J2 = reshape(varargin{2}, 1, []);
J2 = repmat(J2, sz1, 1); % now sz1 ones will be the fastest
J2 = reshape(J2, sz1*sz2, []);
J2 = repmat(J2, sz3, 1);
J = [J1,J2];
if (nargin>2)
    % J3 goes to the slowest
    J3 = reshape(varargin{3}, 1, []);
    J3 = repmat(J3, sz1*sz2, 1); % now sz1 ones will be the fastest
    J3 = reshape(J3, sz1*sz2*sz3, []);
    J = [J,J3];
end;
end

author: oseledets

repo: TT-Toolbox


text:
%This directory contains the files which can be used for various ways of
%visualizing data and analyses.
%
%PLOT_CHANNEL - Plot the classwise averages of one channel of 1D (time OR
% frequency) or 2D data (time x frequency)
%PLOT_SCOREMATRIX - Visualizes score matrix of e.g. r-values
%
%GRID_PLOT - Classwise averaged epochs in a grid layout
%GRID_ADDBARS - Add a colorbar to the bottom of each subplot of a grid plot
%GRID_MARKBYBOX - Draw error or percentiles box inside grid plots
%GRID_MARKINTERVAL - Draw time interval inside grid plots as colored patches
%GRID_MARKRANGE -  Mark range on y axis  inside grid plots as horizontal lines
%GRID_MARKTIMEBOX - Draw time intervals inside grid plots with vertical lines
%GRID_MARKTIMEPOINT -  Draw a time point inside grid plots as a vertical line
%
%MNT_ADAPTMONTAGE - Adapts an electrode montage to another electrode set
%MNT_RESTRICTMONTAGE - Restrict an electrode montage to a subset of channels
%MNT_RESTRICTNIRSMONTAGE - restrict a NIRS montage to a subset of channels
%MNT_GETNIRSMONTAGE - Get montage with NIRS channel positions
%MNT_SCALPTOGRID - Montage for grid plot with boxes at scalp locations
%MNT_SETELECTRODEPOSITIONS - Electrode positions of standard named channels
%MNT_SETGRID - Define a new eletrode grid layout for an electrode montage
%
%PLOT_SCALP - Display a vector as scalp topography
%PLOT_SCALPEVOLUTION - Draws scalp topographies for specified intervals
%PLOT_SCALPEVOLUTIONPLUSCHANNEL - Display evolution of scalp topographies
%PLOT_SCALPPATTERNSPLUSCHANNEL - Display Classwise topographies and one channel
%PLOT_SCALPLOADING - Visualize single channel loadings of a weight vector
%PLOT_SCALPOUTLINE - Schematic sketch of the scalp
%PLOT_SCALPPATTERN - Display average trials as scalp topography
%
%PLOT_CSPANALYSIS - Show CSP filters and patterns as topographies
%
%IMAGE_LOCALCOLORMAP - Display image without using the figure's colormap
%
%PLOT_BARWITHANTENNA - Bar plot with antennas, e.g., to show the SEM

author: bbci


text:
typeOfMaterial = 'ISOTROPIC';
interpolation = 'SIMP_P3';
dim = '2D';

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [y,yp] = fatan(x)
% FATAN  Arctangent function with optional derivative.
%        [Y,YP] = FATAN(X) returns Y=atan(X) and
%        (optionally) YP = 1/(1+X^2)
%
%
y = atan(x);
if nargout == 2
    yp = 1/(1+x^2);
end

repo: Aero-matlab


text:
function y = powcone(x,y,z,alpha)
%POWCONE Defines a power cone x^alpha y^(1-alpha) > norm(z)
%
% Input
%    x,y,z     : SDPVAR objects or doubles
%    alpha     : scalar double 0<=alpha<=1
%
% Example
%    F = pcone(x,y,z,alpha)
%
% An alternative syntax with only one argument is also possible
%    F = pcone(Z) # where Z = [x;y;z;alpha]
%
% A vectorized version is also possible
%    F = pcone(Z) # where Z(:,i) = [x;y;z;alpha]
%
% See also @sdpvar/CONE, @sdpvar/EXPCONE

if nargin == 1
    [n,m] = size(x);
    if n < 4
        error('x must be a vector or matrix of height 4 or more ')
    end
    y = x;
    if m > 1
        y.typeflag = 58;
    else
        y.typeflag = 20;
    end
    y = lmi(y);
else
    y = powcone([x;y;reshape(z,[],1);alpha]);
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function masks = char2img( strings, h, pad )
% Convert ascii text to a binary image using pre-computed templates.
%
% Input strings can only contain standard characters (ascii character
% 32-126). First time char2img() is ever called with a given height h,
% txt2img() is used to create a template for each ascii character. All
% subsequent calls to to char2img() with the given height are very fast as
% the pre-computed templates are used and no display/screen capture is
% needed.
%
% USAGE
%  masks = char2img( strings, h, [pad] )
%
% INPUTS
%  strings  - {1xn} text string or cell array of text strings to convert
%  h        - font height in pixels
%  pad      - [0] amount of extra padding between chars
%
% OUTPUTS
%  masks    - {1xn} binary image masks of height h for each string
%
% EXAMPLE
%  mask=char2img('hello world',50); im(mask{1})
%  mask=char2img(num2str(pi),35); im(mask{1})
%
% See also TXT2IMG
%
% Piotr's Computer Vision Matlab Toolbox      Version 2.65
% Copyright 2014 Piotr Dollar.  [pdollar-at-gmail.com]
% Licensed under the Simplified BSD License [see external/bsd.txt]

% load or create character templates (or simply store persistently)
persistent chars h0
if(isempty(h0) || h0~=h)
  pth=fileparts(mfilename('fullpath'));
  fName=sprintf('%s/private/char2img_h%03i.mat',pth,h); h0=h;
  if(exist(fName,'file')), load(fName); else
    chars=char(32:126); chars=num2cell(chars);
    chars=txt2img(chars,h,{'Interpreter','none'});
    save(fName,'chars');
  end
end
% add padding to chars
if(nargin<3 || isempty(pad)), pad=0; end
charsPad=chars; if(pad), pad=ones(h,pad,'uint8');
  for i=1:length(chars), charsPad{i}=[pad chars{i} pad]; end; end
% create actual string using templates
if(~iscell(strings)), strings={strings}; end
n=length(strings); masks=cell(1,n);
for s=1:n, str=strings{s};
  str(str<32 | str>126)=32; str=uint8(str-31);
  M=[charsPad{str}]; masks{s}=M;
end

author: pdollar

repo: toolbox

save_path: github-repos/MATLAB/pdollar-toolbox


text:
function PrintFeature( input_A, input_B, all, back, outPath)
%PrintFeature print the feature correspondence.
%   Detailed explanation goes here
    fileListA = dir(input_A);
    fileListA = fileListA(3:length(fileListA));
    fileListB = dir(input_B);
    fileListB = fileListB(3:length(fileListB));
    nFrames = length(fileListA);
    if nFrames ~= size(all, 1)
        error('frame number not match');
    end
    for frameIndex = 1:nFrames
        IA = imread([input_A fileListA(frameIndex).name]);
        IB = imread([input_B fileListB(frameIndex).name]);
        pa_all = squeeze(all(frameIndex, :, 1:2));
        pb_all = squeeze(all(frameIndex, :, 3:4));
        pa_back = squeeze(back(frameIndex, :, 1:2));
        pb_back = squeeze(back(frameIndex, :, 3:4));

        IA = DrawFeature(IA, pa_all, 'red', 'o');
        IB = DrawFeature(IB, pb_all, 'red', 'o');

        IA = DrawFeature(IA, pa_back, 'green', 's');
        IB = DrawFeature(IB, pb_back, 'green', 's');
        if ~exist(outPath, 'dir')
            mkdir(outPath);
        end
        imwrite(IA, [outPath '/A' int2str(frameIndex) '.jpg']);
        imwrite(IB, [outPath '/B' int2str(frameIndex) '.jpg']);
    end
end

author: SuTanTank

repo: VideoStitchingViaShakinessRemoving


text:
function dicomrt_contourslice(inputdose,dose_xmesh,dose_ymesh,dose_zmesh,axis,plane,type,VOI,voi2plot)
% dicomrt_contourslice(inputdose,dose_xmesh,dose_ymesh,dose_zmesh,axis,plane,type,VOI,voi2plot)
%
% Show 3D dose distributions.
%
% inputdose contains the 3D dose distribution
% inputdose can be an rtplan generated or a monte carlo generated 3D dose distribution or a dose difference.
%
% dose_xmesh,dose_ymesh,dose_zmesh are x-y-z coordinates of the center of the dose-pixel
% axis is the plane to be used for slicing the 3D dose matrix (can be x,y or z. Key insensitive).
% plane is a vector containing the slice number of the slices to be displayed
% type OPTIONAL parameter which determines whether display the slices in the same
%      plot (type==0) or in independent subplots (type~=0)
% VOI is an OPTIONAL cell array which contain the patients VOIs as read by dicomrt_loadVOI
% voi2plot is an OPTIONAL vector pointing to the number of VOIs to be displayed
%
% Example:
%
% dicomrt_contourslice(A,xmesh,ymesh,zmesh,'z',[10 15 25], A_voi,[1 2 3])
%
% plot a figure that contains: 3 plots with inputdose sliced at voxel number 10 15 and 25
% All plots show also VOIs # 1 2 and 3.
%
% See also dicomrt_loadmcdose
%
% Copyright (C) 2002 Emiliano Spezi (emiliano.spezi@physics.org)

% Check number of argument and set-up some parameters and variables
error(nargchk(6,9,nargin))

% Check type
if exist('type')~=1
    type=0;
end

% Check case and set-up some parameters and variables
[dose_temp,type_dose,doselabel,PatientPosition]=dicomrt_checkinput(inputdose);
dose=dicomrt_varfilter(dose_temp);

% Check axis
if ischar(axis) ~=1
    error('dicomrt_contourslice: Axis is not a character. Exit now!')
elseif axis=='x' | axis=='X'
    dir=1;
elseif axis=='y' | axis=='Y'
    dir=2;
elseif axis=='z' | axis=='Z'
    dir=3;
else
    error('dicomrt_contourslice: Axis can only be X Y or Z. Exit now!')
end

% rebuilding mesh matrices
[rtdose_xmesh,rtdose_ymesh,rtdose_zmesh]=dicomrt_rebuildmatrix(dose_xmesh,dose_ymesh,dose_zmesh);

% Initialize plot parameters:
color = char('r','g','b','c','m','y','k','w');
line = char('-','--',':','-.');
marker = char('+','o','*','x','s','d','^');
width = [0.5,1.0,1.5,2.0,2.5,3.0];

% Rendering 3D dose
figure
set(gcf,'Name',['dicomrt_contourslice: ',inputname(1)]);
hold

if length(plane)==1
    if dir==1
        p = slice(rtdose_xmesh,rtdose_ymesh,rtdose_zmesh,dose,dose_xmesh(plane),[],[]);
    elseif dir==2
        p = slice(rtdose_xmesh,rtdose_ymesh,rtdose_zmesh,dose,[],dose_ymesh(plane),[]);
    elseif dir==3
        p = slice(rtdose_xmesh,rtdose_ymesh,rtdose_zmesh,dose,[],[],dose_zmesh(plane));
    else
        error('dicomrt_contourslice: Error while rendering. Exit now');
    end
    view(30,30);
    set(p,'FaceColor','interp','EdgeColor','none','DiffuseStrength',.8);
    xlabel('X axis (cm)')
    ylabel('Y axis (cm)')
    zlabel('Z axis (cm)')
    if exist('VOI')==1 % VOI is defined
        dicomrt_rendervoi(VOI,voi2plot,1,1,1);
    end
    grid on
else
    if type~=0
        for i=1:length(plane)
            subplot(2,round((length(plane)/2)),i);
            if dir==1
                p = slice(rtdose_xmesh,rtdose_ymesh,rtdose_zmesh,dose,dose_xmesh(plane(i)),[],[]);
                set(p,'FaceColor','interp','EdgeColor','none','DiffuseStrength',.8);
            elseif dir==2
                p = slice(rtdose_xmesh,rtdose_ymesh,rtdose_zmesh,dose,[],dose_ymesh(plane(i)),[]);
                set(p,'FaceColor','interp','EdgeColor','none','DiffuseStrength',.8);
            elseif dir==3
                p = slice(rtdose_xmesh,rtdose_ymesh,rtdose_zmesh,dose,[],[],dose_zmesh(plane(i)));
                set(p,'FaceColor','interp','EdgeColor','none','DiffuseStrength',.8);
            else
                error('dicomrt_contourslice: Error while rendering. Exit now');
            end
            view(30,30);
            set(p,'FaceColor','interp','EdgeColor','none','DiffuseStrength',.8);
            xlabel('X axis (cm)')
            ylabel('Y axis (cm)')
            zlabel('Z axis (cm)')
            if exist('VOI')==1 % VOI is defined
                hold
                dicomrt_rendervoi(VOI,voi2plot,1,1,1);
            end
        end
    else
        for i=1:length(plane)
            if dir==1
                p = slice(rtdose_xmesh,rtdose_ymesh,rtdose_zmesh,dose,dose_xmesh(plane(i)),[],[]);
                set(p,'FaceColor','interp','EdgeColor','none','DiffuseStrength',.8);
            elseif dir==2
                p = slice(rtdose_xmesh,rtdose_ymesh,rtdose_zmesh,dose,[],dose_ymesh(plane(i)),[]);
                set(p,'FaceColor','interp','EdgeColor','none','DiffuseStrength',.8);
            elseif dir==3
                p = slice(rtdose_xmesh,rtdose_ymesh,rtdose_zmesh,dose,[],[],dose_zmesh(plane(i)));
                set(p,'FaceColor','interp','EdgeColor','none','DiffuseStrength',.8);
            else
                error('dicomrt_contourslice: Error while rendering. Exit now');
            end
            view(30,30);
            set(p,'FaceColor','interp','EdgeColor','none','DiffuseStrength',.8);
            xlabel('X axis (cm)')
            ylabel('Y axis (cm)')
            zlabel('Z axis (cm)')
        end
        if exist('VOI')==1 % VOI is defined
            dicomrt_rendervoi(VOI,voi2plot,1,1,1);
        end
        grid on
    end
end

author: cerr

repo: CERR