text:
function [ x, w ] = spquad ( dim, ord, bpt )
%SPQUAD
%   Computes the sparse grid quadrature abscissae and weights
%   on an orthotope/hyperrectangle using the Clenshaw-Curtis rule.
%
%   [X,W]=SPQUAD(DIM,ORD,AB)
%
%   Input Parameters:
%   DIM - number of dimensions
%   ORD - order of the integration rule
%   BPT - boundary points (Optional. Defaults to [-1,1]^dim)
%
%   If used, BPT should be a 2 by DIM matrix containing the
%   endpoints of the hyperrectangle in each column.
%
%   Example usage:
%
%   f=@(x) (1+x(:,1)).*exp(x(:,2).*x(:,3))+(1+x(:,2)).*exp(x(:,3).*x(:,1))
%   [x,w]=spquad(3,4,[-1 0 2; 1 1 3]);
%   Q=w'*f(x);
%
%   written by Greg von Winckel - March 3, 2008
%   Contact: gregory(dot)von-winckel(at)uni-graz(dot)at

%
% Check that the number of dimensions agree.
%
if nargin>2
    if size(bpt,2)~=dim
        error('dimension mismatch');
    end
else
    bpt=repmat([-1,1]',1,dim);
end
%
% Length and midpoint in each dimension.
%
len=diff(bpt); mpt=mean(bpt);
%
% zeroth order case is just the midpoint rule
%
if ord==0
    x=mpt; w=prod(len);
%
%  1D case is special.
%
elseif ord~=0 && dim==1
    x=mpt+len*cos(pi*2^(-ord)*(0:(2^ord))')/2;
    w=clencurt(2^(ord)+1)*len/2;
else
    x0=mpt;
    w0=prod(len);
    for j=1:ord
        [x,w]=sparsegridnd(dim,j,mpt,len);
        [C,I]=intersect(x,x0,'rows');
        w(I)=w(I)+w0; x0=x; w0=w;
    end

end

%% Multidimensional nodes and difference weights.

function [ x, w ] = sparsegridnd ( n, ord, mpt, len )

% Generates the n-dimensional sparse grid of order ord on the
% hyperrectangle with centroid mpt and side lengths len, based
% on Chebyshev-Gauss-Lobatto points.

%
% Define 1D grid points
%
p=@(i) unique((i~=0)*cos(pi*2^(-i)*(0:(2^i)))');
%
% Get configurations of all possible subgrids
%
v=genindex(n,ord);
vmax=1+v(1,n);
%
% Compute all orders of one-dimensional quadrature rules
%
P=uaf(@(j) p(j), (0:vmax-1));
Q=uaf(@(j) diffweight(j), (0:vmax));
%
% Take the union of all possible subgrids
%
m=size(v,1);
xw=uaf(@(k) getpts(P,Q,v(k,:),mpt,len), (1:m)',1);
x=xw(:,1:n);
w=xw(:,n+1);
%
% Kludge to deal with small errors introduced by ndgrid
% need a better way to do this
%
roundn=@(a,n) round(a*10^n)/10^n;
x=roundn(x,15);
%
% Get unique points and weights
%
[x,ii,jj]=unique(x,'rows');
rows=length(ii);
cols=length(jj);
%
% Do node condensation for combining weights
%
w=sparse(jj,(1:cols)',ones(cols,1),rows,cols)*w;

%% One dimensional difference weights

function dw = diffweight ( ii )

if ii==0
    dw=2;
elseif ii==1
    dw=[1;-2;1]/3;
else
    q=@(i) clencurt(2^(i)+1);
    dw=q(ii); dw(1:2:end)=dw(1:2:end)-q(ii-1);
end

%% Compute the subset grid points and difference weights.

function [XW] = getpts ( P, Q, v, mpt, len )

n=size(v,2);
[x{1:n}]=ndgrid(P{1+v});
[w{1:n}]=ndgrid(Q{1+v});
d=size(v,2);

X=uaf(@(k) mpt(k)+x{k}(:)*len(k)/2, (1:d),1);
W=prod(uaf(@(k) w{k}(:), (1:d),1),2);
XW=[X,W];

%% Find all possible combinations of bases

function v = genindex ( n, L1, head )

if n==1
    v = L1;
else
    v = uaf(@(j) genindex(n-1, L1-j, j),(0:L1)',1);
end

if nargin>=3
    v = [head+zeros(size(v,1),1) v];
end

%% Compute 1D Clenshaw-Curtis weights
% Reference: Jörg Waldvogel,
%  "Fast construction of the Fejér and Clenshaw-Curtis quadrature rules,"
%  BIT Numerical Mathematics 43 (1), p. 001-018 (2004).

function w = clencurt ( N1 )

if N1==1
     w=2;
else
    N=N1-1;  c=zeros(N1,1);
    c(1:2:N1,1)=(2./[1 1-(2:2:N).^2 ])';
    f=real(ifft([c(1:N1,:);c(N:-1:2,:)]));
    w=2*([f(1,1); 2*f(2:N,1); f(N1,1)])/2;
end

%% Shorthand array function with uniform output
% Optional third argument converts output to matrix

function result = uaf ( arg, dex, varargin )

result=arrayfun(arg,dex,'UniformOutPut',false);

if nargin>2
    result=cell2mat(result);
end


text:
% TEST_CONV_SUB_1D Test case for CONV_SUB_1D
%
% See also
%   CONV_SUB_2D, WAVELET_1D

classdef test_conv_sub_1d < matlab.unittest.TestCase
    methods (Test)
        function testBasic(testcase)
            jsig = 14;
            sig_length = 2^jsig;
            load handel;
            x = y(1:sig_length);
            xf = fft(x);
            ds = randi(jsig)-1;

            white_list = {'filter_format'};
            type = {'d'};
            values = ...
                {{'fourier','fourier_multires','fourier_truncated'}};

            filt_opt = generate_random_options(white_list,type,values)
            filters = morlet_filter_bank_1d(sig_length,filt_opt);
            filter = filters.psi.filter{randi(jsig-ds)};

            for j=1:jsig
                y_old = old_conv_sub_1d(xf,filter,ds); % see below
                y = conv_sub_1d(xf,filter,ds);
                difference = norm(y_old-y)
                is_different = round(10^13*difference);
                assert(is_different==false);
            end
        end
    end
end

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function hammersley_seed_set ( seed )

%*****************************************************************************80
%
%% HAMMERSLEY_SEED_SET sets the seed of the leaped Hammersley subsequence.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 January 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer SEED(1:DIM_NUM), the seed of the leaped Hammersley subsequence.
%    Each base should be greater than 1.
%
  global hammersley_BASE
  global hammersley_LEAP
  global hammersley_DIM_NUM
  global hammersley_SEED
  global hammersley_STEP

  dim_num = hammersley_DIM_NUM;

  seed(1:dim_num) = floor ( seed(1:dim_num) );

  if ( ~halham_seed_check ( dim_num, seed ) )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'HAMMERSLEY_SEED_SET - Fatal error!\n' );
    error ( 'HAMMERSLEY_SEED_SET - Fatal error!' );
  end

  hammersley_SEED(1:dim_num) = seed(1:dim_num);

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function Callbacks_glottal_pulse_GUI25(f,C)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x=C{1,1};
y=C{1,2};
a=C{1,3};
b=C{1,4};
u=C{1,5};
v=C{1,6};
m=C{1,7};
n=C{1,8};
lengthbutton=C{1,9};
widthbutton=C{1,10};
enterType=C{1,11};
enterString=C{1,12};
enterLabel=C{1,13};
noPanels=C{1,14};
noGraphicPanels=C{1,15};
noButtons=C{1,16};
labelDist=C{1,17};%distance that the label is below the button
noTitles=C{1,18};
buttonTextSize=C{1,19};
labelTextSize=C{1,20};
textboxFont=C{1,21};
textboxString=C{1,22};
textboxWeight=C{1,23};
textboxAngle=C{1,24};
labelHeight=C{1,25};
fileName=C{1,26};
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%PANELS
for j=0:noPanels-1
uipanel('Parent',f,...
'Units','Normalized',...
'Position',[x(1+4*j) y(1+4*j) x(2+4*j)-x(1+4*j) y(3+4*j)-y(2+4*j)]);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%GRAPHIC PANELS
for i=0:noGraphicPanels-1
switch (i+1)
case 1
graphicPanel1 = axes('parent',f,...
'Units','Normalized',...
'Position',[a(1+4*i) b(1+4*i) a(2+4*i)-a(1+4*i) b(3+4*i)-b(2+4*i)],...
'GridLineStyle','--');
case 2
graphicPanel2 = axes('parent',f,...
'Units','Normalized',...
'Position',[a(1+4*i) b(1+4*i) a(2+4*i)-a(1+4*i) b(3+4*i)-b(2+4*i)],...
'GridLineStyle','--');
case 3
graphicPanel3 = axes('parent',f,...
'Units','Normalized',...
'Position',[a(1+4*i) b(1+4*i) a(2+4*i)-a(1+4*i) b(3+4*i)-b(2+4*i)],...
'GridLineStyle','--');
case 4
graphicPanel4 = axes('parent',f,...
'Units','Normalized',...
'Position',[a(1+4*i) b(1+4*i) a(2+4*i)-a(1+4*i) b(3+4*i)-b(2+4*i)],...
'GridLineStyle','--');
case 5
graphicPanel5 = axes('parent',f,...
'Units','Normalized',...
'Position',[a(1+4*i) b(1+4*i) a(2+4*i)-a(1+4*i) b(3+4*i)-b(2+4*i)],...
'GridLineStyle','--');
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%TITLE BOXES
for k=0:noTitles-1
switch (k+1)
case 1
titleBox1 = uicontrol('parent',f,...
'Units','Normalized',...
'Position',[u(1+4*k) v(1+4*k) u(2+4*k)-u(1+4*k) v(3+4*k)-v(2+4*k)],...
'Style','text',...
'FontSize',textboxFont{k+1},...
'String',textboxString(k+1),...
'FontWeight',textboxWeight{k+1},...
'FontAngle',textboxAngle{k+1});
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%BUTTONS
for i=0:(noButtons-1)
enterColor='w';
if strcmp(enterType{i+1},'pushbutton')==1 ||strcmp(enterType{i+1},'text')==1
enterColor='default';
end
if (strcmp(enterLabel{1,(i+1)},'')==0 &&...
        strcmp(enterLabel{1,(i+1)},'...')==0) %i.e. there is a label
%creating a label for some buttons
uicontrol('Parent',f,...
'Units','Normalized',...
'Position',[m(1+2*i) n(1+2*i)-labelDist-labelHeight(i+1) ...
(m(2+2*i)-m(1+2*i)) labelHeight(i+1)],...
'Style','text',...
'String',enterLabel{i+1},...
'FontSize', labelTextSize(i+1),...
'HorizontalAlignment','center');
end
switch (i+1)
case 1
button1=uicontrol('Parent',f,...
'Units','Normalized',...
'Position',[m(1+2*i) n(1+2*i) (m(2+2*i)-m(1+2*i)) (n(2+2*i)-n(1+2*i))],...
'Style',enterType{i+1},...
'String',enterString{i+1},...
'FontSize', buttonTextSize(1+i),...
'BackgroundColor',enterColor,...
'HorizontalAlignment','center',...
'Callback',@button1Callback);
case 2
button2=uicontrol('Parent',f,...
'Units','Normalized',...
'Position',[m(1+2*i) n(1+2*i) (m(2+2*i)-m(1+2*i)) (n(2+2*i)-n(1+2*i))],...
'Style',enterType{i+1},...
'String',enterString{i+1},...
'FontSize', buttonTextSize(1+i),...
'BackgroundColor',enterColor,...
'HorizontalAlignment','center',...
'Callback',@button2Callback);
case 3
button3=uicontrol('Parent',f,...
'Units','Normalized',...
'Position',[m(1+2*i) n(1+2*i) (m(2+2*i)-m(1+2*i)) (n(2+2*i)-n(1+2*i))],...
'Style',enterType{i+1},...
'String',enterString{i+1},...
'FontSize', buttonTextSize(1+i),...
'BackgroundColor',enterColor,...
'HorizontalAlignment','center',...
'Callback',@button3Callback);
case 4
button4=uicontrol('Parent',f,...
'Units','Normalized',...
'Position',[m(1+2*i) n(1+2*i) (m(2+2*i)-m(1+2*i)) (n(2+2*i)-n(1+2*i))],...
'Style',enterType{i+1},...
'String',enterString{i+1},...
'FontSize', buttonTextSize(1+i),...
'BackgroundColor',enterColor,...
'HorizontalAlignment','center',...
'Callback',@button4Callback);
case 5
button5=uicontrol('Parent',f,...
'Units','Normalized',...
'Position',[m(1+2*i) n(1+2*i) (m(2+2*i)-m(1+2*i)) (n(2+2*i)-n(1+2*i))],...
'Style',enterType{i+1},...
'String',enterString{i+1},...
'FontSize', buttonTextSize(1+i),...
'BackgroundColor',enterColor,...
'HorizontalAlignment','center',...
'Callback',@button5Callback);
case 6
button6=uicontrol('Parent',f,...
'Units','Normalized',...
'Position',[m(1+2*i) n(1+2*i) (m(2+2*i)-m(1+2*i)) (n(2+2*i)-n(1+2*i))],...
'Style',enterType{i+1},...
'String',enterString{i+1},...
'FontSize', buttonTextSize(1+i),...
'BackgroundColor',enterColor,...
'HorizontalAlignment','center',...
'Callback',@button6Callback);
case 7
button7=uicontrol('Parent',f,...
'Units','Normalized',...
'Position',[m(1+2*i) n(1+2*i) (m(2+2*i)-m(1+2*i)) (n(2+2*i)-n(1+2*i))],...
'Style',enterType{i+1},...
'String',enterString{i+1},...
'FontSize', buttonTextSize(1+i),...
'BackgroundColor',enterColor,...
'HorizontalAlignment','center',...
'Callback',@button7Callback);
case 8
button8=uicontrol('Parent',f,...
'Units','Normalized',...
'Position',[m(1+2*i) n(1+2*i) (m(2+2*i)-m(1+2*i)) (n(2+2*i)-n(1+2*i))],...
'Style',enterType{i+1},...
'String',enterString{i+1},...
'FontSize', buttonTextSize(1+i),...
'BackgroundColor',enterColor,...
'HorizontalAlignment','center',...
'Callback',@button8Callback);
case 9
button9=uicontrol('Parent',f,...
'Units','Normalized',...
'Position',[m(1+2*i) n(1+2*i) (m(2+2*i)-m(1+2*i)) (n(2+2*i)-n(1+2*i))],...
'Style',enterType{i+1},...
'String',enterString{i+1},...
'FontSize', buttonTextSize(1+i),...
'BackgroundColor',enterColor,...
'HorizontalAlignment','center',...
'Callback',@button9Callback);
end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%USER CODE FOR THE VARIABLES AND CALLBACKS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialize Variables
alpha1=25;
alpha2=10;
period=100;
nfft=1024;
fs=10000;
ivowel=1;
fsd=10000;
imp_train=[];
pulse_train=[];
zout=[];
g=[];
vowels='IY';

% Name the GUI
    set(f,'Name','glottal_pulse_GUI');

% CALLBACKS
% Callback for button1 -- alpha1, editable button value for glottal pulse
% opening cycle (in %)
 function button1Callback(h,eventdata)
     alpha1=str2num(get(button1,'string'))/100;
     if alpha1 < 0
         waitfor(msgbox('alpha1 must lie between 0 and 100'));
     elseif alpha1 > 1
         waitfor(msgbox('alpha1 must lie between 0 and 100'));
     end
 end

% Callback for button2 -- alpha2, editable button value for glottal pulse
% closing cycle (in %)
 function button2Callback(h,eventdata)
     alpha2=str2num(get(button2,'string'))/100;
     if alpha2 < 0
         waitfor(msgbox('alpha2 must lie between 0 and 100'));
     elseif (alpha2 > 1)
         waitfor(msgbox('alpha2 must lie between 0 and 100'));
     end
 end

% Callback for button3 -- period, editable button duration of pitch
% period (in samples) at fs=10000 Hz
 function button3Callback(h,eventdata)
     period=str2num(get(button3,'string'));
     if period < 10
         waitfor(msgbox('pitch period must lie between 10 and 250 samples'));
     elseif period > 250
         waitfor(msgbox('pitch period must lie between 10 and 250 samples'));
     end
 end

% Callback for button4 -- ivowel: editable button for vowel sound impulse
% response to be convolved with glottal pulse
 function button4Callback(h,eventdata)
     ivowel=get(button4,'val');
 end

% Callback for button9 -- Generate Excitation -- both Rosenberg pulse and
% periodic extension at designated pitch period
    function button9Callback(h,eventdata)
        button3Callback(h,eventdata);
        button4Callback(h,eventdata);

     alpha1=str2num(get(button1,'string'))/100;
     alpha2=str2num(get(button2,'string'))/100;

% form rosenberg glottal pulse, g
    n1=round(period*alpha1);
    n2=round(period*alpha2);
    g=[];
    x1=0:n1;
    g=[g 0.5*(1-cos(pi*x1/n1))];
    x2=1:n2;
    g=[g cos(pi*x2/(2*n2))];
    g=[g zeros(1,period-length(g))];

% plot g as a time domain signal over a single period in graphics Panel 1
    reset(graphicPanel2);
    axes(graphicPanel2);

    x=0:length(g)-1;
    plot(x,g,'b','LineWidth',2);
    xpp=['Time in Samples; fs=',num2str(fs),' samples/second'];
    xlabel(xpp);
    ylabel('Pulse Amplitude');axis tight; grid on;
    s1=sprintf('Rosenberg pulse -- period: %d,',period);
    s2=sprintf(' alpha1, alpha2: %6.3f, %6.3f ',alpha1,alpha2);
    stitle=strcat(s1,s2);
    legend('Rosenberg Pulse');

% display glottal pulse parameters in titleBox1
    set(titleBox1,'String',stitle);
    set(titleBox1,'FontSize',25);

% plot log magnitude spectrum of glottal pulse, g, in graphics Panel 1
    ge=[g zeros(1,nfft-length(g))];
    gs=20*log10(abs(fft(ge)));
    freq=0:fs/nfft:fs/2;
    x2=[0 fs/2];
    s2=[0 0];

    reset(graphicPanel1);
    axes(graphicPanel1);
    plot(freq,gs(1:nfft/2+1),'r','LineWidth',2);
    hold on;
    plot(x2,s2,'g','LineWidth',2);
    xlabel('Frequency in Hz'),ylabel('Log Magnitude (dB)');
    axis tight; grid on;legend('Rosenberg Pulse Spectrum');

% play out periodic pulse train of impulses spaced period samples apart,
% and of total length of fsd=10000 samples (fsd defined in list of
% variables at beginning of Callbacks)
    imp_train=zeros(1,fsd);
    imp_train(1:period:fsd)=1;

 % display impulse train for 1000 samples on Graphics Panel 5
    reset(graphicPanel5);
    axes(graphicPanel5);
    plot(0:1000,imp_train(1:1001),'r','LineWidth',2);
    xpp=['Time in Samples; fs=',num2str(fs),' samples/second'];
    xlabel(xpp),ylabel('Value'),axis tight, grid on; legend('Periodic Impulses');

% generate excitation pulses using Rosenberg pulse
    pulse_train=conv(g,imp_train);

% display glottal pulse train for 1000 samples on Graphics Panel 4
    reset(graphicPanel4);
    axes(graphicPanel4);
    plot(0:1000,pulse_train(1:1001),'b','LineWidth',2);
    xpp=['Time in Samples; fs=',num2str(fs),' samples/second'];
    xlabel(xpp),ylabel('Value'),axis tight, grid on; legend('Glottal Pulses');

% load vowel info (formants and bandwidths)
% load 'vowels_fmts_bw.mat'
    str=load('vowels_fmts_bw.mat');
    vowels=str.vowels;
    formants=str.formants;
    bandwidths=str.bandwidths;

% generate vowel impulse response
    fmts=formants(ivowel,:);
    fmts=[fmts 4500];
    yout=vowel_ir(fmts,bandwidths,fsd);

% convolve vowel impulse response with excitation
    zout=conv(yout,pulse_train);

% display periodic vowel sequence for 1000 samples on Graphics Panel 3
    reset(graphicPanel3);
    axes(graphicPanel3);
    plot(0:1000,zout(1:1001),'k','LineWidth',2);
    xpp=['Time in Samples; fs=',num2str(fs),' samples/second'];
    xlabel(xpp),ylabel('Value'),axis tight, grid on;

    xleg=['Vowel Response:',vowels(ivowel,:),'  Period:',num2str(period)];
    legend(xleg);
    end

% Callback for button5 -- Play Periodic Impulses
    function button5Callback(h,eventdata)
    soundsc(imp_train,fsd);
    end

% Callback for button6 -- Play Glottal Pulses
    function button6Callback(h,eventdata)
    soundsc(pulse_train,fsd);
    end

% Callback for button7 -- Play Vowel Excited by Glottal Pulses
    function button7Callback(h,eventdata)
    soundsc(zout,fsd);
    end

% Callback for button8 -- close GUI
 function button8Callback(h,eventdata)
     close(gcf);
 end
end

author: Sable

repo: mcbench-benchmarks


text:
classdef LHSintegrator_StiffnessColumn < LHSintegrator

    properties (Access = public)
        geometry
        stiffnessMatrix
    end

    properties (Access = private)
        freeNodes
    end

    methods (Access = public)

        function obj = LHSintegrator_StiffnessColumn(cParams)
            obj.init(cParams)
            obj.initStiffnessColumn(cParams);
            obj.createQuadrature();
            obj.createInterpolation();
            obj.createGeometry();
        end

        function LHS = compute(obj)
            lhs = obj.computeElementalLHS();
            LHS = obj.assembleMatrix(lhs);
            obj.stiffnessMatrix = LHS;
        end

        function [Kfree,free] = provideFreeStiffnessMatrix(obj)
            free = obj.freeNodes;
            K = obj.stiffnessMatrix;
            Kfree  = K(free,free);
        end

    end

    methods (Access = protected)

        function lhs = computeElementalLHS(obj)
            d = obj.dim;
            nElem = obj.mesh.nelem;
            Edof = d.ndofPerElement;
            Ke = zeros(Edof,Edof,nElem);
            l = obj.computeLength();
            [c1,c2,c3,c4,c5] = obj.coeffsStiffness(l);
            Ke(1,1,:) = c1.*c2;
            Ke(1,2,:) = c1.*c3;
            Ke(1,3,:) = -c1.*c2;
            Ke(1,4,:) = c1.*c3;
            Ke(2,1,:) = c1.*c3;
            Ke(2,2,:) = c1.*c4;
            Ke(2,3,:) = -c1.*c3;
            Ke(2,4,:) = -c1.*c5;
            Ke(3,1,:) = -c1.*c2;
            Ke(3,2,:) = -c1.*c3;
            Ke(3,3,:) = c1.*c2;
            Ke(3,4,:) = -c1.*c3;
            Ke(4,1,:) = c1.*c3;
            Ke(4,2,:) = -c1.*c5;
            Ke(4,3,:) = -c1.*c3;
            Ke(4,4,:) = c1.*c4;
            lhs = Ke;
        end

    end

    methods (Access = private)

        function initStiffnessColumn(obj,cParams)
            obj.freeNodes = cParams.freeNodes;
        end

        function createGeometry(obj)
            q   = obj.quadrature;
            int = obj.interpolation;
            int.computeShapeDeriv(q.posgp);
            s.mesh = obj.mesh;
            g = Geometry.create(s);
            g.computeGeometry(q,int);
            obj.geometry = g;
        end

        function l = computeLength(obj)
            g = obj.geometry;
            l = sum(g.dvolu,2);
        end

        function [c1,c2,c3,c4,c5] = coeffsStiffness(obj,l)
            c1 = (1./(30*l))';
            c2 = 36*ones(1,length(l));
            c3 = (3*l)';
            c4 = (4*l.^2)';
            c5 = (l.^2)';
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
% compute recall and viewpoint accuracy
% modified from PASCAL3D+ VDPM/compute_recall_precision_accuracy.m
% added an input param: prediction filename
% changed a few relative path names
function [recall, precision, accuracy, ap, aa] = compute_recall_precision_accuracy_3dview(cls, prediction_filename, show_curve, ARP)

BASE_DIR = fullfile(mfilename('fullpath'),'../');
addpath(fullfile(BASE_DIR, '../'));
global_variables;

% backward compatibility
if nargin < 3
    show_curve = 1;
end
if nargin < 4
    ARP = 0;
end

%azimuth_interval = [0 (360/(vnum_test*2)):(360/vnum_test):360-(360/(vnum_test*2))];

% viewpoint annotation path
%path_ann_view = '../Annotations';
path_ann_view = fullfile(g_pascal3d_root_folder, 'Annotations');

% read ids of validation images
addpath(fullfile(g_pascal3d_root_folder, 'VDPM'));
addpath(fullfile(g_pascal3d_root_folder, 'PASCAL/VOCdevkit/VOCcode'));
VOCinit;
%pascal_init;
ids = textread(sprintf(VOCopts.imgsetpath, 'val'), '%s');
M = numel(ids);

% open prediction file
%filename = sprintf('data/%s_%d_test_flip.mat', cls, vnum_train);
filename = prediction_filename;
object = load(filename);
dets_all = object.dets;

% load pre-packaged records (from extract_records.m)
% this step greatly accelerates testing.
try
    object = load(fullfile(BASE_DIR, 'voc12val_records.mat'));
catch
    run(fullfile(BASE_DIR, 'extract_records.m'));
    object = load(fullfile(BASE_DIR, 'voc12val_records.mat'));
end
voc12val_records = object.voc12val_records;

energy = [];
correct = [];
correct_view = [];
overlap = [];
count = zeros(M,1);
num = zeros(M,1);
num_pr = 0;
for i = 1:M
    fprintf('%s %d/%d\n', cls, i, M);
    % read ground truth bounding box
    rec = voc12val_records{i}; %PASreadrecord(sprintf(VOCopts.annopath, ids{i}));
    clsinds = strmatch(cls, {rec.objects(:).class}, 'exact');
    diff = [rec.objects(clsinds).difficult];
    clsinds(diff == 1) = [];
    n = numel(clsinds);
    bbox = zeros(n, 4);
    for j = 1:n
        bbox(j,:) = rec.objects(clsinds(j)).bbox;
    end
    count(i) = size(bbox, 1);
    det = zeros(count(i), 1);

    % read ground truth viewpoint
    if isempty(clsinds) == 0
        filename = fullfile(path_ann_view, sprintf('%s_pascal/%s.mat', cls, ids{i}));
        object = load(filename);
        record = object.record;
        view_gt = zeros(n, 3);
        for j = 1:n
            viewpoint = record.objects(clsinds(j)).viewpoint;
            if record.objects(clsinds(j)).viewpoint.distance == 0
                azimuth = viewpoint.azimuth_coarse;
                elevation = viewpoint.elevation_coarse;
                tilt = viewpoint.theta;
            else
                azimuth = viewpoint.azimuth;
                elevation = viewpoint.elevation;
                tilt = viewpoint.theta;
            end
            %view_gt(j) = find_interval(azimuth, azimuth_interval);
            view_gt(j,:) = [azimuth, elevation, tilt];
        end
    else
        view_gt = [];
    end

    % get predicted bounding box
    dets = dets_all{i};
    num(i) = size(dets, 1);
    % for each predicted bounding box
    for j = 1:num(i)
        num_pr = num_pr + 1;
        energy(num_pr) = dets(j, 8); %dets(j, 6);
        bbox_pr = dets(j, 1:4);
        view_pr = dets(j, 5:7); %find_interval((dets(j, 5) - 1) * (360 / vnum_train), azimuth_interval);

        % compute box overlap
        if isempty(bbox) == 0
            o = box_overlap(bbox, bbox_pr);
            [maxo, index] = max(o);
            if maxo >= 0.5 && det(index) == 0
                overlap{num_pr} = index;
                correct(num_pr) = 1;
                det(index) = 1;
                % check viewpoint
                R_pr = angle2dcm(view_pr(1)/180*pi, view_pr(2)/180*pi, view_pr(3)/180*pi);
                R_gt = angle2dcm(view_gt(index,1)/180*pi, view_gt(index,2)/180*pi, view_gt(index,3)/180*pi);
                R_angle = norm(logm(R_pr' * R_gt)) / sqrt(2); % in [0,2pi]
                azimuth_angle = min(abs(view_pr(1) - view_gt(index,1)), 360-abs(view_pr(1) - view_gt(index,1))); % in degree
                if ARP && R_angle < pi/6  % ARP_pi/6
                    correct_view(num_pr) = 1;
                elseif ARP == 0 && azimuth_angle < 30 % AVP_pi/6
                    correct_view(num_pr) = 1;
                else
                    correct_view(num_pr) = 0;
                end
            else
                overlap{num_pr} = [];
                correct(num_pr) = 0;
                correct_view(num_pr) = 0;
            end
        else
            overlap{num_pr} = [];
            correct(num_pr) = 0;
            correct_view(num_pr) = 0;
        end
    end
end
overlap = overlap';

[threshold, index] = sort(energy, 'descend');
correct = correct(index);
correct_view = correct_view(index);
n = numel(threshold);
recall = zeros(n,1);
precision = zeros(n,1);
accuracy = zeros(n,1);
num_correct = 0;
num_correct_view = 0;
for i = 1:n
    % compute precision
    num_positive = i;
    num_correct = num_correct + correct(i);
    if num_positive ~= 0
        precision(i) = num_correct / num_positive;
    else
        precision(i) = 0;
    end

    % compute accuracy
    num_correct_view = num_correct_view + correct_view(i);
    if num_correct ~= 0
        accuracy(i) = num_correct_view / num_positive;
    else
        accuracy(i) = 0;
    end

    % compute recall
    recall(i) = num_correct / sum(count);
end

ap = VOCap(recall, precision);
fprintf('AP = %.4f\n', ap);

aa = VOCap(recall, accuracy);
fprintf('AA = %.4f\n', aa);

% draw recall-precision and accuracy curve
if show_curve
figure;
hold on;
plot(recall, precision, 'r', 'LineWidth',3);
plot(recall, accuracy, 'g', 'LineWidth',3);
xlabel('Recall');
ylabel('Precision/Accuracy');
tit = sprintf('Average Precision = %.1f / Average Accuracy = %.1f', 100*ap, 100*aa);
title(tit);
hold off;
end

function ind = find_interval(azimuth, a)

for i = 1:numel(a)
    if azimuth < a(i)
        break;
    end
end
ind = i - 1;
if azimuth > a(end)
    ind = 1;
end

author: ShapeNet

repo: RenderForCNN


text:
function h = plotVectors(Y, varargin);
% Plot a matrix of amplitude vectors in a series of subplots.
%
% h = plotVectors(Y);
%
% Each column in Y is plotted in a separate subplot.
%
% If you enter additional input arguments, e.g.:
%
% h = plotVectors(X, Y, Z);
%
% where X, Y, and Z are all the same size, corresponding columns of each
% input will be plotted in the same subplot.
%
% ras, 12/2007.
if length(varargin) >= 1
	% check for additional data vectors
	try
		for n = 1:length(varargin)
			Y(:,:,n+1) = varargin{n};
		end
	catch
		size(Y)
		for n = 1:length(varargin)
			size(varargin{n})
		end
		error('Additional input vectors must be the same size.')
	end
end

% get color order for figure
co = plotVectors_colorOrder(size(Y, 3));

% open the figure
h(1) = figure('Color', 'w', 'Position', [100 100 600 500]);

N = size(Y, 2);
nrows = ceil( sqrt(N) );
ncols = ceil( N / nrows );

% want the position of the lower left-hand subplot,
% which is the only one that will have visible axes
lowerLeft = (nrows-1) * ncols + 1;

for n = 1:N
	h(n+1) = subplot(nrows, ncols, n);
	hold on

	for z = 1:size(Y, 3)
		plot(Y(:,n,z), 'LineWidth', 1);
	end

	axis tight
	set(gca, 'Box', 'off', 'TickDir', 'out');

	if n ~= lowerLeft
		axis off
	end

	setLineColors(co);
end

normAxes(h(2:end));

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

% /----------------------------------------------------------/ %
function co = plotVectors_colorOrder(N);
% create a nice color order for plotAmps (the default one isn't
% so great).

co = [0 0 0; ...
	  1 0 0; ...
	  0 0 1; ...
	  .5 .5 0; ...
	  0 .5 .5; ...
	  .4 .4 .4; ...
	  0 .6 0; ...
	  0 1 1; ...
	  1 0 1];

while size(co, 1) < N
	co = [co; co];
end

co = co(1:N,:);

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function[]=makefigs_psi2fields
%MAKEFIGS_PSI2FIELDS  Makes a sample figure for PSI2FIELDS.

load qgsnapshot,use qgsnapshot
psi=qgsnapshot.psi;

[cv,zeta,N,S,P]=psi2fields(x(2)-x(1),psi);
x=x/1000;
y=y/1000;
figure
subplot(2,2,1),jpcolor(x,y,psi),title('Streamfunction')
subplot(2,2,2),jpcolor(x,y,zeta),title('Vorticity')
subplot(2,2,3),jpcolor(x,y,P),title('Okubo-Weiss')
subplot(2,2,4),ii=1:10:length(x);
quiver(x(ii),y(ii),real(cv(ii,ii)),imag(cv(ii,ii))),title('Velocity')

for i=1:4
    subplot(2,2,i)
    axis equal, axis tight,xtick([-3:3]),ytick([-3:1:3])
    fontsize 16 14 14 14
end

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
function varargout = process_ssp2( varargin )
% PROCESS_SSP2: Artifact rejection for a group of recordings file (calculates SSP from FilesA and applies them to FilesB)
%
% USAGE:  OutputFiles = process_ssp2('Run', sProcess, sInputsA, sInputsB)
%                proj = process_ssp2('Compute', F, chanmask)
%           Projector = process_ssp2('BuildProjector', ListSsp, ProjStatus)
%           Projector = process_ssp2('ConvertOldFormat', OldProj)

% @=============================================================================
% 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, 2011-2022
%          Elizabeth Bock, 2011-2018

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription()
    % Description the process
    sProcess.Comment     = 'SSP: Generic';
    sProcess.Category    = 'Custom';
    sProcess.SubGroup    = 'Artifacts';
    sProcess.Index       = 302;
    sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/Tutorials/ArtifactsSsp?highlight=%28Process2%29#Troubleshooting';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'raw', 'data'};
    sProcess.OutputTypes = {'raw', 'raw'};
    sProcess.nInputs     = 2;
    sProcess.nMinFiles   = 1;
    sProcess.isSeparator = 0;

    % Notice
    sProcess.options.label1.Comment = 'Files A = Artifacts samples (raw or epoched)
 Files B = Files to clean (raw)

';
    sProcess.options.label1.Type    = 'label';
    % Time window
    sProcess.options.timewindow.Comment   = 'Time window: ';
    sProcess.options.timewindow.Type      = 'timewindow';
    sProcess.options.timewindow.Value     = [];
    sProcess.options.timewindow.InputTypes = {'raw'};
    % Event name
    sProcess.options.eventname.Comment = 'Event name (empty=continuous): ';
    sProcess.options.eventname.Type    = 'text';
    sProcess.options.eventname.Value   = 'blink';
    sProcess.options.eventname.InputTypes = {'raw'};
    % Event window
    sProcess.options.eventtime.Comment = 'Event window (ignore if no event): ';
    sProcess.options.eventtime.Type    = 'range';
    sProcess.options.eventtime.Value   = {[-.200, .200], 'ms', []};
    sProcess.options.eventtime.InputTypes = {'raw'};
    % Filter
    sProcess.options.bandpass.Comment = 'Frequency band: ';
    sProcess.options.bandpass.Type    = 'range';
    sProcess.options.bandpass.Value   = {[1.5, 15], 'Hz', 2};
    %  Number of components
    sProcess.options.nicacomp.Comment = 'Number of ICA components (0=default): ';
    sProcess.options.nicacomp.Type    = 'value';
    sProcess.options.nicacomp.Value   = {0, '', 0};
    sProcess.options.nicacomp.Hidden  = 1;
    % Sensor types
    sProcess.options.sensortypes.Comment = 'Sensor types or names (empty=all): ';
    sProcess.options.sensortypes.Type    = 'text';
    sProcess.options.sensortypes.Value   = 'MEG, EEG';
    % Use existing SSPs
    sProcess.options.usessp.Comment = 'Compute using existing SSP/ICA projectors';
    sProcess.options.usessp.Type    = 'checkbox';
    sProcess.options.usessp.Value   = 1;
    sProcess.options.usessp.InputTypes = {'raw'};
    % Ignore bad segments
    sProcess.options.ignorebad.Comment = 'Ignore bad segments';
    sProcess.options.ignorebad.Type    = 'checkbox';
    sProcess.options.ignorebad.Value   = 1;
    sProcess.options.ignorebad.InputTypes = {'raw'};
    sProcess.options.ignorebad.Hidden  = 1;
    % Save ERP
    sProcess.options.saveerp.Comment = 'Save averaged artifact in the database';
    sProcess.options.saveerp.Type    = 'checkbox';
    sProcess.options.saveerp.Value   = 0;
    % Method: Average or PCA
    sProcess.options.label2.Comment = '
Method to calculate the projectors:';
    sProcess.options.label2.Type    = 'label';
    sProcess.options.method.Comment = {'PCA: One component per sensor', 'Average: One component only'};
    sProcess.options.method.Type    = 'radio';
    sProcess.options.method.Value   = 1;
    % Examples: EOG, ECG
    sProcess.options.example.Comment = ['
Examples:
' ...
                                        '   - EOG: [-200,+200] ms, [1.5-15] Hz
' ...
                                        '   - ECG: [-40,+40] ms, [13-40] Hz

'];
    sProcess.options.example.Type    = 'label';
    % Default selection of components
    sProcess.options.select.Comment = 'Selected components:';
    sProcess.options.select.Type    = 'value';
    sProcess.options.select.Value   = {1, 'list', 0};
    sProcess.options.select.Hidden  = 1;
end

%% ===== FORMAT COMMENT =====
function Comment = FormatComment(sProcess)
    if isfield(sProcess.options, 'eventname') && ~isempty(sProcess.options.eventname.Value)
        Comment = ['SSP: ' sProcess.options.eventname.Value];
    else
        Comment = 'SSP';
    end
end

%% ===== RUN =====
function OutputFiles = Run(sProcess, sInputsA, sInputsB)
    OutputFiles = {};
    strOptions = '=> ';

    % ===== GET OPTIONS =====
    % Get time window
    if isfield(sProcess.options, 'timewindow') && isfield(sProcess.options.timewindow, 'Value') && iscell(sProcess.options.timewindow.Value) && ~isempty(sProcess.options.timewindow.Value)
        TimeWindow = sProcess.options.timewindow.Value{1};
    else
        TimeWindow = [];
    end
    % Get options values
    if isfield(sProcess.options, 'eventname') && isfield(sProcess.options.eventname, 'Value') && ~isempty(sProcess.options.eventname.Value)
        % Event name
        evtName = strtrim(sProcess.options.eventname.Value);
        if isempty(evtName)
            bst_report('Warning', sProcess, [], 'Event name is not specified: starting from the beginning of the file.');
        end
        strOptions = [strOptions, 'Event=' evtName ', '];
    else
        evtName = [];
    end
    % Event time window (only used if event is a point in time, not a window)
    if isfield(sProcess.options, 'eventtime') && isfield(sProcess.options.eventtime, 'Value') && ~isempty(sProcess.options.eventtime.Value) && ~isempty(sProcess.options.eventtime.Value{1})
        evtTimeWindow = sProcess.options.eventtime.Value{1};
        strOptions = [strOptions, sprintf('Epoch=[%1.3f,%1.3f]s, ', evtTimeWindow(1), evtTimeWindow(2))];
    else
        evtTimeWindow = [];
    end
    % Get downsampling value
    if isfield(sProcess.options, 'resample') && isfield(sProcess.options.resample, 'Value') && ~isempty(sProcess.options.resample.Value)
        resample = sProcess.options.resample.Value{1};
        strOptions = [strOptions, 'resample=' num2str(resample) ', '];
    else
        resample = 0;
    end
    % Get bandpass range
    BandPass = sProcess.options.bandpass.Value{1};
    if isempty(BandPass) || all(BandPass <= 0)
        bst_report('Info', sProcess, sInputsA, 'The bandpass filter was disabled.');
        BandPass = [];
        strOptions = [strOptions, 'NoBandPass, '];
    else
        strOptions = [strOptions, 'BandPass=[' num2str(BandPass(1)) ',' num2str(BandPass(2)), ']Hz, '];
    end
    % WARNING: Channel file should be the same for all the files
    if ~all(strcmpi(sInputsA(1).ChannelFile, {sInputsA.ChannelFile, sInputsB.ChannelFile}))
        bst_report('Warning', sProcess, sInputsA, [...
            'All the files should relate to the same channel file. ' ...
            'If the head position is different between two files, the values recorded for one channel do not mean ' ...
            'the same thing, therefore the calculation or the application of spatial filters (SSP) cannot be reliable.' 10 ...
            'If you decide to use the results of this process, you should understand this problem and check ' ...
            'your analysis for possible errors.']);
    end
    % Check input file types
    isRawA = strcmpi(sInputsA(1).FileType, 'raw');
    isRawB = strcmpi(sInputsB(1).FileType, 'raw');
    if ~isRawB
        bst_report('Error', sProcess, sInputsB, ['You can apply the calculated projectors only to a continuous file (Link to raw file).' 10 'Files B cannot contain any imported data blocks.']);
        return;
    end
    % Use existing SSPs
    if isfield(sProcess.options, 'usessp') && isfield(sProcess.options.usessp, 'Value') && ~isempty(sProcess.options.usessp.Value)
        UseSsp = sProcess.options.usessp.Value;
    else
        UseSsp = 0;
    end
    strOptions = [strOptions, 'UseSsp=' num2str(UseSsp) ', '];
    % Force selection of specific components
    if isfield(sProcess.options, 'select') && isfield(sProcess.options.select, 'Value') && ~isempty(sProcess.options.select.Value)
        if iscell(sProcess.options.select.Value)
            SelectComp = sProcess.options.select.Value{1};
        elseif isnumeric(sProcess.options.select.Value)
            SelectComp = sProcess.options.select.Value;
        else
            error('Invalid value for option "select".');
        end
    else
        SelectComp = [];
    end
    % Force selection of specific components
    if isfield(sProcess.options, 'nicacomp') && isfield(sProcess.options.nicacomp, 'Value') && ~isempty(sProcess.options.nicacomp.Value)
        nIcaComp = sProcess.options.nicacomp.Value{1};
        strOptions = [strOptions, 'nIcaComp=' num2str(nIcaComp) ', '];
    else
        nIcaComp = 0;
    end
    % Find components correlated to reference signals
    if isfield(sProcess.options, 'icasort')
        icaSort = sProcess.options.icasort.Value;
    else
        icaSort = [];
    end
    % Ignore bad segments
    if panel_record('IsEventBad', evtName)
        % If the event name contains the tag "bad", we need to include the bad segments
        bst_report('Info', sProcess, sInputsA, 'The selected event contains the tag "bad": including all the bad segments...');
        isIgnoreBad = 0;
    elseif isfield(sProcess.options, 'ignorebad') && isfield(sProcess.options.ignorebad, 'Value') && ~isempty(sProcess.options.ignorebad.Value)
        isIgnoreBad = sProcess.options.ignorebad.Value;
    else
        isIgnoreBad = 1;
    end
    strOptions = [strOptions, 'IgnoreBadSegments=' num2str(isIgnoreBad) 10];
    % Save ERP
    if isfield(sProcess.options, 'saveerp') && isfield(sProcess.options.saveerp, 'Value') && ~isempty(sProcess.options.saveerp.Value)
        SaveErp = sProcess.options.saveerp.Value;
    else
        SaveErp = 0;
    end
    % Method
    if isfield(sProcess.options, 'method') && isfield(sProcess.options.method, 'Value') && ~isempty(sProcess.options.method.Value)
        if ischar(sProcess.options.method.Value)
            Method = sProcess.options.method.Value;
        else
            switch (sProcess.options.method.Value)
                case 1,   Method = 'SSP_pca';
                case 2,   Method = 'SSP_mean';
                otherwise, error('Invalid method.');
            end
        end
    else
        Method = 'SSP_pca';
    end
    strOptions = [strOptions, '=> Method=' Method ', '];
    isICA = ~isempty(strfind(Method, 'ICA'));
    % Check method/inputs incompabilities
    if strcmpi(Method, 'SSP_mean') && isRawA && isempty(evtName)
        bst_report('Error', sProcess, sInputsA, 'The method SSP_mean cannot be applied to continuous recordings.');
        return;
    end

    % ===== COMPUTE BANDPASS FILTER =====
    % Get the time vector for the first file
    DataMat = in_bst_data(sInputsA(1).FileName, 'Time');
    sfreq = 1 ./ (DataMat.Time(2) - DataMat.Time(1));
    % Design band-pass filter
    if ~isempty(BandPass) && ~all(BandPass == 0)
        % If we need to resample the recordings
        if (resample > 0)
            filterFreq = resample;
        % Use the original sampling frequency
        else
            filterFreq = sfreq;
        end
        % Design the filter
        isMirror = 0;
        [tmp, FiltSpec] = process_bandpass('Compute', [], filterFreq, BandPass(1), BandPass(2), 'bst-hfilter-2019', isMirror);
        % Estimate transient period (before and after resampling)
        nTransientLoad    = round(FiltSpec.transient * sfreq);
        nTransientDiscard = round(FiltSpec.transient * filterFreq);
        % Show warning when computing epoched files
        if ~isRawA
            bst_report('Warning', sProcess, sInputsA, sprintf('Removing %1.3fs at the beginning and the end of each input for filtering.', FiltSpec.transient));
        end
    else
        nTransientLoad    = [];
        nTransientDiscard = [];
    end

    % ===== READ ARTIFACTS (FILES A) =====
    % Initialize progress bar
    bst_progress('text', 'Reading recordings...');
    progressPos = bst_progress('get');
    % Initialize concatenated data matrix
    F = {};
    Fref = {}; % for holding reference signals (EOG, ECG)
    iBad = [];
    iTimeZero = [];
    nSamples = 0;
    nMaxSamples = 100000;
    projChan = [];
    projChanNames = [];
    % Read together all the files in group A
    for iFile = 1:length(sInputsA)
        % ===== GET CHANNEL FILE =====
        if isempty(projChan) || ~isequal(projChan, sInputsA(iFile).ChannelFile)
            % Load channel file
            ChannelMat = in_bst_channel(sInputsA(iFile).ChannelFile);
            % Check for different channel structures in FilesA
            if isempty(projChanNames)
                projChanNames = {ChannelMat.Channel.Name};
            elseif ~isequal(projChanNames, {ChannelMat.Channel.Name})
                bst_report('Error', sProcess, sInputsA(iFile), 'All the files A must have the exact same list of channels.');
                return;
            end
            % Get channels to process
            iChannels = channel_find(ChannelMat.Channel, sProcess.options.sensortypes.Value);
            if isempty(iChannels)
                bst_report('Error', sProcess, sInputsA(iFile), 'No channels to process.');
                return;
            end
            % Warning if more than one channel type
            allTypes = unique({ChannelMat.Channel(iChannels).Type});
            if (length(allTypes) > 1) && ~isequal(allTypes, {'MEG GRAD', 'MEG MAG'})
                for iType = 1:length(allTypes)-1
                    allTypes{iType} = [allTypes{iType}, ', '];
                end
                bst_report('Warning', sProcess, sInputsA, ...
                    ['Mixing different channel types to compute the projector: ' [allTypes{:}], '.' 10 ...
                     'You should compute projectors separately for each sensor type.']);
            end
            if ~isempty(icaSort)
                iRef = channel_find(ChannelMat.Channel, icaSort);
                if isempty(iRef)
                    bst_report('Error', sProcess, sInputsA(iFile), sprintf('Channels %s not found.', icaSort));
                end
                iSensors = iChannels;
                iChannels = [iSensors iRef];
            end
        end

        % ===== GET DATA =====
        % === RAW: EVENTS ===
        if isRawA && ~isempty(evtName)
            % Load the raw file descriptor
            DataMat = in_bst_data(sInputsA(iFile).FileName, 'F', 'ChannelFlag');
            sFile = DataMat.F;
            sfreq_file = sFile.prop.sfreq;
            if isempty(sFile.events)
                bst_report('Error', sProcess, sInputsA(iFile), 'No events in the input file.');
                return;
            end
            % Get list of bad segments in file
            [badSeg, badEpochs] = panel_record('GetBadSegments', sFile);
            % Get the event to process
            events = sFile.events;
            iEvt = find(strcmpi({events.label}, evtName));
            nOcc = size(events(iEvt).times, 2);
            if isempty(iEvt) || (nOcc == 0)
                bst_report('Error', sProcess, sInputsA(iFile), ['Event type "' evtName '" not found, or has no occurrence.']);
                return;
            end
            % Extended / simple event
            isExtended = (size(events(iEvt).times, 1) == 2);
            % Simple events: get the samples to read around each event
            if ~isExtended
                evtSmpRange = round(evtTimeWindow .* sFile.prop.sfreq);
            % Extended: cannot work with "mean" option
            elseif strcmpi(Method, 'SSP_mean')
                bst_report('Error', sProcess, sInputsA, 'Method "SSP_mean" cannot be used for extended events.');
                return;
            % Extended: read the full block
            else
                evtSmpRange = [0 0];
            end
            % Add transients for bandpass
            if ~isempty(nTransientLoad)
                evtSmpRange = evtSmpRange + [-1 1] .* nTransientLoad;
            end
            % Reading options
            % NOTE: FORCE READING CLEAN DATA (Baseline correction + CTF compensators + Previous SSP)
            ImportOptions = db_template('ImportOptions');
            ImportOptions.ImportMode      = 'Event';
            ImportOptions.EventsTimeRange = evtSmpRange ./ sFile.prop.sfreq;
            ImportOptions.UseCtfComp      = 1;
            ImportOptions.UseSsp          = UseSsp;    % ADDED OPTION (FT: 27-Jun-2014) - Before we were always applying the previous SSPs
            ImportOptions.RemoveBaseline  = 'all';
            ImportOptions.DisplayMessages = 0;
            nInfoBad = 0;
            % Loop on each occurrence of the event
            for iOcc = 1:nOcc
                % Progress bar
                bst_progress('set', progressPos + round(iOcc / nOcc * 50));
                % Simple event: read a time window around the marker
                if ~isExtended
                    SamplesBounds = round(events(iEvt).times(1,iOcc) .* sFile.prop.sfreq) + evtSmpRange;
                % Extended event: read the full event
                else
                    SamplesBounds = round(events(iEvt).times(:,iOcc)' .* sFile.prop.sfreq) + evtSmpRange;
                end
                % Check that this epoch is within the segment of file to consider
                TimeBounds = SamplesBounds ./ sFile.prop.sfreq;
                if ~isempty(TimeWindow) && ((TimeBounds(1) < TimeWindow(1)) || (TimeBounds(2) > TimeWindow(2)))
                    continue;
                end
                % Check if this segment is outside of ALL the bad segments (either entirely before or entirely after)
                if isIgnoreBad && ~isempty(badSeg) && (~all((SamplesBounds(2) < badSeg(1,:)) | (SamplesBounds(1) > badSeg(2,:))))
                    nInfoBad = nInfoBad + 1;
                    continue;
                % Check if this this segment is  outside of the file bounds
                elseif (TimeBounds(1) < sFile.prop.times(1)) || (TimeBounds(2) > sFile.prop.times(2))
                    bst_report('Info', sProcess, sInputsA(iFile), sprintf('Event %s #%d is too close to the beginning or end of the file: ignored...', evtName, iOcc));
                    continue;
                end
                % Read block
                [Fevt, TimeVector] = in_fread(sFile, ChannelMat, events(iEvt).epochs(iOcc), SamplesBounds, [], ImportOptions);
                % SSP_mean: Check that we can get a time zero
                if strcmpi(Method, 'SSP_mean')
                    if ((TimeVector(1) > 0) || (TimeVector(end) < 0))
                        bst_report('Warning', sProcess, sInputsA(iFile), ['File #"' num2str(iFile) '" does not have a time t=0, it cannot be used for the method "mean".']);
                        continue;
                    end
                    % Get the sample where t=0
                    iTimeZero{end+1} = bst_closest(0, TimeVector);
                end
                % Concatenate to final matrix
                F{end+1} = Fevt;
                nSamples = nSamples + size(Fevt,2);
                % Check whether we read already all the samples we need
                if ~isICA && (nSamples >= nMaxSamples)
                    bst_report('Info', sProcess, sInputsA(iFile), sprintf('Reached the maximum number of samples at event %d / %d', iOcc, nOcc));
                    break;
                end
            end
            % Display message with number of events ignored
            if (nInfoBad > 0)
                bst_report('Info', sProcess, sInputsA(iFile), sprintf('Event %s: %d/%d events were in bad segments and ignored.', evtName, nInfoBad, nOcc));
            end

        % === RAW: CONTINUOUS ===
        elseif isRawA && isempty(evtName)
            % Load the raw file descriptor
            DataMat = in_bst_data(sInputsA(iFile).FileName, 'F', 'ChannelFlag');
            sFile = DataMat.F;
            sfreq_file = sFile.prop.sfreq;
            % Options for LoadInputFile()
            LoadOptions.IgnoreBad      = isIgnoreBad;  % From raw files: ignore the bad segments
            LoadOptions.ProcessName    = func2str(sProcess.Function);
            LoadOptions.RemoveBaseline = 'all';
            LoadOptions.UseSsp         = UseSsp;
            % Add transients for bandpass filter
            if ~isempty(nTransientLoad) && ~isempty(TimeWindow)
                rawTime = TimeWindow + [-1 1] .* (nTransientLoad / sfreq_file);
            else
                rawTime = TimeWindow;
            end
            % Load input signals
            [sMat, nSignals, iRows] = bst_process('LoadInputFile', sInputsA(iFile).FileName, [], rawTime, LoadOptions);
            if isempty(sMat.Data)
                bst_report('Error', sProcess, [], 'No data could be read from the input file.');
                return
            end
            % Use the loaded block
            F{end+1} = zeros(length(ChannelMat.Channel), size(sMat.Data,2));
            F{end}(iRows,:) = sMat.Data;
            TimeVector = sMat.Time;

        % === IMPORTED DATA ===
        else
            % Progress bar
            bst_progress('set', progressPos + round(iFile / length(sInputsA) * 50));
            % Load file
            DataMat = in_bst_data(sInputsA(iFile).FileName, 'F', 'ChannelFlag', 'Time');
            % Sampling frquency
            TimeVector = DataMat.Time;
            sfreq_file = 1 ./ (TimeVector(2) - TimeVector(1));
            % SSP_mean: Check that we can get a time zero
            if strcmpi(Method, 'SSP_mean')
                if ((TimeVector(1) > 0) || (TimeVector(end) < 0))
                    bst_report('Warning', sProcess, sInputsA(iFile), ['File #"' num2str(iFile) '" does not have a time t=0, it cannot be used for the method "mean".']);
                    continue;
                end
                % Get the sample where t=0
                iTimeZero{end+1} = bst_closest(0, TimeVector);
            end
            % Concatenate to final matrix
            F{end+1} = DataMat.F;
            nSamples = nSamples + size(DataMat.F,2);
            % Check whether we read already all the samples we need
            if ~isICA && (nSamples >= nMaxSamples)
                bst_report('Info', sProcess, sInputsA(iFile), sprintf('Reached the maximum number of samples at file %d / %d', iFile, length(sInputsA)));
                break;
            end
        end
        % Check that the sampling frequency is the same
        if (abs(sfreq_file - sfreq) > 1e-3)
            bst_report('Error', sProcess, sInputsA(iFile), 'Input files have different sampling rates.');
            break;
        end
        % Add bad channels in this file to the global list
        iBad = union(iBad, find(DataMat.ChannelFlag == -1));
    end
    % No data was read...
    if isempty(F)
        bst_report('Error', sProcess, sInputsA, 'No data could be read from the input file.');
        return;
    end

    % ===== DOWNSAMPLE =====
    % Downsample data to specified rate.
    if (resample > 0)
        bst_progress('text', 'Downsampling recordings...');
        bst_progress('set', progressPos + 25);
        for iFile = 1:length(F)
            [F{iFile}, TimeResample] = process_resample('Compute', F{iFile}, TimeVector, resample);
        end
        sfreq = resample;
        TimeVector = TimeResample;
    end

    % ===== COMPUTE AVERAGE =====
    % Compute the average of the data blocks (unfiltered)
    if SaveErp
        % Check that the size is the same for all the blocks
        if ~all(cellfun(@(c)isequal(size(c), size(F{1})), F))
            bst_report('Warning', sProcess, [], 'Input files have different sizes, cannot compute average');
            SaveErp = 0;
        else
            Favg = mean(cat(3, F{:}), 3);
        end
    end

    % ===== PROCESS DATA =====
    % Set the progress bar to 50%
    bst_progress('text', 'Processing recordings...');
    bst_progress('set', progressPos + 50);
    % Initializations for this loop
    Fzero = zeros(size(F{1},1), 1);
    nAvg = length(F);
    % Apply montages and filters
    for iBlock = 1:length(F)
        % Progress bar
        bst_progress('set', progressPos + 50 + round(iBlock / length(F) * 50));
        % Remove the bad channels from the matrix
        [iBad, iChanRemove] = intersect(iChannels, iBad);
        if ~isempty(iBad) && ~isempty(F{iBlock})
            iChannels(iChanRemove) = [];
        end
        % Filter recordings
        if ~isempty(BandPass) && ~all(BandPass == 0)
            F{iBlock}(iChannels,:) = process_bandpass('Compute', F{iBlock}(iChannels,:), sfreq, FiltSpec);
        end
        % Compute the average at time zero (filtered)
        if strcmpi(Method, 'SSP_mean')
            Fzero = Fzero + F{iBlock}(:,iTimeZero{iBlock});
        end
        % Filter recordings: Remove transients
        if ~isempty(BandPass) && ~all(BandPass == 0)
            F{iBlock} = F{iBlock}(:, (nTransientDiscard+1):(end-nTransientDiscard));
        end
        % Keep only the needed channels
        if ~isempty(icaSort)
            Fref{iBlock} = F{iBlock}(iRef,:);
            F{iBlock} = F{iBlock}(iSensors,:);
        else
            F{iBlock} = F{iBlock}(iChannels,:);
        end
    end
    % Comment
    nSamplesFinal = sum(cellfun(@(c)size(c,2), F));
    strOptions = [strOptions, 'Nsamples=' num2str(nSamplesFinal) ' from ' num2str(length(F)) ' blocks'];
    % Concatenate all the loaded data
    F = [F{:}];
    if ~isempty(icaSort)
        Fref = [Fref{:}];
        iChannels = iSensors;
    end
    % Error if nothing was loaded
    if isempty(F)
        bst_report('Error', sProcess, sInputsA, 'No data could be read from the input files.');
        return;
    elseif (size(F,2) < size(F,1))
        bst_report('Error', sProcess, sInputsA, 'Not enough data could be read from the input files.');
        return;
    elseif ~isICA && (size(F,1) < 32)
        bst_report('Warning', sProcess, sInputsA, ['You selected only ' num2str(size(F,1)) ' channels to compute the projectors. ' 10 ...
            'The SSP approach may not work correctly with a low number of sensors. ' ...
            'Make sure you selected correctly all the channels of data you have in your file.']);
    end
    % ICA: Check that the number of time samples is sufficient
    % Recommended number of time samples (from S Makeig):  Ntime/Nchan^2 >> 10
    nTime = size(F,2);
    if (nIcaComp > 0)
        nChan = nIcaComp;
    else
        nChan = size(F,1);
    end
    if isICA && (nTime / nChan^2 < 10)
        % Warning message
        strWarning = ['There is probably not enough data for a correct ICA decomposition.' 10 ...
            'Number time samples in input: '    num2str(nTime)      ' (' num2str(round(nTime / sfreq)) 's)' 10 ...
            'Recommend number of time samples:' num2str(10*nChan^2) ' (' num2str(round(10*nChan^2 / sfreq)) 's)'];
        % Add to the report
        bst_report('Warning', sProcess, [], strWarning);
        % Display on the console
        disp([10 'WARNING: ' strWarning]);
    end

    % ===== COMPUTE PROJECTORS =====
    Y = [];
    bst_progress('text', 'Computing projector...');
    switch (Method)

        % === SSP: PCA ===
        case 'SSP_pca'
            % === CHECK NUMBER OF SAMPLES ===
            % Minimum number of time samples required to estimate covariance
            nMinSmp = 10 * size(F,1);
            if (size(F,2) < nMinSmp)
                nTimePerBlock = length(TimeVector);
                nBlock = ceil(size(F,2) / nTimePerBlock);
                nBlockTotal = ceil(nMinSmp / nTimePerBlock);
                if isRawA && ~isempty(evtName)
                    errMsg = sprintf(' - Add %d events (Total: %d)', nBlockTotal - nBlock, nBlockTotal);
                    if ~isExtended
                        nAddTime = ceil((nMinSmp - size(F,2)) / nBlock / 2);
                        newTimeWin = round([evtSmpRange(1) - nAddTime, evtSmpRange(2) + nAddTime] ./ sFile.prop.sfreq .* 1000);
                        errMsg = sprintf([errMsg, 10, ' - Increase the time window around each event to [%d,%d] ms'], newTimeWin(1), newTimeWin(2));
                    end
                elseif isRawA && strcmpi(sFile.format, 'CTF') && length(sFile.epochs) > 1
                    errMsg = ' - Convert the input files to continuous.';
                else
                    errMsg = sprintf(' - Add %d files in the process list (Total: %d)', nBlockTotal - nBlock, nBlockTotal);
                end
                bst_report('Error', sProcess, sInputsA, ['Not enough time samples to compute projectors. You may:' 10 errMsg]);
                return;
            end

            % === COMPUTE PROJECTOR ===
            % Create channel mask matrix
            chanmask = zeros(length(ChannelMat.Channel), 1);
            chanmask(iChannels) = 1;
            % Call computation function
            proj = Compute(F, chanmask);
%             % Select the components with a singular value > threshold
%             if isempty(ForceSelect)
%                 singThresh = 0.12;
%                 proj.CompMask = double(proj.SingVal ./ sum(proj.SingVal) > singThresh);
%             % Force selection of specific components
%             else
%                 proj.CompMask(ForceSelect) = 1;
%             end
            % By default: first component selected
            if isempty(SelectComp)
                SelectComp = 1;
            end
            % Apply component selection (default: first component only)
            proj.CompMask(SelectComp) = 1;
            % Set the category as active if there is one selected component
            proj.Status = any(proj.CompMask);

        % === SSP: MEAN ===
        case 'SSP_mean'
            % Check number of samples
            if isempty(Fzero)
                bst_report('Error', sProcess, sInputsA, 'No data could be read from the input files.');
                return;
            else
                bst_report('Info', sProcess, sInputsA, sprintf('Computing projectors based on an average of %d events.', nAvg));
            end
            % Get the average at t=0
            Components = zeros(length(ChannelMat.Channel), 1);
            Components(iChannels) = Fzero(iChannels,:) ./ nAvg;
            % Normalize columns of the components
            Components = Components ./ sqrt(sum(Components .^2));
            % Build projector structure
            proj = db_template('projector');
            proj.Components = Components;
            proj.CompMask   = 1;
            proj.Status     = 1;
            proj.SingVal    = [];

        % === ICA: JADE ===
        case 'ICA_jade'
            bst_progress('text', 'Calling external function: jadeR()...');
            % Run decomposition
            if ~isempty(nIcaComp) && (nIcaComp ~= 0)
                W = jadeR(F, nIcaComp);
            else
                W = jadeR(F);
            end
            % Error handling
            if isempty(W)
                bst_report('Error', sProcess, sInputsA, 'Function "jadeR" did not return any results.');
                return;
            end

        % === INFOMAX ===
        case 'ICA_infomax'
            bst_progress('text', 'Calling external function: EEGLAB''s runica()...');
            % Remove the mean
            F = bst_bsxfun(@minus, F, mean(F,2));
            % Run EEGLAB ICA function
            if ~isempty(nIcaComp) && (nIcaComp ~= 0)
                [icaweights, icasphere] = runica(F, 'pca', nIcaComp, 'lrate', 0.001, 'extended', 1, 'interupt', 'off');
            else
                [icaweights, icasphere] = runica(F, 'lrate', 0.001, 'extended', 1, 'interupt', 'off');
            end
            % Error handling
            if isempty(icaweights) || isempty(icasphere)
                bst_report('Error', sProcess, sInputsA, 'Function "runica" did not return any results.');
                return;
            end
            % Reconstruct unmixing matrix
            W = icaweights * icasphere;

        % === ICA: PICARD ===
        case 'ICA_picard'
            bst_progress('text', 'Calling external function: picard()...');
            % Install picard plugin
            [isInstalled, errMsg] = bst_plugin('Install', 'picard');
            if ~isInstalled
                error(errMsg);
            end
            % Run decomposition
            if ~isempty(nIcaComp) && (nIcaComp ~= 0)
                [Y,W] = picard(F, 'pca', nIcaComp);
            else
                [Y,W] = picard(F);
            end

        % === ICA: FASTICA ===
        case 'ICA_fastica'
            bst_progress('text', 'Calling external function: fastica()...');
            % Install fastica plugin
            [isInstalled, errMsg] = bst_plugin('Install', 'fastica');
            if ~isInstalled
                error(errMsg);
            end
            % Scale the values to higher values, so that the pcamat function doesn't complain for small eigenvalues
            F = F ./ mean(abs(F(:)));
            % Run decomposition
            if ~isempty(nIcaComp) && (nIcaComp ~= 0)
                [M,W] = fastica(F, 'numOfIC', nIcaComp);
            else
                [M,W] = fastica(F);
            end

        otherwise
            bst_report('Error', sProcess, sInputsA, ['Invalid method: "' Method '".']);
            return;
    end

    % Finish assembling ICA results
    if ismember(Method, {'ICA_jade', 'ICA_infomax', 'ICA_picard', 'ICA_fastica'})
        % Fill with the missing channels with zeros
        Wall = zeros(length(ChannelMat.Channel), size(W,1));
        Wall(iChannels,:) = W';
        % Build projector structure
        proj = db_template('projector');
        proj.Components = Wall;
        proj.CompMask   = zeros(size(Wall,2), 1);   % No component selected by default
        proj.Status     = 1;
        proj.SingVal    = 'ICA';
        % Apply component selection (if set explicitly)
        if ~isempty(SelectComp)
            proj.CompMask(SelectComp) = 1;
        end
        % Compute IC
        if isempty(Y)
            Y = W * F;
        end
        % Sort ICA components
        if ~isempty(icaSort)
            % By correlation with reference channel
            C = bst_corrn(Fref, Y);
            [corrs, iSort] = sort(max(abs(C),[],1), 'descend');
            proj.Components = proj.Components(:,iSort);
        elseif ismember(Method, {'ICA_picard', 'ICA_fastica'})
            % By explained variance
            if diff(size(W)) == 0
                M = inv(W);
            else
                M = pinv(W);
            end
            var = sum(M.^2, 1) .* sum(Y.^2, 2)';
            [var, iSort] = sort(var, 'descend');
            proj.Components = proj.Components(:,iSort);
        end
    end

    % Modality used in the end
    AllMod = unique({ChannelMat.Channel(iChannels).Type});
    strMod = '';
    for iMod = 1:length(AllMod)
        strMod = [strMod, AllMod{iMod} '+'];
    end
    strMod = strMod(1:end-1);
    % Comment
    if ~isempty(evtName)
        proj.Comment = [evtName ': ' Method ', ' strMod ', ' datestr(clock)];
    else
        proj.Comment = [Method ': ' strMod  ', ' datestr(clock)];
    end

    % ===== APPLY PROJECTORS (FILES B) =====
    bst_progress('text', 'Applying projector to the recordings...');
    % Get all the channel files
    ChannelFiles = unique({sInputsB.ChannelFile});
    % Apply projectors to all the files in input
    for iFile = 1:length(ChannelFiles)
        % Load destination channel file
        ChannelMatB = in_bst_channel(ChannelFiles{iFile});
        destChanNames = {ChannelMatB.Channel.Name};
        % If channel names in projector and destination files do not match at all: ERROR
        [commonNames,I,J] = intersect(projChanNames, destChanNames);
        if isempty(commonNames)
            bst_report('Error', sProcess, sInputsB(iFile), 'List of channels are too different in the source and destination files.');
            return;
        end
        % If the channels list is not the same: re-order them to match the destination channel file
        if ~isequal(projChanNames, destChanNames)
            fixComponents = zeros(length(destChanNames), size(proj.Components,2));
            fixComponents(J,:) = proj.Components(I,:);
            proj.Components = fixComponents;
            bst_report('Warning', sProcess, sInputsB(iFile), sprintf('List of channels differ in FilesA (%d) and FilesB (%d). The common channels (%d) have been re-ordered to match the channels in FilesB.', length(projChanNames), length(destChanNames), length(commonNames)));
        end
        % Add projector to channel file
        [newproj, errMsg] = import_ssp(ChannelFiles{iFile}, proj, 1, 1, strOptions);
        if ~isempty(errMsg)
            bst_report('Error', sProcess, sInputsB(iFile), errMsg);
            return;
        end
    end

    % ===== SAVE THE AVERAGE =====
    if SaveErp && ~isempty(Favg)
        % Divide by number of average
        ChannelFlag = ones(size(Favg,1), 1);
        ChannelFlag(iBad) = -1;
        % Remove transients
        if ~isempty(nTransientDiscard)
            Favg = Favg(:, (nTransientDiscard+1):(end-nTransientDiscard));
            TimeVector = TimeVector((nTransientDiscard+1):(end-nTransientDiscard));
        end

        % === BEFORE ===
        % Create new output structure
        sOutput = db_template('datamat');
        sOutput.F           = Favg;
        sOutput.Comment     = [proj.Comment ' (before)'];
        sOutput.ChannelFlag = ChannelFlag;
        sOutput.Time        = TimeVector;
        sOutput.DataType    = 'recordings';
        sOutput.Device      = 'ArtifactERP';
        sOutput.nAvg        = nAvg;
        sOutput.Leff        = nAvg;
        % Get output study
        [tmp, iOutputStudy] = bst_process('GetOutputStudy', sProcess, sInputsB);
        sOutputStudy = bst_get('Study', iOutputStudy);
        % Output filename
        OutputFileERP = bst_process('GetNewFilename', bst_fileparts(sOutputStudy.FileName), 'data_artifact_before');
        % Save file
        bst_save(OutputFileERP, sOutput, 'v6');
        % Add file to database structure
        db_add_data(iOutputStudy, OutputFileERP, sOutput);

        % === SSP COMPONENTS ===
        if ~isICA
            % Force the first component to be used (allow multiple components by default)
            projErp = proj;
            projErp.Status(1) = 1;
            %projErp.CompMask = 0 * projErp.CompMask;
            projErp.CompMask(1) = 1;
            % Build projector
            Projector = BuildProjector(projErp, 1);
            % Apply to average data
            sOutput.F       = Projector * Favg;
            sOutput.Comment = [proj.Comment ' (after)'];
            % Output filename
            OutputFileERP = bst_process('GetNewFilename', bst_fileparts(sOutputStudy.FileName), 'data_artifact_after');
            % Save file
            bst_save(OutputFileERP, sOutput, 'v6');
            % Add file to database structure
            db_add_data(iOutputStudy, OutputFileERP, sOutput);
        % === ICA COMPONENTS ===
        else
            % Create file structure
            sOutput = db_template('matrixmat');
            sOutput.Value       = proj.Components' * Favg;
            sOutput.Comment     = [proj.Comment ' (ICA)'];
            sOutput.Time        = TimeVector;
            sOutput.ChannelFlag = [];
            sOutput.nAvg        = nAvg;
            sOutput.Leff        = nAvg;
            % Description of the signals: IC*
            sOutput.Description = cell(size(proj.Components,1),1);
            for i = 1:size(proj.Components,1)
                sOutput.Description{i} = sprintf('IC%d', i);
            end
            % Output filename
            OutputFileIC = bst_process('GetNewFilename', bst_fileparts(sOutputStudy.FileName), 'matrix_artifact_ica');
            % Save file
            bst_save(OutputFileIC, sOutput, 'v6');
            % Add file to database structure
            db_add_data(iOutputStudy, OutputFileIC, sOutput);
        end
    end
    % Return all the input files
    OutputFiles = {sInputsB.FileName};
end

%% ===== COMPUTE PROJECTOR =====
function proj = Compute(F, chanmask)
    % SVD decomposition
    [U,S,V] = svd(F, 'econ');
    % Create projector structure
    proj = db_template('projector');
    % Keep all the dimensions
    nChannel = length(chanmask);
    nProj    = size(U,2);
    proj.Components = zeros(nChannel, nProj);
    proj.Components(chanmask == 1,:) = U;
    % Other fields
    proj.SingVal  = diag(S)';
    proj.CompMask = zeros(1,nProj);
    proj.Status   = 1;
end

%% ===== BUILD PROJECTOR ====
% Combine all the projectors in decomposed form to create a [nChan x nChan] matrix
%
% USAGE: Projector = process_ssp2('BuildProjector', ListProj, ProjStatus)
%
% INPUT:
%    - ListProj   : Array of db_template('projector')
%    - ProjStatus : List of the projector status to include
%                   0 = not applied, not used
%                   1 = have to be used, but still have to be applied on the fly
%                   2 = already applied, saved in the file, not revertible
% OUTPUT:
%    - Projector  : [nChannels x nChannels] matrix, projector in the condensed form (I-UUt)
%
% COMMENTS:
%    There are 5 categories of projectors:
%    - SSP_pca:   CompMask=[Ncomp x 1],   SingVal=[Ncomp x 1],   Components=[Nchan x Ncomp]=U
%    - SSP_mean:  CompMask=1,             SingVal=[],            Components=[Nchan x 1]=U
%    - ICA:       CompMask=[Ncomp x 1],   SingVal='ICA',         Components=[Nchan x Ncomp]=W'
%    - REF:       CompMask=[],            SingVal='REF',         Components=[Nchan x Ncomp]=Wmontage
%    - Other:     CompMask=[],            SingVal=[],            Components=[Nchan x Nchan]=Projector=I-UUt
%
%    Description of the notations used here:
%    - W: Unmixing matrix  [Nelectrodes x Ncomponents]
%    - Winv = pinv(W) = [Ncomponents x Nelectrodes]
%    - In EEGLAB:  W = icaweights * icasphere;
%    - Activations_IC = W * Data
%    - CleanData = Winv(:,iComp) * Activations(iComp,:)
%                = Winv(:,iComp) * W(iComp,:) * Data
%
function Projector = BuildProjector(ListProj, ProjStatus) %#ok<*DEFNU>
    % Call on an old form of Projector (I-UUt)
    if ~isstruct(ListProj)
        Projector = ListProj;
        return
    end
    % Initialize matrices
    nChannel = size(ListProj(1).Components,1);

    % === SORT PROJECTORS ===
    % Find projectors to group or remove
    iProjDel = [];
    iProjSsp = [];
    U = [];
    for i = 1:length(ListProj)
        % Is entry not selected: skip
        if ~ismember(ListProj(i).Status, ProjStatus) || (~isempty(ListProj(i).CompMask) && all(ListProj(i).CompMask == 0))
            iProjDel(end+1) = i;
        % New SSP: Stack selected vectors all together
        elseif ~isempty(ListProj(i).CompMask) && ~isequal(ListProj(i).SingVal, 'ICA') && ~isequal(ListProj(i).SingVal, 'REF')
            iProjSsp(end+1) = i;
            U = [U, ListProj(i).Components(:,ListProj(i).CompMask == 1)];
        end
    end

    % === ORTHOGONALIZE SSP VECTORS ===
    if ~isempty(U) && ~all(U(:) == 0)
        % Reorthogonalize the vectors
        [U,S,V] = svd(U,0);
        S = diag(S);
        % Enforce strict zero values (because Matlab SVD function can randomly return small values instead of zeros !!!)
        % If not, sometimes it adds small contributions of some channels that have nothing to do with the projectors (<1e-16).
        % If this added channel contains high values (eg. Stim channel in Volts) it can corrupt significantly the MEG values.
        iZero = find((abs(U) < 1e-15) & (U ~= 0));
        if ~isempty(iZero)
            U(iZero) = 0;
        end
        % Throw away the linearly dependent guys (threshold on singular values: 0.01 * the first one)
        iThresh = find(S < 0.01 * S(1),1);
        if ~isempty(iThresh)
            disp(sprintf('SSP> %d linearly dependent vectors removed...', size(U,2)-iThresh+1));
            U = U(:, 1:iThresh-1);
        end
        % Compute projector in the form: I-UUt
        SspProj = eye(nChannel) - U*U';
        % Copy the projector in the Components field of the first SSP entry
        ListProj(iProjSsp(1)).Components = SspProj;
        % Remove all the other SSP entries
        if (length(iProjSsp) >= 2)
            iProjDel = [iProjDel, iProjSsp(2:end)];
        end
    % Remove all the SSP projectors
    elseif ~isempty(iProjSsp)
        iProjDel = [iProjDel, iProjSsp];
    end

    % === CLEANING LIST ===
    % Remove unwanted projectors
    ListProj(iProjDel) = [];
    % Nothing left: return
    if isempty(ListProj)
        Projector = [];
        return;
    end

    % === BUILD FINAL PROJECTOR ===
    % Initialize returned projector as an identity matrix
    Projector = eye(nChannel);
    % Add the projectors in the order of appearance
    for i = 1:length(ListProj)
        % ICA
        if isequal(ListProj(i).SingVal, 'ICA')
            % Get selected channels (find the non-zero channels)
            iChan = find(any(ListProj(i).Components ~= 0, 2));
            % Get selected components
            iComp = find(ListProj(i).CompMask == 1);
            % Initialize projector
            P = eye(size(ListProj(i).Components,1));
            % Compute projector
            W = ListProj(i).Components(iChan,:)';
            Winv = pinv(W);
            P(iChan,iChan) = eye(size(W,2)) - Winv(:,iComp) * W(iComp,:);
            % Check if there are any complex values in the projector
            if any(~isreal(P))
                warning('WARNING: ICA components contain complex values. Something went wrong in their computation.');
                P = real(P);
            end
        % Other projectors
        else
            P = ListProj(i).Components;
        end
        % Add to the final operator
        Projector = P * Projector;
    end
    % Make sure it is not an identity matrix
    if isequal(Projector, eye(nChannel))
        Projector = [];
    end
end

%% ===== CONVERT OLD FORMAT =====
% Old format: I - UUt
% New format: Structure with decomposed form (U, maskU...)
function proj = ConvertOldFormat(OldProj)
    if isempty(OldProj)
        proj = [];
    elseif ~isstruct(OldProj)
        proj = db_template('projector');
        proj.Components = OldProj;
        proj.Comment    = 'Unnamed';
        proj.Status     = 1;
    else
        proj = OldProj;
    end
end

author: brainstorm-tools


text:
function traj = simulateCannon(init,param)

v0 = init.speed;
th0 = init.angle;

c = param.c;  %Quadratic drag coefficient
nGrid = param.nGrid;

% Set up initial conditions for ode45
x0 = 0;  y0 = 0;  %Start at the origin
dx0 = v0*cos(th0);
dy0 = v0*sin(th0);
if dy0 < 0, error('Cannot point cannon through ground! sin(th0) > 0 required.'); end;

% Set up arguments for ode45
userFun = @(t,z)cannonDynamics(t,z,c);  %Dynamics function
tSpan = [0,100];  %Never plan on reaching final time
z0 = [x0;y0;dx0;dy0];  %Initial condition
options = odeset('Events',@groundEvent,'Vectorized','on');

% Run a simulation
sol = ode45(userFun, tSpan, z0, options);

% Extract the solution on uniform grid:
traj.t = linspace(sol.x(1), sol.x(end), nGrid);
z = deval(sol,traj.t);
traj.x = z(1,:);
traj.y = z(2,:);
traj.dx = z(3,:);
traj.dy = z(4,:);

end

author: MatthewPeterKelly

repo: dscTutorials


text:
function [handles, key_points] = draw_pie_wedge(T1, T2, myradius, varargin)
% Draws a pie wedge with color fill, for polar plots and circular charts
%
% :Usage:
% ::
%
%     handles = draw_pie_wedge(T1, T2, myradius, ['linecolor'], [r g b], ['fillcolor'], [r g b])
%
% ..
%     Dec 2017, Tor Wager
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
% ..
%
% :Inputs:
%
%   **T1:**
%        Starting location, in radians. 0 = to the right on plot; 2*pi is
%        full-circle, also to the right on plot.
%
%   **T2:**
%        Ending location, in radians. 0 = to the right on plot; 2*pi is
%        full-circle, also to the right on plot.
%
% :Optional Inputs:
%   **'linecolor':**
%        Followed by [r g b] triplet specifying line color, or 'none'
%
%   **'fillcolor:**
%        Followed by [r g b] triplet specifying fill color, or 'none'
%
%   **'outside_radius_offset'**
%        Followed by outside_radius_offset value. See below.
%
%   **'linewidth':**
%       Followed by line width
%
%   **'stripes':**
%       Turn on radial striping; can be followed by value true or false
%
%   **'stripedensity':**
%       Number of stripes/lines to draw in wedge
%
% :Outputs:
%
%   **handles:**
%        A structure of line and fill handles
%
%   **key_points:**
%        x and y locations for the mid-point of the arc, inner mid-point of the
%        wedge, and outer mid-point offset by a constant value. This is
%        set by outside_radius_offset, and is useful for text labels
%
% :Examples:
% ::
%
% myradius = 3;
% T1 = 1; % start, in radians
% T2 = 4; % end, in radians
%
% handles = draw_pie_wedge(1, 3, 2); % Draw wedge from radian=1 to 3, radius = 2
%
% Draw reference circle with radius = 2
% handles = draw_pie_wedge(0, 2*pi, 2, 'linecolor', [0 0 0], 'fillcolor', 'none');
%
% Now draw a wedge on top of that:
% handles = draw_pie_wedge(1, 1.5, 2.5, 'linecolor', [.7 .7 0], 'fillcolor', [.3 .7 .7]);
%
% Draw with dense striping pattern:
% handles = draw_pie_wedge(1, 1.5, 2.5, 'linecolor', [.7 .7 0], 'fillcolor', [.3 .7 .7], 'stripes', 'stripedensity', 40);
%
% Draw a series of wedges in different colors:
% n_categories = 7;
% st = linspace(0, 2*pi, n_categories+1); % starting vals in radians
%
% en = st(2:end);
% st = st(1:end-1);
%
% colors = seaborn_colors(n_categories + 1); % add one to make colors more diff
% clear handles
%
% figure; hold on;
%
% for i = 1:n_categories
%
%     handles(i) = draw_pie_wedge(st(i), en(i), i, 'linecolor', colors{i} ./ 2, 'fillcolor', colors{i});
%
% end
%
% axis equal
%
% For another extended example with text labels, see tor_wedge_plot.m

linecolor = [1 0 0];
fillcolor = [.7 .3 .7];
linewidth = 2;
outside_radius_offset = .04; % Value to offset by when calculating outside-wedge midpoint for labels
dostripes = false;
stripedensity = 12;

for i = 1:length(varargin)
    if ischar(varargin{i})
        switch varargin{i}

            case 'linecolor', linecolor = varargin{i+1}; varargin{i+1} = [];
            case 'fillcolor', fillcolor = varargin{i+1}; varargin{i+1} = [];

            case 'linewidth', linewidth = varargin{i+1}; varargin{i+1} = [];

            case 'stripes'
                dostripes = true;
                % Optional argument true or false following keyword
                if length(varargin) > i && ~ischar(varargin{i+1}) % param value entered in next arg
                    dostripes = varargin{i+1}; varargin{i+1} = [];
                end

            case 'stripedensity', stripedensity = varargin{i+1}; varargin{i+1} = [];

            case 'outside_radius_offset', outside_radius_offset = varargin{i+1}; varargin{i+1} = [];
            otherwise, warning(['Unknown input string option:' varargin{i}]);
        end
    end
end

% This would draw a reference circle
%hold on
%h0 = drawCircle(0, 0, myradius);
%set(h0, 'Color', 'k');

% Circle Arc
% ------------------------------------------
%h = drawCircleArc(0, 0, myradius, T1, T2);

[h xt yt] = drawArc(T1, T2, myradius);

% get midpoint of arc
% ------------------------------------------
tmid = (T1+T2)/2;
xmid_arc = 0 + myradius*cos(tmid);
ymid_arc = 0 + myradius*sin(tmid);

% Middle of pie
xmid_pie = 0 + myradius/2*cos(tmid);
ymid_pie = 0 + myradius/2*sin(tmid);

% outside arc
rad2 = myradius + outside_radius_offset; %.1 * myradius;
xmid_outside = 0 + rad2*cos(tmid);
ymid_outside = 0 + rad2*sin(tmid);

% others
radialangle = rad2deg(tmid); % radial angle in degrees
tangentangle = rad2deg(tmid + (pi/2)); % angle of tangent, for text

key_points = struct('tmid_radians', tmid, 'radialangle', radialangle, 'tangentangle', tangentangle, ...
    'xmid_arc', xmid_arc, 'ymid_arc', ymid_arc, 'xmid_outside', xmid_outside, ...
    'ymid_outside', ymid_outside, 'xmid_pie', xmid_pie);

set(h, 'Color', linecolor, 'LineWidth', linewidth)

% Pie wedge
% ------------------------------------------
hold on
[x1,y1] = pol2cart(T1, myradius);
h(2) = plot([0 x1], [0 y1], 'Color', linecolor, 'LineWidth', linewidth);

[x2,y2] = pol2cart(T2, myradius);
h(3) = plot([0 x2], [0 y2], 'Color', linecolor, 'LineWidth', linewidth);

% Fill
% ------------------------------------------

if ischar(fillcolor) && strcmp(fillcolor, 'none')
    % skip it
    hfill = [];
else
    hfill = fill([xt 0], [yt 0], fillcolor, 'FaceAlpha', .65);
end

handles = struct('line_han', h, 'fill_han', hfill);

% Striping
% ------------------------------------------
if dostripes

    handles.stripe_han = drawStripes(T1, T2, myradius, stripedensity);
    set(handles.stripe_han, 'Color', linecolor, 'LineWidth', linewidth)

else

    handles.stripe_han = [];

end

end % main function

function [h xt yt] = drawArc(T1, T2, myradius)

% key bits
t = T1:.01:T2;          % Arc
xt = 0 + myradius*cos(t);
yt = 0 + myradius*sin(t);
xt = [xt 0+myradius*cos(T2)];
yt = [yt 0+myradius*sin(T2)];

h = line(xt, yt);

end

function h = drawStripes(T1, T2, myradius, stripedensity)

r = linspace(0, myradius, stripedensity);

for i = 1:length(r)

    h(i) = drawArc(T1, T2, r(i));

end

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
% Digital Video Stabilization and Rolling Shutter Correction using Gyroscopes
% Copyright (C) 2011 Alexandre Karpenko
%
% 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
% any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

function W = imalign(T, I, p)
% align the template image T to the input image I

if size(T,3) == 3
    T = rgb2gray(T);
end
if size(I,3) == 3
    I = rgb2gray(I);
end
if max(T(:)) > 1
    T = double(T) / 255;
end
if max(I(:)) > 1
    I = double(I) / 255;
end

gradTx = filter2(fspecial('sobel')', T);
gradTy = filter2(fspecial('sobel') , T);
%{
figure(1); showim(reshape(gradTx, size(T)));
figure(2); showim(reshape(gradTy, size(T)));
pause;
%}

x = repmat((1:size(T,2)) , size(T,1), 1);
y = repmat((1:size(T,1))', 1, size(T,2));

gradT_dWdp = [gradTx(:) .* x(:) ...
              gradTy(:) .* x(:) ...
              gradTx(:) .* y(:) ...
              gradTy(:) .* y(:) ...
              gradTx(:)         ...
              gradTy(:)         ];

H = gradT_dWdp' * gradT_dWdp;

p = p(:);

for i=1:100
    W = [1+p(1)  p(3)  p(5); ...
          p(2)  1+p(4) p(6)];

    figure(1); clf;
    imshow(imtransform(T,  maketform('affine', W'), 'XData', [1 size(I,2)], 'YData', [1 size(I,1)]));
    hold on;
    h = imshow(I);
    set(h, 'AlphaData', 0.6);
    pause(0.01);

    WI = imtransform(I,  maketform('affine', W'), 'XData', [1 size(I,2)], 'YData', [1 size(I,1)], 'FillValues', NaN);
    errorIm = WI - T;
    errorIm(isnan(errorIm)) = 0;
    delp = H \ (gradT_dWdp' * errorIm(:));
    p = p + delp;
    display(p);
end

author: alex-golts

repo: Video-Stabilization


text:
%%*****************************************************************************
%% sqlpmisc:
%% unscale and produce infeasibility certificates if appropriate
%%
%% SDPT3: version 3.1
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last Modified: 16 Sep 2004.
%%*****************************************************************************

   function [X,y,Z,termcode,resid,reldist,msg] = sqlpmisc(blk,At,C,b,X,y,Z,permZ,param);

   termcode    = param.termcode;
   iter        = param.iter;
   obj         = param.obj;
   relgap      = param.relgap;
   prim_infeas = param.prim_infeas;
   dual_infeas = param.dual_infeas;
   homRd       = param.homRd;
   homrp       = param.homrp;
   AX          = param.AX;
   ZpATynorm   = param.ZpATynorm;
   m0          = param.m0;
   indeprows   = param.indeprows;
   normX0      = param.normX0;
   normZ0      = param.normZ0;
   inftol      = param.inftol;
   maxit       = param.maxit;
   scale_data  = param.scale_data;
   printlevel  = param.printlevel;
%%
   resid = []; reldist = []; msg = [];
   if (scale_data)
      normA = param.normA; normC = param.normC; normb = param.normb;
   else
      normA = 1; normC = 1; normb = 1;
   end
   Anorm = ops(At,'norm'); xnorm = ops(X,'norm'); ynorm = norm(y);
   infeas = max(prim_infeas,dual_infeas);
%%
   if (iter >= maxit)
      termcode = -6;
      msg = 'sqlp stop: maximum number of iterations reached';
      if (printlevel); fprintf('\n  %s',msg); end
   end
   if (termcode <= 0)
      %%
      %% To detect near-infeasibility when the algorithm provides
      %% a "better" certificate of infeasibility than of optimality.
      %%
      err = max(infeas,relgap);
      iflag = 0;
      if (obj(2) > 0)
         if (homRd < 0.1*sqrt(err*inftol))
            iflag = 1;
            msg = sprintf('prim_inf,dual_inf,relgap = %3.2e, %3.2e, %3.2e',...
		  prim_infeas,dual_infeas,relgap);
            if (printlevel); fprintf('\n  %s',msg); end
            termcode = 1;
         end
      end
      if (obj(1) < 0)
         if (homrp < 0.1*sqrt(err*inftol))
            iflag = 1;
            msg = sprintf('prim_inf,dual_inf,relgap = %3.2e, %3.2e, %3.2e',...
                  prim_infeas,dual_infeas,relgap);
            if (printlevel); fprintf('\n  %s',msg); end
            termcode = 2;
         end
      end
      if (iflag == 0)
         if (scale_data == 1)
            X = ops(ops(X,'./',normA),'*',normb);
            y = y*normC;
            Z = ops(ops(Z,'.*',normA),'*',normC);
         end
      end
   end
   if (termcode == 1) & (iter > 3)
      msg = 'sqlp stop: primal problem is suspected of being infeasible';
      if (printlevel); fprintf('\n  %s',msg); end
      if (scale_data == 1)
         X = ops(X,'./',normA); b = b*normb;
      end
      rby = 1/(b'*y); y = rby*y; Z = ops(Z,'*',rby);
      resid = ZpATynorm * rby;
      reldist = ZpATynorm/(Anorm*ynorm);
   end
   if (termcode == 2) & (iter > 3)
      msg = 'sqlp stop: dual problem is suspected of being infeasible';
      if (printlevel); fprintf('\n  %s',msg); end
      if (scale_data == 1)
         C = ops(C,'.*',normC);
         Z = ops(Z,'.*',normA);
      end
      tCX = blktrace(blk,C,X);
      X = ops(X,'*',1/(-tCX));
      resid = norm(AX)/(-tCX);
      reldist = norm(AX)/(Anorm*xnorm);
   end
   if (termcode == 3)
      maxblowup = max(ops(X,'norm')/normX0,ops(Z,'norm')/normZ0);
      msg = sprintf('sqlp stop: primal or dual is diverging, %3.1e',maxblowup);
      if (printlevel); fprintf('\n  %s',msg); end
   end
   [X,Z] = unperm(blk,permZ,X,Z);
   if ~isempty(indeprows)
      ytmp = zeros(m0,1);
      ytmp(indeprows) = y;
      y = ytmp;
   end
%%*****************************************************************************
%% unperm: undo the permutations applied in validate.
%%
%% [X,Z] = unperm(blk,permZ,X,Z);
%%
%% undoes the permutation introduced in validate.
%%*****************************************************************************

  function [X,Z] = unperm(blk,permZ,X,Z);
%%
  for p = 1:size(blk,1)
     if (strcmp(blk{p,1},'s') & ~isempty(permZ{p}))
        per = permZ{p};
        X{p} = X{p}(per,per);
        Z{p} = Z{p}(per,per);
     end
  end
%%*****************************************************************************

repo: REST


text:
function p=drex_invlogit(a)
%DREX subfunction
%Written by Issam El Naqa 2003-2005
%Extracted for generalized use 2005, AJH
%
% Copyright 2010, Joseph O. Deasy, on behalf of the DREES development team.
%
% This file is part of the Dose Response Explorer System (DREES).
%
% DREES development has been led by:  Issam El Naqa, Aditya Apte, Gita Suneja, and Joseph O. Deasy.
%
% DREES has been financially supported by the US National Institutes of Health under multiple grants.
%
% DREES is distributed under the terms of the Lesser GNU Public License.
%
%     This version of DREES 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.
%
% DREES 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 DREES.  If not, see <http://www.gnu.org/licenses/>.

p = 1 ./ (1 + exp(-a));

return

author: mvallieres

repo: radiomics