text:
function x=windinfo(w,fs)
%V_WINDINFO window information and figures of merit X=(W,FS)
% Usage: (1) v_windinfo(v_windows('hamming',720,'ds'),720);         % plot hamming window info
%
%  Inputs:  W        is a vector containing the window
%           FS       is the sampling frequency (default=1)
%
% Outputs:  X.len         length of the window (samples)
%           X.nw          length of the window (samples)
%           X.ewgdelay    energy centroid delay from first sample (samples)
%           X.dcgain      DC gain (dB)
%           X.sidelobe    maximum sdelobe level in dB relative to DC gain
%           X.falloff     rate at which sidelobes decay (dB/octave)
%           X.enbw        equivalent noise bandwidth (*fs/len Hz)
%           X.scallop     scalloping loss (dB)
%           X.ploss       processing loss (dB)
%           X.wcploss     worst case processing loss (dB)
%           X.band3       3dB bandwidth (Hz)
%           X.band6       6 dB bandwidth (Hz)
%           X.band0       essential bandwidth (to first minimum) (Hz)
%           X.gain0       gain at first minimum (Hz)
%           X.olc50       50% overlap correction
%           X.olc75       75% overlap correction
%           X.cola        overlap factors giving constant overlap add (exlcuding multiples)
%           X.cola2       as X.cola but for squared window
%
% If no output argument is given, the window and frequency response
% will be plotted e.g. v_windinfo(v_windows('hamming',720,'ds'),720);
%
% To obtain the figures of merit listed in Table 1 of [1] set
% fs = length(W), multiply X.olc50 and X.olc75 by 100%. The "coherent gain
% listed in the table is 10^(x.dcgain/20)/(max(w)*length(w)).
%
%  [1]  F. J. Harris. On the use of windows for harmonic analysis with the
%       discrete fourier transform. Proc IEEE, 66 (1): 51-83, Jan. 1978.

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

if nargin<2
    fs=1;
end
w=w(:);
nw=length(w);
x.len=nw/fs;
x.nw=nw;
% energy weighted group delay = centre of energy
x.ewgdelay=((1:nw)*w.^2/sum(w.^2)-1)/fs;
% now calculate spectrum
of=16;      % spectrum oversample factor must be even
nwo=of*nw;
f=rfft(w,nwo);
p=f.*conj(f);
% sidelobe attenuation is maximum peak (note DC peak at p(1) is not found)
[kp,vp]=v_findpeaks(p,'q');
[kt,vt]=v_findpeaks(-p,'q');
% dbpo=10*log10(vp(2:end)./vp(1:end-1))./log2((kp(2:end)-1)./(kp(1:end-1)-1)); % slope in dB/octave
if ~numel(kp)
    x.sidelobe=10*log10(min(p)/p(1));
else
    x.sidelobe=10*log10(max(vp)/p(1));
end
np=length(kp);
ipa=floor(np/4);
if ~ipa
    x.falloff=0;
else
    ipb=floor(np/2);
    x.falloff=10*log10(vp(ipb)/vp(ipa))/log2((ipb-1)/(ipa-1));
end
sumw2=sum(w.^2);
sumw=sum(w);
enbwbin=nw*sumw2/sumw^2;
x.enbw=enbwbin*fs/nw;
x.dcgain=20*log10(sumw);
% do linear interpolation in p() to find 3dB and 6dB points
p3=0.5*p(1);
i3=find(p<p3,1);
if ~numel(i3)
    x.band3=Inf;
    x.band6=Inf;
else
    x.band3=2*(i3-(p3-p(i3))/(p(i3-1)-p(i3))-1)/of*fs/nw;
    p6=0.25*p(1);
    i6=find(p<p6,1);
    x.band6=2*(i6-(p6-p(i6))/(p(i6-1)-p(i6))-1)/of*fs/nw;
end
% do linear interpolation in f() to find closest approach to the origin
if~numel(kt)
    x.band0=Inf;
    x.gain0=0;
else
    i0=floor(kt(1));
    df=f(i0+1)-f(i0);
    j0=-real(f(i0)*conj(df))/abs(df)^2;
    x.band0=2*(i0+j0-1)/of*fs/nw;
    p0=abs(f(i0)+j0*df)^2;
    if p0>0
        x.gain0=10*log10(p0/p(1));
    else
        x.gain0=-Inf;
    end
