text:
function rec = perform_bp_benchmark(P,W, slist, ntrial)

% perform_bp_benchmark - benchmark basis pursuit recovery
%
% rec = perform_bp_benchmark(P,W, slist, ntrial);
%
%   Try to solve many (ntrial) problems of the kind
%       x = argmin_a |a|_1    s.t. U*a=U*x0
%   where x_0 is an s-sparse signal.
%   U = P*W
%   The report the number of coordinate of x0 acuatly recovered in x.
%
%   rec(i) is the average ratio of recoverred dirac for sparsity slist(i).
%   rec=1 means perfect recovery and rec=0 means no recovery (bad).
%
%   Copyright (c) 2007 Gabriel Peyr?

if nargin<2
    slist = 10:18;
end
if nargin<3
    ntrial = 50;
end

% BP options
maxIters = 20; lambda = 0; OptTol = 1e-5;

U = P*W;
m = size(U,2);
rec = [];

% detection threshold
T = 0.1;

for is = 1:length(slist)
    progressbar(is,length(slist));
    % sparsity
    s = slist(is);
    e = 0;
    for k = 1:ntrial
        x = compute_rand_sparse(m,s);
        y = U*x;
        x0 = SolveBP(U, y, m, maxIters, lambda, OptTol);
        % compute the number of coordinates recovered
%        [tmp,I] = sort(x0);
%        x0 = x0 * 0; x0(I(end-s+1:end)) = 1;
%        e = e + sum(x0.*x)/s;
        e = e + norm( (W*x)-(W*x0) );
    end
    rec(end+1) = e/ntrial;
end

author: gpeyre

repo: matlab-toolboxes


text:
data_path = tracker_param.seq_path;
location = tracker_param.init_rect;
s1 = tracker_param.s1;
s2 = tracker_param.s2;
roi_size = tracker_param.roi_size;
pf_param = tracker_param.pf_param ;
gnet_solver_def_file = tracker_param.gnet_solver_def_file;
snet_solver_def_file = tracker_param.snet_solver_def_file;
im1_id = tracker_param.startFrame;
ch_num = tracker_param.ch_num;
fnum = tracker_param.endFrame;
check_num = im1_id + 19;

repo: FCNT


text:
%%*************************************************************************
%% mybicgstab
%%
%% [xx,resnrm,flag] = mybicgstab(A,b,M1,tol,maxit)
%%
%% iterate on  bb - (M1)*AA*x
%%
%% r = b-A*xtrue;
%%
%%*************************************************************************

function [xx,resnrm,flag] = mybicgstab(A,b,M1,tol,maxit,printlevel)

N = length(b);
if (nargin < 6); printlevel = 1; end
if (nargin < 5) || isempty(maxit); maxit = max(30,length(A.mat22)); end;
if (nargin < 4) || isempty(tol); tol = 1e-10; end;
tolb = min(1e-4,tol*norm(b));
flag = 1;

x = zeros(N,1);
if (norm(x))
    if isstruct(A); r = b-matvec(A,x); else r = b-mexMatvec(A,x); end;
else
    r =b;
end
err = norm(r); resnrm(1) = err;  minresnrm = err; xx = x;
%%if (err < 1e-3*tolb); return; end

