text:
%% Analyzing Neural Time Series Data
% Matlab code for Chapter 34
% Mike X Cohen
%
% This code accompanies the book, titled "Analyzing Neural Time Series Data"
% (MIT Press). Using the code without following the book may lead to confusion,
% incorrect data analyses, and misinterpretations of results.
% Mike X Cohen assumes no responsibility for inappropriate or incorrect use of this code.

%% load sample data

load sampleEEGdata

% note: Most of these figures take a while to generate. Have patience!

%% extract TF power (create data that are used for the rest of this chapter)

% definitions, selections...
chan2use = 'fcz';

min_freq = 3;
max_freq = 30;
num_frex = 20;

% define wavelet parameters
time = -1:1/EEG.srate:1;
frex = logspace(log10(min_freq),log10(max_freq),num_frex);
s    = logspace(log10(3),log10(10),num_frex)./(2*pi*frex);

% definte convolution parameters
n_wavelet            = length(time);
n_data               = EEG.pnts*EEG.trials;
n_convolution        = n_wavelet+n_data-1;
n_conv_pow2          = pow2(nextpow2(n_convolution));
half_of_wavelet_size = (n_wavelet-1)/2;

% note that you don't need the wavelet itself, you need the FFT of the wavelet
wavelets = zeros(num_frex,n_conv_pow2);
for fi = 1:num_frex
    wavelets(fi,:)  = fft( sqrt(1/(s(fi)*sqrt(pi))) * exp(2*1i*pi*frex(fi).*time) .* exp(-time.^2./(2*(s(fi)^2))) , n_conv_pow2 );
end

% get FFT of data
eegfft = fft(reshape(EEG.data(strcmpi(chan2use,{EEG.chanlocs.labels}),:,:),1,EEG.pnts*EEG.trials),n_conv_pow2);

% initialize
eegpower = zeros(num_frex,EEG.pnts,EEG.trials); % frequencies X time X trials
eegphase = zeros(num_frex,EEG.pnts,EEG.trials); % frequencies X time X trials

% loop through frequencies and compute synchronization
for fi=1:num_frex

    % convolution
    eegconv = ifft(wavelets(fi,:).*eegfft);
    eegconv = eegconv(1:n_convolution);
    eegconv = eegconv(half_of_wavelet_size+1:end-half_of_wavelet_size);

    % reshape to time X trials
    eegpower(fi,:,:) = abs(reshape(eegconv,EEG.pnts,EEG.trials)).^2;
    eegphase(fi,:,:) = exp(1i*angle(reshape(eegconv,EEG.pnts,EEG.trials)));
end