end
% overlap factors
i50=round(nw*0.5);
x.olc50=sum(w(1:nw-i50).*w(1+i50:nw))/sumw2;
i75=round(nw*0.25);
x.olc75=sum(w(1:nw-i75).*w(1+i75:nw))/sumw2;
% processing loss and scalloping loss
x.scallop=10*log10(p(1)/p(1+of/2));
x.ploss=10*log10(enbwbin);
x.wcploss=x.ploss+x.scallop;
co=zeros(1,nw);
co2=co;
w2=w.^2;
for i=1:nw
    if i*fix(nw/i)==nw  % check i is a factor of the window length
        if co(i)==0
            co(i)=all(abs(sum(reshape(w',nw/i,i),2)-i*mean(w))<max(abs(w))*1e-4);
            if co(i)
                co(2*i:i:nw)=-1; % ignore multiples
            end
        end
        if co2(i)==0
            co2(i)=all(abs(sum(reshape(w2',nw/i,i),2)-i*mean(w2))<max(w2)*1e-4);
            if co2(i)
                co2(2*i:i:nw)=-1; % ignore multiples
            end
        end
    end
end
co(co<0)=0;
co2(co2<0)=0;
x.cola=find(co);
x.cola2=find(co2);
%
% now plot it if no output arguments given
%
if ~nargout
    clf;
    subplot(212);
    nf=min(max(floor(2*max(x.band6,x.band0)*of*nw/fs)+1,of*8),length(p));
    ff=(0:nf-1)*fs/(of*nw);
    fqi=[x.enbw x.band3 x.band6]/2;
    if ff(end)>2000
        ff=ff/1000;
        fqi=fqi/1000;
        xlab='kcyc/L';
    else
        xlab='cyc/L';
    end
    dbn=20*log10(x.nw); % window width in dB
    dbrange=min(100,-1.5*x.sidelobe);
    dd=10*log10(max(p(1:nf),p(1)*0.1^(dbrange/10)));
    ffs=[0 ff(end)];
    dbs=repmat(x.dcgain+x.sidelobe,1,2);
    ffb=[0 fqi(1) fqi(1)];
    dbb=[dd(1) dd(1) dd(1)-dbrange];
    ff3=[0 fqi(2) fqi(2)];
    db3=[dd(1)+db(0.5)/2 dd(1)+db(0.5)/2 dd(1)-dbrange];
    ff6=[0 fqi(3) fqi(3)];
    db6=[dd(1)+db(0.5) dd(1)+db(0.5) dd(1)-dbrange];
    area(ffb,dbb-dbn,max(dd)-dbrange-dbn,'facecolor',[1 0.7 0.7]);
    hold on
    plot(ffs,dbs-dbn,':k',ff3,db3-dbn,':k',ff6,db6-dbn,':k',ffb,dbb-dbn,'r',ff,dd-dbn,'b');
    legend(['Equiv Noise BW = ' sprintsi(x.enbw,-2) 'cyc/L'],['Max sidelobe = ' sprintf('%.0f',x.sidelobe) ' dB'],['-3 & -6dB BW = ' sprintf('%.2g',(x.band3)) ' & ' sprintf('%.2g',(x.band6)) ' cyc/L']);
    hold off
    axis([0 ff(end) max(dd)-dbrange-dbn max(dd)+2-dbn]);
    ylabel('Gain/N (dB)');
    xlabel(sprintf('Freq (%s)',xlab));
    %
    % Now plot the window itself
    %
    subplot(211);
    tax=(0:nw-1)/fs-x.ewgdelay;
    area(tax,w,'FaceColor',[0.7 0.7 1]);
    ylabel('Window');
    xlabel('Time/L');
    dtax=(tax(end)-tax(1))*0.02;
    axv=[tax(1)-dtax tax(end)+dtax min(0,min(w)) max(w)*1.05];
    texthvc(tax(end),max(w),sprintf('N=%d',nw),'rtk');
    if length(x.cola)>3
        tcola=sprintf(',%d',x.cola(1:3));
        tcola=[tcola ',...'];
    else
        tcola=sprintf(',%d',x.cola);
    end
    if length(x.cola2)>3
        tcola2=sprintf(',%d',x.cola2(1:3));
        tcola2=[tcola2 ',...'];
    else
        tcola2=sprintf(',%d',x.cola2);
    end
    texthvc(tax(1),max(w),sprintf('COLA=%s\nCOLA^2=%s',tcola(2:end),tcola2(2:end)),'ltk');
    axis(axv);
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function vis_logreg(pred,acc,idxPos,idxNeg)

p1_pos =  pred(idxPos,1);
p1_neg =  pred(idxNeg,1);
cost_pos = log((1-p1_pos)./p1_pos);
cost_neg = log((1-p1_neg)./p1_neg);

% vis_violate_pairwise(cost_positive, cost_negative);

% visualize training examples VS pairwise cost
thr_cost = 5;
x_range_cost = -thr_cost:0.1:thr_cost;
x_range_prob =  0:0.01:1;

% positive
subplot(2,3,1);
[hist_cost_pos,x] = hist(cost_pos, x_range_cost);
bar(x,hist_cost_pos);
xlim([x_range_cost(1), x_range_cost(end)]);
title(['Accuracy: ' num2str(acc(1))]);

hist_cost_pos = hist_cost_pos./sum(hist_cost_pos);
subplot(2,3,2);
plot(x_range_cost,hist_cost_pos,'Linewidth',2);
title('norm cost');

subplot(2,3,3);
plot(x_range_cost,cumsum(hist_cost_pos) ,'Linewidth',2);
title('cum norm cost');ylim([0, 1]); % plot cumulative normalized histogram

% negative
subplot(2,3,1);hold on;
[hist_cost_neg,x] =  hist(cost_neg, x_range_cost);
bar(x,hist_cost_neg, 'r');
xlim([x_range_cost(1), x_range_cost(end)]);
title('cost');
legend('pos ex costs', 'neg ex costs');
hist_cost_neg = hist_cost_neg./sum(hist_cost_neg);

subplot(2,3,2); hold on;
plot(x_range_cost,hist_cost_neg, 'Color','r');
title('norm cost');
legend('pos ex costs', 'neg ex costs');
subplot(2,3,3); hold on;
plot(x_range_cost,cumsum(hist_cost_neg), 'Color','r', 'Linewidth',2);
title('cum norm cost');ylim([0, 1]); % plot cumulative normalized histogram
legend('pos ex costs', 'neg ex costs');

% visualize training example VS pairwise probability
% positive
subplot(2,3,4);[g1_pos,x] =  hist(p1_pos,x_range_prob);
bar(x,g1_pos);
xlim([x_range_prob(1), x_range_prob(end)]);
title('hist of probability');  % plot histogram of distance

g1_pos = g1_pos./sum(g1_pos);
subplot(2,3,5); plot(x_range_prob, g1_pos); title('norm prob'); % plot normalized histogram

subplot(2,3,6);
plot(x_range_prob,cumsum(g1_pos), 'Linewidth',2);
title('cum norm prob'); ylim([0, 1]); % plot cumulative normalized histogram

% negative
subplot(2,3,4); hold on;
[g1_neg,x] = hist(p1_neg,x_range_prob);
bar(x,g1_neg,  'r'); xlim([x_range_prob(1), x_range_prob(end)]);
title('prob');% plot histogram of distance
legend('p(i) = 1 pos ', 'p(i) = 1 neg');

g1_neg = g1_neg./sum(g1_neg);
subplot(2,3,5);hold on;
plot(x_range_prob,g1_neg,'Color','r');
title('norm prob'); % plot normalized histogram
legend('p(i) = 1 pos ', 'p(i) = 1 neg');

subplot(2,3,6); hold on;
plot(x_range_prob, cumsum(g1_neg), 'Color','r','Linewidth',2);
title('cum norm prob');ylim([0, 1]);% plot cumulative normalized histogram
legend('p(i) = 1 pos ', 'p(i) = 1 neg');
end

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
% MINIMAL_LSTM
%   Demonstrates LSTMs on a toy binary addition problem.
%   The task is to predict the sequence of bits of the sum of two numbers,
%   given their binary sequences (one bit at a time).

run('../../setup_autonn.m') ;  % add AutoNN to the path
rng(0) ;  % set random seed

T = 8 ;  % number of bits / time steps
d = 16 ;  % dimensionality of the LSTM state
iters = 1500 ;  % number of iterations

% inputs
x = Input() ;
y = Input() ;

% initialize the shared parameters for an LSTM with d hidden units and
% two inputs
[W, b] = vl_nnlstm_params(d, 2) ;

% initial state
h = cell(T+1, 1);
c = cell(T+1, 1);
h{1} = zeros(d, 1, 'single');
c{1} = zeros(d, 1, 'single');

% run LSTM over all time steps
for t = 1:T
  [h{t+1}, c{t+1}] = vl_nnlstm(x(:,t), h{t}, c{t}, W, b);
end

% concatenate output into m x T matrix, ignoring initial state (h{1})
H = [h{2:end}] ;

% final projection (note the same projection is applied at all time steps)
prediction = vl_nnconv(reshape(H, 1, 1, d, T), 'size', [1, 1, d, 1]) ;

% define loss, and classification error
loss = vl_nnloss(prediction, y, 'loss', 'logistic') ;
err = vl_nnloss(prediction, y, 'loss','binaryerror') ;

% use workspace variables' names as the layers' names, and compile net
Layer.workspaceNames() ;
net = Net(loss, err) ;

% initialize solver
solver = solvers.Adam() ;
solver.learningRate = 1e-2 ;

losses = zeros(1, iters) ;
errors = zeros(1, iters) ;

for iter = 1:iters
  % generate a simple binary addition problem
  A = randi([0, 2^(T-1) - 1]) ;  % last bit is always 0 to prevent overflow
  B = randi([0, 2^(T-1) - 1]) ;
  C = A + B ;  % true answer

  % convert to vectors of binary digits
  data_x = [dec2bin(A, T); dec2bin(B, T)] ;  % concatenate binary numbers
  data_x = single(fliplr(data_x) == '1') ;  % convert from string to double, and reverse sequence
  data_y = dec2bin(C, T) ;  % same for the true answer
  data_y = single(fliplr(data_y) == '1') * 2 - 1;  % bit classes for prediction will be -1 or 1

  % evaluate network to compute gradients
  net.eval({'x', data_x, 'y', data_y}) ;

  % take one SGD step
  solver.step(net) ;

  % plot loss and error
  losses(iter) = net.getValue(loss) ;
  errors(iter) = net.getValue(err) ;

  if mod(iter, 100) == 0
    % display the current prediction
    fprintf('Iteration %i\n', iter);
    fprintf('True sequence: %s\n', sprintf('%i ', data_y > 0));
    fprintf('Predicted:     %s\n', sprintf('%i ', net.getValue('prediction') > 0));
    fprintf('Errors: %i\n\n', errors(iter) * T);
  end
end

figure(3) ;
subplot(1, 2, 1) ; plot(1 : 10 : iters, losses(1 : 10 : iters)) ;
xlabel('Iteration') ; ylabel('Loss') ;
subplot(1, 2, 2) ; plot(1 : 10 : iters, errors(1 : 10 : iters)) ;
xlabel('Iteration') ; ylabel('Error') ;

loss
net

repo: MFT


text:
function maxzlta()
    % maxzlta calculates the maximum z value for the LTA function.
    % The parameter step (window) can be defined by the user.

    ZG=ZmapGlobal.Data; % used by get_zmap_globals
    report_this_filefun();

    winlen_days = ZG.compare_window_dur_v3 / ZG.bin_dur;

    [len, ncu] = size(cumuall);       % redefine ncu
    len = len -2;
    lta = 1:1:ncu-2;
    var1 = zeros(1,ncu);
    var2 = zeros(1,ncu);
    lta = zeros(1,ncu);
    maxlta = zeros(1,ncu);
    maxlta = maxlta -5;
    cu = [cumuall(1:ti-1,:) ; cumuall(ti+winlen_days+1:len,:)];
    mean1 = mean(cu(:,:));
    wai = waitbar(0,'Please wait...')
    set(wai,'Color',[0.8 0.8 0.8],'NumberTitle','off','Name','Percent done');
    for i = 1:ncu
        var1(i) = cov(cu(:,i));
    end     % for i

    for it = 1:step: len - winlen_days

        waitbar(it/len)

        mean2 = mean(cumuall(it:it+winlen_days,:));
        for i = 1:ncu
            var2(i) = cov(cumuall(it:it+winlen_days,i));
        end     % for i
        lta = (mean1 - mean2)./(sqrt(var1/it+var2/(len-it)));
        maxlta2 = [maxlta ;  lta ];
        maxlta = max(maxlta2);

    end    % for it

    valueMap = reshape(maxlta,length(gy),length(gx));

    close(wai)

    stri = [  'Maximum z  Map of   '  file1];
    stri2 = ['winlen_days = ' char(days(winlen)) ];
    in = 'lta';
    view_max(valueMap,gx,gy,stri,'');

end

author: CelsoReyes


text:
function MW_src = perform_histogram_matching_wavelet(MW_src,MW_tgt, Jmin, options)

% perform_histogram_matching_wavelet - match the histogram of a wavelet transform
%
% Matching of wavelet coefficients only:
%   options.dotransform=0
%   MW_src = perform_histogram_matching_wavelet(MW_src,MW_tgt,Jmin,options);
% Matching of image + wavelet coefficients:
%   options.dotransform=1
%   M_src = perform_histogram_matching_wavelet(M_src,M_tgt,Jmin,options);
%
%   Match the spacial histogram of the image and the histogram of each
%   wavelet sub-band.
%
%   Works also for color images.
%
%   You can set options.use_histomatching=0 if you want to equalize the
%   kurtosis and skewness of subbands and not their histograms
%   (works well for natural images).
%
%   Copyright (c) 2004 Gabriel Peyr?

if nargin>=4
    if ~isstruct(options)
        error('options should be a structure.');
    end
end

options.null = 0;
dotransform = getoptions(options, 'dotransform', 1);
niter_synthesis = getoptions(options, 'niter_synthesis', 1);

if nargin<3
    Jmin = 3;
end

if niter_synthesis>1
    options.niter_synthesis = 1;
    for i=1:niter_synthesis
        MW_src = perform_histogram_matching_wavelet(MW_src,MW_tgt, Jmin, options);
    end
    return;
end

if size(MW_src,3)==3
    options.color_mode = 'pca';
    [MW_src,options.ColorP] = change_color_mode(MW_src,+1,options);
    [MW_tgt,options.ColorP] = change_color_mode(MW_tgt,+1,options);
    for i=1:size(MW_src,3)
        MW_src(:,:,i) = perform_histogram_matching_wavelet(MW_src(:,:,i),MW_tgt(:,:,i),Jmin, options);
    end
    [MW_src,options.ColorP] = change_color_mode(MW_src,-1,options);
    return;
end

if size(MW_src,3)>1
	for i=1:size(MW_src,3)
        MW_src(:,:,i) = perform_histogram_matching_wavelet(MW_src(:,:,i),MW_tgt(:,:,i),Jmin, options);
    end
    return;
end

% for spacial histogram, do not consider absval
options.absval = 0;
options.rows = 0;
options.cols = 0;

if dotransform == 1
    % perform image extension
    n1 = size(MW_src,1);
    n2 = size(MW_tgt,1);
    n = max(n1,n2);
    n = 2^( ceil(log2(n)) );
    MW_src = perform_image_extension(MW_src,n);
    MW_tgt = perform_image_extension(MW_tgt,n);
    options.wavelet_type = 'biorthogonal_swapped';
    options.wavelet_vm = 4;
    % perform spacial matching
    MW_src = perform_matching(MW_src, MW_tgt, options);
    % compute transform
    MW_src = perform_wavelet_transform(MW_src, Jmin, +1, options);
    MW_tgt = perform_wavelet_transform(MW_tgt, Jmin, +1, options);
    % perform coefficients matching
    options.dotransform = 0;
    MW_src = perform_histogram_matching_wavelet(MW_src,MW_tgt,Jmin, options);
    % undo transforms
    MW_src = perform_wavelet_transform(MW_src, Jmin, -1, options);
    MW_tgt = perform_wavelet_transform(MW_tgt, Jmin, -1, options);
    % perform spacial matching
    MW_src = perform_matching(MW_src, MW_tgt, options);
    MW_src = MW_src(1:n1,1:n1);
    return;
end

if size(MW_src,1)~=size(MW_tgt,1)
    error('Wavelets coefficients should be of the same size.');
end

mm = size(MW_src);
Jmax = log2(mm)-1;

for j=Jmax:-1:Jmin
    for q=1:3
        [selx,sely] = compute_quadrant_selection(j,q);
        MW_src(selx,sely) = match_statistics(MW_src(selx,sely), MW_tgt(selx,sely), options);
    end
end

% match low scales
selx = 1:2^Jmin; sely = 1:2^Jmin;
MW_src(selx,sely) = perform_matching(MW_src(selx,sely), MW_tgt(selx,sely), options);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function x = perform_matching(x,y, options)

% x = perform_histogram_matching(x,y, options);
x = perform_histogram_equalization(x,y, options);

function M = perform_image_extension(M,n)

m = size(M,1);
k = n-m;
while k>size(M,1)
    M = perform_image_extension(M,size(M,1)*2);
    k = k - size(M,1)/2;
end
M = [M; M(end:-1:end-k+1,:)];
M = [M, M(:,end:-1:end-k+1)];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function x = match_statistics(x,y, options)

options.null = 0;
if isfield(options, 'use_histomatching')
    use_histomatching = options.use_histomatching;
else
    use_histomatching = 1;
end

if isfield(options, 'nb_bins')
    nb_bins = options.nb_bins;
else
    nb_bins = 100;
end

if use_histomatching
    options.absval = 1;
    x = perform_matching(x, y, options);
else
    x = perform_kurtosis_equalization(x,y);
end

author: gpeyre

repo: matlab-toolboxes


text:
function transforms3d(varargin)
%TRANSFORMS3D  Conventions for manipulating 3D affine transforms.
%
%   By 'transform' we mean an affine transform. A 3D affine transform
%   is represented by a 4*4 matrix. The last row of the matrix is equal to
%   [0 0 0 1].
%
%
%
%   Example:
%   % create a translation by the vector [10 20 30]:
%   T = createTranslation3d([10 20 30]);
%   % Transform a basic point:
%   PT1 = [4 5 6];
%   PT2 = transformPoint3d(PT1, T)
%   % returns:
%   PT2 =
%       14   25   36
%
%   See also
%   createTranslation3d, createScaling3d, , createBasisTransform3d
%   createRotationOx, createRotationOy, createRotationOz
%   rotation3dAxisAndAngle, rotation3dToEulerAngles,
%   createRotation3dLineAngle, eulerAnglesToRotation3d
%   transformPoint3d, transformVector3d, transformLine3d, transformPlane3d
%   composeTransforms3d, recenterTransform3d
%

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2008-10-13,    using Matlab 7.4.0.287 (R2007a)
% Copyright 2008 INRA - BIA PV Nantes - MIAJ Jouy-en-Josas.

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
 function lists = subset_lists(nsubset, nview, varargin)
%function lists = subset_lists(nsubset, nview [, options])
%
% in
%	nsubset		# of subsets (1 to nview)
%	nview		# of views
% option
%	(none now, hard-wired to "bit-reversal" ordering)
% out
%	lists {cell}	{nsubset} lists of view indices for each subset
%
% Copyright 2000-4-?, Jeff Fessler, The University of Michigan

if nargin < 2, ir_usage, end

starts = subset_start(nsubset);

for ii=1:nsubset
	lists{ii} = [starts(ii):nsubset:nview];
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function m = size(t,idx)
%SIZE Size of ktensor.
%
%   D = SIZE(T) returns the size of the tensor.
%
%   I = SIZE(T,DIM) returns the size of the dimension specified by
%   the scalar DIM.
%
%   See also KTENSOR, KTENSOR/NDIMS.
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2012) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

if isempty(t.lambda)
    m = [];
end

if exist('idx','var')
    m = size(t.u{idx}, 1);
else
    for i = 1 : ndims(t)
	m(i) = size(t.u{i}, 1);
    end
end

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
% An example file for function 'compute_THD'.
% type 'help compute_THD' for more details.

clc;
close all;

% check cos wave:
'cos wave example'
t = 0:0.0001:1;
freq = 75;
x = 20 * cos(2*pi*freq*t + pi/4);
[ THD_cos, ph, amp ] = compute_THD( t,x, freq );
[THD_cos,  ph*180/pi,  amp]

% check square wave:
'square wave example'
freq = 10; % Hz
number_of_cycles = 4;
dt = 0.0001;
t = 0:dt:(number_of_cycles/freq);
x = (mod(t,1/freq) < 0.5/freq);
x = 1*(2*x-1);
[ THD_square, ph, amp ] = compute_THD( t,x, freq );

plot(t,x); xlabel('t [sec]'); title('square wave'); ylim([-1.1 1.1]);
[THD_square,  ph*180/pi,  amp]

%%% double check THD of square wave using the fourier series:
n = 1:10000;
freq_vec = freq*n;  % fourier series frequency vector
amp_vec = (4/pi)  * (1./n).*(mod(n,2)==1);
% compute THD by definition:
THD_square_theoretical = (sum(amp_vec(2:end).^2) / amp_vec(1)^2)^0.5;
THD_square_theoretical

author: Sable

repo: mcbench-benchmarks


text:
function obj = setDriftVector(obj, vf)
    % Set the drift vector field Fvec(x) or Fvec(x,dx) of the
    % system
    %
    % Parameters:
    %  vf: the f(x) @type cell

    % validate the inputs vf(x,dx)
    if ~iscell(vf), vf = {vf}; end
    for i=1:numel(vf)
        %         assert((length(vf{i})==obj.numState) && isvector(vf{i}),...
        %             'The %d-th drift vector field should be a (%d x 1) column vector.',i,obj.numState);
        if isa(vf{i},'SymFunction')
            vars = vf{i}.Vars;
            if strcmp(obj.Type,'SecondOrder') % second order system
                assert(length(vars)==2 && vars{1} == obj.States.x && vars{2}==obj.States.dx,...
                    'The SymFunction (vf{%d}) must be a function of states (x) and (dx).',i);
            else % first order system
                assert(length(vars)==1 && vars{1} == obj.States.x,...
                    'The SymFunction (vf{%d}) must be a function of states (x).',i);
            end
        else
            s_vf = SymExpression(vf{i});
            fvec_name = ['Fvec' num2str(i) '_' obj.Name];
            if strcmp(obj.Type,'SecondOrder') % second order system
                sfun_vf = SymFunction(fvec_name,s_vf,{obj.States.x,obj.States.dx});
            else                         % first order system
                sfun_vf = SymFunction(fvec_name,s_vf,{obj.States.x});
            end
            vf{i} = sfun_vf;
        end
    end

    obj.Fvec = vf;
    obj.FvecName_ = cellfun(@(f)f.Name, vf,'UniformOutput',false);
end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


text:
function y = tan(x)
%TAN          Implements  tan(x)  for intervals
%
%   y = tan(x)
%
%interval standard function implementation
%

% written  12/30/98     S.M. Rump
% modified 09/13/99     S.M. Rump  complex allowed, sparse input, NaN input,
%                                  major revision, improved accuracy
% modified 09/02/00     S.M. Rump  rounding unchanged after use
% modified 04/04/04     S.M. Rump  set round to nearest for safety
%                                    accelaration for sparse input
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 11/20/05     S.M. Rump  fast check for rounding to nearest
% modified 12/04/05     S.M. Rump  extreme values for approximate part
% modified 09/06/07     S.M. Rump  approximate std fcts removed
% modified 08/26/12     S.M. Rump  global variables removed
% modified 10/13/12     S.M. Rump  INTLAB_INTVAL_STDFCTS
%

  if issparse(x)
    [ix,jx,sx] = find(x);
    [m,n] = size(x);
    y = sparse(ix,jx,tan(full(sx)),m,n);
    return
  end

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

  if x.complex
    y = sin(x)./cos(x);
    setround(rndold)
    return
  end

  INTLAB_STDFCTS_PI = getappdata(0,'INTLAB_STDFCTS_PI');

  y = x;

  % transform x.inf and x.sup mod pi/2
  [ xinfinf , xinfsup , Sinf ] = modpi2(x.inf);
  [ xsupinf , xsupsup , Ssup ] = modpi2(x.sup);

  % indices with result +/- inf
  setround(1)
  delta = x.sup-x.inf;
  indexinf = ( delta >= INTLAB_STDFCTS_PI.PIINF ) | ...
    ( floor(Sinf/4) ~= floor(Ssup/4) );
  Sinf(Sinf>3) = Sinf(Sinf>3) - 4;
  Ssup(Ssup>3) = Ssup(Ssup>3) - 4;

  % transformation of input arguments by modpi2:
  %   [ xinf,xsup,s ] = modpi2(y)  ==>  0 <= s <= 7, 0 <= x <= pi/4 and
  %   x = -pi/2 + y + s*pi/4 + 2k*pi        for s even
  %   x =       - y + (s-1)*pi/4 + 2k*pi    for s odd

  y.inf(:) = -inf;
  y.sup(:) = inf;

  % treat non-infinity intervals
  Sinf(indexinf) = -1;
  Ssup(indexinf) = -1;

  % save warning status
  wng = warning;
  warning off

  % treat infimum
  index = ( Sinf==0 );
  if any(index(:))
    y.inf(index) = 1 ./ ( - tan_pos(xinfinf(index),-1) );
  end
  index = ( Sinf==1 );
  if any(index(:))
    y.inf(index) = - tan_pos(xinfsup(index),1);
  end
  index = ( Sinf==2 );
  if any(index(:))
    y.inf(index) = tan_pos(xinfinf(index),-1);
  end
  index = ( Sinf==3 );
  if any(index(:))
    res = tan_pos(xinfsup(index),1);
    setround(-1)
    y.inf(index) = 1 ./ res;
  end

  % treat supremum
  index = ( Ssup==0 );
  if any(index(:))
    y.sup(index) = 1 ./ ( - tan_pos(xsupsup(index),1) );
  end
  index = ( Ssup==1 );
  if any(index(:))
    y.sup(index) = - tan_pos(xsupinf(index),-1);
  end
  index = ( Ssup==2 );
  if any(index(:))
    y.sup(index) = tan_pos(xsupsup(index),1);
  end
  index = ( Ssup==3 );
  if any(index(:))
    res = tan_pos(xsupinf(index),-1);
    setround(1)
    y.sup(index) = 1 ./ res;
  end

  % restore warning status
  warning(wng);

  index = isnan(x.inf);
  if any(index(:))
    y.inf(index) = NaN;
    y.sup(index) = NaN;
  end

  setround(rndold)

repo: OpenMAS


text:
function [common, d1, d2] = comp_struct(s1,s2,prt,pse,tol,n1,n2)
% check two structures for differences - i.e. see if structure s1 == structure s2
% function [common, d1, d2] = comp_struct(s1,s2,prt,pse,tol)
%
% inputs  5 - 4 optional
% s1      structure one                              class structure
% s2      structure two                              class structure - optional
% prt     print test results (0 / 1 / 2 / 3)         class integer - optional
% pse     pause flag (0 / 1 / 2)                     class integer - optional
% tol     tol default tolerance (real numbers)       class integer - optional
%
% outputs 5 - 5 optional
% common  matching fields                            class struct - optional
% d1      non-matching fields for structure one      class struct - optional
% d2      non-matching fields for structure two      class struct - optional
%
% prt:
%	0 --> no print
%	1 --> print class erros
%	2 --> print all errors
% pse:
%	1 --> pause for class erros
%	2 --> pause for all errors
%
% example:	[same, er1, er2] = comp_struct(data1,data2,1,1,1e-6)
% michael arant - may 27, 2013
%
% updated - aug 22, 2015
%
% hint:
% passing just one structure causes the program to copy the structure
% and compare the two.  This is an easy way to list the structure
%
% note:  n1 and n2 are not required at launch - used for recursive purposes

%% default arguent check
if nargin < 2
	help comp_struct; error('I / O error');
end

if nargin < 3 || isempty(prt); prt = 0; end
if nargin < 4 || isempty(pse); pse = 0; elseif pse ~= 1 && prt == 0; pse = 0; end
if nargin < 5 || isempty(tol); tol = 1e-20; end
if pse > prt, pse = prt; end

if ~exist('n1','var'); n1 =  inputname(1); end
if ~exist('n2','var'); n2 = inputname(2); end

%% structure definition
d1 = s1; d2 = s2; common = s1;

%% begin analysis
flag = [0 0];
% test entire structure
if ~isequal(s1,s2)
	% differences noted - parse
	if isstruct(s1) && isstruct(s2)
		% both structures - once sub structures are tested, do not
		% modify the parrent
		flag(2) = 1;
		% both structures - get the field names for each structure
		fn1 = fieldnames(s1);
		fn2 = fieldnames(s2);
		% missing fields?  (common was a copy of s1 - so only s1 needs checking)
		temp = find(~ismember(fn1,fn2));
		for ii = 1:numel(temp)
			% drop unique fields
			common = rmfield(common,fn1{temp(ii)});
		end
		% get the common fields
		fn = fieldnames(common);
		% missing fields in set 1
		for ii = 1:numel(fn)
			% common field - recurse and test
			for jj = 1:min([numel(d1) numel(d2)])
				[common(jj).(fn{ii}) d1(jj).(fn{ii}) d2(jj).(fn{ii})] = ...
					comp_struct(s1(jj).(fn{ii}),s2(jj).(fn{ii}), ...
					prt,pse,tol, ...
					[n1 '(' num2str(jj) ').' fn{ii}], ...
					[n2 '(' num2str(jj) ').' fn{ii}]);
			end
		end
	% one or both not structures evaluate
	elseif isstruct(s1)
		% first variable is structure - second is not
		if prt
			fprintf('Error:	%s is a structure and %s is not a structure\n',n1,n2);
		end
		if pse
			uiwait(msgbox(sprintf('Error:  %s is a structure and %s is not a structure\n',n1,n2)))
		end
		% flag purge
		flag(1) = 1;

	elseif isstruct(s2)
		% second variable is structure - first is not
		if prt
			fprintf('Error:	%s is not a structure and %s is a structure\n',n1,n2);
		end
		if pse; pause; end
		% flag purge
		flag(1) = 1;
	else
		% the same?
		if ~isequal(s1,s2)
			flag(1) = 1;
			% not the same - difference?
			% class error
			if ~strcmp(class(s1),class(s2))
				% different classes
				if prt
					fprintf('Error:	%s is class %s and %s class %s\n',n1,class(s1),n2,class(s2));
				end
				if pse; pause; end
				flag(1) = 1;
			else
				% tolerance error?
				if min(size(s1) == size(s2)) && ...
						(isa(s1,'single') || isa(s1,'double'))
					% tolerance match?
					if ((max(s1)==0 || max(s2) ==0) && max(abs(s1-s2)) < tol) || numel(find(abs(s1-s2)./s1 < tol)) == numel(s1)
						flag(1) = 0;
					end
				else
					% print and pause?
					if prt > 1
						fprintf('%s is ',n1); disp(s1);
						fprintf('\b, %s is ',n2); disp(s2);
					end
					if pse > 1; pause; end
				end
			end
		end
	end
end

%% keep or delete
if flag(1) && ~flag(2)
	% denote error
	common = [];
elseif ~flag(2)
	% remove from error structures
	d1 = [];
	d2 = [];
end

%% purge empty fields
% test common
if isstruct(common)
	% fieldnames
	fn = fieldnames(common);
	for ii = 1:numel(fn)
		temp = 1;
		% test for empty field
		for jj = 1:numel(common)
			if ~isempty(common(jj).(fn{ii})); temp = 0; end
		end
		% purge if field was empty
		if temp
			common = rmfield(common,fn{ii});
		end
	end
end
% test d1
if isstruct(d1)
	% fieldnames
	fn = fieldnames(d1);
	for ii = 1:numel(fn)
		temp = 1;
		% test for empty field
		for jj = 1:numel(d1)
			if ~isempty(d1(jj).(fn{ii})); temp = 0; end
		end
		% purge if field was empty
		if temp
			d1 = rmfield(d1,fn{ii});
		end
	end
end
% test d2
if isstruct(d2)
	% fieldnames
	fn = fieldnames(d2);
	for ii = 1:numel(fn)
		temp = 1;
		% test for empty field
		for jj = 1:numel(d2)
			if ~isempty(d2(jj).(fn{ii})); temp = 0; end
		end
		% purge if field was empty
		if temp
			d2 = rmfield(d2,fn{ii});
		end
	end
end

%% test for null fields
if isstruct(common)
	% fieldnames
	if isempty(fieldnames(common)); common = []; end
end
if isstruct(d1)
	% fieldnames
	if isempty(fieldnames(d1)); d1 = []; end
end
if isstruct(d2)
	% fieldnames
	if isempty(fieldnames(d2)); d2 = []; end
end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function varargout = process_average( varargin )
% PROCESS_AVERAGE: Average files, by subject, by condition, or all at once.
%
% USAGE:                    OutputFiles = process_average('Run', sProcess, sInputs)
%                            OutputFile = process_average('AverageFiles', sProcess, sInputs, KeepEvents, isScaleDspm, isWeighted, isMatchRows, isZeroBad)
%                        [sMat,isFixed] = process_average('FixWarpedSurfaceFile', sMat, sInput, sStudyDest)
%  [iGroups, GroupComments, GroupNames] = process_average('SortFiles', sInputs, avgtype)

% @=============================================================================
% 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: Francois Tadel, 2010-2018

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % Description the process
    sProcess.Comment     = 'Average files';
    sProcess.Category    = 'Custom';
    sProcess.SubGroup    = 'Average';
    sProcess.Index       = 301;
    sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/Tutorials/Averaging#Averaging';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'data', 'results', 'timefreq', 'matrix'};
    sProcess.OutputTypes = {'data', 'results', 'timefreq', 'matrix'};
    sProcess.nInputs     = 1;
    sProcess.nMinFiles   = 2;
    sProcess.isSeparator = 1;
    % Definition of the options
    % === AVERAGE TYPE
    sProcess.options.label1.Comment = 'Group files:';
    sProcess.options.label1.Type    = 'label';
    sProcess.options.avgtype.Comment = {'Everything', 'By subject', 'By folder (subject average)', 'By folder (grand average)', 'By trial group (folder average)', 'By trial group (subject average)', 'By trial group (grand average)'};
    sProcess.options.avgtype.Type    = 'radio';
    sProcess.options.avgtype.Value   = 1;
    % === FUNCTION
    sProcess.options.label2.Comment = 'Function:';
    sProcess.options.label2.Type    = 'label';
    sProcess.options.avg_func.Comment = {'Arithmetic average:  mean(x)', ...
                                         'Average absolute values:  mean(abs(x))', ...
                                         'Root mean square (RMS):  sqrt(sum(x.^2)/N)', ...
                                         'Standard deviation:  sqrt(var(x))', ...
                                         'Standard error:  sqrt(var(x)/N)', ...
                                         'Arithmetic average + Standard deviation', ...
                                         'Arithmetic average + Standard error', ...
                                         'Median:  median(x)'};
    sProcess.options.avg_func.Type    = 'radio';
    sProcess.options.avg_func.Value   = 1;
    % === WEIGHTED AVERAGE
    sProcess.options.weighted.Comment    = 'Weighted average:  mean(x) = sum(Leff_i * x(i)) / sum(Leff_i)';
    sProcess.options.weighted.Type       = 'checkbox';
    sProcess.options.weighted.Value      = 0;
    sProcess.options.weightedlabel.Comment    = '      Leff_i = Effective number of averages for file #i';
    sProcess.options.weightedlabel.Type       = 'label';
    % === KEEP EVENTS
    sProcess.options.keepevents.Comment    = 'Keep all the event markers from the individual epochs';
    sProcess.options.keepevents.Type       = 'checkbox';
    sProcess.options.keepevents.Value      = 0;
    sProcess.options.keepevents.InputTypes = {'data', 'matrix'};
    % === SCALE NORMALIZE SOURCE MAPS (DEPRECATED OPTION AFTER INVERSE 2018)
    sProcess.options.scalenormalized.Comment    = 'Adjust normalized source maps for SNR increase.
Example: dSPM(Average) = sqrt(Navg) * Average(dSPM)';
    sProcess.options.scalenormalized.Type       = 'checkbox';
    sProcess.options.scalenormalized.Value      = 0;
    sProcess.options.scalenormalized.InputTypes = {'results'};
    sProcess.options.scalenormalized.Hidden     = 1;
    % === MATCH ROWS WITH NAMES
    sProcess.options.matchrows.Comment    = 'Match signals between files using their names';
    sProcess.options.matchrows.Type       = 'checkbox';
    sProcess.options.matchrows.Value      = 1;
    sProcess.options.matchrows.InputTypes = {'timefreq', 'matrix'};
    % === EXCLUDE ZEROS FROM THE AVERAGE
    sProcess.options.iszerobad.Comment    = 'Exclude the flat signals from the average (zero at all times)';
    sProcess.options.iszerobad.Type       = 'checkbox';
    sProcess.options.iszerobad.Value      = 1;
    sProcess.options.iszerobad.InputTypes = {'timefreq', 'matrix'};
end

%% ===== FORMAT COMMENT =====
function Comment = FormatComment(sProcess)
    % Function
    if isfield(sProcess.options, 'avg_func')
        switch(sProcess.options.avg_func.Value)
            case 1,  Comment = 'Average: ';
            case 2,  Comment = 'Average/abs: ';
            case 3,  Comment = 'RMS: ';
            case 4,  Comment = 'Standard deviation: ';
            case 5,  Comment = 'Standard error: ';
            case 6,  Comment = 'Average+Std: ';
            case 7,  Comment = 'Average+Stderr: ';
            case 8,  Comment = 'Median: ';
        end
    else
        Comment = 'Average: ';
    end
    % Weighted
    if isfield(sProcess.options, 'weighted') && isfield(sProcess.options.weighted, 'Value') && ~isempty(sProcess.options.weighted.Value) && sProcess.options.weighted.Value
        Comment = ['Weighted ' Comment];
    end
    % Average type
    iAvgType = sProcess.options.avgtype.Value;
    Comment = [Comment, sProcess.options.avgtype.Comment{iAvgType}];
end

%% ===== RUN =====
function OutputFiles = Run(sProcess, sInputs) %#ok<DEFNU>
    % Initialize returned list
    OutputFiles = {};
    % Get current progressbar position
    if bst_progress('isVisible')
        curProgress = bst_progress('get');
        bst_progress('text', 'Grouping files...');
    else
        curProgress = [];
    end
    % Weighted
    if isfield(sProcess.options, 'weighted') && isfield(sProcess.options.weighted, 'Value') && ~isempty(sProcess.options.weighted.Value)
        isWeighted = sProcess.options.weighted.Value;
    else
        isWeighted = 0;
    end
    % Keep events
    if isfield(sProcess.options, 'keepevents') && isfield(sProcess.options.keepevents, 'Value') && ~isempty(sProcess.options.keepevents.Value)
        KeepEvents = sProcess.options.keepevents.Value;
    else
        KeepEvents = 0;
    end
    % Scale normalized source maps (DEPRECATED AFTER INVERSE 2018)
    if isfield(sProcess.options, 'scalenormalized') && isfield(sProcess.options.scalenormalized, 'Value') && ~isempty(sProcess.options.scalenormalized.Value)
        isScaleDspm = sProcess.options.scalenormalized.Value;
    else
        isScaleDspm = 0;
    end
    % Match signals between files using their names
    if isfield(sProcess.options, 'matchrows') && isfield(sProcess.options.matchrows, 'Value') && ~isempty(sProcess.options.matchrows.Value)
        isMatchRows = sProcess.options.matchrows.Value;
    else
        isMatchRows = 1;
    end
    % Exclude zero values
    if isfield(sProcess.options, 'iszerobad') && isfield(sProcess.options.iszerobad, 'Value') && ~isempty(sProcess.options.iszerobad.Value)
        isZeroBad = sProcess.options.iszerobad.Value;
    else
        isZeroBad = 1;
    end

    % Check type of input files
    if ~isempty(strfind(sInputs(1).FileName, 'timefreq_pac')) || ~isempty(strfind(sInputs(1).FileName, 'timefreq_dpac'))
        bst_report('Error', sProcess, [], ['Cannot average PAC maps after their computation.' 10 'Please use the process option "Save average PAC across trials" instead.']);
        return;
    end

    % Group files
    [iGroups, GroupComments] = SortFiles(sInputs, sProcess.options.avgtype.Value);
    % Average each group
    for i = 1:length(iGroups)
        % Set progress bar at the same level for each loop
        if ~isempty(curProgress)
            bst_progress('set', curProgress);
        end
        % Do not process if there is only one input
        if (length(iGroups{i}) == 1)
            bst_report('Warning', sProcess, sInputs(iGroups{i}(1)), 'File is alone in its trial/comment group. Not processed.');
            continue;
        end
        % Set the comment of the output file
        if ~isempty(GroupComments)
            sProcess.options.Comment.Value = GroupComments{i};
        end
        % Compute group average
        OutputFiles{end+1} = AverageFiles(sProcess, sInputs(iGroups{i}), KeepEvents, isScaleDspm, isWeighted, isMatchRows, isZeroBad);
    end
end

%% ===== SORT FILES =====
function [iGroups, GroupComments, GroupNames] = SortFiles(sInputs, avgtype)
    GroupComments = [];
    GroupNames = {};
    % Group files in different ways: by subject, by condition, or all together
    switch (avgtype)
        % === EVERYTHING ===
        case 1
            iGroups = {1:length(sInputs)};
            GroupNames{1} = 'All';

        % === BY SUBJECT ===
        case 2
            [uniqueSubj,I,J] = unique({sInputs.SubjectFile});
            iGroups = cell(1, length(uniqueSubj));
            for i = 1:length(uniqueSubj)
                iGroups{i} = find(J == i)';
                GroupNames{i} = sInputs(iGroups{i}(1)).SubjectName;
            end

        % === BY CONDITION ===
        case {3,4}
            % Subject average
            if (avgtype == 3)
                inputCondPath = cell(1,length(sInputs));
                for iInput = 1:length(sInputs)
                    inputCondPath{iInput} = [sInputs(iInput).SubjectName, '/', sInputs(iInput).Condition];
                end
            % Grand average
            else
                inputCondPath = {sInputs.Condition};
            end
            % Process each condition independently
            [uniqueCond,I,J] = unique(inputCondPath);
            iGroups = cell(1, length(uniqueCond));
            GroupComments = cell(1, length(uniqueCond));
            for i = 1:length(uniqueCond)
                iGroups{i} = find(J == i)';
                GroupComments{i} = sInputs(iGroups{i}(1)).Condition;
            end
            GroupNames = uniqueCond;

        % === BY TRIAL GROUPS ===
        case {5,6,7}
            % Get the condition path (SubjectName/Condition/CommentBase) or (Condition/CommentBase) for each input file
            CondPath = cell(1, length(sInputs));
            trialComment = cell(1, length(sInputs));
            for iInput = 1:length(sInputs)
                % Default comment
                trialComment{iInput} = sInputs(iInput).Comment;
                % If results/timefreq and attached to a data file
                if any(strcmpi(sInputs(iInput).FileType, {'results','timefreq'})) && ~isempty(sInputs(iInput).DataFile)
                    switch (file_gettype(sInputs(iInput).DataFile))
                        case 'data'
                            [sStudyAssoc, iStudyAssoc, iFileAssoc] = bst_get('DataFile', sInputs(iInput).DataFile);
                            if ~isempty(sStudyAssoc)
                                trialComment{iInput} = sStudyAssoc.Data(iFileAssoc).Comment;
                            else
                                bst_report('Warning', 'process_average', sInputs(iInput), ['File skipped, the parent node has been deleted:' 10 sInputs(iInput).DataFile]);
                            end

                        case {'results', 'link'}
                            [sStudyAssoc, iStudyAssoc, iFileAssoc] = bst_get('ResultsFile', sInputs(iInput).DataFile);
                            if ~isempty(sStudyAssoc)
                                [sStudyAssoc2, iStudyAssoc2, iFileAssoc2] = bst_get('DataFile', sStudyAssoc.Result(iFileAssoc).DataFile);
                                if ~isempty(sStudyAssoc2)
                                    trialComment{iInput} = sStudyAssoc2.Data(iFileAssoc2).Comment;
                                else
                                    bst_report('Warning', 'process_average', sInputs(iInput), ['File skipped, the parent node has been deleted:' 10 sStudyAssoc.Result(iFileAssoc).DataFile]);
                                end
                            else
                                bst_report('Warning', 'process_average', sInputs(iInput), ['File skipped, the parent node has been deleted:' 10 sInputs(iInput).DataFile]);
                            end

                        case 'matrix'
                            [sStudyAssoc, iStudyAssoc, iFileAssoc] = bst_get('MatrixFile', sInputs(iInput).DataFile);
                            if ~isempty(sStudyAssoc)
                                trialComment{iInput} = sStudyAssoc.Matrix(iFileAssoc).Comment;
                            else
                                bst_report('Warning', 'process_average', sInputs(iInput), ['File skipped, the parent node has been deleted:' 10 sInputs(iInput).DataFile]);
                            end
                    end
                end
                % Condition average
                if (avgtype == 5)
                    CondPath{iInput} = [sInputs(iInput).SubjectName, '/', sInputs(iInput).Condition, '/', str_remove_parenth(trialComment{iInput})];
                % Subject average
                elseif (avgtype == 6)
                    CondPath{iInput} = [sInputs(iInput).SubjectName, '/', str_remove_parenth(trialComment{iInput})];
                % Grand average
                else
                    CondPath{iInput} = str_remove_parenth(trialComment{iInput});
                end
            end

            % Process each condition independently
            [uniquePath,I,J] = unique(CondPath);
            iGroups = cell(1, length(uniquePath));
            GroupComments = cell(1, length(uniquePath));
            for i = 1:length(uniquePath)
                % Skip empty paths
                if isempty(uniquePath{i})
                    continue;
                end
                % Find files in this condition
                iGroups{i} = find(J == i)';
                GroupComments{i} = str_remove_parenth(trialComment{iGroups{i}(1)});
            end
            GroupNames = uniquePath;
    end
end

%% ===== AVERAGE FILES =====
function OutputFile = AverageFiles(sProcess, sInputs, KeepEvents, isScaleDspm, isWeighted, isMatchRows, isZeroBad)
    OutputFile = [];
    % Parse inputs
    if (nargin < 7) || isempty(isZeroBad)
        isZeroBad = 1;
    end
    if (nargin < 6) || isempty(isMatchRows)
        isMatchRows = 1;
    end
    if (nargin < 5) || isempty(isWeighted)
        isWeighted = 0;
    end
    if (nargin < 4) || isempty(isScaleDspm)
        isScaleDspm = 0;
    end
    if (nargin < 3) || isempty(KeepEvents)
        KeepEvents = 0;
    end

    % === PROCESS AVERAGE ===
    % Get function
    isResults = strcmpi(sInputs(1).FileType, 'results');
    % if isResults && isfield(sProcess.options, 'avg_func')
    if isfield(sProcess.options, 'avg_func') && isfield(sProcess.options.avg_func, 'Value') && ~isempty(sProcess.options.avg_func.Value)
        switch (sProcess.options.avg_func.Value)
            case 1,  Function = 'mean';   isVariance = 0;   strComment = 'Avg';
            case 2,  Function = 'abs';    isVariance = 0;   strComment = 'Avg(abs)';
            case 3,  Function = 'rms';    isVariance = 0;   strComment = 'RMS';
            case 4,  Function = 'mean';   isVariance = 1;   strComment = 'Std';
            case 5,  Function = 'mean';   isVariance = 1;   strComment = 'StdError';
            case 6,  Function = 'mean';   isVariance = 1;   strComment = 'AvgStd';
            case 7,  Function = 'mean';   isVariance = 1;   strComment = 'AvgStderr';
            case 8,  Function = 'median'; isVariance = 0;   strComment = 'Median';
        end
    else
        Function = 'mean';   isVariance = 0;   strComment = 'Avg';
    end
    % Compute average
    [Stat, Messages, iAvgFile, Events] = bst_avg_files({sInputs.FileName}, [], Function, isVariance, isWeighted, isMatchRows, isZeroBad, 1);

    % Apply corrections on the variance value
    if strcmpi(strComment, 'Std') || strcmpi(strComment, 'AvgStd')
        Stat.var = sqrt(Stat.var);
    elseif strcmpi(strComment, 'StdError') ||  strcmpi(strComment, 'AvgStderr')
        Stat.var = sqrt(Stat.var / length(sInputs));
    end
    % Add messages to report
    if ~isempty(Messages)
        if isempty(Stat)
            bst_report('Error', sProcess, sInputs, Messages);
            return;
        else
            bst_report('Warning', sProcess, sInputs, Messages);
        end
    end

    % Load first file of the list
    [sMat, matName] = in_bst(sInputs(iAvgFile(1)).FileName);

    % === SCALE dSPM VALUES (DEPRECATED AFTER INVERSE 2018) ===
    % Apply a scaling to the dSPM functions, to compensate for the fact that the scaling applied to the NoiseCov was not correct
    if isScaleDspm && isResults && isfield(sMat, 'Function') && ismember(sMat.Function, {'dspm','mnp','glsp','lcmvp'}) && isfield(sMat, 'ImageGridAmp')
        if ~isWeighted
            bst_report('Warning', sProcess, [], 'You cannot scale the normalized maps if you do not compute a weighted average. Select the option "Weighted" to enable this option.');
        elseif (sMat.nAvg ~= Stat.nAvg)
            Factor = sqrt(Stat.nAvg) / sqrt(sMat.nAvg);
            bst_report('Warning', sProcess, [], sprintf('Averaging normalized maps (%s): scaling the values by %1.3f to match the number of trials averaged (%d => %d)', sMat.Function, Factor, sMat.nAvg, Stat.nAvg));
            % Apply on both .var and .mean fields
            if isfield(Stat, 'mean') && ~isempty(Stat.mean)
                Stat.mean = Factor * Stat.mean;
            end
            if isfield(Stat, 'var') && ~isempty(Stat.var)
                Stat.var = (Factor ^ 2) * Stat.var;
            end
        end
        % disp('BST> Warning: Averaging dSPM maps is different from computing the dSPM of the average. ');
        % disp('BST>          => dSPM(Average) = sqrt(N) * Average(dSPM)');
    end

    % === CREATE OUTPUT STRUCTURE ===
    % Get output study
    [sStudy, iStudy, Comment, uniqueDataFile] = bst_process('GetOutputStudy', sProcess, sInputs);
    % Comment: forced in the options
    if isfield(sProcess.options, 'Comment') && isfield(sProcess.options.Comment, 'Value') && ~isempty(sProcess.options.Comment.Value)
        Comment = [strComment ': ' sProcess.options.Comment.Value, ' (' num2str(length(sInputs)) ' files)'];
    % Comment: Process default
    else
        % % Replace the number of files
        % Comment = strrep(Comment, [num2str(length(sInputs)) ' '], [num2str(length(iAvgFile)) ' ']);
        % % Add file name
        % Comment = [strComment ': ' Comment, ' (' num2str(length(sInputs)) ' files)'];
        Comment = [strComment ': ' Comment];
    end
    % Weighted
    if isWeighted && ~strcmpi(strComment, 'Median')
        Comment = ['W' Comment];
    end
    % Copy fields from Stat structure
    switch (strComment)
        case {'Avg', 'Avg(abs)', 'RMS', 'Median'}
            sMat.(matName) = Stat.mean;
        case {'Std', 'StdError'}
            sMat.(matName) = Stat.var;
        case {'AvgStd', 'AvgStderr'}
            sMat.(matName) = Stat.mean;
            sMat.Std = Stat.var;
    end
    sMat.ChannelFlag = Stat.ChannelFlag;
    sMat.Time        = Stat.Time;
    sMat.nAvg        = Stat.nAvg;
    sMat.Leff        = Stat.Leff;
    sMat.Comment     = Comment;

    % History: Average
    if isfield(sMat, 'History')
        % Copy the history of the first file (but remove the entries "import_epoch" and "import_time")
        prevHistory = sMat.History;
        if ~isempty(prevHistory)
            % Remove entry 'import_epoch'
            iLineEpoch = find(strcmpi(prevHistory(:,2), 'import_epoch'));
            if ~isempty(iLineEpoch)
                prevHistory(iLineEpoch,:) = [];
            end
            % Remove entry 'import_time'
            iLineTime  = find(strcmpi(prevHistory(:,2), 'import_time'));
            if ~isempty(iLineTime)
                prevHistory(iLineTime,:) = [];
            end
        end
        % History for the new average file
        sMat = bst_history('reset', sMat);
        sMat = bst_history('add', sMat, 'average', FormatComment(sProcess));
        sMat = bst_history('add', sMat, 'average', 'History of the first file:');
        sMat = bst_history('add', sMat, prevHistory, ' - ');
    else
        sMat = bst_history('add', sMat, 'average', Comment);
    end
    % History: Number of files averaged for each channel
    if any(Stat.nGoodSamples(1) ~= Stat.nGoodSamples) && isfield(Stat, 'RowNames') && ~isempty(Stat.RowNames) && iscell(Stat.RowNames) && (length(Stat.RefRowNames) <= 1)
        strInfo = 'Number of files that were averaged, for each signal:';
        sMat = bst_history('add', sMat, 'average', strInfo);
        uniqueGood = unique(Stat.nGoodSamples);
        for i = length(uniqueGood):-1:1
            strHistoryGood = [sprintf(' - %d files: ', uniqueGood(i)), sprintf('%s ', Stat.RowNames{Stat.nGoodSamples == uniqueGood(i)})];
            sMat = bst_history('add', sMat, 'average', strHistoryGood);
            strInfo = [strInfo, 10, strHistoryGood];
        end
        % Generate process warning
        bst_report('Warning', sProcess, sInputs, strInfo);
    end
    % History: List files
    sMat = bst_history('add', sMat, 'average', 'List of averaged files:');
    for i = 1:length(iAvgFile)
        sMat = bst_history('add', sMat, 'average', [' - ' sInputs(iAvgFile(i)).FileName]);
    end

    % Averaging results from the different data file: reset the "DataFile" field
    if isfield(sMat, 'DataFile') && ~isempty(sMat.DataFile) && (length(uniqueDataFile) ~= 1)
        sMat.DataFile = [];
    end
    % Copy all the events found in the input files
    if KeepEvents && ~isempty(Events)
        sMat.Events = Events;
    else
        sMat.Events = [];
    end
    % Rownames
    if isfield(Stat, 'RowNames') && ~isempty(Stat.RowNames)
        if strcmpi(matName, 'TF')
            sMat.RowNames = Stat.RowNames;
        elseif strcmpi(matName, 'Value')
            sMat.Description = Stat.RowNames;
        end
    end

    % === AVERAGE WARPED BRAINS ===
    if isResults
        sMat = FixWarpedSurfaceFile(sMat, sInputs(1), sStudy);
    end

    % === SAVE FILE ===
    % Output filename
    if strcmpi(sInputs(1).FileType, 'data')
        allFiles = {};
        for i = 1:length(sInputs)
            [tmp, allFiles{end+1}, tmp] = bst_fileparts(sInputs(i).FileName);
        end
        fileTag = str_common_path(allFiles, '_');
    else
        fileTag = bst_process('GetFileTag', sInputs(1).FileName);
    end
    OutputFile = bst_process('GetNewFilename', bst_fileparts(sStudy.FileName), [fileTag, '_average']);
    % Save on disk
    bst_save(OutputFile, sMat, 'v6');
    % Register in database
    db_add_data(iStudy, OutputFile, sMat);
end

%% ===== FIX SURFACE FOR WARPED BRAINS =====
% Average source files coming from different subjects that are all different deformations of the default brain:
% should re-use the initial cortex surface instead of the first cortex surface available
function [sMat, isFixed] = FixWarpedSurfaceFile(sMat, sInput, sStudyDest)
    isFixed = 0;
    % Not a warped surface: skip
    if ~isfield(sMat, 'SurfaceFile') || isempty(sMat.SurfaceFile) || isempty(strfind(sMat.SurfaceFile, '_warped'))
        return;
    end
    % Must be from non-default to default anatomy
    isDestDefaultSubj = ismember(bst_fileparts(sStudyDest.BrainStormSubject), {bst_get('DirDefaultSubject'), bst_get('NormalizedSubjectName')});
    isSrcDefaultSubj  = ismember(bst_fileparts(sInput.SubjectFile),           {bst_get('DirDefaultSubject'), bst_get('NormalizedSubjectName')});
    if isSrcDefaultSubj || ~isDestDefaultSubj
        return;
    end
    % Rebuild possible target surface
    [tmp, fBase, fExt] = bst_fileparts(sMat.SurfaceFile);
    SurfaceFileDest = [bst_get('DirDefaultSubject'), '/', strrep([fBase, fExt], '_warped', '')];
    % Find destination file in database
    [sSubjectDest, iSubjectDes, iSurfDest] = bst_get('SurfaceFile', SurfaceFileDest);
    if isempty(iSurfDest)
        return;
    end
    % If this surface exists: use it if it has the same number of vertices as the source surface
    % Get the vertices number from source and destination surface files
    VarInfoSrc  = whos('-file', file_fullpath(sMat.SurfaceFile), 'Vertices');
    VarInfoDest = whos('-file', file_fullpath(SurfaceFileDest), 'Vertices');
    % Number of vertices match: change the surface
    if (VarInfoSrc.size(1) == VarInfoDest.size(1))
        sMat.SurfaceFile = SurfaceFileDest;
        isFixed = 1;
    end
end

author: brainstorm-tools


text:
function [lambda_vec, error_train, error_val] = ...
    validationCurve(X, y, Xval, yval)
%VALIDATIONCURVE Generate the train and validation errors needed to
%plot a validation curve that we can use to select lambda
%   [lambda_vec, error_train, error_val] = ...
%       VALIDATIONCURVE(X, y, Xval, yval) returns the train
%       and validation errors (in error_train, error_val)
%       for different values of lambda. You are given the training set (X,
%       y) and validation set (Xval, yval).
%

% Selected values of lambda (you should not change this)
lambda_vec = [0 0.001 0.003 0.01 0.03 0.1 0.3 1 3 10]';

% You need to return these variables correctly.
error_train = zeros(length(lambda_vec), 1);
error_val = zeros(length(lambda_vec), 1);

% ====================== YOUR CODE HERE ======================
% Instructions: Fill in this function to return training errors in
%               error_train and the validation errors in error_val. The
%               vector lambda_vec contains the different lambda parameters
%               to use for each calculation of the errors, i.e,
%               error_train(i), and error_val(i) should give
%               you the errors obtained after training with
%               lambda = lambda_vec(i)
%
% Note: You can loop over lambda_vec with the following:
%
%       for i = 1:length(lambda_vec)
%           lambda = lambda_vec(i);
%           % Compute train / val errors when training linear
%           % regression with regularization parameter lambda
%           % You should store the result in error_train(i)
%           % and error_val(i)
%           ....
%
%       end
%
%
for i = 1:length(lambda_vec)
	lambda = lambda_vec(i);
	theta = trainLinearReg(X, y, lambda);
    error_train(i) = linearRegCostFunction(X, y, theta, 0);
	error_val(i) = linearRegCostFunction(Xval,yval,theta,0);
end

% =========================================================================

end

author: AvaisP