omega  = 1.0;
r_tld = r;
%%
%%
%%
breakyes = 0;
smtol = 1e-40;
for iter = 1:maxit,

    rho   = (r_tld'*r);
    if (abs(rho) < smtol)
        flag = 2;
        if (printlevel); fprintf('*'); end;
        breakyes = 1;
        break;
    end
    if (iter > 1)
        beta  = (rho/rho_1)* (alp/omega);
        p = r + beta*(p - omega*v);
    else
        p = r;
    end
    p_hat = precond(A,M1,p);
    if isstruct(A); v = matvec(A,p_hat); else v = mexMatvec(A,p_hat); end;
    alp = rho / (r_tld'*v);
    s = r - alp*v;
    %%
    s_hat = precond(A,M1,s);
    if isstruct(A); t = matvec(A,s_hat); else t = mexMatvec(A,s_hat); end;
    omega = (t'*s) / (t'*t);
    x = x + alp*p_hat + omega*s_hat;
    r = s - omega*t;
    rho_1 = rho;
    %%
    %% check convergence
    %%
    err = norm(r); resnrm(iter+1) = err; %#ok
    if (err < minresnrm);
        xx = x; minresnrm = err;
    end
    if (err < tolb)
        break;
    end
    if (err > 10*minresnrm)
        if (printlevel); fprintf('^'); end
        breakyes = 2;
        break;
    end
    if (abs(omega) < smtol)
        flag = 2;
        if (printlevel); fprintf('*'); end
        breakyes = 1;
        break;
    end
end
if (~breakyes) && (printlevel >=3); fprintf(' '); end
%%
%%*************************************************************************
%%*************************************************************************
%% precond:
%%*************************************************************************

function Mx = precond(A,L,x)

m = L.matdim; m2 = length(x)-m;
Mx = zeros(length(x),1);

for iter = 1
    if norm(Mx); r = x - matvec(A,Mx); else r = x; end
    if (m2 > 0)
        r1 = full(r(1:m));
    else
        r1 = full(r);
    end
    if (m2 > 0)
        r2 = r(m+1:m+m2);
        w = linsysolvefun(L,r1);
        z = mexMatvec(A.mat12,w,1) - r2;
        z = L.Mu \ (L.Ml \ (L.Mp*z));
        r1 = r1 - mexMatvec(A.mat12,z);
    end
    d = linsysolvefun(L,r1);
    if (m2 > 0)
        d = [d; z]; %#ok
    end
    Mx = Mx + d;
end
%%*************************************************************************
%%*************************************************************************
%% matvec: matrix-vector multiply.
%% matrix = [A.mat11, A.mat12; A.mat12', A.mat22]
%%*************************************************************************

function Ax = matvec(A,x)

m = length(A.mat11); m2 = length(x)-m;
if issparse(x); x = full(x); end
if (m2 > 0)
    x1 = x(1:m);
else
    x1 = x;
end
Ax = mexMatvec(A.mat11,x1);
if (m2 > 0)
    x2 = x(m+1:m+m2);
    Ax = Ax + mexMatvec(A.mat12,x2);
    Ax2 = mexMatvec(A.mat12,x1,1) + mexMatvec(A.mat22,x2);
    Ax = [full(Ax); full(Ax2)];
end
%%*************************************************************************

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function feat = poly2featuresOnly(u, v, imsize, v0, yc, f)

imh = imsize(1);
imw = imsize(2);

if isempty(f)
    f = 1.38; %S*max(size(im)) / imh;
end

u = u(:)';
v = v(:)';

[tmp, ind] = min(u);
u = [u(ind:end)  u(1:ind-1)];
v = [v(ind:end)  v(1:ind-1)];

u = (u - imw/2) ./ imh;
v = 1 - (v ./ imh);

[v1, ind1] = min(v); % lowest point on object in image

[footx, footz] = computeGroundPosition([min(u) max(u)], [v1 v1], v0, yc, f);
footz(2) = footz(1) + footx(2)-footx(1);

[cx, cz] = computeGroundPosition(u, v, v0, yc, f);

data.u = (u*imh+imw/2)';
data.v = ((1-v)*imh)';
data.x3d = cx;
data.z3d = cz;
data.foot = [footx footz];

feat = contactdata2features(data);

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

function [x, z] = computeGroundPosition(u, v, v0, yc, f)

z = yc*f./max((v0-v), 0.001);
x = u.*z./f;

author: Cloud-CV

repo: object-proposals


text:
%system model for wander mechanization (with dcm) under small angle
%assumption and no wander error. (Based on Phi-model definition)
%dt==0 -> return cont. time model
%Depending on input arguments this function may return (A,N), (Ad,N),
%(Ac,Qc),(Ad,Qd)
%this function is nothing but sys_wander_large with deleted wander states.

%Although this is for a wander-azimuth mechanization, this implementation
%is singular as it calls geoparam_v000.

%For a real non-singular sys-model implementation see sys_wander_psi_v000

%(For an example usage of this see: example_largeheading)

function [RetA RetB]=sys_wander_dcm(llh, vel_n, Cbn, wander, acc, dt, Aimu, Qimu, Cimu, Rimu)
% position   (1-3) (not it llh but in northing, easting, height)
% velocity   (4-6)
% attitude   (7-9)
% wander     (10-11) ([del_sa, del_ca])

%system disturbance coefs
N=zeros(9,6);
N(7:9,4:6)=-Cbn; %attitude
N(4:6,1:3)=Cbn; %velocity

%system matrix
A=zeros(9);

[Rn, Re, g, sL, cL, WIE_E]=geoparam_v000(llh);
Cgn=[wander(2) wander(1) 0;-wander(1) wander(2) 0; 0 0 1];
Rn_h=Rn+llh(3);
Re_h=Re+llh(3);
R=sqrt(Rn_h*Re_h);
acc_n=Cbn*acc;

mx_b=[wander(2)*(-WIE_E)*sL/Rn_h 0;wander(1)*(-WIE_E)*sL/Rn_h 0; (-WIE_E)*cL/Rn_h 0];
mx_c=[0 1/R 0;-1/R 0 0;0 0 0];
mx_d=[-vel_n(2)/R^2; vel_n(1)/R^2;0];

mx_e=Cgn*[0 1/(Re_h) 0; -1/(Rn_h) 0 0;0 0 0]*Cgn'*vel_n+2*Cgn*[WIE_E*cL; 0; -WIE_E*sL]; %2wie_n+wen_n
mx_f=Cgn*[0 1/(Re_h) 0; -1/(Rn_h) 0 0;0 0 0]*Cgn'*vel_n+Cgn*[WIE_E*cL; 0; -WIE_E*sL]; %wie_n+wen_n

%%position errors
A(1:2,4:5)=Cgn(1:2,1:2)';
A(3,6)=-1;

%%Velocity errors
A(4:6,1:2)=2*skew(vel_n)*mx_b;
A(4:6,3)=skew(vel_n)*mx_d+[0;0;-2*g/R];
A(4:6,4:6)=skew(vel_n)*mx_c-skew(mx_e);
A(4:6,7:9)=skew(acc_n);

%%Attitude errors
A(7:9,1:2)=mx_b;
A(7:9,3)=mx_d;
A(7:9,4:6)=mx_c;
A(7:9,7:9)=-skew(mx_f);

%%%%%%%%
%%Combine navigation models with imu models and return
if ~(isempty(Aimu) || isempty(Qimu) || isempty(Rimu) || isempty(Cimu))
    nst_imu=size(Aimu,1);
    nst_nav=size(A,1);
    Ac=[A N*Cimu;zeros(nst_imu,nst_nav) Aimu];
    Qc=[N*Rimu*N' zeros(nst_nav,nst_imu);zeros(nst_imu,nst_nav) Qimu];

    RetA=Ac;
    RetB=Qc;

    if (dt>0)
        %Discretize the model
        nst=nst_nav+nst_imu;
        mx_a = dt*[-Ac,Qc;zeros(nst),Ac'];
        mx_b = expm(mx_a);
        STM = mx_b(nst+1:2*nst,nst+1:2*nst)';
        Qd = STM*mx_b(1:nst,nst+1:2*nst);
        RetA=STM;
        RetB=Qd;
    end
else
    RetA=A;
    RetB=N;
    if (dt>0)
        RetA=expm(A*dt);
    end
end

repo: ekfmonoslam


text:
clear; clc;
%% desired signal
% n = (1:1000)';
%#codegen
% s = sin(0.075*pi*n);

%% noise signal
% v = 0.8*randn(1000,1); % Random noise part.
% ar = [1,1/2];          % Autoregression coefficients.
% v1 = filter(1,ar,v);   % Noise signal. Applies a 1-D digital filter.

%% primary input (noise corrupted signal)
% x = s + v1;
[x,fs1] = audioread('x14.wav');

%% reference input
% ma = [1, -0.8, 0.4 , -0.2];
% v2 = filter(ma,1,v);
[v2,fs1] = audioread('x15.wav');

%% adaptive filter
% initialization
L = 512;
% lms = dsp.LMSFilter(L,'Method','LMS');
nlms = dsp.LMSFilter(L,'Method','Normalized LMS','LeakageFactor',1,'AdaptInputPort',false,...
'WeightsResetInputPort',false,'WeightsOutput','Last');
% filter step size
% [mumaxlms,mumaxmselms]   = maxstep(lms,x);
[mumaxnlms,mumaxmsenlms] = maxstep(nlms,x); % maxstep function of dsp.LMSFilter
% lms.StepSize  = mumaxmselms/30;
nlms.StepSize = mumaxmsenlms/6;

%% filter with the adaptive filter
% [ylms,elms,wlms] = lms(v2,x);
[ynlms,enlms,wnlms] = nlms(v2,x);

%% compute the optimal solution (FIR Wiener filter)
% bw = firwiener(L-1,v2,x); % Optimal FIR Wiener filter
% yw = filter(bw,1,v2);   % Estimate of x using Wiener filter
% ew = x - yw;            % Estimate of actual sinusoid

%% plot
% primary input
% plot(n(900:end),x(900:end),'k:')
% xlabel('Time index (n)');
% ylabel('Amplitude');
plot([1:length(x)], x, 'm:')
xlabel('Time index (n)');
ylabel('Amplitude');

% denoised result
hold on;
% plot(n(900:end),[ew(900:end), elms(900:end),enlms(900:end)]);
% legend('Wiener filter denoised sinusoid',...
%     'LMS denoised sinusoid','NLMS denoised sinusoid');
% xlabel('Time index (n)');
% ylabel('Amplitude');
plot([1:length(x)], enlms);
xlabel('Time index (n)');
ylabel('Amplitude');
legend('primary input','denoised result');
hold off;

%% reset nlms
reset(nlms);

%% convergence investigation through learning curve
% M = 10; % Decimation factor
% msenlms = msesim(nlms,v2,x,M);
% figure;
% plot(1:M:x(end),msenlms)
% legend('LMS learning curve','NLMS learning curve')
% xlabel('Time index (n)');
% ylabel('MSE');

%% theorectical learning curves
% reset(nlms);
% figure;
% [mmselms,emselms,meanwlms,pmselms] = msepred(nlms,v2,x,M);
% plot(1:M:x(end),[mmselms*ones(500,1),emselms*ones(500,1),...
%         pmselms,mselms])
% legend('MMSE','EMSE','predicted LMS learning curve',...
%     'LMS learning curve')
% xlabel('Time index (n)');
% ylabel('MSE');

author: CharlesThaCat

repo: acoustic-interference-cancellation


text:
function [S, mets, badRxns, reversible]=constructS(equations,mets,rxns)
% constructS
%   Constructs a stoichiometric matrix from a cell array of equations
%
%   equations   cell array of equations on the form 'A + 2 B <=> 3 C',
%               where <=> indicates reversible and => irreversible reactions
%   mets        cell array of metabolites. All metabolites in the equations
%               must be present in the list (opt, default generated from
%               the equations)
%   rxns        cell array of reaction ids. This is only used for printing
%               reaction ids instead of equations in warnings/errors (opt,
%               default [])
%
%   S           the resulting stoichiometric matrix mets cell array with
%               metabolites that corresponds to the order in the S matrix
%   badRxns     boolean vector with the reactions that have one or more
%               metabolites as both substrate and product. An example would
%               be the phosphotransferase ATP + ADP <=> ADP + ATP. In the
%               stoichiometric matrix this equals to an empty reaction
%               which can be problematic
%   reversible  boolean vector with true if the equation is reversible
%
%   Usage: [S, mets, badRxns, reversible]=constructS(equations,mets)

equations=convertCharArray(equations);
switch nargin
    case 2
        mets=convertCharArray(mets);
    case 3
        rxns=convertCharArray(rxns);
end

badRxns=false(numel(equations),1);

%Check that no equations are too short to have reversibility data
I=cellfun(@numel,equations);
I=find(I<4,1);
if any(I)
    if isempty(rxns)
        EM=['The following equation does not have reversibility data: ' equations{I} ];
        dispEM(EM);
    else
        EM=['The reaction ' rxns{I} ' does not have reversibility data'];
        dispEM(EM);
    end
end

%Makes life a little easier
equations=strtrim(equations);
equations=fixEquations(equations);

if nargin<2
    mets=parseRxnEqu(equations);
end
if nargin<3
    rxns=[];
end

%Get which reactions are reversible
reversible=cellfun(@any,strfind(equations,' <=> '));

%Make them all reversible. This is not all that neat, but nevermind
equations=strrep(equations,' => ',' <=> ');

%Replace the the plus signs with some weird character that will be used for
%parsing
equations=strrep(equations,' + ', '€');

%Generate the stoichiometric matrix
S=zeros(numel(mets),numel(equations));

%Keep track of coefficients to be added to S-matrix
metsToS = cell(100000,1);
rxnsToS = zeros(100000,1);
coefToS = zeros(100000,1);
newEntry = 1;
%Loop through the equations and add the info to the S matrix
for i=1:numel(equations)
    %Start by finding the position of the (=> or <=>)
    arrowIndex=strfind(equations{i},' <=> ');

    if numel(arrowIndex)~=1
        if isempty(rxns)
            EM=['The following equation does not have reversibility data: ' equations{i} ];
            dispEM(EM);
        else
            EM=['The reaction ' rxns{i} ' does not have reversibility data'];
            dispEM(EM);
        end
    end

    reactants=regexp(equations{i}(1:arrowIndex-1),'€','split');
    products=regexp(equations{i}(arrowIndex+5:end),'€','split');

    %If the splitting character is at the end (if exchange rxns), then an
    %empty string will exist together with the real ones. Remove it
    reactants(cellfun(@isempty,reactants))=[];
    products(cellfun(@isempty,products))=[];

    %A vector where an element is -1 is the corresponding metabolite is a
    %reactant and 1 if it's a product
    multiplyWith=[ones(numel(reactants),1)*-1; ones(numel(products),1)];

    metabolites=strtrim([reactants products]);

    %Now loop through the reactants and see if the metabolite has a
    %coefficient (it will look as 'number name')
    for j=1:numel(metabolites)
        space=strfind(metabolites{j},' ');

        if isempty(space)
            %No coefficient
            coeff=1;
            name=metabolites{j};
        else
            coeff=str2double(metabolites{j}(1:space(1)));

            %If it was not a coefficiant
            if isnan(coeff)
                coeff=1;
                name=strtrim(metabolites{j});
            else
                name=strtrim(metabolites{j}(space+1:end));
            end
        end

        %Find the name in the mets list [a b]=ismember(name,mets);
        metsToS{newEntry}=name;
        rxnsToS(newEntry)=i;
        coefToS(newEntry)=coeff*multiplyWith(j);
        newEntry=newEntry+1;
    end
end
%Remove unused fields
metsToS(newEntry:end)=[];
rxnsToS(newEntry:end)=[];
coefToS(newEntry:end)=[];

%Match to mets array
[metsPresent,metsLoc]=ismember(metsToS,mets);

%Find badRxns
[~,I]=unique([rxnsToS,metsLoc],'rows','stable');
x=1:length(rxnsToS);
x(I)=[];
x=unique(rxnsToS(x));
badRxns(x)=true;

if any(~metsPresent)
    if isempty(rxns)
        error(['Could not find the following metabolites in the metabolite list: ',...
        strjoin(unique(metsToS(~metsPresent)),', ')],'')
    else
        missingMet = find(~metsPresent);
        missingMet = char(strcat(metsToS(missingMet),' (reaction:',rxns(rxnsToS(missingMet)),')\n'));
        error(['Could not find the following metabolites (reaction indicated) in the metabolite list: \n' ...
            missingMet '%s'],'');
    end
end
linearIndices=sub2ind(size(S),metsLoc,rxnsToS);
S(linearIndices)=coefToS;
S=sparse(S);
end

function equ=fixEquations(equ)
%If the equation starts with "=>" or "<=>" then add a space again. This is
%an alternative way to represent uptake reactions. The opposite way for
%producing reactions
equ=equ(:);
for i=1:numel(equ)
    if strcmp(equ{i}(1:2),'=>') || strcmp(equ{i}(1:3),'<=>')
        equ{i}=[' ' equ{i}];
    else
        if strcmp(equ{i}(end-1:end),'=>') || strcmp(equ{i}(end-2:end),'<=>')
            equ{i}=[equ{i} ' '];
        end
    end
end
end

author: SysBioChalmers

repo: RAVEN


text:
function c = tapas_hgf_binary_config
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Contains the configuration for the Hierarchical Gaussian Filter (HGF)
% for binary inputs in the absence of perceptual uncertainty.
%
% The HGF is the model introduced in
%
% Mathys C, Daunizeau J, Friston, KJ, and Stephan KE. (2011). A Bayesian foundation
% for individual learning under uncertainty. Frontiers in Human Neuroscience, 5:39.
%
% The binary HGF model has since been augmented with a positive factor kappa1 which
% scales the second level with respect to the first, i.e., the relation between the
% first and second level is
%
% p(x1=1|x2) = s(kappa1*x2), where s(.) is the logistic sigmoid.
%
% By default, kappa1 is fixed to 1, leading exactly to the model introduced in
% Mathys et al. (2011).
%
% This file refers to BINARY inputs (Eqs 1-3 in Mathys et al., (2011));
% for continuous inputs, refer to tapas_hgf_config.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% The HGF configuration consists of the priors of parameters and initial values. All priors are
% Gaussian in the space where the quantity they refer to is estimated. They are specified by their
% sufficient statistics: mean and variance (NOT standard deviation).
%
% Quantities are estimated in their native space if they are unbounded (e.g., the omegas). They are
% estimated in log-space if they have a natural lower bound at zero (e.g., the sigmas).
%
% Parameters can be fixed (i.e., set to a fixed value) by setting the variance of their prior to
% zero. Aside from being useful for model comparison, the need for this arises whenever the scale
% and origin at the j-th level are arbitrary. This is the case if the observation model does not
% contain the representations mu_j and sigma_j. A choice of scale and origin is then implied by
% fixing the initial value mu_j_0 of mu_j and either kappa_j-1 or omega_j-1.
%
% Fitted trajectories can be plotted by using the command
%
% >> tapas_hgf_binary_plotTraj(est)
%
% where est is the stucture returned by tapas_fitModel. This structure contains the estimated
% perceptual parameters in est.p_prc and the estimated trajectories of the agent's
% representations (cf. Mathys et al., 2011). Their meanings are:
%
%         est.p_prc.mu_0       row vector of initial values of mu (in ascending order of levels)
%         est.p_prc.sa_0       row vector of initial values of sigma (in ascending order of levels)
%         est.p_prc.rho        row vector of rhos (representing drift; in ascending order of levels)
%         est.p_prc.ka         row vector of kappas (in ascending order of levels)
%         est.p_prc.om         row vector of omegas (in ascending order of levels)
%
% Note that the first entry in all of the row vectors will be NaN because, at the first level,
% these parameters are either determined by the second level (mu_0 and sa_0) or undefined (rho,
% kappa, and omega).
%
%         est.traj.mu          mu (rows: trials, columns: levels)
%         est.traj.sa          sigma (rows: trials, columns: levels)
%         est.traj.muhat       prediction of mu (rows: trials, columns: levels)
%         est.traj.sahat       precisions of predictions (rows: trials, columns: levels)
%         est.traj.v           inferred variance of random walk (rows: trials, columns: levels)
%         est.traj.w           weighting factors (rows: trials, columns: levels)
%         est.traj.da          volatility prediction errors  (rows: trials, columns: levels)
%         est.traj.ud          updates with respect to prediction  (rows: trials, columns: levels)
%         est.traj.psi         precision weights on prediction errors  (rows: trials, columns: levels)
%         est.traj.epsi        precision-weighted prediction errors  (rows: trials, columns: levels)
%         est.traj.wt          full weights on prediction errors (at the first level,
%                                  this is the learning rate) (rows: trials, columns: levels)
%
% Note that in the absence of sensory uncertainty (which is the assumption here), the first
% column of mu, corresponding to the first level, will be equal to the inputs. Likewise, the
% first column of sa will be 0 always.
%
% Tips:
% - When analyzing a new dataset, take your inputs u and use
%
%   >> est = tapas_fitModel([], u, 'tapas_hgf_binary_config', 'tapas_bayes_optimal_binary_config');
%
%   to determine the Bayes optimal perceptual parameters (given your current priors as defined in
%   this file here, so choose them wide and loose to let the inputs influence the result). You can
%   then use the optimal parameters as your new prior means for the perceptual parameters.
%
% - If you get an error saying that the prior means are in a region where model assumptions are
%   violated, lower the prior means of the omegas, starting with the highest level and proceeding
%   downwards.
%
% - Alternatives are lowering the prior means of the kappas, if they are not fixed, or adjusting
%   the values of the kappas or omegas, if any of them are fixed.
%
% - If the log-model evidence cannot be calculated because the Hessian poses problems, look at
%   est.optim.H and fix the parameters that lead to NaNs.
%
% - Your guide to all these adjustments is the log-model evidence (LME). Whenever the LME increases
%   by at least 3 across datasets, the adjustment was a good idea and can be justified by just this:
%   the LME increased, so you had a better model.
%
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2012-2017 Christoph Mathys, TNU, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

% Config structure
c = struct;

% Model name
c.model = 'hgf_binary';

% Number of levels (minimum: 3)
c.n_levels = 3;

% Input intervals
% If input intervals are irregular, the last column of the input
% matrix u has to contain the interval between inputs k-1 and k
% in the k-th row, and this flag has to be set to true
c.irregular_intervals = false;

% Sufficient statistics of Gaussian parameter priors

% Initial mus and sigmas
% Format: row vectors of length n_levels
% For all but the first two levels, this is usually best
% kept fixed to 1 (determines origin on x_i-scale). The
% first level is NaN because it is determined by the second,
% and the second implies neutrality between outcomes when it
% is centered at 0.
c.mu_0mu = [NaN, 0, 1];
c.mu_0sa = [NaN, 0, 0];

c.logsa_0mu = [NaN,   log(0.1), log(1)];
c.logsa_0sa = [NaN,          0,      0];

% Rhos
% Format: row vector of length n_levels.
% Undefined (therefore NaN) at the first level.
% Fix this to zero to turn off drift.
c.rhomu = [NaN, 0, 0];
c.rhosa = [NaN, 0, 0];

% Kappas
% Format: row vector of length n_levels-1.
% Fixing log(kappa1) to log(1) leads to the original HGF model.
% Higher log(kappas) should be fixed (preferably to log(1)) if the
% observation model does not use mu_i+1 (kappa then determines the
% scaling of x_i+1).
c.logkamu = [log(1), log(1)];
c.logkasa = [     0,      0];

% Omegas
% Format: row vector of length n_levels.
% Undefined (therefore NaN) at the first level.
c.ommu = [NaN,  -3,  -6];
c.omsa = [NaN, 4^2, 4^2];

% Gather prior settings in vectors
c.priormus = [
    c.mu_0mu,...
    c.logsa_0mu,...
    c.rhomu,...
    c.logkamu,...
    c.ommu,...
         ];

c.priorsas = [
    c.mu_0sa,...
    c.logsa_0sa,...
    c.rhosa,...
    c.logkasa,...
    c.omsa,...
         ];

% Check whether we have the right number of priors
expectedLength = 3*c.n_levels+2*(c.n_levels-1)+1;
if length([c.priormus, c.priorsas]) ~= 2*expectedLength;
    error('tapas:hgf:PriorDefNotMatchingLevels', 'Prior definition does not match number of levels.')
end

% Model function handle
c.prc_fun = @tapas_hgf_binary;

% Handle to function that transforms perceptual parameters to their native space
% from the space they are estimated in
c.transp_prc_fun = @tapas_hgf_binary_transp;

return;

author: translationalneuromodeling

repo: tapas


text:
function sFile = in_fopen_manscan(DataFile)
% IN_FOPEN_MANSCAN: Open a MANSCAN file (continuous recordings)
%
% USAGE:  sFile = in_fopen_manscan(DataFile)

% @=============================================================================
% This software is part of the Brainstorm software:
% http://neuroimage.usc.edu/brainstorm
%
% Copyright (c)2000-2013 Brainstorm by the University of Southern California
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPL
% 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, 2012

%% ===== READ HEADER =====
% Get text file (.mbi)
MbiFile = strrep(DataFile, '.mb2', '.mbi');
% If doesn't exist: error
if ~exist(MbiFile, 'file')
    ft_error('Cannot open file: missing text file .mbi');
end
% Initialize header
iEpoch = 1;
hdr.epoch(iEpoch).Comment = {};
hdr.epoch(iEpoch).Channel = [];
hdr.Events = [];
curBlock = '';
% Read file line by line
fid = fopen_or_error(MbiFile,'r');
while(1)
    % Reached the end of the file: exit the loop
    if feof(fid)
        break;
    end;
    % Read one line
    buf = fgetl(fid);
    if isempty(buf)
        curBlock = '';
        continue;
    end
    % Split line based on space characters
    splitBuf = str_split(buf, [0 9 32]);
    if isempty(splitBuf)
        continue;
    end
    % Current block
    if ~isempty(curBlock)
        switch curBlock
            case 'event'
                % Ignore MANSCAN auto events
                if any(strcmpi(evtProp, 'Channel'))
                    continue;
                end
                % New event
                iEvt = length(hdr.Events) + 1;
                % Read event name and epoch
                hdr.Events(iEvt).Name   = splitBuf{1};
                hdr.Events(iEvt).iEpoch = iEpoch;
                % Read each property
                for iProp = 2:length(evtProp)
                    if strcmpi(evtProp{iProp}, 'Channel')
                        hdr.Events(iEvt).Name = [hdr.Events(iEvt).Name '_' splitBuf{iProp}];
                    elseif strcmpi(evtType{iProp}, 'String')
                        hdr.Events(iEvt).(evtProp{iProp}) = splitBuf{iProp};
                    else
                        hdr.Events(iEvt).(evtProp{iProp}) = str2num(splitBuf{iProp});
                    end
                end

            case 'channel'
                % New channel
                iChan = length(hdr.epoch(iEpoch).Channel) + 1;
                % Read channel name
                hdr.epoch(iEpoch).Channel(iChan).Name = splitBuf{1};
                % Read each property
                for iProp = 2:length(chanProp)
                    if strcmpi(chanProp{iProp}, 'UnitDisplayToFile')
                        propValue = str2num(splitBuf{iProp});
                    else
                        propValue = splitBuf{iProp};
                    end
                    hdr.epoch(iEpoch).Channel(iChan).(chanProp{iProp}) = propValue;
                end
        end
        % Next line
        continue;
    end
    % Detect keyword
    switch lower(splitBuf{1})
        case 'datatypeid'
            iEpoch = iEpoch + 1;
            hdr.epoch(iEpoch).Comment = {};
            hdr.epoch(iEpoch).Channel = [];
        case 'comment'
            hdr.epoch(iEpoch).Comment{end+1} = strtrim(strrep(buf, splitBuf{1}, ''));
        case 'event'
            curBlock = 'event';
            evtProp = splitBuf;
            evtType = str_split(fgetl(fid), [0 9 32]);
        case 'channel'
            curBlock = 'channel';
            chanProp = splitBuf;
        case 'offsetdisplayexperiment'
            hdr.epoch(iEpoch).OffsetDisplay = str2num(splitBuf{2});
        case 'experimenttime'
            hdr.epoch(iEpoch).ExperimentTime = strtrim(strrep(buf, splitBuf{1}, ''));
        case 'worddatafile'
            % Read channel order
            hdr.epoch(iEpoch).ChannelOrder = str_split(fgetl(fid), [0 9 32]);
            % Read record info
            datainfo = str_split(fgetl(fid), [0 9 32]);
            hdr.epoch(iEpoch).StartData1 = str2num(datainfo{1});
            hdr.epoch(iEpoch).StartData2 = str2num(datainfo{2});
            hdr.epoch(iEpoch).Sweeps     = str2num(datainfo{3});
            hdr.epoch(iEpoch).BinFile    = datainfo{4};
    end
end

%% ===== CREATE BRAINSTORM SFILE STRUCTURE =====
% Initialize returned file structure
sFile = []; %db_template('sfile');
% Add information read from header
sFile.byteorder  = 'l';
sFile.filename   = DataFile;
sFile.format = 'EEG-MANSCAN';
sFile.device = 'MANSCAN';
sFile.channelmat = [];
% Comment: short filename
[tmp__, sFile.comment, tmp__] = fileparts(DataFile);
% Consider that the sampling rate of the file is the sampling rate of the first signal
sFile.prop.sfreq = 256;
sFile.prop.nAvg  = 1;
% No info on bad channels
sFile.channelflag = ones(length(hdr.epoch(iEpoch).ChannelOrder), 1);

%% ===== EPOCHS =====
if (length(hdr.epoch) <= 1)
    sFile.prop.samples = [0, hdr.epoch(1).Sweeps-1];
    sFile.prop.times   = sFile.prop.samples ./ sFile.prop.sfreq;
else
    % Build epochs structure
    for iEpoch = 1:length(hdr.epoch)
        sFile.epochs(iEpoch).label   = sprintf('Epoch #%02d', iEpoch);
        sFile.epochs(iEpoch).samples = [0, hdr.epoch(iEpoch).Sweeps-1];
        sFile.epochs(iEpoch).times   = sFile.epochs(iEpoch).samples ./ sFile.prop.sfreq;
        sFile.epochs(iEpoch).nAvg    = 1;
        sFile.epochs(iEpoch).select  = 1;
        sFile.epochs(iEpoch).bad         = 0;
        sFile.epochs(iEpoch).channelflag = [];

        % Check if all the epochs have the same channel list
        if (iEpoch > 1) && ~isequal(hdr.epoch(iEpoch).ChannelOrder, hdr.epoch(1).ChannelOrder)
            ft_error('Channel list must remain constant across epochs.');
        end
    end
    % Extract global min/max for time and samples indices
    sFile.prop.samples = [min([sFile.epochs.samples]), max([sFile.epochs.samples])];
    sFile.prop.times   = [min([sFile.epochs.times]),   max([sFile.epochs.times])];
end

%% ===== CREATE EMPTY CHANNEL FILE =====
nChannels = length(hdr.epoch(1).ChannelOrder);
ChannelMat.Comment = [sFile.device ' channels'];
%ChannelMat.Channel = repmat(db_template('channeldesc'), [1, nChannels]);
hdr.Gains = [];
% For each channel
for iChan = 1:nChannels
    % Find channel in description list
    chName = hdr.epoch(1).ChannelOrder{iChan};
    iDesc = find(strcmpi({hdr.epoch(1).Channel.Name}, chName));
    sDesc = hdr.epoch(1).Channel(iDesc);
    % Type
    if isfield(sDesc, 'IsEEG') && ~isempty(sDesc.IsEEG)
        if strcmpi(sDesc.IsEEG, 'Yes')
            ChannelMat.Channel(iChan).Type = 'EEG';
        else
            ChannelMat.Channel(iChan).Type = 'Misc';
        end
    elseif isfield(sDesc, 'IsEEG')
        ChannelMat.Channel(iChan).Type = 'EEG REF';
    else
        ChannelMat.Channel(iChan).Type = 'EEG';
    end
    % Create structure
    ChannelMat.Channel(iChan).Name    = chName;
    ChannelMat.Channel(iChan).Loc     = [0; 0; 0];
    ChannelMat.Channel(iChan).Orient  = [];
    ChannelMat.Channel(iChan).Weight  = 1;
    ChannelMat.Channel(iChan).Comment = '';
    % Check that this channel has a gain
    if isempty(hdr.epoch(1).Channel(iDesc).UnitDisplayToFile)
        hdr.Gains(iChan,1) = 1;
    else
        hdr.Gains(iChan,1) = hdr.epoch(1).Channel(iDesc).UnitDisplayToFile;
    end
end
% Return channel structure
sFile.channelmat = ChannelMat;

%% ===== FORMAT EVENTS =====
if ~isempty(hdr.Events)
    % Get all the epochs names
    uniqueEvt = unique({hdr.Events.Name});
    % Create one category for each event
    for iEvt = 1:length(uniqueEvt)
        % Get all the occurrences
        iOcc = find(strcmpi({hdr.Events.Name}, uniqueEvt{iEvt}));
        % Get the samples for all the occurrences
        sample   = [hdr.Events(iOcc).BeginSample];
        duration = [hdr.Events(iOcc).DurationSamples];
        % Extended event = duration is not 1 for all the markers
        if ~all(duration == 1)
            sample = [sample; sample + duration];
        end
        % Create event structure
        sFile.events(iEvt).label   = uniqueEvt{iEvt};
        sFile.events(iEvt).samples = sample;
        sFile.events(iEvt).times   = sample ./ sFile.prop.sfreq;
        sFile.events(iEvt).epochs  = [hdr.Events(iOcc).iEpoch];
        sFile.events(iEvt).select  = 1;
    end
end
% Save file header
sFile.header = hdr;

function splStr = str_split( str, delimiters, isCollapse )
% STR_SPLIT: Split string.
%
% USAGE:  str_split( str, delimiters, isCollapse=1 ) : delimiters in an array of char delimiters
%         str_split( str )             : default are file delimiters ('\' and '/')
%
% INPUT:
%    - str        : String to split
%    - delimiters : String that contains all the characters used to split, default = '/\'
%    - isCollapse : If 1, remove all the empty entries
%
% OUTPUT:
%    - splStr : cell-array of blocks found between separators

% @=============================================================================
% This software is part of the Brainstorm software:
% http://neuroimage.usc.edu/brainstorm
%
% Copyright (c)2000-2013 Brainstorm by the University of Southern California
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPL
% 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, 2008

% Default delimiters: file delimiters ('\', '/')
if (nargin < 3) || isempty(isCollapse)
    isCollapse = 1;
end
if (nargin < 2) || isempty(delimiters)
    delimiters = '/\';
end
% Empty input
if isempty(str)
    splStr = {};
    return
end

% Find all delimiters in string
iDelim = [];
for i=1:length(delimiters)
    iDelim = [iDelim strfind(str, delimiters(i))];
end
iDelim = unique(iDelim);

% If no delimiter: return the whole string
if isempty(iDelim)
    splStr = {str};
    return
end

% Allocates the split array
splStr = cell(1, length(iDelim)+1);

% First part (before first delimiter)
if (iDelim(1) ~= 1)
    iSplitStr = 1;
    splStr{iSplitStr} = str(1:iDelim(1)-1);
else
    iSplitStr = 0;
end

% Loop over all other delimiters
for i = 2:length(iDelim)
    if (isCollapse  && (iDelim(i) - iDelim(i-1) > 1)) || ...
       (~isCollapse && (iDelim(i) - iDelim(i-1) >= 1))
        iSplitStr = iSplitStr + 1;
        splStr{iSplitStr} = str(iDelim(i-1)+1:iDelim(i)-1);
    end
end

% Last part (after last delimiter)
if (iDelim(end) ~= length(str))
    iSplitStr = iSplitStr + 1;
    splStr{iSplitStr} = str(iDelim(end)+1:end);
end

% Remove all the unused entries
if (iSplitStr < length(splStr))
    splStr(iSplitStr+1:end) = [];
end

author: spm


text:
function [d,d1] = util_norm_ratio(f,g,displ)
if nargin < 3 || isempty(displ)
    displ = false;
end
 d=norm(f - g)/norm(g + f);
d1=norm(f - g,1)/norm(g + f,1); %%
if displ
    fprintf(1,' Norm1 difference: %d\n Norm2 difference: %d\n',d1,d);
end

author: SheffieldML

repo: vargplvm


text:
function [varargout] = likGamma(link, hyp, y, mu, s2, inf, i)

% likGamma - Gamma likelihood function for strictly positive data y. The
% expression for the likelihood is
%   likGamma(f) = al^al*y^(al-1)/gamma(al) * exp(-y*al/mu) / mu^al with
% mean=mu and variance=mu^2/al where mu = g(f) is the Gamma intensity, f is a
% Gaussian process, y is the strictly positive data. Hence, we have -- with
% log(Zy) = log(gamma(al)) - al*log(al) + (1-al)*log(y)
%   llik(f) = log(likGamma(f)) = -al*( log(g(f)) + y/g(f) ) - log(Zy).
% The larger one chooses al, the stronger the likelihood resembles a Gaussian
% since skewness = 2/sqrt(al) and kurtosis = 6/al.
%
% We provide two inverse link functions 'exp' and 'logistic':
%   g(f) = exp(f) and g(f) = log(1+exp(f))).
% The link functions are located at util/glm_invlink_*.m.
%
% Note that for neither link function the likelihood lik(f) is log concave.
%
% The hyperparameters are:
%
% hyp = [  log(al)  ]
%
% Several modes are provided, for computing likelihoods, derivatives and moments
% respectively, see likFunctions.m for the details. In general, care is taken
% to avoid numerical issues when the arguments are extreme.
%
% See also LIKFUNCTIONS.M.
%
% Copyright (c) by Hannes Nickisch, 2013-10-16.

if nargin<4, varargout = {'1'}; return; end   % report number of hyperparameters

al = exp(hyp);

if nargin<6                              % prediction mode if inf is not present
  if numel(y)==0,  y = zeros(size(mu)); end
  s2zero = 1; if nargin>4, if norm(s2)>0, s2zero = 0; end, end         % s2==0 ?
  if s2zero                                                    % log probability
    lg = g(mu,link);
    lZy = gammaln(al) - al*log(al) + (1-al)*log(y);     % normalisation constant
    lp = -al*(lg+y./exp(lg)) - lZy;
  else
    lp = likGamma(link, hyp, y, mu, s2, 'infEP');
  end
  ymu = {}; ys2 = {};
  if nargout>1                                 % compute y moments by quadrature
    n = max([length(y),length(mu),length(s2)]); on = ones(n,1);
    N = 20; [t,w] = gauher(N); oN = ones(1,N); lw = ones(n,1)*log(w');
    mu = mu(:).*on; sig = sqrt(s2(:)).*on;                        % vectors only
    lg = g(sig*t'+mu*oN,link);
    ymu = exp(logsumexp2(lg+lw));     % first moment using Gaussian-Hermite quad
    if nargout>2
      elg = exp(lg);
      yv = elg.^2/al;                  % second y moment from Gamma distribution
      ys2 = (yv+(elg-ymu*oN).^2)*w;
    end
  end
  varargout = {lp,ymu,ys2};
else
  switch inf
  case 'infLaplace'
    [lg,dlg,d2lg,d3lg] = g(mu,link); elg = exp(lg);
    if nargin<7                                             % no derivative mode
      lZy = gammaln(al) - al*log(al) + (1-al)*log(y);   % normalisation constant
      lp = -al*(lg+y./elg) - lZy;
      dlp = {}; d2lp = {}; d3lp = {};                         % return arguments
      if nargout>1
        dlp = -al*dlg.*(1-y./elg);           % dlp, derivative of log likelihood
        if nargout>2                    % d2lp, 2nd derivative of log likelihood
          d2lp = -al*d2lg.*(1-y./elg) - al*dlg.*dlg.*y./elg;
          if nargout>3                  % d3lp, 3rd derivative of log likelihood
            d3lp = -al*d3lg.*(1-y./elg) + al*dlg.*(dlg.*dlg-3*d2lg).*y./elg;
          end
        end
      end
      varargout = {lp,dlp,d2lp,d3lp};
    else                                                       % derivative mode
      dlZy = al*psi(0,al) - al*(log(al) + 1 + log(y));
      lp_dhyp = -al*(lg+y./elg) - dlZy; % derivative of log likelihood w.r.t. al
      dlp_dhyp = -al*dlg.*(1-y./elg);                         % first derivative
      d2lp_dhyp = -al*d2lg.*(1-y./elg) - al*dlg.*dlg.*y./elg;  % and also second
      varargout = {lp_dhyp,dlp_dhyp,d2lp_dhyp};
    end

  case 'infEP'
    if nargin<7                                             % no derivative mode
      % Since we are not aware of an analytical expression of the integral,
      % we use quadrature.
      varargout = cell(1,nargout);
      [varargout{:}] = lik_epquad({@likGamma,link},hyp,y,mu,s2);
    else                                                       % derivative mode
      varargout = {[]};                                     % deriv. wrt hyp.lik
    end

  case 'infVB'
    error('infVB not supported')
  end
end

% compute the log intensity using the inverse link function
function varargout = g(f,link)
  varargout = cell(nargout, 1);  % allocate the right number of output arguments
  if strcmp(link,'exp')
    [varargout{:}] = glm_invlink_exp(f);
  else
    [varargout{:}] = glm_invlink_logistic(f);
  end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function [bbci_cls, data_adapt]= bbci_adaptation_pcovmean(bbci_cls, data_adapt, marker, varargin)
%BBCI_ADAPTATIOIN_PCOVMEAN - Supervised adapt class means and pooled covariance matrix of an LDA classifier
%
%Technique by Carmen Vidaurre, see
%  Viduarre C et al, J Neural Eng, 2011
%  http://dx.doi.org/10.1088/1741-2560/8/2/025009
%
%Synopsis:
%  [BBCI_CLS, DATA_ADAPT]= ...
%      bbci_adaptation_pcovmean(BBCI_CLS, DATA_ADAPT, 'init', PARAMs, ...)
%  [BBCI_CLS, DATA_ADAPT]= ...
%      bbci_adaptation_pcovmean(BBCI_CLS, DATA_ADAPT, MARKER, FEATURE)
%
%This function is called internally by bbci_apply_adaptation.
%
%The selectable parameters are
%  'UC_mean' - [DOUBLE, default 0.05] Update coefficient for supervised
%              adaptation of the class means.
%  'UC_pcov' - [DOUBLE, default 0.03] Update coefficient for unsupervised
%              adaptation of the pooled covariance matrix (inverse of the
%              extended ..., to be exact).
%  'ival' - [1x2 DOUBLE] Time interval relative to the start marker, during
%              which features are averaged and used for adaptation;
%              no default - needs to be set.
%              Note: When choosing the start of the adaptation interval,
%              you have to take into account bbci.feature.ival (meaning
%              that feature are based on retrospective signals).
%  'mrk_start' - {1xnClasses CELL} defines the markers that trigger the
%              adaptation for each class; no default - needs to be set.
%  'mrk_end' - [1 nMarkers INT, default []] specifies markers that indicate
%              the end of a trial. Gathering feature information for the
%              adaptation of one trial ends, when one of those end-markers
%              is received AND the end of the adaptation interval 'ival'
%              (see above) is reached. If 'mrk_end' is empty, the first part
%              of the condition is ignored, i.e., only 'ival' is taken into
%              account.
%  'scaling' - [BOOL, default true] Defines whether the updated classifier
%              should be (re-) scaled.
%  'log_mean_limit' - [INT, default 100].

% 01-2012 Benjamin Blankertz
% 04-2011 Claudia Sannelli

if ischar(marker) && strcmp(marker, 'init'),
  props= {'UC_mean'          0.05    '!DOUBLE[1]'
          'UC_pcov'          0.03    '!DOUBLE[1]'
          'ival'             []      '!DOUBLE[2]'
          'scaling'          1       '!BOOL'
          'mrk_start'        {}      'CELL'
          'mrk_end'          []      'INT'
          'log_mean_limit'   100     '!INT'
          'tag'              ''      'CHAR'
         };
  opt= opt_proplistToStruct(varargin{:});
  opt= opt_setDefaults(opt, props, 2);
  data_adapt.opt= opt;
  data_adapt.feature= zeros(size(bbci_cls.C.w));
  data_adapt.trial_start= NaN;
  data_adapt.lastcheck= -inf;
  bbci_log_write(data_adapt.log.fid, ...
                 '# %s started with bias=%g and options %s.', ...
                 data_adapt.opt.tag, bbci_cls.C.b, util_toString(opt));
  return;
else
  feature= varargin{1};
end

time= marker.current_time;
check_ival= [data_adapt.lastcheck time];
events= bbci_apply_queryMarker(marker, check_ival);
data_adapt.lastcheck= time;

if ~isempty(events) && isnan(data_adapt.trial_start),
  data_adapt.clidx= 0;
  for k= 1:length(data_adapt.opt.mrk_start),
    marker_idx= find(ismember([events.desc], data_adapt.opt.mrk_start{k},'legacy'));
    if ~isempty(marker_idx),
      data_adapt.clidx= k;
      break;
    end
  end
  if data_adapt.clidx>0,
    data_adapt.trial_start= events(marker_idx(1)).time;
    data_adapt.end_marker_received= isempty(data_adapt.opt.mrk_end);
    data_adapt.counter= 0;
    bbci_log_write(data_adapt.log.fid, ...
        ['# %s at ' data_adapt.log.time_fmt ...
         ' trial started with marker %d -> class %d.'], ...
        data_adapt.opt.tag, data_adapt.trial_start/1000, ...
        events(marker_idx(1)).desc, data_adapt.clidx);
  end
end

if isnan(data_adapt.trial_start),
  return;
end

% within the adaptation interval, add up the features and count them
adapt_ival= data_adapt.trial_start + data_adapt.opt.ival;
if time >= adapt_ival(1) && time <= adapt_ival(2),
  data_adapt.feature= data_adapt.feature + feature.x(:);
  data_adapt.counter= data_adapt.counter + 1;
end

if ~isempty(events) && any(ismember([events.desc], data_adapt.opt.mrk_end,'legacy'));
  data_adapt.end_marker_received= 1;
end
if data_adapt.end_marker_received && time >= adapt_ival(2),
  if data_adapt.counter==0,
    return;
  end
  % handy definitions
  UCm= data_adapt.opt.UC_mean;
  UCp= data_adapt.opt.UC_pcov;
  feat= data_adapt.feature/data_adapt.counter;
  % adapt class means (supervised)
  bbci_cls.C.mean(:,data_adapt.clidx)= UCm * feat + ...
      (1 - UCm) * bbci_cls.C.mean(:,data_adapt.clidx);
  % adapt (inverse of extended) pooled covariance matrix (unsupervised)
  extfeat= [1; feat];
  v= bbci_cls.C.extinvcov*extfeat;
  bbci_cls.C.extinvcov= (1/(1 - UCp)) * ...
    (bbci_cls.C.extinvcov - UCp/(1 - UCp + UCp * extfeat' * v) * v*v');
  bbci_cls.C.extinvcov= 0.5 * (bbci_cls.C.extinvcov + bbci_cls.C.extinvcov');
  % recalculate classifier
  diff_mean= diff(bbci_cls.C.mean, 1, 2);
  bbci_cls.C.w= bbci_cls.C.extinvcov(2:end, 2:end) * diff_mean;

  % rescale projection vector
  if data_adapt.opt.scaling,
    bbci_cls.C.w= bbci_cls.C.w/(bbci_cls.C.w' * diff_mean)*2;
  end
  bbci_cls.C.b= -bbci_cls.C.w' * mean(bbci_cls.C.mean, 2);
  if numel(bbci_cls.C.mean) < data_adapt.opt.log_mean_limit,
    bbci_log_write(data_adapt.log.fid, ...
                   ['# %s at ' data_adapt.log.time_fmt ...
                   ' bias adapted to %g, mean adapted to %s' ...
                   ' with %d features.'], ...
                   data_adapt.opt.tag, marker.current_time/1000, ...
                   bbci_cls.C.b, util_toString(bbci_cls.C.mean), ...
                   data_adapt.counter);
  else
    bbci_log_write(data_adapt.log.fid, ...
                   ['# %s at ' data_adapt.log.time_fmt ...
                   ' bias adapted to %g with %d features.'], ...
                   data_adapt.opt.tag, marker.current_time/1000, ...
                   bbci_cls.C.b, data_adapt.counter);
  end
  data_adapt.feature(:)= 0;
  data_adapt.trial_start= NaN;
end

author: bbci


text:
function lik = ml_gauss_lik(y,X,w,e_var)
%ML_GAUSS_LIK Gaussian Likelihood
%
%   input -----------------------------------------------------------------
%
%       o y     : (N x 1), output.
%
%       o X     : (N x D), input, N samples of dimension D.
%
%       o w     : (D+1 x 1), linear regressor weights.
%
%       o e_var : (1 x 1), noise variance.
%
%   output ----------------------------------------------------------------
%
%       o lik   : (N x 1),
%

[N,D] = size(X);

% (N x 1) - (N x D+1) (D+1 x 1)
% (N x 1) - (N x 1)

diff = y - [X,ones(N,1)] * w';
z     = sum(diff .* diff);

lik  = 1.0/( (2 * pi * e_var)^(D/2)) * exp(-1/(2*e_var) *  z);

end

author: epfl-lasa


text:
function [ a, kpvt, info ] = dsifa ( a, lda, n )

%*****************************************************************************80
%
%% DSIFA factors a real symmetric matrix.
%
%  Discussion:
%
%    To solve A*X = B, follow DSIFA by DSISL.
%
%    To compute inverse(A)*C, follow DSIFA by DSISL.
%
%    To compute determinant(A), follow DSIFA by DSIDI.
%
%    To compute inertia(A), follow DSIFA by DSIDI.
%
%    To compute inverse(A), follow DSIFA by DSIDI.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 June 2005
%
%  Author:
%
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Dongarra, Moler, Bunch and Stewart,
%    LINPACK User's Guide,
%    SIAM, (Society for Industrial and Applied Mathematics),
%    3600 University City Science Center,
%    Philadelphia, PA, 19104-2688.
%    ISBN 0-89871-172-X
%
%  Parameters:
%
%    Input, real A(LDA,N), the symmetric matrix to be factored.  Only
%    the diagonal and upper triangle are used.
%
%    Input, integer LDA, the leading dimension of the array A.
%
%    Input, integer N, the order of the matrix.
%
%    Output, real A(LDA,N), a block diagonal matrix and the multipliers which
%    were used to obtain it.  The factorization can be written A = U*D*U'
%    where U is a product of permutation and unit upper triangular
%    matrices, U' is the transpose of U, and D is block diagonal
%    with 1 by 1 and 2 by 2 blocks.
%
%    Output, integer KPVT(N), the pivot indices.
%
%    Output, integer INFO, error flag.
%    0, normal value.
%    K, if the K-th pivot block is singular.  This is not an error
%    condition for this subroutine, but it does indicate that DSISL
%    or DSIDI may divide by zero if called.
%

%
%  ALPHA is used in choosing pivot block size.
%
  alpha = ( 1.0 + sqrt ( 17.0 ) ) / 8.0;

  info = 0;
%
%  Main loop on K, which goes from N to 1.
%
  k = n;

  while ( 0 < k )

    if ( k == 1 )
      kpvt(1) = 1;
      if ( a(1,1) == 0.0 )
        info = 1;
      end
      return
    end
%
%  This section of code determines the kind of
%  elimination to be performed.  When it is completed,
%  KSTEP will be set to the size of the pivot block, and
%  SWAP will be set to .true. if an interchange is required.
%
    km1 = k - 1;
    absakk = abs ( a(k,k) );
%
%  Determine the largest off-diagonal element in column K.
%
    imax = idamax ( k-1, a(1:k-1,k), 1 );
    colmax = abs ( a(imax,k) );

    if ( alpha * colmax <= absakk )

      kstep = 1;
      swap = 0;
%
%  Determine the largest off-diagonal element in row IMAX.
%
    else

      rowmax = 0.0;
      imaxp1 = imax + 1;
      for j = imaxp1 : k
        rowmax = max ( rowmax, abs ( a(imax,j) ) );
      end

      if ( imax ~= 1 )
        jmax = idamax ( imax-1, a(1:imax-1,imax), 1 );
        rowmax = max ( rowmax, abs ( a(jmax,imax) ) );
      end

      if ( alpha * rowmax <= abs ( a(imax,imax) ) )
        kstep = 1;
        swap = 1;
      elseif ( alpha * colmax * ( colmax / rowmax ) <= absakk )
        kstep = 1;
        swap = 0;
      else
        kstep = 2;
        swap = ( imax ~= k-1 )
      end

    end
%
%  Column K is zero.
%  Set INFO and iterate the loop.
%
    if ( max ( absakk, colmax ) == 0.0 )

      kpvt(k) = k;
      info = k;
%
%  1 x 1 pivot block.
%
%  Perform an interchange.
%
    elseif ( kstep ~= 2 )

      if ( swap )

        [ a(1:imax,imax), a(1:imax,k) ] = ...
          dswap ( imax, a(1:imax,imax), 1, a(1:imax,k), 1 );

        for jj = imax : k
          j = k + imax - jj;
          t = a(j,k);
          a(j,k) = a(imax,j);
          a(imax,j) = t;
        end

      end
%
%  Perform the elimination.
%
      for jj = 1 : k-1
        j = k - jj;
        mulk = -a(j,k) / a(k,k);
        t = mulk;
        a(1:j,j) = daxpy ( j, t, a(1:j,k), 1, a(1:j,j), 1 );
        a(j,k) = mulk;
      end
%
%  Set the pivot array.
%
      if ( swap )
        kpvt(k) = imax;
      else
        kpvt(k) = k;
      end
%
%  2 x 2 pivot block.
%
%  Perform an interchange.
%
    else

      if ( swap )

        [ a(1:imax,imax), a(1:imax,k-1) ] = ...
          dswap ( imax, a(1:imax,imax), 1, a(1:imax,k-1), 1 );

        for jj = imax : k-1
          j = k-1 + imax - jj;
          t = a(j,k-1);
          a(j,k-1) = a(imax,j);
          a(imax,j) = t;
        end

        t = a(k-1,k);
        a(k-1,k) = a(imax,k);
        a(imax,k) = t;

      end
%
%  Perform the elimination.
%
      if ( k-2 ~= 0 )

        ak = a(k,k) / a(k-1,k);
        akm1 = a(k-1,k-1) / a(k-1,k);
        denom = 1.0 - ak * akm1;

        for jj = 1 : k-2

          j = k - 1 - jj;
          bk = a(j,k) / a(k-1,k);
          bkm1 = a(j,k-1) / a(k-1,k);
          mulk = ( akm1 * bk - bkm1 ) / denom;
          mulkm1 = ( ak * bkm1 - bk ) / denom;
          t = mulk;
          a(1:j,j) = daxpy ( j, t, a(1:j,k), 1, a(1:j,j), 1 );
          t = mulkm1;
          a(1:j,j) = daxpy ( j, t, a(1:j,k-1), 1, a(1:j,j), 1 );
          a(j,k) = mulk;
          a(j,k-1) = mulkm1;

        end

      end
%
%  Set the pivot array.
%
      if ( swap )
        kpvt(k) = -imax;
      else
        kpvt(k) = 1 - k;
      end

      kpvt(k-1) = kpvt(k);

    end

    k = k - kstep;

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% STD_ERSP - Compute ERSP and/or ITC transforms for ICA components
%              or data channels of a dataset. Save results into Matlab
%              float files.
%
% Function description:
%              The function computes the mean ERSP or ITC for the selected
%              dataset ICA components or data channels in the requested
%              frequency range and time window (the two are dependent).
%              Frequencies are equally log spaced. Options specify component
%              numbers, desired frequency range,  time window length,
%              frequency resolution, significance level, and wavelet
%              cycles. See >> help newtimef and >> timef details
%
%              Two Matlab files are saved (for ERSP and ITC). These contain
%              the ERSP|ITC image, plus the transform parameters
%              used to compute them. Saves the computed dataset mean images
%              in dataset-name files with extensions '.icaersp' and '.icaitc'
%              for ICA components or '.datersp', '.datitc' for data channels.
% Usage:
%              >> [X times logfreqs ] = std_ersp(EEG, 'key', 'val', ...);
% Inputs:
%   EEG          - a loaded epoched EEG dataset structure. May be an array
%                  of such structure containing several datasets.
%
% Other inputs:
%   'trialindices' - [cell array] indices of trials for each dataset.
%                  Default is EMPTY (no trials). NEEDS TO BE SET.
%   'components' - [numeric vector] components of the EEG structure for which
%                  activation spectrum will be computed. Note that because
%                  computation of ERP is so fast, all components spectrum are
%                  computed and saved. Only selected component
%                  are returned by the function to Matlab
%                  {default|[] -> all}
%   'channels'   - [cell array] channels of the EEG structure for which
%                  activation spectrum will be computed. Note that because
%                  computation of ERP is so fast, all channels spectrum are
%                  computed and saved. Only selected channels
%                  are returned by the function to Matlab
%                  {default|[] -> none}
%   'recompute'  - ['on'|'off'] force recomputing data file even if it is
%                  already on disk.
%   'rmcomps'    - [integer array] remove artifactual components (this entry
%                  is ignored when plotting components). This entry contains
%                  the indices of the components to be removed. Default is none.
%   'interp'     - [struct] channel location structure containing electrode
%                  to interpolate ((this entry is ignored when plotting
%                  components). Default is no interpolation.
%   'fileout'    - [string] name of the file to save on disk. The default
%                  is the same name (with a different extension) as the
%                  dataset given as input.
%  'savetrials'  - ['on'|'off'] save single-trials ERSP. Requires a lot of disk
%                  space (dataset space on disk times 10) but allow for refined
%                  single-trial statistics. This option is obsolete as
%                  trials are now always saved.
%  'savefile'    - ['on'|'off'] save file or simply return measures.
%                  Default is to save files ('on').
%  'getparams'   - ['on'|'off'] return optional parameters for the newtimef
%                  function (and do not compute anything). This argument is
%                  obsolete (default is 'off').
%
% ERSP optional inputs:
%   'type'       - ['ersp'|'itc'|'ersp&itc'] save ERSP, ITC, or both data
%                  types to disk {default: 'ersp'}
%   'freqs'      - [minHz maxHz] the ERSP/ITC frequency range to compute
%                  and return. {default: 3 to EEG sampling rate divided by 3}
%   'timelimits' - [minms maxms] time window (in ms) to compute.
%                  {default: whole input epoch}.
%   'cycles'     - [wavecycles (factor)]. If 0 -> DFT (constant window length
%                  across frequencies).
%                  If >0 -> the number of cycles in each analysis wavelet.
%                  If [wavecycles factor], wavelet cycles increase with
%                  frequency, beginning at wavecyles. (0 < factor < 1)
%                  factor = 0 -> fixed epoch length (DFT, as in FFT).
%                  factor = 1 -> no increase (standard wavelets)
%                  {default: [0]}
%   'padratio'   - (power of 2). Multiply the number of output frequencies
%                  by dividing their frequency spacing through 0-padding.
%                  Output frequency spacing is (low_freq/padratio).
%   'alpha'      - If in (0, 1), compute two-tailed permutation-based
%                  probability thresholds and use these to mask the output
%                  ERSP/ITC images {default: NaN}
%   'powbase'    - Deprecated. Note that baseline can be readjusted after
%                  computation as single trial spectral decompositions are stored.
%
% Other optional inputs:
%   This function will take any of the NEWTIMEF optional inputs (for instance
%   to compute log-space frequencies)...
%
% Outputs:
%   X         - the masked log ERSP/ITC of the requested ICA components/channels
%               in the selected frequency and time range. Note that for
%               optimization reasons, this parameter is now empty or 0. X
%               thus must be read from the datafile saved on disk.
%   times     - vector of time points for which the ERSPs/ITCs were computed.
%   logfreqs  - vector of (equally log spaced) frequencies (in Hz) at which the
%               log ERSP/ITC was evaluated.
%   parameters - parameters given as input to the newtimef function.
%
% Files written or modified:
%              [dataset_filename].icaersp   <-- saved component ERSPs
%              [dataset_filename].icaitc    <-- saved component ITCs
%              [dataset_filename].icatimef  <-- saved component single
%                                               trial decompositions.
%  OR for channels
%              [dataset_filename].datersp   <-- saved channel ERSPs
%              [dataset_filename].datitc    <-- saved channel ITCs
%              [dataset_filename].dattimef  <-- saved channel single
%                                               trial decompositions.
% Example:
%            % Create mean ERSP and ITC images on disk for all comps from
%            % dataset EEG use three-cycle wavelets (at 3 Hz) to more than
%            % three-cycle wavelets at 50 Hz. See >> help newtimef
%            % Return the (equally log-freq spaced, probability-masked) ERSP.
%            >> [Xersp, times, logfreqs] = std_ersp(EEG, ...
%                       'type', 'ersp', 'freqs', [3 50], 'cycles', [3 0.5]);
%
% See also: TIMEF, STD_ITC, STD_ERP, STD_SPEC, STD_TOPO, STD_PRECLUST
%
% Authors: Arnaud Delorme, Hilit Serby, SCCN, INC, UCSD, January, 2005-

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

function [X, times, logfreqs, parameters] = std_ersp(EEG, varargin)

if nargin < 1
    help std_ersp;
    return;
end

X = [];
options = {};
if length(varargin) > 1
    if ~ischar(varargin{1})
        if length(varargin) > 0, options = { options{:} 'components' varargin{1} }; end
        if length(varargin) > 1, options = { options{:} 'freqs'      varargin{2} }; end
        if length(varargin) > 2, options = { options{:} 'timewindow' varargin{3} }; end
        if length(varargin) > 3, options = { options{:} 'cycles'     varargin{4} }; end
        if length(varargin) > 4, options = { options{:} 'padratio'   varargin{5} }; end
        if length(varargin) > 5, options = { options{:} 'alpha'      varargin{6} }; end
        if length(varargin) > 6, options = { options{:} 'type'       varargin{7} }; end
        if length(varargin) > 7, options = { options{:} 'powbase'    varargin{8} }; end
    else
        options = varargin;
    end
end

[g, timefargs] = finputcheck(options, { ...
                        'components'    'integer'               []          [];
                        'channels'      { 'cell','integer' }    { [] [] }   {};
                        'powbase'       'real'                  []          [];
                        'trialindices' { 'integer','cell' }     []          [];
                        'savetrials'    'string'      { 'on','off' }        'off';
                        'plot'          'string'      { 'on','off' }        'off'; % not documented for debugging purpose
                        'recompute'     'string'      { 'on','off' }        'off';
                        'getparams'     'string'      { 'on','off' }        'off';
                        'savefile'      'string'      { 'on','off' }        'on';
                        'parallel'      'string'      { 'on','off' }        'off';
                        'timewindow'    'real'                  []          [];    % ignored, deprecated
                        'fileout'       'string'                []          '';
                        'timelimits'    'real'                  []          [EEG(1).xmin EEG(1).xmax]*1000;
                        'cycles'        'real'                  []          [3 .5];
                        'padratio'      'real'                  []          1;
                        'trialinfo'     'struct'                []          struct([]);
                        'freqs'         'real'                  []          [0 EEG(1).srate/2];
                        'rmcomps'       'cell'                  []          cell(1,length(EEG));
                        'interp'        'struct'                { }         struct([]);
                        'freqscale'     'string'                []         'log';
                        'alpha'         'real'                  []          NaN;
                        'baseline'      'real'                  []          0;
                        'type'          'string'      { 'ersp','itc','both','ersp&itc' }  'both'}, 'std_ersp', 'ignore');
if ischar(g), error(g); end
if isempty(g.trialindices), g.trialindices = cell(length(EEG)); end
if ~iscell(g.trialindices), g.trialindices = { g.trialindices }; end
if ~isempty(g.powbase), disp('''powbase'' parameter is no longer supported at computation time'); end

% checking input parameters
% -------------------------
if isempty(g.components) && isempty(g.channels)
    if isempty(EEG(1).icaweights)
        error('EEG.icaweights not found');
    end
    g.components = 1:size(EEG(1).icaweights,1);
    disp('Computing ERSP with default values for all components of the dataset');
end

% select ICA components or data channels
% --------------------------------------
if isempty(g.fileout), g.fileout = fullfile(EEG(1).filepath, EEG(1).filename(1:end-4)); end
if ~isempty(g.components)
    g.indices = g.components;
    prefix = 'comp';
    filenameersp   = [ g.fileout '.icaersp'  ];
    filenameitc    = [ g.fileout '.icaitc'   ];
    filenametrials = [ g.fileout '.icatimef' ];
    if ~isempty(g.channels)
        error('Cannot compute ERSP/ITC for components and channels at the same time');
    end
elseif ~isempty(g.channels)
    if iscell(g.channels)
        if ~isempty(g.interp)
            g.indices = eeg_chaninds(g.interp, g.channels, 0);
        else
            g.indices = eeg_chaninds(EEG(1), g.channels, 0);
            for ind = 2:length(EEG)
                if ~isequal(eeg_chaninds(EEG(ind), g.channels, 0), g.indices)
                    error([ 'Channel information must be consistent when ' 10 'several datasets are merged for a specific design' ]);
                end
            end
        end
    else
        g.indices = g.channels;
    end
    prefix = 'chan';
    filenameersp   = [ g.fileout '.datersp'  ];
    filenameitc    = [ g.fileout '.datitc'   ];
    filenametrials = [ g.fileout '.dattimef' ];
end

% Check if ERSP/ITC information found in datasets and if fits requested parameters
% ----------------------------------------------------------------------------
if exist( filenametrials ) && strcmpi(g.recompute, 'off')
    fprintf('Use existing file for ERSP: %s; check the ''recompute checkbox'' to force recomputing.\n', filenameersp);
    return;
end

% Compute ERSP parameters
% -----------------------
parameters = { 'cycles', g.cycles, 'padratio', g.padratio, ...
               'alpha', g.alpha, 'freqscale', g.freqscale, timefargs{:} };
defaultlowfreq = 3;
[time_range] = compute_ersp_times(g.cycles,  EEG(1).srate, ...
                                 [EEG(1).xmin EEG(1).xmax]*1000 , defaultlowfreq, g.padratio);
if time_range(1) < time_range(2) && g.freqs(1) == 0
     g.freqs(1) = defaultlowfreq; % for backward compatibility
end
parameters = { parameters{:} 'freqs' g.freqs };
if strcmpi(g.plot, 'off')
    parameters = { parameters{:} 'plotersp', 'off', 'plotitc', 'off', 'plotphase', 'off' };
end
parameters{end+1} = 'baseline';
parameters{end+1} = g.baseline;

% return parameters
% -----------------
if strcmpi(g.getparams, 'on')
    X = []; times = []; logfreqs = [];
    if strcmpi(g.savetrials, 'on')
        parameters = { parameters{:} 'savetrials', g.savetrials };
    end
    return;
end

options = {};
if ~isempty(g.rmcomps), options = { options{:} 'rmcomps' g.rmcomps }; end
if ~isempty(g.interp),  options = { options{:} 'interp' g.interp }; end
if isempty(g.channels)
     X = eeg_getdatact(EEG, 'component', g.indices, 'trialindices', g.trialindices );
else X = eeg_getdatact(EEG, 'channel'  , g.indices, 'trialindices', g.trialindices, 'rmcomps', g.rmcomps, 'interp', g.interp);
end
if size(X, 3) == 1
    error('The data is continuous for one of the dataset. ERSP can only be computed when data trials are present');
end

% frame range
% -----------
pointrange1 = round(max((g.timelimits(1)/1000-EEG(1).xmin)*EEG(1).srate, 1));
pointrange2 = round(min(((g.timelimits(2)+1000/EEG(1).srate)/1000-EEG(1).xmin)*EEG(1).srate, EEG(1).pnts));
pointrange = [pointrange1:pointrange2];

% Compute ERSP && ITC
% ------------------
allTrialsTmp   = cell(1,length(g.indices));
allTrialsTime  = cell(1,length(g.indices));
allTrialsFreqs = cell(1,length(g.indices));
eeglab_options;
usesingle = option_single;

disp('Computing time/frequency decomposition...');
parfor k = 1:length(g.indices)
    tmpparams = parameters;
    if length(g.indices) > 1
        tmpparams{end+1} = 'verbose';
        tmpparams{end+1} = 'off';
    end

    % Run timef() to get ERSP
    % ------------------------
    timefdata  = reshape(X(k,pointrange,:), 1, length(pointrange)*size(X,3));
    mytimes = [];
    mylogfreqs = [];
    alltfX   = [];
    if ~isempty(timefdata)
        [logersp,logitc,logbase,mytimes,mylogfreqs,logeboot,logiboot,alltfX] ...
              = newtimef( timefdata, length(pointrange), g.timelimits, EEG(1).srate, tmpparams{2:end});
        %figure; newtimef( TMP.data(32,:), EEG.pnts, [EEG.xmin EEG.xmax]*1000, EEG.srate, cycles, 'freqs', freqs);
        %figure; newtimef( timefdata, length(pointrange), g.timelimits, EEG.srate, cycles, 'freqs', freqs);
    end
    %if strcmpi(g.plot, 'on'), return; end
    if usesingle
        alltfX = single(alltfX);
    end

    allTrialsTmp{k}   = single( alltfX );
    allTrialsTime{k}  = mytimes;
    allTrialsFreqs{k} = mylogfreqs;
end
all_trials = [];
for k = 1:length(g.indices)  % for each (specified) component/channel
    all_trials = setfield( all_trials, [ prefix int2str(g.indices(k)) ], allTrialsTmp{k});
end
X = allTrialsTmp{1};

% Save ERSP into file
% -------------------
logfreqs             = allTrialsFreqs{1};
times                = allTrialsTime{1};
all_trials.freqs     = allTrialsFreqs{1};
all_trials.times     = allTrialsTime{1};
all_trials.parameters = { options{:} parameters{:} };
all_trials.datatype   = 'TIMEF';
all_trials.datafiles  = computeFullFileName( { EEG.filepath }, { EEG.filename });
all_trials.datatrials = g.trialindices;

all_trials.parameters = parameters;
if ~isempty(g.channels)
    if ~isempty(g.interp)
        all_trials.labels = { g.interp(g.indices).labels };
    elseif ~isempty(EEG(1).chanlocs)
        tmpchanlocs = EEG(1).chanlocs;
        all_trials.labels = { tmpchanlocs(g.indices).labels };
    end
end
all_trials.trialinfo = g.trialinfo;

if strcmpi(g.savefile, 'on')
    std_savedat( filenametrials , all_trials );
end

% compute full file names
% -----------------------
function res = computeFullFileName(filePaths, fileNames);
for index = 1:length(fileNames)
    res{index} = fullfile(filePaths{index}, fileNames{index});
end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
% LORENZ system
% System identification: NARX

clear all, close all, clc
figpath = '../FIGURES/LORENZ/';
datapath = '../DATA/LORENZ/';
addpath('../utils');

SystemModel = 'LORENZ';

%% Generate Data
ONLY_TRAINING_LENGTH = 1;
InputSignalType = 'sphs'; % prbs; chirp; noise; sine2; sphs; mixed
getTrainingData

%% NARX: Training
ModelName = 'NARX';
Nvar = 3;
rng(2,'twister')
SUBSTRACT_MEAN = 0;
xmean = mean(x)';

if SUBSTRACT_MEAN == 1
    xtrain = x' - repmat(xmean',[1 size(x',2)]);
    utrain = u;
else
    xtrain = x';
    utrain = u;
end

% prepare training data
yt = con2seq(xtrain);
yi = con2seq(utrain);

% Neural network
stateDelays = 1;        % state delay vector
inputDelays = 1;        % input delay vector
hiddenSizes = [10];     % network structure (number of neurons per layer)

% Nonlinear autoregressive neural network
net = narxnet(inputDelays,stateDelays, hiddenSizes);

% Training parameters %nnstart
net.trainFcn = 'trainlm';
net.trainParam.min_grad = 1e-10;
net.trainParam.showCommandLine = 1;

% Prepares training data (shifting, copying feedback targets into inputs as needed, etc.)
[Us,Ui,Si,Ss] = preparets(net,yi,{},yt);

% Train net with prepared training data in open-loop
tic
net = train(net,Us,Ss,Ui,Si);
toc
% view(net)

% Close loop for recursive prediction
netc = closeloop(net);

%% Prediction over training phase
% Prepare validation data / Get initial state from training data
[Us,Ui,Si,So] = preparets(netc,yi,{},yt);

% Predict on validation data
predict = netc(Us,Ui,Si);
xNARX = cell2mat(predict)';

if SUBSTRACT_MEAN == 1
    xNARX = xNARX + repmat(xmean,[size(xNARX,1) 1]);
end

% Error
e = cell2mat(gsubtract(So,predict));

%% Show validation
clear ph
figure,box on,
ccolors = get(gca,'colororder');
ccolors_valid = [ccolors(1,:)-[0 0.2 0.2];ccolors(2,:)-[0.1 0.2 0.09];ccolors(3,:)-[0.1 0.2 0.09]];
for i = 1:Nvar
    ph(i) = plot(tspan,x(:,i),'-','Color',ccolors(i,:),'LineWidth',1); hold on
end
for i = 1:Nvar
    ph(Nvar+i) = plot(tspan(2:end),xNARX(:,i),'--','Color',ccolors_valid(i,:),'LineWidth',2);
end
xlim([0 (length(tspan)-1)*dt]), ylim([-25 50])
xlabel('Time')
ylabel('xi')
legend(ph([1,4]),'True',ModelName)
set(gca,'LineWidth',1, 'FontSize',14)
set(gcf,'Position',[100 100 300 200])
set(gcf,'PaperPositionMode','auto')
if SUBSTRACT_MEAN == 1
    print('-depsc2', '-loose', '-cmyk', [figpath,'EX_',SystemModel,'_SI_',ModelName,'_',InputSignalType,'_SUBSTRACT_MEAN.eps']);
else
    print('-depsc2', '-loose', '-cmyk', [figpath,'EX_',SystemModel,'_SI_',ModelName,'_',InputSignalType,'.eps']);
end

%% Validation 3D
filename = ['EX_',SystemModel,'_SI_',ModelName,'_',InputSignalType,'_train'];
xModel = xNARX;
iModel = 3;
xTRUTH = x;
color_type = 'models';
VIZ_3D_MODELvsTRUTH

%% Prediction
% prepare validation data
if SUBSTRACT_MEAN == 1
    xvalid = xv' - repmat(xmean',[1 size(xv',2)]);
    uvalid = uv;
else
    xvalid = xv';
    uvalid = uv;
end

yt_valid = con2seq(xvalid);
yi_valid = con2seq(uvalid);
[Us,Ui,Si,So] = preparets(netc,yi_valid,{},yt_valid);

% Reference
xA      = xv;
tA      = tv;

% Predict on validation data
predict = netc(Us,Ui,Si);
xB = cell2mat(predict)';
tB = tA(2:end);

if SUBSTRACT_MEAN == 1
    xB = xB + repmat(xmean,[size(xB,1) 1]);
end
%% Show training and prediction
u = u'; uv = uv';
VIZ_SI_Validation

%% Validation 3D
filename = ['EX_',SystemModel,'_SI_',ModelName,'_',InputSignalType,'_valid'];
xModel = xB;
xTRUTH = xA;
iModel = 3;
color_type = 'models';
VIZ_3D_MODELvsTRUTH

%% Save Data
Model.name = 'NARX';
Model.net = netc;
Model.xmean = xmean;
Model.stateDelays = stateDelays;
Model.inputDelays = inputDelays;
Model.hiddenSizes = hiddenSizes;
Model.SUBSTRACT_MEAN = SUBSTRACT_MEAN;
Model.dt = dt;
save(fullfile(datapath,['EX_',SystemModel,'_SI_',ModelName,'_',InputSignalType,'.mat']),'Model')

author: eurika-kaiser

repo: SINDY-MPC


text:
%  INTERNAL FUNCTION: multikronecker
%
%  ::
%
%    C=kronall(A1,A2,...,An)
%
%  Args:
%
%     - **Ai** [matrix]: input for the kronecker product
%
%  Returns:
%     :
%
%     - **C** [matrix]: kron(A1,kron(A2,kron(A3,...)))
%
%  See also:
%     tensorperm
%
%

author: jmaih


text:
%% Edge Element Discretization of Maxwell Equations
% We test Maxwell solvers in iFEM.

clear all; close all
rowNames ={'h=1/2';'h=1/4';'h=1/8'};
colHeaders = {'H^curl Error','L^2 Error'};

%% The data of the pde
%
% * pde = Maxwelldata1; % zero Neumann boundary condition and curl u = 0
% * pde = Maxwelldata2; % non-homogenous Neumann boundary condition
% * pde = Maxwelldata3; % polynomial data and curl u = 0
% * pde = Maxwelldata4; % zero Dirichlet boundary condition
% * pde = Maxwelldata5; % linear polynomial data
% * pde = planewavedataC; % plane wave with complex coefficients
% * pde = planewavedata1; % plane wave with real coefficients

%% Positive Definite Case
% curl curl E + E = f.

help Maxwelldata2
%%
% Dirichlet boundary condition
[node,elem,HB] = cubemesh([-1,1,-1,1,-1,1],0.25);
pde = Maxwelldata2;
bdFlag = setboundary3(node,elem,'Dirichlet');
option.solver = 'cg';
cubeMaxwell2;
%%
makeHtmlTable([energyErr L2Err],[],rowNames,colHeaders);

%%
% Neumann boundary condition
[node,elem,HB] = cubemesh([-1,1,-1,1,-1,1],1);
pde = Maxwelldata2;
bdFlag = setboundary3(node,elem,'Neumann');
option.solver = 'cg';
cubeMaxwell2;
%%
makeHtmlTable([energyErr L2Err],[],rowNames,colHeaders);

%%
% Optimal first order of convergence is achieved. HX preconditioned CG
% converges around 20 steps.

%% Indefinite with real coefficients
% curl curl E - E = f.

help planewavedata1
%%
% Dirichlet boundary condition
[node,elem,HB] = cubemesh([-1,1,-1,1,-1,1],1);
pde = planewavedata1;
bdFlag = setboundary3(node,elem,'Dirichlet');
option.solver = 'cg';
cubeMaxwell2;
%%
makeHtmlTable([energyErr L2Err],[],rowNames,colHeaders);

%%
% Neumann boundary condition
[node,elem,HB] = cubemesh([-1,1,-1,1,-1,1],1);
pde = planewavedata1;
bdFlag = setboundary3(node,elem,'Neumann');
option.solver = 'cg';
cubeMaxwell2;
%%
makeHtmlTable([energyErr L2Err],[],rowNames,colHeaders);

%%
% Optimal first order of convergence is achieved. HX preconditioned CG
% converges around 40 steps although the system is indefinite.

%% Indefinite: complex coefficients, real solution
% curl curl E - (1-i)E = f.

help planewavedataC
%%
% Dirichlet boundary condition
[node,elem,HB] = cubemesh([-1,1,-1,1,-1,1],1);
pde = planewavedataC;
bdFlag = setboundary3(node,elem,'Dirichlet');
option.solver = 'gmres';
cubeMaxwell2;
%%
makeHtmlTable([energyErr L2Err],[],rowNames,colHeaders);

%%
% Neumann boundary condition
[node,elem,HB] = cubemesh([-1,1,-1,1,-1,1],1);
pde = planewavedataC;
bdFlag = setboundary3(node,elem,'Neumann');
option.solver = 'gmres';
cubeMaxwell2;
%%
makeHtmlTable([energyErr L2Err],[],rowNames,colHeaders);

%%
% Optimal first order of convergence is achieved. HX preconditioned GMRES
% converges around 90 steps although the system is indefinite.
%
% Bug: For Neumann boundary condition, the rate of L2 error is not quite right.

%% Indefinite: real coefficents, complex solution
help planewavedata
%%
% Dirichlet boundary condition
[node,elem,HB] = cubemesh([-1,1,-1,1,-1,1],1);
bdFlag = setboundary3(node,elem,'Dirichlet');
option.solver = 'cg';
planewaveMaxwell2;
%%
makeHtmlTable([energyErr L2Err],[],rowNames,colHeaders);

%%
% Neumann boundary condition
[node,elem,HB] = cubemesh([-1,1,-1,1,-1,1],1);
bdFlag = setboundary3(node,elem,'Neumann');
option.solver = 'cg';
planewaveMaxwell2;
%%
makeHtmlTable([energyErr L2Err],[],rowNames,colHeaders);

%%
% The computation of error can't handle complex functions. So in
% planewaveMaxwell the error between uI and uh is computed. Therefore
% slightly better rate of convergence is observed. The rate of L2 error is
% almost second order.

repo: Aero-matlab


text:
function output = trim_signals(input, offset)
    % Trim signals according to offsets
    N = length(offset);
    output = cell(size(input));
    frame_index = floor(offset/N);
%     display(frame_index);
    trim_index = max(frame_index) - frame_index;
%     display(trim_index);
    for i=1:N
        output{i} = input{i}((trim_index(i)+1):end);
    end

    % trim end so all signals have same length
    % required by non-uniform reconstruction algorithm
    minlen = min(cellfun('length', output));
    for i = 1:N
       output{i} = output{i}(1:minlen);
    end
end

author: yueyuzhao

repo: gyrophone

save_path: github-repos/MATLAB/yueyuzhao-gyrophone


text:
%--------------------------------------------------------------------------
%   [dataout] = steeringvector(array_pos,sig_pos)
%--------------------------------------------------------------------------
%   功能:
%   空间传播方法计算导向矢量
%--------------------------------------------------------------------------
%   输入:
%           array_pos               接收信号阵列坐标，按照列排列
%           sig_pos                 发射信号阵列坐标 按照列排列
%           lambda                  载频波长
%   输出:
%           A                       导向矢量
%--------------------------------------------------------------------------
%  例子:
% c = 3e8;
% fc = 20e9;
% lambda = c/fc;
% tgt_angle = [50 90 140];
% R = [1000 5600 9000];
%
% tgt_pos = [R.*cosd(tgt_angle);R.*sind(tgt_angle);zeros(1,3)];
% radar_pos = [zeros(1,8);lambda/2*(1:8)-4*lambda/2 - lambda/4;zeros(1,8)];
% steering_vector(radar_pos,tgt_pos(:,1),lambda)
%--------------------------------------------------------------------------
function A = steering_vector(array_pos,sig_pos,lambda)
for idx = 1:size(array_pos,2)
    R(idx,:) = norm(array_pos(:,idx)-sig_pos);
end
A = exp(1j.*2*pi*R./lambda);
A = conj(A);
end

repo: SP


text:
report_this_filefun(mfilename('fullpath'));

bz2 = [];
p1 = [];
p99 = [];
p50 = [];
pmab = [];
pma = [];
pmib = [];
pmi = [];
si = [];

tdiff = round((teb - t0b)*365/par1);
def = {'100','300','100','2'};
tit ='Random Zmax calculation';
prompt={'Number of bins in each time series?', 'Number of random samples drawn? (Grid-size)',...
    'Number of repeats?','Window length in years ?'};

ni2 = inputdlg(prompt,tit,1,def);
l = ni2{4};
iwl = str2double(l);
l = ni2{3};
nr = str2double(l);
l = ni2{2};
no = str2double(l);
l = ni2{1};
ni = str2double(l);

iwl0 = iwl;
iwl = floor(iwl*365/par1);
zr = (-15:0.1:15)*0;
n0 = no;
na = no;

titStr ='Warning!                                        ';

messtext= ...
    ['                                                '
    ' This rotine sometimes takes a long time!       '
    '  and may run out of memory. You can interupt   '
    ' the calculation with a ^C. The results         '
    ' calculated so far are stored in the variable   '
    ' pmab and save to a *.mat file                  '];

zmap_message_center.set_message(titStr,messtext);
figure_w_normalized_uicontrolunits(mess)

[newmatfile, newpath] = uiputfile([hodi '/*.mat'], 'Filename for saving of intermediate results ');

wai = waitbar(0,' Please Wait ...  ');
set(wai,'NumberTitle','off','Name','Percent done');;
watchon
think

con = 0;
rng('shuffle')
for k=1:nr
    na = n0; no = n0;
    while na+100 > 100;
        if na  > 100 ; na == 100; end
        con = con+1;
        l = ceil(rand([ni na])*a.Count);
        [cumu, xt] = hist(reshape(a(l,3),ni,na),(t0b:par1/365:teb));
        for j = 2:tdiff-iwl
            cu = [cumu(1:j-1,:) ; cumu(j+iwl+1:length(cumu(:,1)),:)];
            mean1 = mean(cu);
            mean2 = mean(cumu(j:j+iwl,:));
            %var1 = diag(cov(cu));
            var1 = (std(cu)).^2;;
            var2 = (std(cumu(j:j+iwl,:)).^2);
            as = (mean1 - mean2)./(sqrt(var1/(length(cumu(:,1))-iwl)+var2/iwl));
            pma = [pma max(as)];
            pmi = [pmi min(as)];
        end     % for j
        pmab = [pmab max(pma) ];
        pmib = [pmib min(pmi) ];
        pma = [];pmi =[];
        no = no - 100;
        na = no;
    end   % while na
    pmab = [pmab max(pma) ];
    pmib = [pmib min(pmi) ];
    waitbar(k/nr);
    do = ['save '  newpath newmatfile  ' pmab pma pmib pmi iwl0  con n0 nr ni '];
    eval(do,'disp(''Error while trying to save intermediate results! Permission? '') ')
end   % for k

close(wai)

figure

histogram(pmab)
set(gca,'box','on',...
    'SortMethod','childorder','TickDir','out','FontWeight',...
    'bold','FontSize',ZmapGlobal.Data.fontsz.m,'Linewidth',1.2)
grid
xlabel('Zmax')
ylabel('Number ')
clear title
title(['ni =' num2str(ni) ', #samples = ' num2str(n0) ', #repeats=' num2str(con*nr), ', Tw= ' num2str(iwl0)]) ;

matdraw

figure

histogram(pmib)
set(gca,'box','on',...
    'SortMethod','childorder','TickDir','out','FontWeight',...
    'bold','FontSize',ZmapGlobal.Data.fontsz.m,'Linewidth',1.2)
grid
xlabel('Zmin')
ylabel('Number ')
clear title
title(['ni =' num2str(ni) ', #samples = ' num2str(n0) ', #repeats=' num2str(con*nr), ', Tw= ' num2str(iwl0)]) ;

matdraw

author: CelsoReyes


text:
function pick = nmsMe(boxes, overlap)
% Non-maximum suppression. (FAST VERSION)
% Greedily select high-scoring detections and skip detections
% that are significantly covered by a previously selected
% detection.
% NOTE: This is adapted from Pedro Felzenszwalb's version (nms.m),
% but an inner loop has been eliminated to significantly speed it
% up in the case of a large number of boxes

% modified based on Tomasz Malisiewicz's esvm code

if isempty(boxes)
  pick = [];
  return;
end

x1 = boxes(:,1);
y1 = boxes(:,2);
x2 = boxes(:,3);
y2 = boxes(:,4);
s = boxes(:,end);

area = (x2-x1+1) .* (y2-y1+1);
[~, I] = sort(s);

pick = s*0;
counter = 1;
while ~isempty(I)

  last = length(I);
  i = I(last);
  pick(counter) = i;
  counter = counter + 1;

  xx1 = max(x1(i), x1(I(1:last-1)));
  yy1 = max(y1(i), y1(I(1:last-1)));
  xx2 = min(x2(i), x2(I(1:last-1)));
  yy2 = min(y2(i), y2(I(1:last-1)));

  w = max(0.0, xx2-xx1+1);
  h = max(0.0, yy2-yy1+1);

  inter = w.*h;
  o = inter ./ (area(i) + area(I(1:last-1)) - inter);

  I([last; find(o>overlap)]) = [];
end

pick = pick(1:(counter-1));

author: thusiyuan


text:
function I = sum3(f)
%SUM3   Triple definite integral of a CHEBFUN3T over [a, b] x [c, d] x [e g].
%   This is a direct generalization of Clenshaw-Curtis quadrature as
%   implemented in Page 77 of Battles' PhD thesis.

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

if ( isempty(f) )
    I = [];
    return;
end

coeffs = f.coeffs;
% Rescaling factors:
rescaleFactor = 0.5*diff(f.domain(1:2));                 % (b - a)/2
rescaleFactor = rescaleFactor * 0.5*diff(f.domain(3:4)); % (d - c)/2
rescaleFactor = rescaleFactor * 0.5*diff(f.domain(5:6)); % (g - e)/2

if ( nargin == 1 )
    [nx, ny, nz] = size(coeffs);
    I1 = zeros(nx,ny);
    for i = 1:nx
        for j = 1:ny
            I1(i,j) = squeeze(coeffs(i,j,1:2:nz)).'*(2./(1-(0:2:nz-1)'.^2));
        end
    end
    I2 = zeros(nx,1);
    for i = 1:nx
        I2(i) = I1(i,1:2:ny)*(2./(1-(0:2:ny-1)'.^2));
    end
    I = I2(1:2:nx).'*(2./(1-(0:2:nx-1)'.^2));
end

% Rescale the output:
I = I*rescaleFactor;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [ x, seed ] = normal_ms_sample ( mu, sigma, seed )

%*****************************************************************************80
%
%% NORMAL_MS_SAMPLE samples the Normal MS PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real MU, SIGMA, the parameters of the PDF.
%    0.0 < SIGMA.
%
%    Input, integer SEED, a seed for the random number generator.
%
%    Output, real X, a sample of the PDF.
%
%    Output, integer SEED, an updated seed for the random number generator.
%
  [ y, seed ] = normal_01_sample ( seed );

  x = mu + sigma * y;

  return
end