% remove edge artifacts
time_s = dsearchn(EEG.times',-500);
time_e = dsearchn(EEG.times',1200);

eegpower = eegpower(:,time_s:time_e,:);
tftimes  = EEG.times(time_s:time_e);
nTimepoints = numel(tftimes);

%% Figure 34.1

voxel_pval   = 0.01;
cluster_pval = 0.05;

% note: try to use 1000 or more permutations for real data
n_permutes = 1000;

baseidx(1) = dsearchn(tftimes',-500);
baseidx(2) = dsearchn(tftimes',-100);

% compute actual t-test of difference
realbaselines = squeeze(mean(eegpower(:,baseidx(1):baseidx(2),:),2));
realmean      = 10*log10(bsxfun(@rdivide, mean(eegpower,3), mean(realbaselines,2)));

% initialize null hypothesis matrices
permuted_maxvals = zeros(n_permutes,2,num_frex);
permuted_vals    = zeros(n_permutes,num_frex,numel(tftimes));
max_clust_info   = zeros(n_permutes,1);

for permi=1:n_permutes
    cutpoint = randsample(2:nTimepoints-diff(baseidx)-2,1);
    permuted_vals(permi,:,:) = 10*log10(bsxfun(@rdivide,mean(eegpower(:,[cutpoint:end 1:cutpoint-1],:),3),mean(realbaselines,2)) );
    % btw, using bsxfun instead of repmat increases the speed
    % of this loop by a factor of ~5.
end

zmap = (realmean-squeeze(mean(permuted_vals))) ./ squeeze(std(permuted_vals));
threshmean = realmean;
threshmean(abs(zmap)<norminv(1-voxel_pval))=0;

figure
subplot(221)
contourf(tftimes,frex,realmean,40,'linecolor','none')
axis square
set(gca,'clim',[-3 3],'xlim',[-500 1200])
title('power map')
xlabel('Time (ms)'), ylabel('Frequency (Hz)')

subplot(222)
contourf(tftimes,frex,zmap,40,'linecolor','none')
axis square
set(gca,'clim',[-3 3],'xlim',[-500 1200])
title('unthresholded Z map')
xlabel('Time (ms)'), ylabel('Frequency (Hz)')

subplot(223)
contourf(tftimes,frex,threshmean,40,'linecolor','none')
axis square
set(gca,'clim',[-3 3],'xlim',[-500 1200])
title('Uncorrected power map')
xlabel('Time (ms)'), ylabel('Frequency (Hz)')

% this time, the cluster correction will be done on the permuted data, thus
% making no assumptions about parameters for p-values
for permi = 1:n_permutes

    % for cluster correction, apply uncorrected threshold and get maximum cluster sizes
    fakecorrsz = squeeze((permuted_vals(permi,:,:)-mean(permuted_vals,1)) ./ std(permuted_vals,[],1) );
    fakecorrsz(abs(fakecorrsz)<norminv(1-voxel_pval))=0;

    % get number of elements in largest supra-threshold cluster
    clustinfo = bwconncomp(fakecorrsz);
    max_clust_info(permi) = max([ 0 cellfun(@numel,clustinfo.PixelIdxList) ]); % the zero accounts for empty maps
    % using cellfun here eliminates the need for a slower loop over cells
end

% apply cluster-level corrected threshold
zmapthresh = zmap;
% uncorrected pixel-level threshold
zmapthresh(abs(zmapthresh)<norminv(1-voxel_pval))=0;
% find islands and remove those smaller than cluster size threshold
clustinfo = bwconncomp(zmapthresh);
clust_info = cellfun(@numel,clustinfo.PixelIdxList);
clust_threshold = prctile(max_clust_info,100-cluster_pval*100);

% identify clusters to remove
whichclusters2remove = find(clust_info<clust_threshold);

% remove clusters
for i=1:length(whichclusters2remove)
    zmapthresh(clustinfo.PixelIdxList{whichclusters2remove(i)})=0;
end

subplot(224)
contourf(tftimes,frex,zmapthresh,40,'linecolor','none')
axis square
set(gca,'clim',[-3 3],'xlim',[-500 1200])
title('Cluster-corrected Z map')
xlabel('Time (ms)'), ylabel('Frequency (Hz)')

%% Figure 34.3

voxel_pval = 0.05;
mcc_voxel_pval = 0.05; % mcc = multiple comparisons correction
mcc_cluster_pval = 0.05;

% note: try to use 1000 or more permutations for real data
n_permutes = 1000;

real_condition_mapping = [ -ones(1,floor(EEG.trials/2)) ones(1,ceil(EEG.trials/2)) ];

% compute actual t-test of difference (using unequal N and std)
tnum   = squeeze(mean(eegpower(:,:,real_condition_mapping==-1),3) - mean(eegpower(:,:,real_condition_mapping==1),3));
tdenom = sqrt( (std(eegpower(:,:,real_condition_mapping==-1),0,3).^2)./sum(real_condition_mapping==-1) + (std(eegpower(:,:,real_condition_mapping==1),0,3).^2)./sum(real_condition_mapping==1) );
real_t = tnum./tdenom;

% initialize null hypothesis matrices
permuted_tvals  = zeros(n_permutes,num_frex,nTimepoints);
max_pixel_pvals = zeros(n_permutes,2);
max_clust_info  = zeros(n_permutes,1);

% generate pixel-specific null hypothesis parameter distributions
for permi = 1:n_permutes
    fake_condition_mapping = sign(randn(EEG.trials,1));

    % compute t-map of null hypothesis
    tnum   = squeeze(mean(eegpower(:,:,fake_condition_mapping==-1),3)-mean(eegpower(:,:,fake_condition_mapping==1),3));
    tdenom = sqrt( (std(eegpower(:,:,fake_condition_mapping==-1),0,3).^2)./sum(fake_condition_mapping==-1) + (std(eegpower(:,:,fake_condition_mapping==1),0,3).^2)./sum(fake_condition_mapping==1) );
    tmap   = tnum./tdenom;

    % save all permuted values
    permuted_tvals(permi,:,:) = tmap;

    % save maximum pixel values
    max_pixel_pvals(permi,:) = [ min(tmap(:)) max(tmap(:)) ];

    % for cluster correction, apply uncorrected threshold and get maximum cluster sizes
    % note that here, clusters were obtained by parametrically thresholding
    % the t-maps
    tmap(abs(tmap)<tinv(1-voxel_pval,EEG.trials-1))=0;

    % get number of elements in largest supra-threshold cluster
    clustinfo = bwconncomp(tmap);
    max_clust_info(permi) = max([ 0 cellfun(@numel,clustinfo.PixelIdxList) ]); % notes: cellfun is superfast, and the zero accounts for empty maps
end

% now compute Z-map
zmap = (real_t-squeeze(mean(permuted_tvals,1)))./squeeze(std(permuted_tvals));

figure
subplot(221)
contourf(tftimes,frex,zmap,40,'linecolor','none')
axis square
set(gca,'clim',[-3 3],'xlim',[-500 1200])
title('Unthresholded Z map')
xlabel('Time (ms)'), ylabel('Frequency (Hz)')

% apply uncorrected threshold
subplot(222)
contourf(tftimes,frex,zmap,40,'linecolor','none')
zmapthresh = zmap;
zmapthresh(abs(zmapthresh)<norminv(1-voxel_pval))=false;
zmapthresh=logical(zmapthresh);
hold on
contour(tftimes,frex,zmapthresh,1,'linecolor','k')

axis square
set(gca,'clim',[-3 3],'xlim',[-500 1200])
title('Unthresholded Z map')
xlabel('Time (ms)'), ylabel('Frequency (Hz)')

% apply pixel-level corrected threshold
lower_threshold = prctile(max_pixel_pvals(:,1),    mcc_voxel_pval*100/2);
upper_threshold = prctile(max_pixel_pvals(:,2),100-mcc_voxel_pval*100/2);

zmapthresh = zmap;
zmapthresh(zmapthresh>lower_threshold & zmapthresh<upper_threshold)=0;
subplot(223)
contourf(tftimes,frex,zmapthresh,40,'linecolor','none')
axis square
set(gca,'clim',[-3 3],'xlim',[-500 1200])
title('Pixel-corrected Z map')
xlabel('Time (ms)'), ylabel('Frequency (Hz)')

% apply cluster-level corrected threshold
zmapthresh = zmap;
% uncorrected pixel-level threshold
zmapthresh(abs(zmapthresh)<norminv(1-voxel_pval))=0;
% find islands and remove those smaller than cluster size threshold
clustinfo = bwconncomp(zmapthresh);
clust_info = cellfun(@numel,clustinfo.PixelIdxList);
clust_threshold = prctile(max_clust_info,100-mcc_cluster_pval*100);

% identify clusters to remove
whichclusters2remove = find(clust_info<clust_threshold);

% remove clusters
for i=1:length(whichclusters2remove)
    zmapthresh(clustinfo.PixelIdxList{whichclusters2remove(i)})=0;
end

subplot(224)
contourf(tftimes,frex,zmapthresh,40,'linecolor','none')
axis square
set(gca,'clim',[-3 3],'xlim',[-500 1200])
title('Cluster-corrected Z map')
xlabel('Time (ms)'), ylabel('Frequency (Hz)')

%% Figure 34.4

voxel_pval = 0.01;
mcc_voxel_pval = 0.05; % mcc = multiple comparisons correction
mcc_cluster_pval = 0.05;

% note: try to use 1000 or more permutations for real data
n_permutes = 1000;

rts=zeros(1,EEG.trials);
for ei=1:EEG.trials
    % In this task, the button press always followed the stimulus at
    % time=0. Thus, finding the RT involves finding the latency of the
    % event that occurs after the time=0 event.
    % If you follow a procedure like this in your data, you may need to
    % include special exceptions, e.g., if there was no response or if
    % a non-response marker could have occurred between stimulus and response.
    time0event = find(cell2mat(EEG.epoch(ei).eventlatency)==0);
    rts(ei)    = EEG.epoch(ei).eventlatency{time0event+1};
end

% rank-transform RTs
rtsrank = tiedrank(rts);

% rank-transform power data (must be transformed)
eegpowerreshaped = reshape(eegpower,num_frex*nTimepoints,EEG.trials)';
eegpowerrank = tiedrank(eegpowerreshaped);

% technically, you want to perform a correlation, but a linear least-squares fit provides the same conceptual results
% while being ~15 times faster, and we don't care here about the actual scale of the data. For completeness,
% the following line shows you how to compute the Spearman correlation coefficient
% realcorrs = 1-6*sum((eegpowerrank-repmat(rtsrank',1,size(eegpowerrank,2))).^2)/(EEG.trials*(EEG.trials^2-1));
realcorrs = (rtsrank*rtsrank')\rtsrank*eegpowerrank;
realcorrs = reshape(realcorrs,num_frex,nTimepoints);

% initialize null hypothesis matrices
permuted_rvals  = zeros(n_permutes,num_frex,nTimepoints);
max_pixel_rvals = zeros(n_permutes,2);
max_clust_info  = zeros(n_permutes,1);

% generate pixel-specific null hypothesis parameter distributions
for permi = 1:n_permutes
    fake_rt_mapping = rtsrank(randperm(EEG.trials));

    % compute t-map of null hypothesis
    fakecorrs = (fake_rt_mapping*fake_rt_mapping')\fake_rt_mapping*eegpowerrank;

    % reshape to 2D map for cluster-correction
    fakecorrs = reshape(fakecorrs,num_frex,nTimepoints);

    % save all permuted values
    permuted_rvals(permi,:,:) = fakecorrs;

    % save maximum pixel values
    max_pixel_rvals(permi,:) = [ min(fakecorrs(:)) max(fakecorrs(:)) ];
end

% this time, the cluster correction will be done on the permuted data, thus
% making no assumptions about parameters for p-values
for permi = 1:n_permutes

    % indices of permutations to include in thresholding at this iteration
    perms2use4distribution = true(1,n_permutes);
    perms2use4distribution(permi) = 0;

    % for cluster correction, apply uncorrected threshold and get maximum cluster sizes
    fakecorrsz = squeeze((permuted_rvals(permi,:,:)-mean(permuted_rvals(perms2use4distribution,:,:),1)) ./ std(permuted_rvals(perms2use4distribution,:,:),[],1) );
    fakecorrsz(abs(fakecorrsz)<norminv(1-voxel_pval))=0;

    % get number of elements in largest supra-threshold cluster
    clustinfo = bwconncomp(fakecorrsz);
    max_clust_info(permi) = max([ 0 cellfun(@numel,clustinfo.PixelIdxList) ]); % the zero accounts for empty maps
end

% now compute Z-map
zmap = (realcorrs-squeeze(mean(permuted_rvals,1)))./squeeze(std(permuted_rvals));

figure
subplot(221)
contourf(tftimes,frex,zmap,40,'linecolor','none')
axis square
set(gca,'clim',[-4 4],'xlim',[-500 1200])
title('Unthresholded Z map')
xlabel('Time (ms)'), ylabel('Frequency (Hz)')

% apply uncorrected threshold
zmapthresh = zmap;
zmapthresh(abs(zmapthresh)<norminv(1-voxel_pval))=0;
zmapthresh=logical(zmapthresh);
subplot(222)
contourf(tftimes,frex,zmap,40,'linecolor','none')
hold on
contour(tftimes,frex,zmapthresh,1,'linecolor','k')
axis square
set(gca,'clim',[-4 4],'xlim',[-500 1200])
title('Uncorrected thresholded Z map')
xlabel('Time (ms)'), ylabel('Frequency (Hz)')

% apply pixel-level corrected threshold
lower_threshold = prctile(max_pixel_rvals(:,1),    mcc_voxel_pval*100/2);
upper_threshold = prctile(max_pixel_rvals(:,2),100-mcc_voxel_pval*100/2);

zmapthresh = zmap;
zmapthresh(realcorrs>lower_threshold & realcorrs<upper_threshold)=0;
subplot(223)
contourf(tftimes,frex,zmapthresh,40,'linecolor','none')
axis square
set(gca,'clim',[-4 4],'xlim',[-500 1200])
title('Pixel-corrected Z map')
xlabel('Time (ms)'), ylabel('Frequency (Hz)')

% apply cluster-level corrected threshold
zmapthresh = zmap;
% uncorrected pixel-level threshold
zmapthresh(abs(zmapthresh)<norminv(1-voxel_pval))=0;
% find islands and remove those smaller than cluster size threshold
clustinfo = bwconncomp(zmapthresh);
clust_info = cellfun(@numel,clustinfo.PixelIdxList);
clust_threshold = prctile(max_clust_info,100-mcc_cluster_pval*100);

% identify clusters to remove
whichclusters2remove = find(clust_info<clust_threshold);

% remove clusters
for i=1:length(whichclusters2remove)
    zmapthresh(clustinfo.PixelIdxList{whichclusters2remove(i)})=0;
end

subplot(224)
contourf(tftimes,frex,zmapthresh,40,'linecolor','none')
axis square
set(gca,'clim',[-4 4],'xlim',[-500 1200])
title('Cluster-corrected Z map')
xlabel('Time (ms)'), ylabel('Frequency (Hz)')

%% Figure 34.5

chan2use = 'o1';
time2use = dsearchn(EEG.times',[0 250]');
freq2use = dsearchn(frex',10);

eegfft = fft(reshape(EEG.data(strcmpi(chan2use,{EEG.chanlocs.labels}),:,:),1,EEG.pnts*EEG.trials),n_conv_pow2);
eegconv = ifft(wavelets(freq2use,:).*eegfft);
eegconv = eegconv(1:n_convolution);
eegconv = eegconv(half_of_wavelet_size+1:end-half_of_wavelet_size);

% reshape to time X trials
temp  = abs(reshape(eegconv,EEG.pnts,EEG.trials)).^2;
o1power = zscore(mean(temp(time2use(1):time2use(2),:),1));

% define covariates (RT and trial number)
X = [ zscore(rts') o1power' ]';
eegpowerrank = tiedrank(eegpowerreshaped)';

figure
subplot(311)
imagesc(X)

subplot(212)
imagesc(eegpowerrank)

%% Figure 34.6

voxel_pval = 0.01;
mcc_cluster_pval = 0.05;

% note: try to use 1000 or more permutations for real data
n_permutes = 1000;

realbeta = (X*X')\X*eegpowerrank';
realbeta = reshape(realbeta,[2 num_frex nTimepoints]);

% initialize null hypothesis matrices
permuted_bvals = zeros(n_permutes,2,num_frex,nTimepoints);
max_clust_info = zeros(n_permutes,2);

% generate pixel-specific null hypothesis parameter distributions
for permi = 1:n_permutes

    % randomly shuffle trial order
    fakeX = X(:,randperm(EEG.trials));

    % compute beta-map of null hypothesis
    fakebeta = (fakeX*fakeX')\fakeX*eegpowerrank';

    % reshape to 2D map for cluster-correction
    fakebeta = reshape(fakebeta,[2 num_frex nTimepoints ]);

    % save all permuted values
    permuted_bvals(permi,:,:,:) = fakebeta;
end

% this time, the cluster correction will be done on the permuted data, thus
% making no assumptions about parameters for p-values
for permi = 1:n_permutes

    for testi=1:2
        % for cluster correction, apply uncorrected threshold and get maximum cluster sizes
        fakecorrsz = squeeze((permuted_bvals(permi,testi,:,:)-mean(permuted_bvals(:,testi,:,:),1)) ./ std(permuted_bvals(:,testi,:,:),[],1) );
        fakecorrsz(abs(fakecorrsz)<norminv(1-voxel_pval))=0;
        % get number of elements in largest supra-threshold cluster
        clustinfo = bwconncomp(fakecorrsz);
        max_clust_info(permi,testi) = max([ 0 cellfun(@numel,clustinfo.PixelIdxList) ]); % the zero accounts for empty maps
    end
end

figure
for testi=1:2

    % now compute Z-map
    zmap = (squeeze(realbeta(testi,:,:))-squeeze(mean(permuted_bvals(:,testi,:,:),1))) ./ squeeze(std(permuted_bvals(:,testi,:,:),[],1));

    subplot(2,3,1+(testi-1)*3)
    contourf(tftimes,frex,zmap,40,'linecolor','none')
    axis square
    set(gca,'clim',[-3 3],'xlim',[-500 1200])
    title('Unthresholded Z map')
    xlabel('Time (ms)'), ylabel('Frequency (Hz)')

    % apply uncorrected threshold
    zmapthresh = zmap;
    zmapthresh(abs(zmapthresh)<norminv(1-voxel_pval))=0;
    subplot(2,3,2+(testi-1)*3)
    contourf(tftimes,frex,zmapthresh,40,'linecolor','none')
    axis square
    set(gca,'clim',[-3 3],'xlim',[-500 1200])
    title('Uncorrected thresholded Z map')
    xlabel('Time (ms)'), ylabel('Frequency (Hz)')

    % apply cluster-level corrected threshold
    zmapthresh = zmap;
    % uncorrected pixel-level threshold
    zmapthresh(abs(zmapthresh)<norminv(1-voxel_pval))=0;
    % find islands and remove those smaller than cluster size threshold
    clustinfo = bwconncomp(zmapthresh);
    clust_info = cellfun(@numel,clustinfo.PixelIdxList);
    clust_threshold = prctile(max_clust_info(:,testi),100-mcc_cluster_pval*100);

    % identify clusters to remove
    whichclusters2remove = find(clust_info<clust_threshold);

    % remove clusters
    for i=1:length(whichclusters2remove)
        zmapthresh(clustinfo.PixelIdxList{whichclusters2remove(i)})=0;
    end

    subplot(2,3,3+(testi-1)*3)
    contourf(tftimes,frex,zmapthresh,40,'linecolor','none')
    axis square
    set(gca,'clim',[-3 3],'xlim',[-500 1200])
    title('Cluster-corrected Z map')
    xlabel('Time (ms)'), ylabel('Frequency (Hz)')
end

%% Figure 34.7

a = rand(10000,1);
b = rand(10000,1);

figure
clear h

subplot(221)
[y,x] = hist(a,50);
h(1)=bar(x,y,'histc');
set(gca,'xlim',[-.05 1.05])

subplot(222)
[y,x] = hist(b,50);
h(2)=bar(x,y,'histc');
set(gca,'xlim',[-.05 1.05])

subplot(212)
[y,x] = hist(atanh(a-b),50);
h(3)=bar(x,y,'histc');
set(gca,'xlim',[-2 2])
title('ITPC differences')
xlabel('difference value'), ylabel('Count')

set(h,'linestyle','none','facecolor','k')

%% Figure 34.8

% The code to produce this figure is presented in chapter19.m, between the
% cells for figures 19.6 and 19.7. To generate figure 34.8, you will need to run
% the code for figures 19.2-6.

%% Figure 34.9

voxel_pval   = 0.01;
cluster_pval = 0.05;

% note: try to use 1000 or more permutations for real data
n_permutes = 1000;

% compute actual t-test of difference
realitpc = squeeze(abs(mean(eegphase,3)));

% initialize null hypothesis matrices
permuted_maxvals = zeros(n_permutes,2,num_frex);
permuted_vals    = zeros(n_permutes,num_frex,EEG.pnts);
max_clust_info   = zeros(n_permutes,1);
eegtemp          = zeros(size(realitpc));

for permi=1:n_permutes
    for triali=1:EEG.trials
        cutpoint = randsample(2:nTimepoints-2,1);
        eegtemp(:,:,triali) = eegphase(:,[cutpoint:end 1:cutpoint-1],triali);
    end
    permuted_vals(permi,:,:) = squeeze(abs(mean(eegtemp,3)));

    % note: the following lines produce fairly similar results as the loop above
    % cutpoint = randsample(2:nTimepoints-2,1);
    %permuted_vals(permi,:,:) = squeeze(abs(mean(eegphase(:,[cutpoint:end 1:cutpoint-1],:),3)));
end

zmap = (realitpc-squeeze(mean(permuted_vals))) ./ squeeze(std(permuted_vals));
threshmean = realitpc;
threshmean(abs(zmap)<norminv(1-voxel_pval))=0;

figure
subplot(221)
contourf(EEG.times,frex,realitpc,40,'linecolor','none')
axis square
set(gca,'clim',[0 .5],'xlim',[-200 1000])
title('power map')
xlabel('Time (ms)'), ylabel('Frequency (Hz)')

subplot(222)
contourf(EEG.times,frex,zmap,40,'linecolor','none')
axis square
set(gca,'clim',[-5 5],'xlim',[-200 1000])
title('unthresholded Z map')
xlabel('Time (ms)'), ylabel('Frequency (Hz)')

subplot(223)
contourf(EEG.times,frex,threshmean,40,'linecolor','none')
axis square
set(gca,'clim',[0 .5],'xlim',[-200 1000])
title('Uncorrected power map')
xlabel('Time (ms)'), ylabel('Frequency (Hz)')

% this time, the cluster correction will be done on the permuted data, thus
% making no assumptions about parameters for p-values
for permi = 1:n_permutes

    % for cluster correction, apply uncorrected threshold and get maximum cluster sizes
    fakecorrsz = squeeze((permuted_vals(permi,:,:)-mean(permuted_vals,1)) ./ std(permuted_vals,[],1) );
    fakecorrsz(abs(fakecorrsz)<norminv(1-voxel_pval))=0;

    % get number of elements in largest supra-threshold cluster
    clustinfo = bwconncomp(fakecorrsz);
    max_clust_info(permi) = max([ 0 cellfun(@numel,clustinfo.PixelIdxList) ]); % the zero accounts for empty maps
    % using cellfun here eliminates the need for a slower loop over cells
end

% apply cluster-level corrected threshold
zmapthresh = realitpc;
% uncorrected pixel-level threshold
zmapthresh(abs(zmap)<norminv(1-voxel_pval))=0;
% find islands and remove those smaller than cluster size threshold
clustinfo = bwconncomp(zmapthresh);
clust_info = cellfun(@numel,clustinfo.PixelIdxList);
clust_threshold = prctile(max_clust_info,100-cluster_pval*100);

% identify clusters to remove
whichclusters2remove = find(clust_info<clust_threshold);

% remove clusters
for i=1:length(whichclusters2remove)
    zmapthresh(clustinfo.PixelIdxList{whichclusters2remove(i)})=0;
end

subplot(224)
contourf(EEG.times,frex,zmapthresh,40,'linecolor','none')
axis square
set(gca,'clim',[0 .5],'xlim',[-200 1000])
title('Cluster-corrected Z map')
xlabel('Time (ms)'), ylabel('Frequency (Hz)')

%% Figure 33.5

% The code for figures 33.5/6 are presented here and in the next cell. You
% will need first to run the code for figure 34.3 to run this code.

% compute actual t-test of difference (using unequal N and std)
tnum   = squeeze(mean(eegpower(4,400,real_condition_mapping==-1),3) - mean(eegpower(4,400,real_condition_mapping==1),3));
tdenom = sqrt( (std(eegpower(4,400,real_condition_mapping==-1),0,3).^2)./sum(real_condition_mapping==-1) + (std(eegpower(4,400,real_condition_mapping==1),0,3).^2)./sum(real_condition_mapping==1) );
real_t = tnum./tdenom;

n_permutes = round(linspace(100,3000,200));
zvals   = zeros(size(n_permutes));

for grandpermi=1:length(n_permutes)

    % initialize null hypothesis matrices
    permuted_tvals  = zeros(n_permutes(grandpermi),1);

    % generate pixel-specific null hypothesis parameter distributions
    for permi = 1:n_permutes(grandpermi)
        fake_condition_mapping = sign(randn(EEG.trials,1));
        % compute t-map of null hypothesis
        tnum   = squeeze(mean(eegpower(4,400,fake_condition_mapping==-1),3)-mean(eegpower(4,400,fake_condition_mapping==1),3));
        tdenom = sqrt( (std(eegpower(4,400,fake_condition_mapping==-1),0,3).^2)./sum(fake_condition_mapping==-1) + (std(eegpower(4,400,fake_condition_mapping==1),0,3).^2)./sum(fake_condition_mapping==1) );
        % save all permuted values
        permuted_tvals(permi) = tnum./tdenom;
    end

    zvals(grandpermi) = (real_t-mean(permuted_tvals))/std(permuted_tvals);

    % display progress
    if mod(grandpermi,20)==0, disp([ 'Metapermutation #' num2str(grandpermi) ]); end
end

figure
subplot(211)
plot(n_permutes,zvals)
xlabel('Number of iterations'), ylabel('Z-value')

subplot(212)
hist(zvals,30)
xlabel('Z-value at different runs of permutation test'), ylabel('Count')

%% Figure 33.6

% compute actual t-test of difference (using unequal N and std)
tnum   = squeeze(mean(eegpower(:,:,real_condition_mapping==-1),3) - mean(eegpower(:,:,real_condition_mapping==1),3));
tdenom = sqrt( (std(eegpower(:,:,real_condition_mapping==-1),0,3).^2)./sum(real_condition_mapping==-1) + (std(eegpower(:,:,real_condition_mapping==1),0,3).^2)./sum(real_condition_mapping==1) );
real_t = tnum./tdenom;

n_permutes = round(linspace(100,3000,200));
zvals = zeros(length(n_permutes),size(tnum,1),size(tnum,2));

for grandpermi=1:length(n_permutes)

    % initialize null hypothesis matrices
    permuted_tvals  = zeros(n_permutes(grandpermi),size(tnum,1),size(tnum,2));

    % generate pixel-specific null hypothesis parameter distributions
    for permi = 1:n_permutes(grandpermi)
        fake_condition_mapping = sign(randn(EEG.trials,1));
        % compute t-map of null hypothesis
        tnum   = squeeze(mean(eegpower(:,:,fake_condition_mapping==-1),3)-mean(eegpower(:,:,fake_condition_mapping==1),3));
        tdenom = sqrt( (std(eegpower(:,:,fake_condition_mapping==-1),0,3).^2)./sum(fake_condition_mapping==-1) + (std(eegpower(:,:,fake_condition_mapping==1),0,3).^2)./sum(fake_condition_mapping==1) );
        % save all permuted values
        permuted_tvals(permi,:,:) = tnum./tdenom;
    end

    zvals(grandpermi,:,:) = (real_t-squeeze(mean(permuted_tvals)))./squeeze(std(permuted_tvals));

    % display progress
    if mod(grandpermi,20)==0, disp([ 'Metapermutation # ' num2str(grandpermi) ]); end
end

figure
subplot(211)
plot(n_permutes,squeeze(zvals(:,4,400)))
subplot(212)
hist(squeeze(zvals(:,4,400)),30)

zvalsall = reshape(zvals,length(n_permutes),size(tnum,1)*size(tnum,2));
figure
plot(mean(zvalsall,1),std(zvalsall,[],1),'.')
set(gca,'xlim',[-3.5 3.5],'ylim',[0 .12])
xlabel('Average Z-statistic')
ylabel('Standard deviation of Z-statistics')

figure
clear h

[~,z0]=min(abs(mean(zvalsall,1)));
[x0,y0]=ind2sub(size(tnum),z0);
[yy,xx]=hist(squeeze(zvals(:,x0,y0)),30);
h(1) = bar(xx,yy,'histc');
hold on
plot([0 0],get(gca,'ylim'),'k:')

[~,z2]=min(abs(mean(zvalsall,1)-2));
[x2,y2]=ind2sub(size(tnum),z2);
[yy,xx]=hist(squeeze(zvals(:,x2,y2)),30);
h(2) = bar(xx,yy,'histc');
plot([2 2],get(gca,'ylim'),'k:')

[~,z3]=min(abs(mean(zvalsall,1)--3));
[x3,y3]=ind2sub(size(tnum),z3);
[yy,xx]=hist(squeeze(zvals(:,x3,y3)),30);
hold on
h(3) = bar(xx,yy,'histc');
plot([-3 -3],get(gca,'ylim'),'k:')

set(h,'linestyle','none')
set(gca,'xlim',[-3.5 3.5])
xlabel('Z-value')
ylabel('Count of possible z-values')

%% end.

author: mikexcohen

repo: AnalyzingNeuralTimeSeries


text:
%% [A3_] = tensor_slices_horizontal(A)
%
function [A3_] = tensor_slices_horizontal(A)
  [n1, n2, n3] = size(A);
  for i = 1:n1
    a = A(i,:,:);
    A3_{i} = reshape(a,[n2 n3]);
  end
end

author: andrewssobral

repo: lrslibrary


text:
function a = imag(a)
%IMAG         imaginary part of hessians
%
%   c = imag(a)
%

% written  04/04/04     S.M. Rump
% modified 04/06/05     S.M. Rump  rounding unchanged
%

  a.x = imag(a.x);
  a.dx = imag(a.dx);
  a.hx = imag(a.hx);

repo: OpenMAS


text:
function [x,y]=henon(n,level,a,b,x0,y0)
%Syntax: [x,y]=henon(n,level,a,b,x0,y0)
%______________________________________
%
% Simulation of the Henon map.
%    x'=1-a*x^2+y
%    y'=b*x
%
% x and y are the simulated time series.
% n is the number of the simulated points.
% level is the noise standard deviation divided by the standard deviation of the
%   noise-free time series. We assume Gaussian noise with zero mean.
% a is the a parameter.
% b is the b parameter.
% x0 is the initial value for x.
% y0 is the initial value for y.
%
%
% Reference:
%
% Henon M (1976):A two-dimensional mapping with a strange attractor.
% Communications in Mathematical Physics 50: 69-77
%
%
% Alexandros Leontitsis
% Department of Education
% University of Ioannina
% 45110 - Dourouti
% Ioannina
% Greece
%
% University e-mail: me00743@cc.uoi.gr
% Lifetime e-mail: leoaleq@yahoo.com
% Homepage: http://www.geocities.com/CapeCanaveral/Lab/1421
%
% 5 Nov 2001

if nargin<1 | isempty(n)==1
    n=500;
else
    % n must be scalar
    if sum(size(n))>2
        error('n must be scalar.');
    end
    % n must be positive
    if n<0
        error('n must be positive.');
    end
    % n must be an integer
    if round(n)-n~=0
        error('n must be an integer.');
    end
end

if nargin<2 | isempty(level)==1
    level=0;
else
    % level must be scalar
    if sum(size(level))>2
        error('level must be scalar.');
    end
    % level must be positive
    if level<0
        error('level must be positive.');
    end
end

if nargin<3 | isempty(a)==1
    a=1.4;
else
    % a must be scalar
    if sum(size(a))>2
        error('a must be scalar.');
    end
end

if nargin<4 | isempty(b)==1
    b=0.3;
else
    % b must be scalar
    if sum(size(b))>2
        error('s must be a scalar.');
    end
end

if nargin<5 | isempty(x0)==1
    x0=0.1;
else
    % x0 must be scalar
    if sum(size(x0))>2
        error('x0 must be scalar.');
    end
end

if nargin<6 | isempty(y0)==1
    y0=0.1;
else
    % y0 must be scalar
    if sum(size(y0))>2
        error('y0 must be scalar.');
    end
end

% Initialize
x(1,1)=1-a*x0^2+b*y0;
y(1,1)=b*x0;

% Simulate
for i=2:n
    x(i,1)=1-a*x(i-1,1)^2+y(i-1,1);
    y(i,1)=b*x(i-1,1);
end

% Add normal white noise
x=x+randn(n,1)*level*std(x);
y=y+randn(n,1)*level*std(y);

author: Sable

repo: mcbench-benchmarks


text:
function [V,Beta,R] = qr_left (A)
%QR_LEFT left-looking Householder QR factorization.
% Example:
%  [V,Beta,R] = qr_left (A)
% See also: cs_demo

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

[m n] = size (A) ;
V = zeros (m,n) ;
Beta = zeros (1,n) ;
R = zeros (m,n) ;
for k = 1:n
    x = A (:,k) ;
    for i = 1:k-1
        v = V (i:m,i) ;
        beta = Beta (i) ;
        x (i:m) = x (i:m) - v * (beta * (v' * x (i:m))) ;
    end
    [v,beta,s] = gallery ('house', x (k:m), 2) ;
    V (k:m,k) = v ;
    Beta (k) = beta ;
    R (1:(k-1),k) = x (1:(k-1)) ;
    R (k,k) = s ;
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
classdef MW13 < PROBLEM
% <multi> <real> <large/none> <constrained>
% Constrained benchmark MOP proposed by Ma and Wang

%------------------------------- Reference --------------------------------
% Z. Ma and Y. Wang, Evolutionary constrained multiobjective optimization:
% Test suite construction and performance comparisons. IEEE Transactions on
% Evolutionary Computation, 2019, 23(6): 972-986.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        %% Default settings of the problem
        function Setting(obj)
            obj.M = 2;
            if isempty(obj.D); obj.D = 15; end
            obj.lower    = zeros(1,obj.D);
            obj.upper    = ones(1,obj.D);
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values and constraint violations
        function Population = Evaluation(obj,varargin)
            X = varargin{1};
            X = max(min(X,repmat(obj.upper,size(X,1),1)),repmat(obj.lower,size(X,1),1));
            z = 1 - exp(-10*(X(:,obj.M:end) - (repmat(obj.M:obj.D,size(X,1),1) - 1)/obj.D).^2);
            g = 1 + sum((1.5 + (0.1/obj.D)*z.^2 - 1.5*cos(2*pi*z)),2);
            PopObj(:,1) = g.*X(:,1)*1.5;
            PopObj(:,2) = g.*(5 - exp(PopObj(:,1)./g) - abs(0.5*sin(3*pi*PopObj(:,1)./g)));
            PopCon(:,1) = (5 - exp(PopObj(:,1)) - 0.5*sin(3*pi*PopObj(:,1)) - PopObj(:,2)).*(5 - (1 + 0.4*PopObj(:,1)) - 0.5*sin(3*pi*PopObj(:,1)) - PopObj(:,2));
            PopCon(:,2) = -(5 - (1 + PopObj(:,1) + 0.5*PopObj(:,1).^2) - 0.5*sin(3*pi*PopObj(:,1)) - PopObj(:,2)).*(5 - (1 + 0.7*PopObj(:,1)) - 0.5*sin(3*pi*PopObj(:,1)) - PopObj(:,2));
            Population  = SOLUTION(X,PopObj,PopCon,varargin{2:end});
            obj.FE      = obj.FE + length(Population);
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(obj,N)
            R(:,1)  = (0:1.5/(N-1):1.5)';
            R(:,2)  = 5 - exp(R(:,1)) - 0.5*abs(sin(3*pi*R(:,1)));
            c1      = (5-exp(R(:,1))-0.5*sin(3*pi*R(:,1))-R(:,2)).*(5-(1+0.4*R(:,1))-0.5*sin(3*pi*R(:,1))-R(:,2));
            invalid = c1>0;
            while any(invalid)
                R(invalid,:) = R(invalid,:).*1.001;
                c1      = (5-exp(R(:,1))-0.5*sin(3*pi*R(:,1))-R(:,2)).*(5-(1+0.4*R(:,1))-0.5*sin(3*pi*R(:,1))-R(:,2));
                invalid = c1>0;
            end
            R = R(NDSort(R,1)==1,:);
        end
        %% Generate the feasible region
        function R = GetPF(obj)
            [x,y] = meshgrid(linspace(0,2,400),linspace(0,4.5,400));
            z     = nan(size(x));
            fes1  = (5-exp(x)-0.5*sin(3*pi*x)-y).*(5-(1+0.4*x)-0.5*sin(3*pi*x)-y) <= 0;
            fes2  = -(5-(1+x+0.5*x.^2)-0.5*sin(3*pi*x)-y).*(5-(1+0.7*x)-0.5*sin(3*pi*x)-y) <= 0;
            z(fes1 & fes2 & exp(x)+abs(0.5*sin(3*pi*x))+y>=5) = 0;
            R = {x,y,z};
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [dPriors,dMu,dSigma,res ] = log_gmm_div(X,Priors,Mu,Sigma )
%LOG_GMM_DIV Derivative of the Log of a Guassian Mixture Model
%
%
%   input -----------------------------------------------------------------
%
%       o X         : (D x N),      Samples
%
%       o Priors    : (1 x K),      Weights
%
%       o Mu        : (D x K),      Means
%
%       o Sigma     : (D x D x K),  Covariance
%
%   output ----------------------------------------------------------------
%
%
%
%

[D,N]   = size(X);

K       = size(Priors,2);

% (N x 1)
L       = gmm_pdf(X,Priors,Mu,Sigma);

dPriors = zeros(1,K);
dMu     = zeros(D,K);
dSigma  = zeros(D,D,K);

% for each parameter k, compute the gradient for each point and sum

k_prior = 1./K;

res = zeros(K,N);
% compute gradient for Priors = responsability factor
for k=1:K
    res(k,:) =  Priors(k) .* gaussPDF(X,Mu(:,k),Sigma(:,:,k)) ./ L;
end

for k=1:K
    %invSigma_k    = inv(Sigma(:,:,k));

   % dPriors(k)    =  (1/N) .* sum(gaussPDF(X,Mu(:,k),Sigma(:,:,k)) ./ L);
    Nk            =  sum(res(k,:));
    dPriors(k)    =  Nk/N;
    dMu(:,k)      =  (1/Nk) .* sum( repmat(res(k,:)',1,D) .* X' ,1)';
    diff          =  X - repmat(Mu(:,K),1,N); % (D x N)
    %                                   (N x N) * (N x D)

    for n=1:N
         dSigma(:,:,k) = dSigma(:,:,k) +  res(k,n) * (X(:,n) - Mu(:,K)) * (X(:,n) - Mu(:,K))';
    end
     dSigma(:,:,k)  =  sqrt((1/Nk) .* dSigma(:,:,k));
   %  (1/Nk) .*  diff * diag(res(k,:)) * diff';

   % dSigma(:,:,k) =  dPriors(k) * 0.5 *(-invSigma_k + invSigma_k * Cov * invSigma_k );
    %dSigma(:,:,k) = 0.5 * ( dSigma(:,:,k) +  dSigma(:,:,k)');
end

%dPriors = dPriors ./ sum(dPriors);

end

author: epfl-lasa


text:
% DESCRIPTION:
%       subscript reorder the sensor points if a binary sensor mask was
%       used for Cartesian sensor mask nearest neighbour interpolation
%
% ABOUT:
%       author      - Bradley Treeby
%       date        - 2nd September 2012
%       last update - 18th August 2014
%
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009-2014 Bradley Treeby and Ben Cox

% This file is part of k-Wave. k-Wave is free software: you can
% redistribute it and/or modify it under the terms of the GNU Lesser
% General Public License as published by the Free Software Foundation,
% either version 3 of the License, or (at your option) any later version.
%
% k-Wave 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 Lesser General Public License for
% more details.
%
% You should have received a copy of the GNU Lesser General Public License
% along with k-Wave. If not, see <http://www.gnu.org/licenses/>.

% update command line status
disp('  reordering Cartesian measurement data...');

if record.p
    % calculate the position for reordering data
    new_col_pos = length(sensor_data.p(1,:)) + 1;

    % append the reordering data
    sensor_data.p(:, new_col_pos) = reorder_index;

    % reorder based on the order_index
    sensor_data.p = sortrows(sensor_data.p, new_col_pos);

    % remove the reordering data
    sensor_data.p = sensor_data.p(:, 1:new_col_pos - 1);
end

if record.p_max
    % append the reordering data
    sensor_data.p_max(:, 2) = reorder_index;

    % reorder based on the order_index
    sensor_data.p_max = sortrows(sensor_data.p_max, 2);

    % remove the reordering data
    sensor_data.p_max = sensor_data.p_max(:, 1);
end

if record.p_min
    % append the reordering data
    sensor_data.p_min(:, 2) = reorder_index;

    % reorder based on the order_index
    sensor_data.p_min = sortrows(sensor_data.p_min, 2);

    % remove the reordering data
    sensor_data.p_min = sensor_data.p_min(:, 1);
end

if record.p_rms
    % append the reordering data
    sensor_data.p_rms(:, 2) = reorder_index;

    % reorder based on the order_index
    sensor_data.p_rms = sortrows(sensor_data.p_rms, 2);

    % remove the reordering data
    sensor_data.p_rms = sensor_data.p_rms(:, 1);
end

if record.u
    % calculate the position for reordering data
    new_col_pos = length(sensor_data.ux(1,:)) + 1;

    % append the reordering data
    sensor_data.ux(:, new_col_pos) = reorder_index;

    % reorder based on the order_index
    sensor_data.ux = sortrows(sensor_data.ux, new_col_pos);

    % remove the reordering data
    sensor_data.ux = sensor_data.ux(:, 1:new_col_pos - 1);

    % repeat for y-dimension if 2D or 3D
    if kgrid.dim > 1
        sensor_data.uy(:, new_col_pos) = reorder_index;
        sensor_data.uy = sortrows(sensor_data.uy, new_col_pos);
        sensor_data.uy = sensor_data.uy(:, 1:new_col_pos - 1);
    end

    % repeat for z-dimension if 3D
    if kgrid.dim > 2
        sensor_data.uz(:, new_col_pos) = reorder_index;
        sensor_data.uz = sortrows(sensor_data.uz, new_col_pos);
        sensor_data.uz = sensor_data.uz(:, 1:new_col_pos - 1);
    end
end

if record.u_max
    % append the reordering data
    sensor_data.ux_max(:, 2) = reorder_index;

    % reorder based on the order_index
    sensor_data.ux_max = sortrows(sensor_data.ux_max, 2);

    % remove the reordering data
    sensor_data.ux_max = sensor_data.ux_max(:, 1);

    % repeat for y-dimension if 2D or 3D
    if kgrid.dim > 1
        sensor_data.uy_max(:, 2) = reorder_index;
        sensor_data.uy_max = sortrows(sensor_data.uy_max, 2);
        sensor_data.uy_max = sensor_data.uy_max(:, 1);
    end

    % repeat for z-dimension if 3D
    if kgrid.dim > 2
        sensor_data.uz_max(:, 2) = reorder_index;
        sensor_data.uz_max = sortrows(sensor_data.uz_max, 2);
        sensor_data.uz_max = sensor_data.uz_max(:, 1);
    end
end

if record.u_min
    % append the reordering data
    sensor_data.ux_min(:, 2) = reorder_index;

    % reorder based on the order_index
    sensor_data.ux_min = sortrows(sensor_data.ux_min, 2);

    % remove the reordering data
    sensor_data.ux_min = sensor_data.ux_min(:, 1);

    % repeat for y-dimension if 2D or 3D
    if kgrid.dim > 1
        sensor_data.uy_min(:, 2) = reorder_index;
        sensor_data.uy_min = sortrows(sensor_data.uy_min, 2);
        sensor_data.uy_min = sensor_data.uy_min(:, 1);
    end

    % repeat for z-dimension if 3D
    if kgrid.dim > 2
        sensor_data.uz_min(:, 2) = reorder_index;
        sensor_data.uz_min = sortrows(sensor_data.uz_min, 2);
        sensor_data.uz_min = sensor_data.uz_min(:, 1);
    end
end

if record.u_rms
    % append the reordering data
    sensor_data.ux_rms(:, 2) = reorder_index;

    % reorder based on the order_index
    sensor_data.ux_rms = sortrows(sensor_data.ux_rms, 2);

    % remove the reordering data
    sensor_data.ux_rms = sensor_data.ux_rms(:, 1);

    % repeat for y-dimension if 2D or 3D
    if kgrid.dim > 1
        sensor_data.uy_rms(:, 2) = reorder_index;
        sensor_data.uy_rms = sortrows(sensor_data.uy_rms, 2);
        sensor_data.uy_rms = sensor_data.uy_rms(:, 1);
    end

    % repeat for z-dimension if 3D
    if kgrid.dim > 2
        sensor_data.uz_rms(:, 2) = reorder_index;
        sensor_data.uz_rms = sortrows(sensor_data.uz_rms, 2);
        sensor_data.uz_rms = sensor_data.uz_rms(:, 1);
    end
end

if record.I
    % calculate the position for reordering data
    new_col_pos = length(sensor_data.Ix(1,:)) + 1;

    % append the reordering data
    sensor_data.Ix(:, new_col_pos) = reorder_index;

    % reorder based on the order_index
    sensor_data.Ix = sortrows(sensor_data.Ix, new_col_pos);

    % remove the reordering data
    sensor_data.Ix = sensor_data.Ix(:, 1:new_col_pos - 1);

    % repeat for y-dimension if 2D or 3D
    if kgrid.dim > 1
        sensor_data.Iy(:, new_col_pos) = reorder_index;
        sensor_data.Iy = sortrows(sensor_data.Iy, new_col_pos);
        sensor_data.Iy = sensor_data.Iy(:, 1:new_col_pos - 1);
    end

    % repeat for z-dimension if 3D
    if kgrid.dim > 2
        sensor_data.Iz(:, new_col_pos) = reorder_index;
        sensor_data.Iz = sortrows(sensor_data.Iz, new_col_pos);
        sensor_data.Iz = sensor_data.Iz(:, 1:new_col_pos - 1);
    end
end

if record.I_avg
    % append the reordering data
    sensor_data.Ix_avg(:, 2) = reorder_index;

    % reorder based on the order_index
    sensor_data.Ix_avg = sortrows(sensor_data.Ix_avg, 2);

    % remove the reordering data
    sensor_data.Ix_avg = sensor_data.Ix_avg(:, 1);

    % repeat for y-dimension if 2D or 3D
    if kgrid.dim > 1
        sensor_data.Iy_avg(:, 2) = reorder_index;
        sensor_data.Iy_avg = sortrows(sensor_data.Iy_avg, 2);
        sensor_data.Iy_avg = sensor_data.Iy_avg(:, 1);
    end

    % repeat for z-dimension if 3D
    if kgrid.dim > 2
        sensor_data.Iz_avg(:, 2) = reorder_index;
        sensor_data.Iz_avg = sortrows(sensor_data.Iz_avg, 2);
        sensor_data.Iz_avg = sensor_data.Iz_avg(:, 1);
    end
end

repo: Aero-matlab


text:
function [hs,probs] = forestApply( data, forest, maxDepth, minCount )
% Apply learned forest classifier.
%
% USAGE
%  [hs,probs] = forestApply( data, forest, [maxDepth], [minCount] )
%
% INPUTS
%  data     - [NxF] N length F feature vectors
%  forest   - learned forest classification model
%  maxDepth - [] maximum depth of tree
%  minCount - [] minimum number of data points to allow split
%
% OUTPUTS
%  hs       - [Nx1] predicted output labels
%  probs    - [NxH] predicted output label probabilities
%
% EXAMPLE
%
% See also forestTrain
%
% Piotr's Image&Video Toolbox      Version 3.01
% Copyright 2012 Piotr Dollar.  [pdollar-at-caltech.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the Simplified BSD License [see external/bsd.txt]
if(nargin<3 || isempty(maxDepth)), maxDepth=0; end
if(nargin<4 || isempty(minCount)), minCount=0; end
assert(isa(data,'single')); M=length(forest);
for i=1:M, tree=forest(i);
  if(maxDepth>0), tree.child(tree.depth>=maxDepth) = 0; end
  if(minCount>0), tree.child(tree.count<=minCount) = 0; end
  ids = forestInds(data,tree.thrs,tree.fids,tree.child);
  p=tree.distr(ids,:); if(i==1), probs=p; else probs=probs+p; end
end
[~,hs] = max(probs,[],2);
end

author: zouchuhang

repo: LayoutNet


text:
classdef TiGE2 < ALGORITHM
% <many> <real/integer/label/binary/permutation> <constrained/none>
% Tri-Goal Evolution Framework for CMaOPs

%------------------------------- Reference --------------------------------
% Y. Zhou, Z. Min, J. Wang, Z. Zhang, and J.Zhang, Tri-goal evolution
% framework for constrained many-objective optimization, IEEE Transactions
% on Systems Man and Cybernetics Systems, 2020, 50(8): 3086-3099.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        function main(Algorithm,Problem)
            %% Parameter setting
            [Epsilon0,row] = Algorithm.ParameterSet(0.05,1.01);

            %% Generate random population
            Population = Problem.Initialization();
            [fpr,fcd]  = Estimation(Population.objs,1/Problem.N^(1/Problem.M));
            fcv        = Calculate_fcv(Population);
            Epsilon    = Epsilon0;
            PopObj_1   = [fpr,fcd];
            [fm,~]     = NDSort(PopObj_1,Problem.N);
            PopObj     = [fm' + Epsilon * fcv,fcv];
            [frank,~]  = NDSort(PopObj,Problem.N);
            fitness    = frank' + fcv./(fcv+1);

            %% Optimization
            while Algorithm.NotTerminated(Population)
                MatingPool = TournamentSelection(2,Problem.N,fitness);
                Offspring  = OperatorGA(Problem,Population(MatingPool));
                [fpr,fcd]  = Estimation(Offspring.objs,1/Problem.N^(1/Problem.M));
                fcv = Calculate_fcv(Offspring);
                OffObj_1   = [fpr,fcd];
                [fm,~] = NDSort(OffObj_1,Problem.N);
                OffObj = [fm' + Epsilon * fcv,fcv];
                [Population,fitness] = EnvironmentalSelection([Population,Offspring],PopObj,OffObj,Problem.N);
                [fpr,fcd] = Estimation(Population.objs,1/Problem.N^(1/Problem.M));
                fcv = Calculate_fcv(Population);
                PopObj_1 = [fpr,fcd];
                [fm,~]   = NDSort(PopObj_1,Problem.N);
                PopObj   = [fm' + Epsilon * fcv,fcv];
                Epsilon  = row * Epsilon;
            end
        end
    end
end


text:
function [tfr,t,f] = tfrsp(x,t,N,h,trace);
%TFRSP	Spectrogram time-frequency distribution.
%	[TFR,T,F]=TFRSP(X,T,N,H,TRACE) computes the Spectrogram
%	distribution of a discrete-time signal X.
%
%	X     : signal.
%	T     : time instant(s)          (default : 1:length(X)).
%	N     : number of frequency bins (default : length(X)).
%	H     : analysis window, H being normalized so as to
%	        be  of unit energy.      (default : Hamming(N/4)).
%	TRACE : if nonzero, the progression of the algorithm is shown
%	                                 (default : 0).
%	TFR   : time-frequency representation. When called without
%	        output arguments, TFRSP runs TFRQVIEW.
%	F     : vector of normalized frequencies.
%
%	Example :
%	 sig=fmlin(128,0.1,0.4);
%	 h=tftb_window(17,'Kaiser'); tfrsp(sig,1:128,64,h,1);
%
%        [tfr,t,freq]=tfrsp(sig,1:128,64,h,1); plot(fftshift(freq),fftshift(tfr(:,100)))
%
%	See also all the time-frequency representations listed in
%	 the file CONTENTS (TFR*)

%	F. Auger, May-August 1994, July 1995.
%	Copyright (c) 1996 by CNRS (France).
%
%  This program is free software; you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published by
%  the Free Software Foundation; either version 2 of the License, or
%  (at your option) any later version.
%
%  This program is distributed in the hope that it will be useful,
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
%
%  You should have received a copy of the GNU General Public License
%  along with this program; if not, write to the Free Software
%  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA

if (nargin == 0),
 error('At least 1 parameter required');
end;
[xrow,xcol] = size(x);
if (nargin < 1),
 error('At least 1 parameter is required');
elseif (nargin <= 2),
 N=xrow;
end;

hlength=floor(N/4);
hlength=hlength+1-rem(hlength,2);

if (nargin == 1),
 t=1:xrow; h = tftb_window(hlength); trace=0;
elseif (nargin == 2) | (nargin == 3),
 h = tftb_window(hlength); trace=0;
elseif (nargin == 4),
 trace = 0;
end;

if (N<0),
 error('N must be greater than zero');
end;
[trow,tcol] = size(t);
if (xcol==0) | (xcol>2),
 error('X must have one or two columns');
elseif (trow~=1),
 error('T must only have one row');
elseif (2^nextpow2(N)~=N),
 fprintf('For a faster computation, N should be a power of two\n');
end;

[hrow,hcol]=size(h); Lh=(hrow-1)/2;
if (hcol~=1)|(rem(hrow,2)==0),
 error('H must be a smoothing window with odd length');
end;

tfr= zeros (N,tcol) ;
if trace, disp('Spectrogram'); end;
for icol=1:tcol,
 ti= t(icol); tau=-min([round(N/2)-1,Lh,ti-1]):min([round(N/2)-1,Lh,xrow-ti]);
 indices= rem(N+tau,N)+1;
 if trace, disprog(icol,tcol,10); end;
 tfr(indices,icol)=x(ti+tau).*conj(h(Lh+1+tau))/norm(h(Lh+1+tau));
end;
tfr=abs(fft(tfr)).^2;
if trace, fprintf('\n'); end;

if (nargout==0),
 tfrqview(tfr,x,t,'tfrsp',h);
elseif (nargout==3),
 if rem(N,2)==0,
  f=[0:N/2-1 -N/2:-1]'/N;
 else
  f=[0:(N-1)/2 -(N-1)/2:-1]'/N;
 end;
end;

author: HeLiangHIT


text:
function [s, cfg] = ft_statfun_indepsamplesT(cfg, dat, design)

% FT_STATFUN_INDEPSAMPLEST calculates the independent samples T-statistic on the
% biological data in dat (the dependent variable), using the information on the
% independent variable (ivar) in design.
%
% Use this function by calling one of the high-level statistics functions as
%   [stat] = ft_timelockstatistics(cfg, timelock1, timelock2, ...)
%   [stat] = ft_freqstatistics(cfg, freq1, freq2, ...)
%   [stat] = ft_sourcestatistics(cfg, source1, source2, ...)
% with the following configuration option
%   cfg.statistic = 'ft_statfun_indepsamplesT'
%
% Configuration options
%   cfg.computestat    = 'yes' or 'no', calculate the statistic (default='yes')
%   cfg.computecritval = 'yes' or 'no', calculate the critical values of the test statistics (default='no')
%   cfg.computeprob    = 'yes' or 'no', calculate the p-values (default='no')
%
% The following options are relevant if cfg.computecritval='yes' and/or
% cfg.computeprob='yes'.
%   cfg.alpha = critical alpha-level of the statistical test (default=0.05)
%   cfg.tail  = -1, 0, or 1, left, two-sided, or right (default=1)
%               cfg.tail in combination with cfg.computecritval='yes'
%               determines whether the critical value is computed at
%               quantile cfg.alpha (with cfg.tail=-1), at quantiles
%               cfg.alpha/2 and (1-cfg.alpha/2) (with cfg.tail=0), or at
%               quantile (1-cfg.alpha) (with cfg.tail=1).
%
% Design specification
%   cfg.ivar  = row number of the design that contains the labels of the conditions that must be
%               compared (default=1). The labels are the numbers 1 and 2.
%
% See also FT_TIMELOCKSTATISTICS, FT_FREQSTATISTICS or FT_SOURCESTATISTICS

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

% set the defaults
cfg.computestat    = ft_getopt(cfg, 'computestat', 'yes');
cfg.computecritval = ft_getopt(cfg, 'computecritval', 'no');
cfg.computeprob    = ft_getopt(cfg, 'computeprob', 'no');
cfg.alpha          = ft_getopt(cfg, 'alpha', 0.05);
cfg.tail           = ft_getopt(cfg, 'tail', 1);

% perform some checks on the configuration
if strcmp(cfg.computeprob,'yes') && strcmp(cfg.computestat,'no')
  % probabilities can only be calculated if the test statistics are calculated
  cfg.computestat = 'yes';
end
if isfield(cfg,'uvar') && ~isempty(cfg.uvar)
  ft_error('cfg.uvar should not exist for an independent samples statistic');
end

% perform some checks on the design and data
sel1     = design(cfg.ivar,:)==1;
sel2     = design(cfg.ivar,:)==2;
nreplc1  = sum(~isnan(dat(:,sel1)), 2);
nreplc2  = sum(~isnan(dat(:,sel2)), 2);
nrepl    = nreplc1 + nreplc2;
hasnans1 = any(nreplc1<sum(sel1));
hasnans2 = any(nreplc2<sum(sel2));

if any(nrepl<size(design,2))
  ft_warning('Not all replications are used for the computation of the statistic.');
end
df = nrepl - 2;

if strcmp(cfg.computestat, 'yes')
  % compute the statistic use nanmean only if necessary
  if hasnans1
    avg1 = nanmean(dat(:,sel1), 2);
    var1 = nanvar(dat(:,sel1), 0, 2);
  else
    avg1 = mean(dat(:,sel1), 2);
    % var1 = var(dat(:,sel1), 0, 2);
    % the following achieves the same as the line above, but faster
    var1 = (sum(dat(:,sel1).^2,2)-(avg1.^2).*nreplc1)./(nreplc1-1);
  end
  if hasnans2
    avg2 = nanmean(dat(:,sel2), 2);
    var2 = nanvar(dat(:,sel2), 0, 2);
  else
    avg2 = mean(dat(:,sel2), 2);
    % var2 = var(dat(:,sel2), 0, 2);
    % the following achieves the same as the line above, but faster
    var2 = (sum(dat(:,sel2).^2,2)-(avg2.^2).*nreplc2)./(nreplc2-1);
  end

  varc = (1./nreplc1 + 1./nreplc2).*((nreplc1-1).*var1 + (nreplc2-1).*var2)./df;

  % in the case of non-equal triallengths, and tfrs as input-data nreplc are
  % vectors with different values. when the triallengths are equal, and the
  % input is a tfr, nreplc are vectors with either zeros (all trials contain nan
  % meaning that t_ftimwin did not fit around data), or the number of trials
  s.stat = (avg1 - avg2)./sqrt(varc);
end

if strcmp(cfg.computecritval,'yes')
  % also compute the critical values
  s.df      = df;
  if all(df==df(1)), df = df(1); end
  if cfg.tail==-1
    s.critval = tinv(cfg.alpha,df);
  elseif  cfg.tail==0
    s.critval = [tinv(cfg.alpha/2,df),tinv(1-cfg.alpha/2,df)];
  elseif cfg.tail==1
    s.critval = tinv(1-cfg.alpha,df);
  end
end

if strcmp(cfg.computeprob,'yes')
  % also compute the p-values
  s.df      = df;
  if cfg.tail==-1
    s.prob = tcdf(s.stat,s.df);
  elseif  cfg.tail==0
    s.prob = 2*tcdf(-abs(s.stat),s.df);
  elseif cfg.tail==1
    s.prob = 1-tcdf(s.stat,s.df);
  end
end

author: spm


text:
function p = abs(p);
%ABS          Polynomial with absolute values of coefficients
%
%   r = abs(p);
%
% r_i = abs(p_i), i.e. result is interval polynomial for interval polynomial input
%

% written  08/28/00     S.M. Rump
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
%

  p.c = abs(p.c);


text:
% Author: Ricardo Baptista and Matthias Poloczek
% Date:   June 2018
%
% See LICENSE.md for copyright information
%

function output = local_search(objective, inputs)
% LOCAL_SEARCH: Function runs binary optimization by searching over the neighborhood
% of single model flips at each iteration

% Extract inputs
nVars = inputs.n_vars;
nEval = inputs.evalBudget;

% Generate initial condition and evaluate objective
model 	  = sample_models(1,nVars);
model_val = objective(model);

% determine the total number of iterations
nIter = ceil(nEval/nVars);

% Setup cells to store model, objective, and runtime
model_iter = zeros(nIter, nVars);
obj_iter   = zeros(nIter,1);
time_iter  = zeros(nIter,1);

% Declare counter
counter = 0;

for i=1:nIter

	% Update counter
	counter = counter + 1;

	ls_iter = tic;

	% Setup vector to store new objective values and difference
	new_obj  = zeros(nVars,1);
	diff_obj = zeros(nVars,1);

	for j=1:nVars

		% Setup new_model with one flipped variable
		new_model = model;
		new_model(j) = 1-new_model(j);

		% Evaluate objective
		new_obj(j)  = objective(new_model);
		diff_obj(j) = model_val - new_obj(j);

	end

	% Check if diff_obj is positive - improvement can be made
	if any(diff_obj > 0)

		% Choose optimal index to flip
    	[~, opt_idx] = max(diff_obj);
    	model(opt_idx) = 1 - model(opt_idx);
    	model_val = new_obj(opt_idx);

    end

	% Save models, model_obj, and runtime
	model_iter(counter,:) = model;
	obj_iter(counter)  = model_val;
	time_iter(counter) = toc(ls_iter);

end

% extend results
model_iter_new = zeros(nVars*nIter, nVars);
time_iter_new  = zeros(nVars*nIter, 1);
for i=1:nIter
	idx = nVars*(i-1)+1:nVars*i;
	model_iter_new(idx,:) = repmat(model_iter(i,:), nVars, 1);
	time_iter_new(idx) = interp1([0,1],[0,time_iter(i)],linspace(1/nVars,1,nVars));
end
obj_iter   = reshape(repmat(obj_iter, 1, nVars)', nVars*nIter, 1);

% save outputs
output = struct;
output.objVals  = obj_iter;
output.optModel = model_iter_new;
output.runTime  = time_iter_new;

end

author: baptistar

repo: BOCS

save_path: github-repos/MATLAB/baptistar-BOCS


text:
clear all;
% -------------------------------------------------------------------------
%   Description:
%       Script to demo LapSRN for one image
%
%   Citation:
%       Deep Laplacian Pyramid Networks for Fast and Accurate Super-Resolution
%       Wei-Sheng Lai, Jia-Bin Huang, Narendra Ahuja, and Ming-Hsuan Yang
%       IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2017
%
%   Contact:
%       Wei-Sheng Lai
%       wlai24@ucmerced.edu
%       University of California, Merced
% -------------------------------------------------------------------------

img_filename = 'emma.jpg';

%% parameters
scale  = 4; % SR upsampling scale
gpu    = 1; % GPU ID

%% setup paths
addpath(genpath('utils'));
addpath(fullfile(pwd, 'matconvnet/matlab'));
vl_setupnn;

%% Load pretrained odel
model_filename = fullfile('pretrained_models', sprintf('LapSRN_x%d.mat', scale));

fprintf('Load %s\n', model_filename);
net = load(model_filename);
net = dagnn.DagNN.loadobj(net.net);
net.mode = 'test' ;

if( gpu ~= 0 )
    gpuDevice(gpu)
    net.move('gpu');
end

%% Load GT image
fprintf('Load %s\n', img_filename);
img_GT = im2double(imread(img_filename));
img_GT = mod_crop(img_GT, scale);

%% Generate LR image
img_LR = imresize(img_GT, 1/scale);

%% apply LapSRN
fprintf('Apply LapSRN for %dx SR\n', scale);
img_HR = SR_LapSRN(img_LR, net, scale, gpu);

%% show results
img_LR = imresize(img_LR, scale);
figure, imshow(cat(2, img_LR, img_HR, img_GT));
title(sprintf('Bicubic %dx    |    LapSRN %dx    |    Ground Truth', scale, scale));

repo: LapSRN


text:
T_10 = readtable('/Users/lbasek/image-denoising-benchmark/results/final_results_sigma10.csv');
T_10(:,'WNNM') = [];
T_10(:,'NCSR') = [];

T_15 = readtable('/Users/lbasek/image-denoising-benchmark/results/final_results_sigma15.csv');
T_15(:,'WNNM') = [];
T_15(:,'NCSR') = [];

T_20 = readtable('/Users/lbasek/image-denoising-benchmark/results/final_results_sigma20.csv');
T_20(:,'WNNM') = [];
T_20(:,'NCSR') = [];

T_25 = readtable('/Users/lbasek/image-denoising-benchmark/results/final_results_sigma25.csv');
T_25(:,'WNNM') = [];
T_25(:,'NCSR') = [];

T_50 = readtable('/Users/lbasek/image-denoising-benchmark/results/final_results_sigma50.csv');
T_50(:,'WNNM') = [];
T_50(:,'NCSR') = [];

time_row_10 = T_10(4,:);
time_row_15 = T_15(4,:);
time_row_20 = T_20(4,:);
time_row_25 = T_25(4,:);
time_row_50 = T_50(4,:);

x = [0,10,15,20,25,50];

bm3d_10 = time_row_10{1,3};
bm3d_15 = time_row_15{1,3};
bm3d_20 = time_row_20{1,3};
bm3d_25 = time_row_25{1,3};
bm3d_50 = time_row_50{1,3};

ksvd_10 = time_row_10{1,4};
ksvd_15 = time_row_15{1,4};
ksvd_20 = time_row_20{1,4};
ksvd_25 = time_row_25{1,4};
ksvd_50 = time_row_50{1,4};

foe_10 = time_row_10{1,5};
foe_15 = ssim_row_15{1,5};
foe_20 = time_row_20{1,5};
foe_25 = time_row_25{1,5};
foe_50 = time_row_50{1,5};

epll_10 = time_row_10{1,6};
epll_15 = time_row_20{1,6};
epll_20 = ssim_row_20{1,6};
epll_25 = time_row_25{1,6};
epll_50 = time_row_50{1,6};

y1 = [time_row_10{1,2}, bm3d_10,bm3d_15,bm3d_20,bm3d_25,bm3d_50];
y2 = [time_row_10{1,2}, ksvd_10,ksvd_15,ksvd_20,ksvd_25,ksvd_50];
y3 = [time_row_10{1,2}, foe_10,foe_15,foe_20,foe_25,foe_50];
y4 = [time_row_10{1,2}, epll_10,epll_15,epll_20,epll_25,epll_50];

h = figure;
hold on
h1 = plot(x,y1,'-s','MarkerSize',10);
h2 = plot(x,y2,'-s','MarkerSize',10);
h3 = plot(x,y3,'-s','MarkerSize',10);
h4 = plot(x,y4,'-s','MarkerSize',10);
hold off
set([h1,h2,h3,h4],'LineWidth',2)
ylabel('Time (sec)')
xlabel('Noise Power')
legend('BM3D','KSVD','FOE','EPLL');
grid on
saveas(h,'./graphs/time','png');

author: lbasek

repo: image-denoising-benchmark


text:
function out=imag(x)

precision=x(1).precision;
out_rval=cell(size(x));
out_ival=cell(size(x));

for ii=1:numel(x)
 [temp,out_rval{ii}]=getVals(x,ii);
 out_ival{ii}='0';
end % for ii=1:max(ex,
out=class(struct('rval',out_rval,...
                 'ival',out_ival,...
                 'precision',precision),'mp');

author: opencobra

repo: cobratoolbox


text:
%%
% WARNING: challenge data from FAUST challenge contain outlier vertices.
% This demo filters them out and works with cleaned-up meshes; the
% resulting matches must then be mapped back to the original (corrupted)
% raw domains before submitting to the challenge!

clear all
close all
clc

addpath('./tools/')
addpath('./tools/flann/')
addpath(genpath('../../matlab/manopt/manopt/'))

options = struct;
options.k = 100;
options.icp_iters = 0;     % 0 for nearest neighbors
options.use_svd   = true;  % false for basic least squares
options.refine_iters = 0;  % 0 for no refinement

%% Load raw FAUST scans

% full resolution shapes
M = load_ply('./data/tr_scan_000.ply');
N = load_ply('./data/tr_scan_001.ply');

% NOTE: raw faust data contains outliers
[M, M.is_outlier] = cleanup(M);
[N, N.is_outlier] = cleanup(N);

%% Load sparse matches (these may come from some matching pipeline)

load('./data/sparse_matches.mat')
[i,j,~] = find(P);
sparse_matches = [i j];
n_matches = size(sparse_matches,1);

colors = create_colormap(N,N);
figure, colormap([1 1 1])
subplot(121)
plot_scalar_map(N, ones(N.n,1)); hold on
plot_cloud_color(N.VERT(sparse_matches(:,2),:), colors(sparse_matches(:,2),:), 5)
axis off; view([0 90])
subplot(122)
plot_scalar_map(M, ones(M.n,1)); hold on
plot_cloud_color(M.VERT(sparse_matches(:,1),:), colors(sparse_matches(:,2),:), 5)
axis off; view([0 90])

%% Compute LBO eigenfunctions

[M.W, ~, M.S] = calc_LB_FEM(M);
[M.evecs, M.evals] = eigs(M.W, M.S, options.k, -1e-5);
M.evals = diag(M.evals);
[M.evals, idx] = sort(M.evals);
M.evecs = M.evecs(:,idx);

[N.W, ~, N.S] = calc_LB_FEM(N);
[N.evecs, N.evals] = eigs(N.W, N.S, options.k, -1e-5);
N.evals = diag(N.evals);
[N.evals, idx] = sort(N.evals);
N.evecs = N.evecs(:,idx);

%% Refine and upscale matches

F = sparse(sparse_matches(:,1), 1:n_matches, 1, M.n, n_matches);
G = sparse(sparse_matches(:,2), 1:n_matches, 1, N.n, n_matches);

if options.refine_iters > 0

    A_init = M.evecs'*(M.S*F);
    B_init = N.evecs'*(N.S*G);
    [u,~,v] = svd(A_init*B_init');
    C_init = u*v';
    C_init = C_init';

    % fps among the input sparse matches
    fps = fps_euclidean(M.VERT(sparse_matches(:,1),:), 1e3, 1);

    matches_upscaled = refine_matches(...
        M, N, F(:,fps), G(:,fps), C_init, options);

    % do a final svd step
    G_svd = sparse(matches_upscaled, 1:M.n, 1, N.n, M.n);
    B_svd = M.evecs'*M.S;
    A_svd = N.evecs'*(N.S*G_svd);
    [u,~,v] = svd(A_svd*B_svd');
    [~, matches_upscaled_svd] = run_icp_fixed(N, M, v*u', options.icp_iters);

else

    B = M.evecs'*(M.S*F);
    A = N.evecs'*(N.S*G);

    if ~options.use_svd
        C_upscaled = A'\B';
        C_upscaled = C_upscaled';
    else
        [u,~,v] = svd(A*B');
        C_upscaled = u*v';
        C_upscaled = C_upscaled';
    end

    [~, matches_upscaled] = run_icp_fixed(N, M, C_upscaled, options.icp_iters);

end

figure
colors_hires = create_colormap(N,N);
subplot(141), colormap(colors_hires), plot_scalar_map(N, 1:N.n); axis off, view([0 90]), freeze_colors
subplot(142), colormap(colors_hires(matches_upscaled,:)), plot_scalar_map(M, 1:M.n); axis off, view([0 90]), freeze_colors
subplot(143), colormap(colors_hires), plot_scalar_map(N, 1:N.n); axis off, view([180 -90]), freeze_colors
subplot(144), colormap(colors_hires(matches_upscaled,:)), plot_scalar_map(M, 1:M.n); axis off, view([180 -90]), freeze_colors

author: OshriHalimi


text:
function [ap, prec, rec] = ml_ap(confidence, gt, draw)
% function [ap, prec, rec] = ml_ap(confidence, gt, draw)
% Average precision, adapted from VOCevaluation
% gt: a vector of 1 and -1. 1 is for positive, -1 is for negative.
% confidence: confidence for belonging to the positive class
% By: Minh Hoai Nguyen (minhhoai@robots.ox.ac.uk)
% Last modified: 23-Nov-2012

    if length(confidence) ~= length(gt)
        error('mismatch');
    end;
    confidence = confidence(:);
    gt = gt(:);
    [~,si]=sort(confidence, 'descend');
    tp=gt(si)>0;
    fp=gt(si)<0;

    fp=cumsum(fp);
    tp=cumsum(tp);
    rec=tp/sum(gt>0);
    prec=tp./(fp+tp);
    ap=VOCap(rec,prec);

    if draw
        % plot precision/recall
        plot(rec,prec,'-');
        grid;
        xlabel 'recall'
        ylabel 'precision'
    end

function ap = VOCap(rec,prec)
    mrec=[0 ; rec ; 1];
    mpre=[0 ; prec ; 0];
    for i=numel(mpre)-1:-1:1
        mpre(i)=max(mpre(i),mpre(i+1));
    end
    i=find(mrec(2:end)~=mrec(1:end-1))+1;
    ap=sum((mrec(i)-mrec(i-1)).*mpre(i));

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [a, tri] = trifacet_signed_area(tri, x)
% TRIFACET_SIGNED_AREA  Signed area of the facets in a 2D triangulation.
%
% A = trifacet_signed_area(TRI, X)
%
%   TRI is a 3-column matrix. Each row contains the 3 nodes that form one
%   triangular facet in the mesh.
%
%   X is a 2-column matrix. X(i, :) contains the xy-coordinates of the
%   i-th node in the mesh.
%
%   A is a vector where A(i) is the signed area of the i-th triangle.
%   Positive areas correspond to counter-clockwise triangles, and negative
%   areas, to clockwise triangles.
%
% [..., TRI2] = trifacet_signed_area(...)
%
%   TRI2 returns TRI reoriented when necessary so that all triangles have
%   a positive signed area.

% Author: Ramon Casero <rcasero@gmail.com>
% Copyright © 2013 University of Oxford
% Version: 0.2.0
%
% University of Oxford means the Chancellor, Masters and Scholars of
% the University of Oxford, having an administrative office at
% Wellington Square, Oxford OX1 2JD, UK.
%
% This file is part of Gerardus.
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details. The offer of this
% program under the terms of the License is subject to the License
% being interpreted in accordance with English Law and subject to any
% action against the University of Oxford being under the jurisdiction
% of the English Courts.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see
% <http://www.gnu.org/licenses/>.

% check arguments
narginchk(2, 2);
nargoutchk(0, 2);

if (size(tri, 2) ~= 3)
    error('TRI must have 3 columns')
end
if (size(x, 2) ~= 2)
    error('X must have 2 columns')
end

% signed area of a triangle with vertices  V1, V2, V3 in counter-clockwise
% orientation
%
% a = 1/2*((x2-x1)(y3-y1) - (x3-x1)(y2-y1))

x2x1 = x(tri(:, 2), 1) - x(tri(:, 1), 1);
y3y1 = x(tri(:, 3), 2) - x(tri(:, 1), 2);
x3x1 = x(tri(:, 3), 1) - x(tri(:, 1), 1);
y2y1 = x(tri(:, 2), 2) - x(tri(:, 1), 2);

a = 0.5 * (x2x1 .* y3y1 - x3x1 .* y2y1);

% reorient negative area triangles
if (nargout > 1)

    % triangles with negative area
    idx = a < 0;

    % reorient triangles with negative area so that they become positive
    tri(idx, 1:2) = tri(idx, [2 1]);

end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function view = resetSlider(view,slider,minVal,maxVal,val)
%
%     view = resetSlider(view,slider,minVal,maxVal,val)
%
% Resets slider to have specified min, max, and value.
% Sets the sliderHandle itself as well as the text label.
% Called by setParameterMap to reset the mapWindow sliders.
%
% djh, 4/99
% ras, 04/04:   made parallel changes with makeSlider and setSlider
%               (now has edit fields w/ label)

% Set slider limits
set(slider.sliderHandle,'min',minVal);
set(slider.sliderHandle,'max',maxVal);

% Set slider value
set(slider.sliderHandle,'Value',val);

% Update slider text labels
str = sprintf('%3.3g',val);
set(slider.labelHandle,'string',str);

% Return the current axes to the main image
% BW:  I don't understand this.  It didn't work when called from open3dwindow,
% so I changed gcf to view.ui.figNum in the hopes that the change was good
% enough.  It all seemed to work and it fixed my problem.
set(view.ui.figNum,'CurrentAxes',view.ui.mainAxisHandle);

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [name, data] = dispData(f)
%DISPDATA   Useful information for DISPLAY.
%   [NAME, DATA] = DISPDATA(F) collects extra information from the given CHEBFUN
%   F and the output will be used by DISPLAY.
%
% See also DISPLAY.

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

% Number of pieces (i.e., funs) information:
numFuns = numel(f.funs);

%% Collect information about the exponents:

% Preallocation:
exps = zeros(numFuns, 2);

% Loop over each FUN:
for j = 1:numFuns
    infoJ = dispData(f.funs{j});
    if ( ~isempty(infoJ) )
        numInfo = numel(infoJ);
        for k = 1:numInfo
            infoJK = infoJ{k};
            if ( strcmpi('exponents', infoJK.name) )
                exps(j,:) = infoJK.data;
            end
        end
    end
end

%% Organize the output which will be passed to DISPLAY:

% Preallocate the extra information:
name = ' ';
data = cell(j,1);

% So far, the only extra information is exponents:
if ( any( exps(:) ) )
    name = '  endpoint exponents';
    for j = 1:numFuns
        data{j} = ['        ' '[' num2str(exps(j,1), '%2.2g') '      ' ...
            num2str(exps(j,2), '%2.2g') ']' '  '];
    end
elseif ( ~isempty(infoJ) )
    name = infoJ{1}.name;
end

% More information for F can be appended to INFO:

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% confusion_matrix(1,1) = fraction of true 1 classified as 1
% confusion_matrix(1,2) = fraction of true 1 classified as 2
% confusion_matrix(1,3) = fraction of true 1 classified as unknown
% confusion_matrix(2,1) = fraction of true 2 classified as 1
% confusion_matrix(2,2) = fraction of true 2 classified as 2
% confusion_matrix(2,3) = fraction of true 2 classified as unknown
function confusion_matrix = ComputeConfusionMatrix(hs,diff_logprob,maxdiff_logprob1,mindiff_logprob2)

confusion_matrix = nan(2,3);
hsPr1 = diff_logprob < maxdiff_logprob1;
hsPr2 = diff_logprob > mindiff_logprob2;
n1 = nnz(hs==1);
n2 = nnz(hs==2);

confusion_matrix(1,1) = nnz(hsPr1(hs==1)) / n1;
confusion_matrix(1,2) = nnz(hsPr2(hs==1)) / n1;
confusion_matrix(1,3) = 1 - confusion_matrix(1,1) - confusion_matrix(1,2);

confusion_matrix(2,1) = nnz(hsPr1(hs==2)) / n2;
confusion_matrix(2,2) = nnz(hsPr2(hs==2)) / n2;
confusion_matrix(2,3) = 1 - confusion_matrix(2,1) - confusion_matrix(2,2);

author: kristinbranson

repo: JAABA


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

% Copyright 2002-2003 The MathWorks, Inc.

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

% Last Modified by GUIDE v2.5 07-Dec-2007 09:40:52

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

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

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

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

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes Robust_sys wait for user response (see UIRESUME)
% uiwait(handles.figure1);

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

% Get default command line output from handles structure
varargout{1} = handles.output;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

% --- Executes on button press in pushbutton1.
function pushbutton1_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

x1 = str2double(get(handles.edit1,'String'));
x2= str2double(get(handles.edit2,'String'));
x3= str2double(get(handles.edit3,'String'));
x4= str2double(get(handles.edit4,'String'));

x5 = str2double(get(handles.edit5,'String'));
x6 = str2double(get(handles.edit6,'String'));
x7 = str2double(get(handles.edit7,'String'));
x8 = str2double(get(handles.edit8,'String'));

x9 = str2double(get(handles.edit9,'String'));
x10 = str2double(get(handles.edit10,'String'));
x11 = str2double(get(handles.edit11,'String'));
x12 = str2double(get(handles.edit12,'String'));

A=[x1,x2;x3,x4];
B=[x5,x6;x7,x8];
C=[x9,x10;x11,x12];
D=0;
[h2n] = normh2(A,B,C,D);

set(handles.text3,'String',num2str(h2n))
% --- Executes on button press in pushbutton2.
function pushbutton2_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
x1 = str2double(get(handles.edit1,'String'));
x2= str2double(get(handles.edit2,'String'));
x3= str2double(get(handles.edit3,'String'));
x4= str2double(get(handles.edit4,'String'));

x5 = str2double(get(handles.edit5,'String'));
x6 = str2double(get(handles.edit6,'String'));
x7 = str2double(get(handles.edit7,'String'));
x8 = str2double(get(handles.edit8,'String'));

x9 = str2double(get(handles.edit9,'String'));
x10 = str2double(get(handles.edit10,'String'));
x11 = str2double(get(handles.edit11,'String'));
x12 = str2double(get(handles.edit12,'String'));

A=[x1,x2;x3,x4]
B=[x5,x6;x7,x8]
C=[x9,x10;x11,x12]
D=0;
[hinfn] = normhinf(A,B,C,D);
set(handles.text4,'String',num2str(hinfn))


function r1 = compute_r1_multiCam(good_image, junk_image, index, queryCam, testCam)
good_cam = testCam(good_image);
good_cam_uni = unique(good_cam);
r1 = ones(1, 6)-3;

% on the same camera
good_cam_now = queryCam;
ngood = length(junk_image)-1;
junk_image_now = [good_image; index(1)];
good_image_now = setdiff(junk_image, index(1));
good_now = 0;
for n = 1:length(index)
    flag = 0;
    if ngood == 0
        r1(good_cam_now) = -1;
        break;
    end
    if ~isempty(find(good_image_now == index(n), 1))
        flag = 1; % good image
        good_now = good_now+1;
    end
    if ~isempty(find(junk_image_now == index(n), 1))
        continue; % junk image
    end
    if flag == 0
        r1(good_cam_now) = 0;
        break;
    end
    if flag == 1%good
        r1(good_cam_now) = 1;
        break;
    end
end

for k = 1:length(good_cam_uni)
    good_cam_now = good_cam_uni(k);
    ngood = length(find(good_cam == good_cam_now));
    pos_junk = find(good_cam ~= good_cam_now);
    junk_image_now = [junk_image; good_image(pos_junk)];
    pos_good = find(good_cam == good_cam_now);
    good_image_now = good_image(pos_good);
    for n = 1:length(index)
        flag = 0;
        if ngood == 0
            r1(good_cam_now) = -1;
            break;
        end
        if ~isempty(find(good_image_now == index(n), 1))
            flag = 1; % good image
        end
        if ~isempty(find(junk_image_now == index(n), 1))
            continue; % junk image
        end

        if flag == 0
            r1(good_cam_now) = 0;
            break;
        end
        if flag == 1%good
            r1(good_cam_now) = 1;
            break;
        end
    end
end

end


text:
function varargout = process_swap_headcoils(varargin)
% PROCESS_SWAP_HEADCOILS: Correct channel names for swapped head localization coils.

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Marc Lalancette, 2019

eval(macro_method);
end

function sProcess = GetDescription() %#ok<DEFNU>
    % Description of the process
    sProcess.Comment     = 'Swap head coils (CTF)';
    sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/Tutorials/HeadMotion#Fixing_swapped_head_coils';
    sProcess.Category    = 'Custom';
    sProcess.SubGroup    = {'Import', 'Channel file'};
    sProcess.Index       = 52;
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'raw', 'data'};
    sProcess.OutputTypes = {'raw', 'data'};
    sProcess.nInputs     = 1;
    sProcess.nMinFiles   = 1;
    sProcess.isSeparator = 1;
    % Options
    sProcess.options.Na.Type    = 'checkbox';
    sProcess.options.Na.Comment = 'Nasion';
    sProcess.options.Na.Value   = 0;
    sProcess.options.Le.Type    = 'checkbox';
    sProcess.options.Le.Comment = 'Left ear';
    sProcess.options.Le.Value   = 0;
    sProcess.options.Re.Type    = 'checkbox';
    sProcess.options.Re.Comment = 'Right ear';
    sProcess.options.Re.Value   = 0;
    sProcess.options.reverse.Type    = 'checkbox';
    sProcess.options.reverse.Comment = 'Reverse order (when swapping all 3 coils)';
    sProcess.options.reverse.Value   = 0;
    %     sProcess.options.display.Type    = 'checkbox';
    %     sProcess.options.display.Comment = 'Display "before" and "after" alignment figures.';
    %     sProcess.options.display.Value   = 0;

end

function Comment = FormatComment(sProcess)
    Comment = sProcess.Comment;
    Na = sProcess.options.Na.Value;
    Le = sProcess.options.Le.Value;
    Re = sProcess.options.Re.Value;
    if Na && Le && Re
        if sProcess.options.reverse.Value
            Comment = [Comment, ': Na<-Le<-Re'];
        else
            Comment = [Comment, ': Na->Le->Re'];
        end
    elseif Na && Le
        Comment = [Comment, ': Na<->Le'];
    elseif Na && Re
        Comment = [Comment, ': Na<->Re'];
    elseif Le && Re
        Comment = [Comment, ': Le<->Re'];
    end

end

function OutputFiles = Run(sProcess, sInputs)
    OutputFiles = {sInputs.FileName};

    Na = sProcess.options.Na.Value;
    Le = sProcess.options.Le.Value;
    Re = sProcess.options.Re.Value;
    if Na + Le + Re < 2
        bst_report('Warning', sProcess, sInputs, 'Fewer than 2 coils selected, nothing to do.');
        return;
    end

    nInFiles = length(sInputs);
    [UniqueChan, iUniqFiles] = unique({sInputs.ChannelFile});
    nFiles = numel(iUniqFiles);
    if nFiles < nInFiles
        bst_report('Warning', sProcess, sInputs, ...
            'Multiple inputs were found for a single channel file. Only the first one will be used for swapping head coils.');
    end

    if Na && Le && Re
        if sProcess.options.reverse.Value
            Comment = 'Na<-Le<-Re';
            iBef = 1:12;
            iAft = [5:12, 1:4];
        else
            Comment = 'Na->Le->Re';
            iBef = [5:12, 1:4];
            iAft = 1:12;
        end
    elseif Na && Le
        Comment = 'Na<->Le';
        iBef = 1:8;
        iAft = [5:8, 1:4];
    elseif Na && Re
        Comment = 'Na<->Re';
        iBef = [1:4, 9:12];
        iAft = [9:12, 1:4];
    elseif Le && Re
        Comment = 'Le<->Re';
        iBef = 5:12;
        iAft = [9:12, 5:8];
    end

    bst_progress('start', 'Swap head coils', ' ', 0, nFiles);
    for iFile = iUniqFiles(:)' % no need to repeat on same channel file.

        ChannelMat = in_bst_channel(sInputs(iFile).ChannelFile);
        bst_progress('inc', 1);

        % Check the input is CTF.
        DataMat = in_bst_data(sInputs(iFile).FileName, 'Device');
        if ~strcmp(DataMat.Device, 'CTF')
            bst_report('Error', sProcess, sInputs(iFile), ...
                'Swap head coils is currently only available for CTF data.');
            continue;
        end

        HluChan = {'HLC0011', 'HLC0012', 'HLC0013', 'HLC0018', ...
            'HLC0021', 'HLC0022', 'HLC0023', 'HLC0028', ...
            'HLC0031', 'HLC0032', 'HLC0033', 'HLC0038'}; % usual order
        [Unused, iHlu] = ismember(HluChan, {ChannelMat.Channel.Name});
        if numel(iHlu) < numel(HluChan)
            bst_report('Error', sProcess, sInputs(iFile), 'Head coil channels not found.');
            continue;
        end

        for i = 1:numel(iBef)
            ChannelMat.Channel(iHlu(iBef(i))).Name = HluChan{iAft(i)};
        end

        ChannelMat = bst_history('add', ChannelMat, 'edit', ...
            ['Swapped head coils: ', Comment]);

        % Save channel file. Need to save before adjusting because LoadHlu
        % gets the channel info from file.
        bst_save(file_fullpath(sInputs(iFile).ChannelFile), ChannelMat, 'v7');

        %         % Attempt to adjust the initial/reference head position.
        %         [ChannelMat, Failed] = process_adjust_coordinates('AdjustHeadPosition', ...
        %             ChannelMat, sInputs(iFile), sProcess);
        %         if ~Failed
        %             % Save channel file.
        %             bst_save(file_fullpath(sInputs(iFile).ChannelFile), ChannelMat, 'v7');
        %             % Show new alignment.
        %             channel_align_manual(sInputs(iFile).ChannelFile, 'MEG', 0);
        %         % else Already noted in report.
        %         end

    end % file loop
    bst_progress('stop');

    OutputFiles = bst_process('CallProcess', 'process_adjust_coordinates', OutputFiles, [], ...
        'reset', 0, 'head', 1, 'bad', 1, 'points', 0, 'remove', 0, 'display', 1);

end

author: brainstorm-tools