text:
% STD_RMDAT - remove datasets from STUDY
%
% Usage:
%   >> [STUDY,ALLEEG] = std_rmdat(STUDY, ALLEEG, 'key', val);
%
% Inputs:
%   STUDY      - EEGLAB STUDY set
%   ALLEEG     - vector of the EEG datasets included in the STUDY structure
%
% Optional inputs:
%   'datinds'     - [integer array] indices of datasets to remove
%   'pntsrange'   - [min max] minimum and maximum of samples. Default is
%                   [0 Inf] (no constraint)
%   'chanrange'   - [min max] minimum and maximum of channels. Default is
%                   [0 Inf] (no constraint)
%   'sraterange'  - [min max] minimum and maximum for sampling rate. Default is
%                   [0 Inf] (no constraint)
%   'trialrange'  - [min max] minimum and maximum of trials. Default is
%                   [1 Inf]
%   'rmvarvalues' - {'string' range} remove datasets having variable value
%                   in the selected range. May also be {'string' 'value'}
%                   for non-numerical variables.
%   'keepvarvalues' - {'string' range} keep datasets having variable value
%                   in the selected range. May also be {'string' 'value'}
%                   for non-numerical variables.
%   'subjectind'  - [integer array] keep only specific subject indices in
%                   STUDY.subject. To remove or keep specificy subjects,
%                   use 'rmvarvalues' and 'keepvarvalues'
%   'checkeventtype' - [cell|array|string] check event type are present.
%   'numeventrange' - [min max] range for number of event of type above.
%                    Default is [1 Inf].
%
% Inputs:
%   STUDY      - EEGLAB STUDY set updated. The fields which is created or
%                updated is STUDY.datasetinfo.trialinfo
%
% Authors: Arnaud Delorme, SCCN/INC/UCSD, July 2022

% Copyright (C) Arnaud Delorme arno@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 [STUDY, ALLEEG, rmDats] = std_rmdat(STUDY, ALLEEG, varargin)

if nargin < 3
    help std_rmdat;
    return;
end

g = finputcheck( varargin, { ...
    'datinds'         'integer'  {} []; ...
    'chanrange'       'integer'  {} [0 Inf]; ...
    'pntsrange'       'integer'  {} [0 Inf]; ...
    'sraterange'      'float'    {} [0 Inf]; ...
    'trialrange'      'float'    {} [1 Inf]; ...
    'checkeventtype'  ''         {} []; ...
    'numeventrange'   'integer'  {} 1; ...
    'subjectind'   'integer'     {} []; ...
    'rmvarvalues'     'cell'     {} {}; ...
    'keepvarvalues'   'cell'     {} {} }, 'std_rmdat');
if isstr(g)
    error(g);
end

allPnts   = [ALLEEG.pnts];
allSrate  = [ALLEEG.srate];
allChans  = [ALLEEG.nbchan];
allTrials = [ALLEEG.trials];

% check pnts range
rmDats = zeros(1, length(allPnts));
if ~isempty(g.datinds)
    rmDats(g.datInds) = 1;
end
rmDats = g.pntsrange(1)  > allPnts   | allPnts   > g.pntsrange(2)  | rmDats;
rmDats = g.sraterange(1) > allSrate  | allSrate  > g.sraterange(2) | rmDats;
rmDats = g.chanrange(1)  > allChans  | allChans  > g.chanrange(2)  | rmDats;
rmDats = g.trialrange(1) > allTrials | allTrials > g.trialrange(2) | rmDats;

% check subject indices
% ---------------------
if ~isempty(g.subjectind)
    if ~isempty(g.keepvarvalues)
        error('Cannot use ''subjectind'' and ''keepvarvalues'' at the same time')
    end
    g.keepvarvalues = { 'subject' STUDY.subject(g.subjectind) };
end

% check variable name values
% --------------------------
if ~isempty(g.rmvarvalues)
    varName   = g.rmvarvalues{1};
    varValues = g.rmvarvalues{2};
    allValues = { STUDY.datasetinfo.(varName) };
    if ischar(varValues)
        allValues = cellfun(@num2str, allValues, 'uniformoutput', false);
        rmDats = rmDats | cellfun(@(x)isequal(x, varValues),  allValues);
    elseif iscell(varValues)
        for iCell = 1:length(varValues)
            allValues = cellfun(@num2str, allValues, 'uniformoutput', false);
            rmDats = rmDats | cellfun(@(x)isequal(x, varValues{iCell}),  allValues);
        end
    elseif length(varValues) ~= 2
        error('When providing nmumerical input for variable selection, there must be 2 values - min and max');
    else
        allValues = cellfun(@(x)fastif(ischar(x), str2double(x), x), allValues);
        rmVals = varValues(1) <= allValues & allValues <= varValues(2); % to deal with NaNs
        rmDats = rmDats | rmVals;
    end
end
if ~isempty(g.keepvarvalues)
    varName   = g.keepvarvalues{1};
    varValues = g.keepvarvalues{2};
    allValues = { STUDY.datasetinfo.(varName) };
    if ischar(varValues)
        allValues = cellfun(@num2str, allValues, 'uniformoutput', false);
        rmDats = rmDats | ~cellfun(@(x)isequal(x, varValues),  allValues);
    elseif iscell(varValues)
        keepDats = zeros(1, length(allPnts));
        for iCell = 1:length(varValues)
            allValues = cellfun(@num2str, allValues, 'uniformoutput', false);
            keepDats = keepDats | cellfun(@(x)isequal(x, varValues{iCell}),  allValues);
        end
        rmDats = rmDats | ~keepDats;
    elseif length(varValues) ~= 2
        error('When providing nmumerical input for variable selection, there must be 2 values - min and max');
    else
        allValues = cellfun(@(x)fastif(ischar(x), str2double(x), x), allValues);
        keepVals = varValues(1) <= allValues & allValues <= varValues(2); % to deal with NaNs
        rmDats = rmDats | ~keepVals;
    end
end

% check event type present
% --------------------------
if ~isempty(g.checkeventtype)
    if ischar(g.checkeventtype) g.checkeventtype = { g.checkeventtype }; end
    rmDatEvents = zeros(1, length(ALLEEG));
    for iDat = 1:length(ALLEEG)
        if isfield(ALLEEG(iDat).event, 'type')
            curEventTypes = { ALLEEG(iDat).event.type };
            if ischar(curEventTypes{1})
                for iType = 1:length(g.checkeventtype)
                    inds = strmatch(g.checkeventtype{iType}, curEventTypes, 'exact');
                    if length(inds) < g.numeventrange(1) || length(inds) > g.numeventrange(2)
                        rmDatEvents(iDat) = true;
                    end
                end
            else
                % numerical or string event types
                curEventTypes = [ ALLEEG(iDat).event.type ];
                for iType = 1:length(g.checkeventtype)
                    inds = find(g.checkeventtype(iType) == curEventTypes);
                    if length(inds) < g.numeventrange(1) || length(inds) > g.numeventrange(2)
                        rmDatEvents(iDat) = true;
                    end
                end
            end
        end
    end
    rmDats = rmDats | rmDatEvents;
end

if sum(rmDats) > 0
    fprintf('%d dataset meet the criteria for removal and have been removed\n', sum(rmDats));
    ALLEEG(rmDats) = [];
    STUDY.datasetinfo(rmDats) = [];
    for iDat = 1:length(STUDY.datasetinfo)
        STUDY.datasetinfo(iDat).index = iDat;
    end
else
    fprintf('No dataset meet the criteria for removal and have been removed\n');
end
STUDY.subject = {};
STUDY = std_checkset(STUDY, ALLEEG);

for iDesign = 1:length(STUDY.design)
    STUDY.design(iDesign).cases.value = intersect(STUDY.design(iDesign).cases.value, STUDY.subject);
end


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

re = []
for i = min(newt2.Magnitude):0.1:max(newt2.Magnitude)-0.5
    i
    l = newt2.Magnitude >= i;
    [bv magco stan av me,  mer] =  bvalcalc(newt2(l,:));
    re = [re ; i bv me];
end

figure
pl = plot(re(:,1),re(:,3),'b')
set(pl,'LineWidth',2.5)
hold on
pl = plot(re(:,1),re(:,2),':r')
set(pl,'LineWidth',2.5)
set(gca,'box','on',...
    'SortMethod','childorder','TickDir','out','FontWeight',...
    'bold','FontSize',ZmapGlobal.Data.fontsz.m,'Linewidth',1.2)

le =legend('-b','b(L2)','ro','b(maxli)');
grid
xlabel('Magnitude')
ylabel('b-value')

author: CelsoReyes


text:
function pass = test_fevalt(pref)
% Test fevalt

if ( nargin == 0 )
    pref = chebfunpref;
end
tol = 1e4*pref.cheb3Prefs.chebfun3eps;

ff = @(x,y,z) cos(x) + sin(x.*y) + sin(z.*x);
f = chebfun3(ff);
seedRNG(42);
x = rand(10,1);
y = rand(10,1);
z = rand(10,1);
F = fevalt(f,x,y,z);

[xx, yy, zz] = ndgrid(x,y,z);
Fexact = ff(xx,yy,zz);
pass(1) = norm(Fexact(:) - F(:)) < tol*vscale(f);

dom = [-1 1 -pi/3 pi/3 -2 0];
f = chebfun3(ff, dom);
x = chebpts(20,dom(1:2));
y = chebpts(20,dom(3:4));
z = chebpts(20,dom(5:6));
F = fevalt(f,x,y,z);

[xx, yy, zz] = ndgrid(x,y,z);
Fexact = ff(xx,yy,zz);
pass(2) = norm(Fexact(:) - F(:)) < tol*vscale(f);

% 'trig' flag + vector inputs
ff = @(x,y,z) sin(pi*(x+y+z));
f = chebfun3(ff, 'trig');
x = linspace(-1, 1, 20)';
y = linspace(-1, 1, 20)';
z = linspace(-1, 1, 20)';
F = fevalt(f,x,y,z);

[xx,yy,zz] = ndgrid(x,y,z);
Fexact = ff(xx,yy,zz);
pass(3) = norm(F(:) - Fexact(:)) < tol*vscale(f);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [inliers, info] = gnc(problem, f, varargin)
%GNC - Graduated Non-Convexity
% Implementation of:
%    - "Graduated Non-Convexity for Robust Spatial Perception: From Non-Minimal Solvers to Global Outlier Rejection"
%      Yang, Antonante, Tzoumas, Carlone (2020). IEEE Robotics and Automation Letters (RA-L), 5(2), 1127–1134.
%      https://arxiv.org/pdf/1909.08605.pdf
%    - "Outlier-Robust Estimation: Hardness, Minimally-Tuned Algorithms, and Applications"
%      Antonante, Tzoumas, Yang, Carlone (2020).
%      https://arxiv.org/pdf/2007.15109.pdf
%
% Syntax:  [inliers, info] = gnc(problem, f, ...)
%
% Inputs:
%    problem - the stracture representing a generic problem (see EmptyProblem)
%    f - a function_handler of a non-minimal solver for problem
%
% Options:
%    - ContinuationFactor: continuation factor for mu
%    - InitialMu: initial mu (default: auto)
%    - NoiseBound: inlier threshold
%    - MaxIterations: maximum number of iterations
%    - FixPriorsWeights: fix priors weights to 1 across all iterations
%    - CostThreshold: if weighted sum of squared residuals is below this value, GNC stops
%    - Debug: whether or not to enable the debug information
%
% Outputs:
%    inliers - the indices of the inliers
%    info - structure containing extended information about the xecution
%
% Example:
%   problem = linearRegressionProblem(100, 0.8);
%   [inliers, info] = gnc(problem, @leastSquareNorm2, ...
%       'NoiseBound', chi2inv(0.99, problem.dof)*problem.MeasurementNoiseStd^2);

% Author: Pasquale Antonante
% email: antonap@mit.edu
% Date: 2021-01-06

params = inputParser;
params.KeepUnmatched = true;
params.addParameter('NoiseBound', 0, @(x) isscalar(x) && x>=0 && isfinite(x));
params.parse(varargin{:});
assert(isProblem(problem), 'The problem doesn''t contain required fields.');

start_t = now;
if params.Results.NoiseBound > 0
    [inliers, info] = gnc_vanilla(problem, f, varargin{:});
else
    error('You need to set the Noise Bound')
end
end_t = now;
info.time = (end_t - start_t) * 24 * 60 * 60; % serial date number to sec
end

function [inliers, info] = gnc_vanilla(problem, f, varargin)
params = inputParser;
params.KeepUnmatched = true;
params.addParameter('ContinuationFactor', 1.4, @(x) isscalar(x) && x>0);
params.addParameter('InitialMu', 'auto', @(x) ischar(x) || (isscalar(x) && x>0));
params.addParameter('NoiseBound', 0, @(x) isscalar(x) && x>=0 && isfinite(x));
params.addParameter('MaxIterations', 1e3, @(x) isscalar(x));
params.addParameter('FixPriorsWeights', true, @(x) islogical(x));
params.addParameter('CostThreshold', 0, @(x) isscalar(x));
params.addParameter('Debug', false, @(x) islogical(x));
params.addParameter('init_', 0);
params.parse(varargin{:});

max_iterations = params.Results.MaxIterations;

if params.Results.Debug
    residuals_history = [];
    weights_history = [];
end

if ismember('init_', params.UsingDefaults)
    try
        [~, f_info] = f(problem);
    catch err
        fprintf('Error message: %s', err.message)
        error("Could not run the global solver")
    end
    assert(isfield(f_info, 'residuals'), 'f should compute residuals');
else
    f_info = params.Results.('init_');
end
barc2 = params.Results.NoiseBound;
weights = ones(1, problem.N);
if ischar(params.Results.InitialMu) && strcmpi(params.Results.InitialMu, 'auto')
    mu = 1 / (2 * max(f_info.residuals) / barc2 - 1 );
elseif isnumeric(params.Results.InitialMu)
    mu = params.Results.InitialMu;
else
    error('Invalid value for InitialMu')
end

prev_f_cost = 0;

i = 1;
info.stopping = 'MaxIterations'; % Worst case stopping
info.status = 1; % worst case status: failure
while i < max_iterations
    if params.Results.Debug
        residuals_history(i,:) = f_info.residuals(:)';
        weights_history(i, :) = weights;
    end
    weights = gncWeightsUpdate(weights, mu, f_info.residuals, barc2);
    if params.Results.FixPriorsWeights
        weights(problem.priors) = 1;
    end
    try
        [~, f_info] = f(problem, 'Weights', weights);
    catch err
        fprintf('Error message: %s', err.message)
        error("Could not run the global solver")
    end
    f_cost = sum(f_info.residuals(:) .* weights(:));
    cost_diff = abs(f_cost - prev_f_cost);
    prev_f_cost = f_cost;

    if (cost_diff < params.Results.CostThreshold) || areBinaryWeights(weights)
        if params.Results.Debug
            residuals_history(i+1,:) = f_info.residuals(:)';
            weights_history(i+1, :) = weights;
        end
        info.residuals = f_info.residuals(:)';
        info.stopping = 'CostThreshold';
        info.status = 1;
        break
    end
    mu = mu * params.Results.ContinuationFactor;
    i = i + 1;
end
inliers = find(weights>1-eps);
info.Iterations = i;
info.params = params.Results;
info.Algorithm = 'GNC';
if params.Results.Debug
    info.mu = mu;
    info.barc2History = repmat(barc2, i+1, 1);
    info.ResidualsHistory = residuals_history;
    info.WeightsHistory = weights_history;
    info.CostDiff = cost_diff;
    info.AreBinaryWeights = areBinaryWeights(weights);
end
end

author: MIT-SPARK

repo: GNC-and-ADAPT


text:
function [F,mptmissing] = filter_enumeration(F_xw,Zmodel,x,w,ops,uncertaintyTypes,separatedZmodel,VariableType)

mptmissing = 0;
if length(F_xw) == 0
    F = [];
    return;
else

    if any(Zmodel.K.q) | any(Zmodel.K.s)
        error('Only polytope uncertainty supported in duality based robustification');
    else
        if isempty(intersect(depends(F_xw),getvariables(w)))
            F = F_xw;
        elseif length(uncertaintyTypes)==1 & isequal(uncertaintyTypes{1},'inf-norm')

            if any(isinf((separatedZmodel{1}.lb))) | any(isinf(separatedZmodel{1}.ub))
                error('You have unbounded uncertain variables')
            else
                n = length(separatedZmodel{1}.lb);
                vertices = [];
                lb = separatedZmodel{1}.lb(:)';
                ub = separatedZmodel{1}.ub(:)';
                E = dec2bin(0:2^n-1,n)';
                E = double(E(:))-48;
                E = reshape(E,n,2^n);
                vertices = (repmat(lb(:),1,2^n) + E.*(repmat(ub(:),1,2^n)-repmat(lb(:),1,2^n)))';
                %for i = 0:2^n-1
                %    vertices = [vertices;lb+dec2decbin(i,n).*(ub-lb)];
                %end
                if ops.verbose
                    disp([' - Enumerated ' num2str(2^n) ' vertices'])
                end
                vertices = unique(vertices,'rows');
                if ops.verbose & 2^n > size(vertices,1)
                    disp([' - Reduced to ' num2str( size(vertices,1)) ' unique vertices'])
                end
                F = replaceVertices(F_xw,w,vertices',VariableType,ops);
            end

        elseif length(uncertaintyTypes)==1 & isequal(uncertaintyTypes{1},'simplex')

            k = abs(Zmodel.F_struc(1,1));
            n = length(w);

            vertices = zeros(n,1);
            for i = 1:n
                v = zeros(n,1);
                v(i) = k;
                vertices = [vertices v];
            end
            if ops.verbose
                disp([' - Enumerated ' num2str(n) ' vertices'])
            end
            vertices = pruneequalities(vertices,Zmodel);
            F = replaceVertices(F_xw,w,vertices,VariableType,ops);

        else
            % FIX : Assumes all uncertainty in all constraints
            K = Zmodel.K;
            A = -Zmodel.F_struc((1+K.f):(K.f + K.l),2:end);
            b =  Zmodel.F_struc((1+K.f):(K.f + K.l),1);

            % Some preprocessing to extract bounds from equality
            % constraints in order to make the uncertainty polytope
            % bounded (required since we are going to run vertex
            % enumeration)
            % We might have x>=0, sum(x)=1, and this code simply extracts
            % the implied bounds x<=1
            [lo,up] = find_lp_bounds(Zmodel.F_struc(1:K.f + K.l,:),K);
            Zmodel.lb = lo;Zmodel.ub = up;
            Zmodel = propagate_bounds_from_equalities(Zmodel);
            up = Zmodel.ub;
            lo = Zmodel.lb;
            upfi = find(~isinf(up));
            lofi = find(~isinf(lo));
            aux = Zmodel;
            aux.F_struc = [aux.F_struc;-lo(lofi) sparse(1:length(lofi),lofi,1,length(lofi),size(A,2))];
            aux.F_struc = [aux.F_struc;up(upfi) -sparse(1:length(upfi),upfi,1,length(upfi),size(A,2))] ;
            aux.K.l = aux.K.l + length(lofi) + length(upfi);
            K = aux.K;
            A = -aux.F_struc((1+K.f):(K.f + K.l),2:end);
            b =  aux.F_struc((1+K.f):(K.f + K.l),1);

            % Extract equalities Ex == f and project
            % Bunch of try-catch definsive as this assumes MPT
            % If fails we try internal crappy enumerator
            P = [];
            if any(K.f)
                f = aux.F_struc(1:K.f,1);
                E = -full(aux.F_struc(1:K.f,2:end));
                En = null(E);
                x0 = (E\f);
                % x = null(E)*z + x0, A*(null(E)*z + x0) <= b
                b = full(b-A*x0);
                A = full(A*En);
                try
                    P = polytope(A,b);
                catch
                end
            else
                x0 = zeros(size(A,2),1);
                En = eye(size(A,2));
                try
                    P = polytope(full(A),full(b));
                catch
                end
            end
            if ~isempty(P)
                try
                    vertices = extreme(P)';
                catch
                    error('The uncertainty space is unbounded (could be an artefact of YALMIPs modelling of nonolinear oeprators).')
                end
            else
                if ~all(b>0)
                    % Crappy enumeration requires feasible point
                    [x_c,R] = chebyball(A*sdpvar(size(A,2),1) <= b);
                else
                    x_c = zeros(size(A,2),1);
                end
                vertices = vertexenumerate(A,b,x_c);
                % This simply gives "nicer" models sometimes
                % First step in cleaning below too
                really_integer = abs(vertices(:)-round(vertices(:)))<eps;
                vertices(really_integer) = round(vertices(really_integer));
            end
            % Map back to original variables
            vertices = repmat(x0,1,size(vertices,2)) + En*vertices;

            % The vertex enumeration was done without any equality constraints.
            % We know check all vertices so see if they satisfy equalities.
            vertices = pruneequalities(vertices,Zmodel);
            if ops.verbose
                disp([' - Enumerated ' num2str(size(vertices,2)) ' vertices'])
            end
            F = replaceVertices(F_xw,w,vertices,VariableType,ops);
        end
    end
end

function F = replaceVertices(F_xw,w,vertices,VariableType,ops)
% Doing LP constraints in a vectorized manner saves a lot of time

F_xw_lp = F_xw(find(is(F_xw,'elementwise')));
F_xw_socp_sdp = F_xw -  F_xw_lp;
F = ([]);

x_Flp = depends(F_xw_lp);
uncAux = yalmip('auxvariablesW');
uncAux = recover(intersect(x_Flp,VariableType.aux_with_w_dependence));
if isequal(ops.robust.auxreduce,'none')
    uncAux = [];
end

if length(F_xw_lp)>0
    rLP = [];
    if ~isempty(uncAux)
        z = sdpvar(repmat(length(uncAux),1,size(vertices,2)),repmat(1,1,size(vertices,2)),'full');
    end
    for i = 1:size(vertices,2)
        temp = replace(sdpvar(F_xw_lp),w,vertices(:,i),0);
        if ~isempty(uncAux)
            temp = replace(temp,uncAux,z{i});
        end
        rLP = [rLP;temp];
    end

    % FIXME: More general detection of silly constraints
    if isa(rLP,'double') & all(rLP>=-eps^0.75)
        F = ([]);
    else
        % Easily generates redundant constraints
        [aux,index] = uniquesafe(getbase(rLP),'rows');
        try
            F = (rLP(index(randperm(length(index)))) >= 0);
        catch
            1
        end
    end
end

% Remaining conic stuff
for j = 1:length(F_xw_socp_sdp)
    for i = 1:size(vertices,2)
        temp = replace(F_xw_socp_sdp(j),w,vertices(:,i),0);
        if ~isempty(uncAux)
            temp = replace(temp,uncAux,z{i});
        end
        F = F + lmi(temp);
    end
end

function vertices = pruneequalities(vertices,Zmodel)
K = Zmodel.K;
% The vertex enumeration was done without any equality constraints.
% We know check all vertices so see if they satisfy equalities.
if any(K.f)
    Aeq = -Zmodel.F_struc(1:K.f,2:end);
    beq =  Zmodel.F_struc(1:K.f,1);
    feasible = sum(abs(Aeq*vertices - repmat(beq,1,size(vertices,2))),1) < 1e-6;
    vertices = vertices(:,feasible);
    if isempty(feasible)
        error('The uncertainty space is infeasible.')
    end
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [struct_irf_record,D_draws,gamma_draws,ETA_record,It,Bu,beta_draws,sigma_draws]=...
    IRFt5_Bayesian(names,betahat,m,n,X,Y,k,p,enddate,startdate,IRFperiods,IRFt,T,arvar,q,It,Bu,lambda1,lambda3,lambda4,pref,strctident)
%%Implementation of a Bayesian Proxy VAR based on the Codes published by Caldara and Herbst (2018)
%%Implementation by Ben Schumann
%% IV identification
%Load IV and make it comparable with the reduced form errors
[EPSIV,IVcut,~,sigmahatIV,sigma_hat,inv_sigma_hat,IV,txt,~,cut1,cut2,cut3,cut4]...
    =bear.loadIV(betahat,k,n,Y,X,T,p,names,startdate,enddate,strctident,pref);

% IV routine
[beta_draws,sigma_draws,~,~,D_draws,gamma_draws,irf_storage,ETA_storage,It,Bu]=...
    bear.irfIV_MH(EPSIV,IVcut,betahat,sigmahatIV,sigma_hat,inv_sigma_hat,IV,txt,cut1,cut2,cut3,cut4,names,It,Bu,n,arvar,lambda1,lambda3,lambda4,m,p,k,q,X,Y,T,startdate,enddate,pref,strctident,IRFperiods,IRFt);

% reorganise
Acc=It-Bu;
% loop over iterations
for ii=1:Acc/strctident.Thin
    for jj=1:IRFperiods
        % loop over variables
        for kk=1:n
            % loop over shocks (only one)
            for ll=1
                %for ll=1:n
                struct_irf_record{kk,ll}(ii,jj)=irf_storage{ii,1}(kk,ll,jj);
            end
        end
    end
end

% loop over variables
for kk=1:n
    % loop over shocks (only one)
    for ll=2:n
        %for ll=1:n
        struct_irf_record{kk,ll}=zeros(Acc/strctident.Thin,IRFperiods);
    end
end

%reorganize Structural Shocks
ETA_record=cell(n,1);
for jj=1:Acc/strctident.Thin
    for kk=1:n
        ETA_record{kk,1}(jj,:)= ETA_storage{jj,1}(kk,:);
    end
end

%finally update It and Bu such that BEAR knows how many draws we kept
It = It/strctident.Thin;
Bu = Bu/strctident.Thin;

author: european-central-bank

repo: BEAR-toolbox


text:
function model  = modelSetOutputWeights(model, W, b)

% MODELSETOUTPUTWEIGHTS Wrapper function to return set output weight and bias matrices.
% FORMAT
% DESC sets the output weight and bias matrices for any mapping model
% that can be interpreted as a generalised linear model (e.g. rbf
% networks, kernel based regressions, multi layer perceptrons, linear).
% RETURN model : the model with updated weights and bias matrices.
% ARG model : the mapping model.
% ARG W : the output weight matrix.
% ARG b : the output biases.
%
% SEEALSO : mlpCreate, rbfCreate, kbrCreate, linearCreate
%
% COPYRIGHT : Neil D. Lawrence, 2008

% MLTOOLS

switch model.type
 case 'mlp'
  model.w2 = W;
  model.b2 = b;
 case 'rbf'
  model.w2 = W;
  model.b2 = b;
 case 'kbr'
  model.A = W;
  model.bias = b;
 case 'linear'
  model.W = W;
  model.b = b;
 otherwise
  error('Model has no implementation of output weights and biases.');
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function varargout = armorf(x,Nr,Nl,p)
%ARMORF   AR parameter estimation via LWR method by Morf modified.
%   x is a matrix whose every row is one variable's time series
%   Nr is the number of realizations, Nl is the length of every realization
%   If the time series are stationary long, just let Nr=1, Nl=length(x)
%   p is the order of AR model
%
%   A = ARMORF(X,NR,NL,P) returns the polynomial coefficients A corresponding to
%     the AR model estimate of matrix X using Morf's method.
%
%   [A,E] = ARMORF(...) returns the final prediction error E (the
%   covariance matrix of the white noise of the AR model).
%
%   [A,E,K] = ARMORF(...) returns the vector K of reflection
%     coefficients (parcor coefficients).
%
%   Ref: M. Morf, etal, Recursive Multichannel Maximum Entropy Spectral Estimation,
%              IEEE trans. GeoSci. Elec., 1978, Vol.GE-16, No.2, pp85-94.
%        S. Haykin, Nonlinear Methods of Spectral Analysis, 2nd Ed.
%              Springer-Verlag, 1983, Chapter 2
%
%   finished on Aug.9, 2002 by Yonghong Chen

% Initialization
[L,N]=size(x);
R0=zeros(L,L);
R0f=R0;
R0b=R0;
pf=R0;
pb=R0;
pfb=R0;
ap(:,:,1)=R0;
bp(:,:,1)=R0;
En=R0;
for i=1:Nr
    En=En+x(:,(i-1)*Nl+1:i*Nl)*x(:,(i-1)*Nl+1:i*Nl)';
    ap(:,:,1)=ap(:,:,1)+x(:,(i-1)*Nl+2:i*Nl)*x(:,(i-1)*Nl+2:i*Nl)';
    bp(:,:,1)=bp(:,:,1)+x(:,(i-1)*Nl+1:i*Nl-1)*x(:,(i-1)*Nl+1:i*Nl-1)';
end
ap(:,:,1) = inv((chol(ap(:,:,1)/Nr*(Nl-1)))');
bp(:,:,1) = inv((chol(bp(:,:,1)/Nr*(Nl-1)))');
for i=1:Nr
    efp = ap(:,:,1)*x(:,(i-1)*Nl+2:i*Nl);
    ebp = bp(:,:,1)*x(:,(i-1)*Nl+1:i*Nl-1);
    pf = pf + efp*efp';
    pb = pb + ebp*ebp';
    pfb = pfb + efp*ebp';
end
En = chol(En/N)'; % Covariance of the noise

% Initial output variables
coeff = []; %  Coefficient matrices of the AR model
kr=[];  % reflection coefficients

for m=1:p
   % Calculate the next order reflection (parcor) coefficient
   ck = inv((chol(pf))')*pfb*inv(chol(pb));
   kr=[kr,ck];
   % Update the forward and backward prediction errors
   ef = eye(L)- ck*ck';
   eb = eye(L)- ck'*ck;

   % Update the prediction error
   En = En*chol(ef)';
   E = (ef+eb)./2;

   % Update the coefficients of the forward and backward prediction errors
   ap(:,:,m+1) = zeros(L);
   bp(:,:,m+1) = zeros(L);
   pf = zeros(L);
   pb = zeros(L);
   pfb = zeros(L);

   for i=1:m+1
       a(:,:,i) = inv((chol(ef))')*(ap(:,:,i)-ck*bp(:,:,m+2-i));
       b(:,:,i) = inv((chol(eb))')*(bp(:,:,i)-ck'*ap(:,:,m+2-i));
   end
   for k=1:Nr
       efp = zeros(L,Nl-m-1);
       ebp = zeros(L,Nl-m-1);
       for i=1:m+1
           k1=m+2-i+(k-1)*Nl+1;
           k2=Nl-i+1+(k-1)*Nl;
           efp = efp+a(:,:,i)*x(:,k1:k2);
           ebp = ebp+b(:,:,m+2-i)*x(:,k1-1:k2-1);
       end
       pf = pf + efp*efp';
       pb = pb + ebp*ebp';
       pfb = pfb + efp*ebp';
   end
   ap = a;
   bp = b;
end
for j=1:p
    coeff = [coeff,inv(a(:,:,1))*a(:,:,j+1)];
end

varargout{1} = coeff;
if nargout >= 2
    varargout{2} = En*En';
end
if nargout >= 3
    varargout{3} = kr;
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function info=hdr_read_header(filename)
% function for reading header of HDR/IMG Analyze ( .hdr ) volume file
%
% info = hdr_read_header(filename);
%
% examples:
% 1,  info=hdr_read_header()
% 2,  info=hdr_read_header('volume.hdr');

if(exist('filename','var')==0)
    [filename, pathname] = uigetfile('*.hdr', 'Read hdr-file');
    filename = [pathname filename];
end
if(exist('analyze75info','file')>0)
    info=analyze75info(filename);
else
    info=get_info_analyze_hdr(filename);
end

function info=get_info_analyze_hdr(filename)
% Strings for data type (compatible names with analyze75info)
strImgData{1}.ImgDataType = 'DT_UNKNOWN';
strImgData{1}.ColorType = 'unknown';
strImgData{2}.ImgDataType = 'DT_BINARY';
strImgData{2}.ColorType = 'grayscale';
strImgData{3}.ImgDataType = 'DT_UNSIGNED_CHAR';
strImgData{3}.ColorType = 'grayscale';
strImgData{5}.ImgDataType = 'DT_SIGNED_SHORT';
strImgData{5}.ColorType = 'grayscale';
strImgData{9}.ImgDataType = 'DT_SIGNED_INT';
strImgData{9}.ColorType = 'grayscale';
strImgData{17}.ImgDataType = 'DT_FLOAT';
strImgData{17}.ColorType = 'grayscale';
strImgData{33}.ImgDataType = 'DT_COMPLEX';
strImgData{33}.ColorType = 'grayscale';
strImgData{65}.ImgDataType = 'DT_DOUBLE';
strImgData{65}.ColorType = 'grayscale';
strImgData{129}.ImgDataType = 'DT_RGB';
strImgData{129}.ColorType = 'truecolor';
strImgData{256}.ImgDataType = 'DT_ALL';
strImgData{256}.ColorType = 'unknown';
strOriData{1} = 'Transverse unflipped';
strOriData{2} = 'Coronal unflipped';
strOriData{3} = 'Sagittal unflipped';
strOriData{4} = 'Transverse flipped';
strOriData{5} = 'Coronal flipped';
strOriData{6} = 'Sagittal flipped';
strOriData{7} = 'Orientation unavailable';

% Open the HDR-file, change MachineFormat if not the right
% header size
info.Filename=filename;
fid = fopen(filename,'rb','l'); fseek(fid,0,'bof');
info.ByteOrder='ieee-le';
info.HdrFileSize = fread(fid, 1,'int32');
if(info.HdrFileSize>2000);
    fclose(fid);
    fid = fopen(filename,'rb','b'); fseek(fid,0,'bof');
    info.ByteOrder='ieee-be';
    info.HdrFileSize = fread(fid, 1,'int32');
end

% Read the Whole Analyze Header
info.Format='Analyze';
info.HdrDataType = fread(fid,10,'char=>char')';
info.DatabaseName = fread(fid,18,'char=>char')';
info.Extents = fread(fid, 1,'int32');
info.SessionError = fread(fid, 1,'int16');
info.Regular = fread(fid, 1,'char=>char')';
unused= fread(fid, 1,'uint8')';
dim = fread(fid,8,'int16')';
if (dim(1) < 3), dim(1) = 4; end
info.Dimensions=dim(2:dim(1)+1);
info.Width=info.Dimensions(1);
info.Height=info.Dimensions(2);
info.VoxelUnits = fread(fid,4,'char=>char')';
info.CalibrationUnits = fread(fid,8,'char=>char')';
unused = fread(fid,1,'int16');
ImgDataType = fread(fid,1,'int16');
info.ImgDataType=strImgData{ImgDataType+1}.ImgDataType;
info.ColorType=strImgData{ImgDataType+1}.ColorType;
info.BitDepth = fread(fid,1,'int16');
unused = fread(fid,1,'int16');
PixelDimensions = fread(fid,8,'float')';
info.PixelDimensions=PixelDimensions(2:length(info.Dimensions)+1);
info.VoxelOffset = fread(fid,1,'float');
info.RoiScale = fread(fid,1,'float');
unused = fread(fid,1,'float');
unused = fread(fid,1,'float');
info.CalibrationMax = fread(fid,1,'float');
info.CalibrationMin = fread(fid,1,'float');
info.Compressed = fread(fid,1,'int32');
info.Verified = fread(fid,1,'int32');
info.GlobalMax = fread(fid,1,'int32');
info.GlobalMin = fread(fid,1,'int32');
info.Descriptor = fread(fid,80,'char=>char')';
info.AuxFile = fread(fid,24,'char=>char')';
Orientationt = fread(fid, 1,'uint8');
if((Orientationt>=48)&&(Orientationt<=53)), Orientationt = Orientationt -48; end
info.Orientationt=strOriData{min(Orientationt,5)+1};
info.Originator = fread(fid,10,'char=>char')';
info.Generated = fread(fid,10,'char=>char')';
info.Scannumber = fread(fid,10,'char=>char')';
info.PatientID = fread(fid,10,'char=>char')';
info.ExposureDate = fread(fid,10,'char=>char')';
info.ExposureTime = fread(fid,10,'char=>char')';
unused = fread(fid, 3,'char=>char')';
info.Views = fread(fid, 1,'int32');
info.VolumesAdded = fread(fid, 1,'int32');
info.StartField = fread(fid, 1,'int32');
info.FieldSkip = fread(fid, 1,'int32');
info.Omax = fread(fid, 1,'int32');
info.Omin = fread(fid, 1,'int32');
info.Smax = fread(fid, 1,'int32');
info.Smin = fread(fid, 1,'int32');
fclose(fid);

% Remove empty string parts from the info struct
names = fieldnames(info);
for i=1:length(names)
    value=info.(names{i});
    if(ischar(value))
        value(uint8(value)==0)=[];
        if(isempty(value)), info.(names{i})=''; else info.(names{i})=value; end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
%
% Smooth point-set registration method using neighboring constraints
% -------------------------------------------------------------------
%
% Authors: Gerard Sanromà, René Alquézar and Francesc Serratosa
%
% Contact: gsanorma@gmail.com
% Date: 15/02/2012
%
% Demo of the non-rigid method
%
% For the execution of the non-rigid version you will need the 'ctps_gen'
% and 'ctps_warp_pts' functions from the TPS-RPM implementation by
% Haili Chui and Anand Rangarajan ("A new algorithm for non-rigid point
% matching",IEEE CVPR 2000) available from:
% http://www.cise.ufl.edu/~anand/students/chui/rpm/TPS-RPM.zip
%
% For the correct execution of this demo you will need the synthetic
% datasets by Dr. Haili Chui and Prof. Anand Rangarajan, available at:
% http://www.umiacs.umd.edu/~zhengyf/PointMatchDemo/DataChui.zip
%

clear variables;

if exist('ctps_gen') ~= 2 || exist('ctps_warp_pts') ~= 2
    disp('The MATLAB m-files (i.e., functions) "ctps_gen" and "ctps_warp_pts" are not available');
    disp('You can get them from the TPS-RPM implementation by Haili Chui and Anand Rangarajan');
    disp('"A new algorithm for non-rigid point matching",IEEE CVPR 2000');
    disp('http://www.cise.ufl.edu/~anand/students/chui/rpm/TPS-RPM.zip');
    return;
end

if exist('DataChui','dir') == 0
    disp('Fish and Chinese character datasets are not available');
    disp('You can download them from:')
    disp('http://www.umiacs.umd.edu/~zhengyf/PointMatchDemo/DataChui.zip');
    return;
end

% Type of experiments (uncomment the appropriate one)
exp_type = 'fish_def';
% exp_type = 'fish_noise';
% exp_type = 'fish_outlier';
% exp_type = 'chinese_def';
% exp_type = 'chinese_noise';
% exp_type = 'chinese_outlier';

% Path to the Chui's dataset
path = 'DataChui';
file_pattern = strcat('save_',exp_type,'_');

sample_range = 1:100;
if ~isempty(strfind(exp_type,'noise'))
    perturb_range = 1:6;
else
    perturb_range = 1:5;
end

E = inf(length(perturb_range),length(sample_range));

for perturb = perturb_range
    for sample = sample_range

        disp(['perturb = ',num2str(perturb),' sample = ',num2str(sample)]);

        % Read sample
        load(fullfile(path,strcat(file_pattern,num2str(perturb),'_',num2str(sample))),...
            'x1','y2a');
        gM = [-x1(:,2) x1(:,1)];
        gD = [-y2a(:,2) y2a(:,1)];
        lD = length(gD);
        lM = length(gM);

        if ~isempty(strfind(exp_type,'outlier'))
            vgM = match_chui_outl(gM,gD);
        else
            vgM = match_chui(gM,gD);
        end

        % gm
        r = gD(1:lM,:) - vgM;
        E(perturb,sample) = mean(sqrt(diag(r*r')));
        disp(['ERROR=',num2str(E(perturb,sample))]);

    end
end

save(strcat('nonrigid_',exp_type),'E');

author: Sable

repo: mcbench-benchmarks


text:
function area = sphere_unit_area_nd ( dim_num )

%*****************************************************************************80
%
%% SPHERE_UNIT_AREA_ND computes the surface area of a unit sphere in ND.
%
%  Discussion:
%
%    The unit sphere in ND satisfies:
%
%      sum ( 1 <= I <= DIM_NUM ) X(I) * X(I) = 1
%
%    Results for the first few values of N are:
%
%    DIM_NUM   Area
%
%     2    2        * PI
%     3    4        * PI
%     4  ( 2 /   1) * PI**2
%     5  ( 8 /   3) * PI**2
%     6  ( 1 /   1) * PI**3
%     7  (16 /  15) * PI**3
%     8  ( 1 /   3) * PI**4
%     9  (32 / 105) * PI**4
%    10  ( 1 /  12) * PI**5
%
%    For the unit sphere, Area(DIM_NUM) = DIM_NUM * Volume(DIM_NUM)
%
%    Sphere_Unit_Area ( DIM_NUM ) = 2 * PI**(DIM_NUM/2) / Gamma ( DIM_NUM / 2 )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 February 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the dimension of the space.
%
%    Output, real SPHERE_UNIT_AREA_ND, the area of the sphere.
%
  if ( mod ( dim_num, 2 ) == 0 )
    m = floor ( dim_num / 2 );
    area = 2.0 * pi^m;
    for i = 1 : m-1
      area = area / i;
    end
  else
    m = floor ( ( dim_num - 1 ) / 2 );
    area = pi^m * 2.0^dim_num;
    for i = m+1 : 2*m
      area = area / i;
    end
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%DEMO_HIERPRIOR  Demonstration of hierarchical prior structures for
%                  Gaussian processes
%
%  Description
%    This demonstration shows how to set hierarchical prior structures for
%    parameters of covariance function and/or likelihood function. This
%    demo is intended for showing how to set hierarchical prior structures
%    and for other testing purposes, not as a demonstration of model
%    selection procedures or comparisons. We test the hierarchical priors
%    for both regression and classification with Laplace and EP.
%
%
%  See also DEMO_REGRESSION_ROBUST, DEMO_SPATIAL1, DEMO_*
%
% Copyright (c) 2014 Ville Tolvanen

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

%% Regression
% load the data. First 100 variables are for training
% and last 100 for test

S = which('demo_hierprior');
L = strrep(S,'demo_hierprior.m','demodata/odata.txt');
x = load(L);
y = x(1:100,2);
x = x(1:100,1);
[n, nin] = size(x);

% Test data
xt = [-2.7:0.01:2.7]';
yt = 0.3+0.4*xt+0.5*sin(2.7*xt)+1.1./(1+xt.^2);

% Create hierarchial prior
pl=prior_t('mu_prior',prior_t);

% Create the model
gpcf = gpcf_sexp('magnSigma2_prior', pl, 'lengthScale_prior', pl);
gpcf2 = gpcf_linear('coeffSigma2_prior', pl);
% gpcf2 = gpcf_neuralnetwork('weightSigma2_prior', pl, 'biasSigma2_prior', pl);
lik = lik_gaussian('sigma2_prior', pl);
gp = gp_set('lik', lik, 'cf', {gpcf gpcf2}, 'jitterSigma2',1e-9);

% Parameters
[w,s]=gp_pak(gp)

% Optimize the parameters
opt=optimset('TolX',1e-5,'TolFun',1e-5,'Display','iter');
gp = gp_optim(gp,x,y,'opt',opt);

derivs={};
derivs{1}=derivativecheck(gp_pak(gp), @(ww) gp_eg(ww, gp, x, y));

[Ef,Varf,lpyt,Ey,Vary]=gp_pred(gp,x,y,xt,'yt',yt);

%% Regression with sparse approximations
% Here we optimize all the model parameters, including the inducing inputs

X_u=linspace(-2,3, 10)';
gp = gp_set('lik', lik, 'cf', {gpcf gpcf2}, 'jitterSigma2',1e-6, 'type', 'FIC', ...
  'X_u', X_u, 'infer_params', 'covariance+likelihood+inducing', ...
  'Xu_prior', prior_gaussian('mu_prior', prior_t()));
gp = gp_optim(gp,x,y,'opt',opt);
derivs{2}(:,1)=derivativecheck(gp_pak(gp), @(ww) gp_eg(ww, gp, x, y));
Ef2=gp_pred(gp,x,y,xt,'yt',yt);

gp = gp_set('lik', lik, 'cf', {gpcf gpcf2}, 'jitterSigma2',1e-6, 'type', 'VAR', ...
  'X_u', X_u, 'infer_params', 'covariance+likelihood+inducing', ...
  'Xu_prior', prior_gaussian('mu_prior', prior_t()));
gp = gp_optim(gp,x,y,'opt',opt);
derivs{2}(:,2)=derivativecheck(gp_pak(gp), @(ww) gp_eg(ww, gp, x, y));
Ef3=gp_pred(gp,x,y,xt,'yt',yt);

gp = gp_set('lik', lik, 'cf', {gpcf gpcf2}, 'jitterSigma2',1e-6, 'type', 'SOR', ...
  'X_u', X_u, 'infer_params', 'covariance+likelihood+inducing', ...
  'Xu_prior', prior_gaussian('mu_prior', prior_t()));
gp = gp_optim(gp,x,y,'opt',opt);
derivs{2}(:,3)=derivativecheck(gp_pak(gp), @(ww) gp_eg(ww, gp, x, y));
Ef4=gp_pred(gp,x,y,xt,'yt',yt);

%% Plot

plot(x,y,'.', xt, Ey, xt, Ef2, xt, Ef3, xt, Ef4);
legend('Data', 'Full-GP', 'FIC', 'VAR', 'SOR');

%% Classification with Laplace & EP (full and sparse approximation)
S = which('demo_spatial1');
data = load(strrep(S,'demo_spatial1.m','demodata/spatial1.txt'));

x = data(:,1:2);
ye = data(:,3);
y = data(:,4);
% Remove some of the data
inds=find(x(:,2)>30);
x=x(inds,:);
y=y(inds,:);
ye=ye(inds,:);

dims = [30    60     1    35];
[trindex, Xu] = set_PIC(x, dims, 7, 'corners', 0);
[n,nin] = size(x);

% Create the covariance functions
pl = prior_t('s2',10, 's2_prior', prior_t());
pm = prior_sqrtunif();
pxu={};
for i1=1:size(Xu,1)
  pxu={pxu{:} pl};
end
gpcf1 = gpcf_matern32('lengthScale', 1, 'magnSigma2', 0.03, ...
  'lengthScale_prior', pl, 'magnSigma2_prior', prior_fixed());
gpcf2 = gpcf_ppcs3('nin',nin,'lengthScale', 5, 'magnSigma2', 0.05,...
  'lengthScale_prior', pl, 'magnSigma2_prior', prior_fixed());

lik = lik_negbin();

% Create the GP structures (Laplace)
gp = gp_set('lik', lik, 'cf', gpcf1,'jitterSigma2', 1e-4, ...
            'infer_params', 'covariance+likelihood');
gp2 = gp_set('type', 'FIC', 'lik', lik, 'cf', gpcf1, 'X_u', Xu, 'Xu_prior', pxu, ...
            'jitterSigma2', 1e-4, 'infer_params', 'covariance+likelihood+inducing');
gp3 = gp_set('type', 'PIC', 'lik', lik, 'cf', gpcf1, 'X_u', Xu, 'Xu_prior', pxu,  ...
           'jitterSigma2', 1e-4, 'infer_params', 'covariance+likelihood+inducing',...
           'tr_index', trindex);
gp4 = gp_set('type', 'CS+FIC', 'lik', lik, 'cf', {gpcf1 gpcf2}, 'X_u', Xu, 'Xu_prior', pxu,  ...
           'jitterSigma2', 1e-4, 'infer_params', 'covariance+likelihood+inducing');
gp4.cf{1}.p.lengthScale=prior_fixed();
gp4.cf{1}.p.magnSigma2=prior_fixed();

opt=optimset('TolX',1e-3,'TolFun',1e-3,'Display','iter');
fprintf('Laplace \n \n');
fprintf('Full GP:\n \n')
derivs{3}(:,1)=derivativecheck(gp_pak(gp), @(ww) gp_eg(ww, gp, x, y, 'z', ye));
fprintf('FIC:\n \n')
derivs{4}(:,1)=derivativecheck(gp_pak(gp2), @(ww) gp_eg(ww, gp2, x, y, 'z', ye));
fprintf('PIC:\n \n')
derivs{4}(:,2)=derivativecheck(gp_pak(gp3), @(ww) gp_eg(ww, gp3, x, y, 'z', ye));
fprintf('CS+FIC:\n \n')
derivs{4}(:,3)=derivativecheck(gp_pak(gp4), @(ww) gp_eg(ww, gp4, x, y, 'z', ye));

% Create the GP structures (EP)
gp = gp_set('lik', lik, 'cf', gpcf1,'jitterSigma2', 1e-4, ...
            'infer_params', 'covariance+likelihood', ...
            'latent_method', 'EP');
gp2 = gp_set('type', 'FIC', 'lik', lik, 'cf', gpcf1, 'X_u', Xu,  'Xu_prior', pxu, ...
            'jitterSigma2', 1e-4, 'infer_params', 'covariance+likelihood+inducing', ...
            'latent_method', 'EP');
gp3 = gp_set('type', 'PIC', 'lik', lik, 'cf', gpcf1, 'X_u', Xu, 'Xu_prior', pxu,  ...
           'jitterSigma2', 1e-4, 'infer_params', 'covariance+likelihood+inducing',...
           'tr_index', trindex, 'latent_method', 'EP');
gp4 = gp_set('type', 'CS+FIC', 'lik', lik, 'cf', {gpcf1 gpcf2}, 'X_u', Xu,  'Xu_prior', pxu, ...
           'jitterSigma2', 1e-4, 'infer_params', 'covariance+likelihood+inducing', ...
            'latent_method', 'EP');
gp4.cf{1}.p.lengthScale=prior_fixed();
gp4.cf{1}.p.magnSigma2=prior_fixed();

fprintf('EP \n \n');
fprintf('Full GP:\n \n')
derivs{5}(:,1)=derivativecheck(gp_pak(gp), @(ww) gp_eg(ww, gp, x, y, 'z', ye));
fprintf('FIC:\n \n')
derivs{6}(:,1)=derivativecheck(gp_pak(gp2), @(ww) gp_eg(ww, gp2, x, y, 'z', ye));
fprintf('PIC:\n \n')
derivs{6}(:,2)=derivativecheck(gp_pak(gp3), @(ww) gp_eg(ww, gp3, x, y, 'z', ye));
fprintf('CS+FIC:\n \n')
derivs{6}(:,3)=derivativecheck(gp_pak(gp4), @(ww) gp_eg(ww, gp4, x, y, 'z', ye));

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
function [scoreF,score0] = argmax_fit_type(Minit,lib,list_sid,verbose)
% ARGMAX_FIT_TYPE... Fit a parse to an image, with type-level parameters
% only
%
% Input
%   Minit : instance of MotorProgram class
%   list_sid: (default=all) list of stroke-ids which we want to optimize
%    while the others are frozen
%   verbose: (true/false) describe algorithm progress
%
% Output:
%  scoreF : final score
%  score0 : initial score
%

    M = Minit.copy();
    if ~exist('list_sid','var')
       list_sid = 1:M.ns;
    end
    if ~exist('verbose','var')
       verbose = false;
    end

    has_relations = M.has_relations(list_sid);
    if has_relations
        fmin = @(theta) myscore_HasRel(theta,M,lib,list_sid);
    else
        all_R = cache_enum_all_relations(lib,M);
        fmin = @(theta) myscore_NoRel(theta,M,lib,list_sid,all_R);
    end

    % Set-up the optimization problem
    [theta0,lb,ub] = model_to_vec_fit_type(M,list_sid);
    options = optimset('Display','off');
    if verbose
        options = optimset(options,'Display','iter');
    end
    options = optimset(options,'TolFun',1e-4,'Algorithm','active-set');

    % Run the optimization
    %tic
    score0 = -fmin(theta0);
    try
        thetaF = fmincon(fmin,theta0,[],[],[],[],lb,ub,[],options);
    catch
        fprintf(1,'*Warning*: optimization failed. Error caught.\n');
        thetaF = theta0;
    end
    scoreF = -fmin(thetaF);
    %t = toc;

    refill(thetaF,Minit,list_sid);

    % make sure we have not added/removed relations
    assert(Minit.has_relations==has_relations);

end

function refill(theta,M,list_sid)
    vec_to_model_fit_type(theta,M,list_sid);
end

% fill the MotorProgram with parameters
% and then score
function minscore = myscore_HasRel(theta,M,lib,list_sid)
    Q = M.copy(); % we don't want to modify the shared MotoProgram base
    refill(theta,Q,list_sid);
    ll = scoreMP(Q,lib,'strokes',list_sid,'type',true,'token',true,'stat',false,'image',true);
    minscore = -ll;
end

% fill the MotorProgram with parameters
% and then optimize the relations.
% Finally, return the score "minscore"
function minscore = myscore_NoRel(theta,M,lib,list_sid,all_R)
    Q = M.copy(); % we don't want to modify the shared MotoProgram base
    refill(theta,Q,list_sid);
    argmax_relations(lib,Q,all_R,list_sid);
    ll = scoreMP(Q,lib,'strokes',list_sid,'type',true,'token',true,'stat',false,'image',true);
    minscore = -ll;
end

author: brendenlake

repo: BPL

save_path: github-repos/MATLAB/brendenlake-BPL


text:
% ts_plot for PS_PLOT function
%
%
%
%   Andy Hooper, June 2006
%
%   ======================================================================
%   11/2010 MMC & MA: plotting time series using 'ts' option
%   03/2011 AH convert radius to m, remove line fitting and add mean
%   08/2016 AH Replace break with return
%   ======================================================================

% Place button for new TS plots
% see ps_plot

mEditBox=findobj('Background',[1 1 1]); % get the handle of editable textbox
%radiusfactor=str2num(get(mEditBox,'String')) % get radius factor from editbox
radiusfactor = str2num(char(get(mEditBox,'String')));					% added by david this is also not a single value but a vector
radiusfactor = radiusfactor(1);								% select only the first value, i checked it for a couple of values and radius is correct

%if radiusfactor > 3600
%    disp('radius factor should be <= 3600')
%    return
%end

momfig_name=['(', get(gcf,'name'), ')']; % inherent fig_name from the velocity plot

% LOAD TS mat file ~ ps_plot_ts_v-d.mat
if ~exist('ph_mm','var')
  %fid = fopen(savetxt);
  fid = fopen('ps_plot_ts_matname.txt');
  tsmat=textscan(fid,'%s'); % get mat filename to load parameters
  fclose(fid);
  clear fid
  eval(['load ' tsmat{1}{1}])% load saved matrix
end
% end of load

%
% GET USER INPUT from MOUSE CLICK
%
 disp('Please select a point on the figure to plot time series (TS)')
 [lon0,lat0] = ginput(1); % turn on when final
 disp(['Selected point coordinates (lon,lat):' num2str(lon0),', ', num2str(lat0) ])

% MAKE A CIRCLE AROUND SELECTED POINT (lon0,lat0)
t = linspace(0,2*pi,114); % 114 pts
h=lon0; k=lat0;  % center cn for the circle
r=radiusfactor;
                   % change radius if you want to include more points
x = r*cos(t);
y = r*sin(t);

% SELECT POINTS based on LONLAT
xy=llh2local(lonlat',[lon0,lat0])'*1000;
in = inpolygon(xy(:,1),xy(:,2),x,y);
n_pts_near=sum(in);  % how many ps found

if sum(in) == 0
    disp(['No points found in radius of ', num2str(r) ] )
    disp('Please make new selection increasing radius factor...')
    return     % no pts selecte
end

% PLOT selection and points
figure
    set(gcf,'name',['Found #pt(s): ', num2str(n_pts_near),...
        '  in radius: ', num2str(r), ' m. ', momfig_name ])
    circlell=local2llh([x;y]/1000,[lon0,lat0]);
    plot(circlell(1,:),circlell(2,:));
    hold on
    plot(lon0,lat0,'*')

    lon2=lonlat(in,1);
    lat2=lonlat(in,2);
    plot(lon2,lat2,'dr')
    axis image; hold off

% if ps>1 than average
% [dist,az] = distance(lat0,lon0,lat2,lon2); % or use llh2local

% plot closest point or avg of multiple points.

% phases
% ph_mm holds corrected phase
% ph_all holds radians to meters
% v_all=-m(2,:)'; % ?

% PLOT TS for given point(s)
ts=ph_mm(in,:);
G=[ones(size(day)),day-master_day] ; % [ 1  a ] --> b + ax

offset=pi*1000*lambda/(4*pi);

ts=nanmean(ts,1);
x_hat=G\double(ts');
ts_hat=G*x_hat;
%tsup_hat=ts_hat+offset;
%tslo_hat=ts_hat-offset;

%%% Typical TS plot - no auxilary subplots
% figure
%     set(gcf,'name',[ ' Times series plot for #point(s): ',...
%         num2str(n_pts_near), ' ', momfig_name])
%     plot(day,ts,'--*'); hold on
%     plot(day,ts_hat,'-*r');
%     plot(day,tsup_hat,'-.g');
%     plot(day,tslo_hat,'-.g');
%     hold off
%     grid on
%     ylabel('mm');
%     xlabel('Time [mmmyy]')
%     %datetick('x','mmmyy')  % keepticksdoc
%     %set(gca, 'XTick',day);
%     set(gca, 'XTickLabel', datestr(day,'mmmyy'));

%%% Enhanced TS plot
 figure % main figure
   orient landscape
   % Bperp
%    subplot(10,1,1)     % Bperp
%      bperp(find(bperp==0))=[]; % drop master.
%      bar(bperp)
%      ylabel('Bperp [m]')
%    grid on
   % TS
%    subplot(10,10,[11 87]) % subplot(10,1,2:9)
       set(gcf,'name',[ ' Times series plot for #point(s): ',...
        num2str(n_pts_near), ' ', momfig_name])
    %h1=plot(day,ts,'--*'); hold on
    %plot(day,ts_hat,'-k','LineSmoothing','on'); % mess up ticks
    plot(day,ts_hat,'-r'); % mess up ticks
    hold on
    h2=plot(day,ts,'ob');
    set(h2,'linewidth',2)
    %h3=plot(day,tsup_hat,'-.g');
    %h4=plot(day,tslo_hat,'-.g');
    set(gca,'fontsize',20)
    hold off
    grid on
    ylabel('LOS (mm)');
    %datetick('x','mmmyy')  % keepticks or see below

    ts_years=unique(datenum(datestr(day,'yyyy'),'yyyy'));
    ts_dates=[ts_years(1):365.25:ceil(ts_years(end)+365.25)];
    %ts_dates=[day(1):365:day(end)+365];
    set(gca, 'XTick',ts_dates);
    set(gca, 'XTickLabel', datestr(ts_dates,'yyyy'));
    % annotate velocit slope in ts plot

   % IFG Dates - excluding master
%    subplot(10,10,[18 90]) % subplot for rectangle
%       putdates(0.05,1,datestr(day,'yyyy-mm-dd'),0.035,9) % putdates(xstart, ystart, labels, labeloffset, fontsize)
   %subplot(10,1,10)
   %bar(Bdop)

%EOF

author: dbekaert

repo: StaMPS

save_path: github-repos/MATLAB/dbekaert-StaMPS


text:
%
% Precompute the rotation matrices to rotate the high-resolution kernels (500 directions)
%
% Parameters
% ----------
% lmax : unsigned int
%   Maximum spherical harmonics order to use for the rotation phase
%
function AMICO_PrecomputeRotationMatrices( lmax )
	if nargin < 1, lmax = 12; end
	global AMICO_data_path

	fprintf( '\n-> Precomputing rotation matrices for l_max=%d:\n', lmax );

	filename = fullfile(AMICO_data_path,sprintf('AUX_matrices__lmax=%d.mat',lmax) );
	if exist( filename, 'file' )
		fprintf( '   [ already computed ]\n' );
		return
	end

	TIME = tic();

	% load file with 500 directions
	grad500 = importdata( '500_dirs.txt' );
	for i = 1:size(grad500,1)
		grad500(i,:) = grad500(i,:) ./ norm( grad500(i,:) );
        if grad500(i,2) < 0
            grad500(i,:) = -grad500(i,:);
        end
	end

	% precompute the matrix to fit the SH coefficients
	[colatitude, longitude] = AMICO_Cart2sphere( grad500(:,1), grad500(:,2), grad500(:,3) );
	Ylm = AMICO_CreateYlm( lmax, colatitude, longitude );
	fit = pinv( Ylm'*Ylm ) * Ylm';

	% precompute the matrices to rotate the functions in SH space
	Ylm_rot = {};
	for ox = 0:180
	for oy = 0:180
		Ylm_rot{ox+1,oy+1} = AMICO_CreateYlm( lmax, ox/180.0*pi, oy/180.0*pi );
	end
	end

	save( filename , 'Ylm_rot', 'fit', 'lmax' );
	fprintf( '   [ %.1f seconds ]\n', toc(TIME) );
end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function [batches, batch_padding] = rcnn_extract_regions(im, boxes, rcnn_model)
% [batches, batch_padding] = rcnn_extract_regions(im, boxes, rcnn_model)
%   Extract image regions and preprocess them for use in Caffe.
%   Output is a cell array of batches.
%   Each batch is a 4-D tensor formatted for input into Caffe:
%     - BGR channel order
%     - single precision
%     - mean subtracted
%     - dimensions from fastest to slowest: width, height, channel, batch_index
%
%   im is an image in RGB order as returned by imread
%   boxes are in [x1 y1 x2 y2] format with one box per row

% AUTORIGHTS
% ---------------------------------------------------------
% Copyright (c) 2014, Ross Girshick
%
% This file is part of the R-CNN code and is available
% under the terms of the Simplified BSD License provided in
% LICENSE. Please retain this notice and LICENSE if you use
% this file (or any portion of it) in your project.
% ---------------------------------------------------------

% convert image to BGR and single
im = single(im(:,:,[3 2 1]));
num_boxes = size(boxes, 1);
batch_size = rcnn_model.cnn.batch_size;
num_batches = ceil(num_boxes / batch_size);
batch_padding = batch_size - mod(num_boxes, batch_size);
if batch_padding == batch_size
  batch_padding = 0;
end

crop_mode = rcnn_model.detectors.crop_mode;
image_mean = rcnn_model.cnn.image_mean;
crop_size = size(image_mean,1);
crop_padding = rcnn_model.detectors.crop_padding;

batches = cell(num_batches, 1);
%for batch = 1:num_batches
for batch = 1:num_batches
  batch_start = (batch-1)*batch_size+1;
  batch_end = min(num_boxes, batch_start+batch_size-1);

  ims = zeros(crop_size, crop_size, 3, batch_size, 'single');
  for j = batch_start:batch_end
    bbox = boxes(j,:);

    crop = rcnn_im_crop(im, bbox, crop_mode, crop_size, ...
        crop_padding, image_mean);

    % swap dims 1 and 2 to make width the fastest dimension (for caffe)
    ims(:,:,:,j-batch_start+1) = permute(crop, [2 1 3]);
  end

  batches{batch} = ims;
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function output = tokenize_cell_array(cellarray, num_headers)

% Takes a cell array, one entry per line and tokenizes it

end

author: faridani

repo: MatlabNLP

save_path: github-repos/MATLAB/faridani-MatlabNLP


text:
function [ a, l, r ] = i4vec_part_quick_a ( n, a )

%*****************************************************************************80
%
%% I4VEC_PART_QUICK_A reorders an I4VEC as part of a quick sort.
%
%  Discussion:
%
%    The routine reorders the entries of A.  Using A(1) as a key,
%    all entries of A that are less than or equal to the key will
%    precede the key which precedes all entries that are greater than the key.
%
%  Example:
%
%    Input:
%
%      N = 8
%
%      A = ( 6, 7, 3, 1, 6, 8, 2, 9 )
%
%    Output:
%
%      L = 3, R = 6
%
%      A = ( 3, 1, 2, 6, 6, 8, 9, 7 )
%            -------        -------
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 February 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of entries of A.
%
%    Input, integer A(N), the array to be checked.
%
%    Output, integer A(N), has been reordered as described above.
%
%    Output, integer L, R, the indices of A that define the three segments.
%    Let KEY = the input value of A(1).  Then
%    I <= L                 A(I) < KEY;
%         L < I < R         A(I) = KEY;
%                 R <= I    KEY < A(I).
%
  if ( n < 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4VEC_PART_QUICK_A - Fatal error!\n' );
    fprintf ( 1, '  N < 1.\n' );
    error ( 'I4VEC_PART_QUICK_A - Fatal error!' );
  elseif ( n == 1 )
    l = 0;
    r = 2;
    return;
  end

  key = a(1);
  m = 1;
%
%  The elements of unknown size have indices between L+1 and R-1.
%
  l = 1;
  r = n + 1;

  for i = 2 : n

    if ( key < a(l+1) )
      r = r - 1;
      [ a(r), a(l+1) ] = i4_swap ( a(r), a(l+1) );
    elseif ( a(l+1) == key )
      m = m + 1;
      [ a(m), a(l+1) ] = i4_swap ( a(m), a(l+1) );
      l = l + 1;
    elseif ( a(l+1) < key )
      l = l + 1;
    end

  end
%
%  Now shift small elements to the left, and KEY elements to center.
%
  for i = 1 : l - m
    a(i) = a(i+m);
  end
%
%  Out of bounds here, occasionally.
%
  l = l - m;

  a(l+1:l+m) = key;

  return
end


text:
%%**********************************************************
%% cheby0:
%%
%%    minimize || p(d) ||_infty
%%    p = polynomial of degree <= m such that p(0) = 1.
%%
%%    Here d = n-vector
%%----------------------------------------------------------
%% [blk,Avec,C,b,X0,y0,Z0,objval,p] = cheby0(d,m,solve);
%%
%% d = a vector.
%% m = degree of polynomial.
%% feas  = 1 if want feasible starting point
%%       = 0 if otherwise.
%% solve = 0 if just want initialization
%%       = 1 if want to solve the problem
%%*****************************************************************
%% SDPT3: version 4.0
%% Copyright (c) 1997 by
%% Kim-Chuan Toh, Michael J. Todd, Reha H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*****************************************************************

function [blk,Avec,C,b,X0,y0,Z0,objval,p] = cheby0(d,m,solve);

  if nargin <= 2; solve = 0; end;
  if (size(d,1) < size(d,2)); d = d.'; end;
  cmp = 1-isreal(d);

  tstart=cputime;
  n = length(d);
  e = ones(n,1);
  V(1:n,1) = e/norm(e); R(1,1) = 1/norm(e);
  for i =1:m
      v = d.*V(:,i);
      for j = 1:i                 %% Arnoldi iterations:
          H(j,i) = (V(:,j))'*v;   %% constructing upper-Hessenberg matrix.
          v = v - H(j,i)*V(:,j);  %% orthonormaliztion of Krylov basis.
      end;
      H(i+1,i) = norm(v);
      V(:,i+1) = v/H(i+1,i);
      R(1:i+1,i+1) = (1/H(i+1,i))*([0; R(1:i,i)] - [R(1:i,1:i)*H(1:i,i); 0]);
  end
  if (cmp)
     blk{1,1} = 'q'; blk{1,2} = 3*ones(1,n);
     C = zeros(3*n,1); C(2:3:3*n) = ones(n,1);
     b = [zeros(2*m,1); -1];
     Atmp = [];
     II = [0:3:3*n-3]'; ee = ones(n,1);
     for k=1:m
         dVk = d.*V(:,k);
         Atmp = [Atmp; [2+II, k*ee, real(dVk)]; [3+II, k*ee, imag(dVk)]];
         Atmp = [Atmp; [2+II, (m+k)*ee, -imag(dVk)]; [3+II, (m+k)*ee, real(dVk)]];
     end
     Atmp = [Atmp;  [1+II, (2*m+1)*ee, -ones(n,1)]];
  else
     blk{1,1} = 'l'; blk{1,2} = 2*ones(1,n);
     b = [zeros(m,1); -1];
     C = [ones(n,1); -ones(n,1)];
     Atmp = [];
     II = [1:n]'; ee = ones(n,1);
     for k=1:m
         dVk = d.*V(:,k);
         Atmp = [Atmp; [II, k*ee, dVk]; [n+II, k*ee, -dVk]];
     end
     Atmp = [Atmp; [II, (m+1)*ee, -ee]; [n+II, (m+1)*ee, -ee]];
  end
  Avec = spconvert(Atmp);
  [X0,y0,Z0] = infeaspt(blk,Avec,C,b);
%%
  if (solve)
     [obj,X,y,Z] = sqlp(blk,Avec,C,b,[],X0,y0,Z0);
     if (cmp)
        y  = y(1:m) + sqrt(-1)*y(m+1:2*m);
     else
        y = y(1:m);
     end
     x1 = R(1:m,1:m)*y(1:m);
     p = [-x1(m:-1:1); 1];
     objval = -mean(obj);
  else
     objval = []; p = [];
  end
%%**********************************************************

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function [mGal] = uGal2mGal(uGal)
% Convert acceleration from microgals to milligals.
% Chad A. Greene 2012
mGal = uGal*1e-3;
end


text:
% CXCORR Circular Cross Correlation function estimates.
% CXCORR(a,b), where a and b represent samples taken over time interval T
% which is assumed to be a common period of two corresponded periodic signals.
% a and b are supposed to be length M row vectors, either real or complex.
%
% [x,c]=CXCORR(a,b) returns the length M-1 circular cross correlation sequence c
% with corresponded lags x.
%
% The circular cross correlation is:
%         c(k) = sum[a(n)*conj(b(n+k))]/[norm(a)*norm(b)];
% where vector b is shifted CIRCULARLY by k samples.
%
% The function doesn't check the format of input vectors a and b!
%
% For circular covariance between a and b look for CXCOV(a,b) in
% http://www.mathworks.com/matlabcentral/fileexchange/loadAuthor.do?objectType=author&objectId=1093734
%
% Reference:
% A. V. Oppenheim, R. W. Schafer and J. R. Buck, Discrete-Time Signal Processing,
% Upper Saddler River, NJ : Prentice Hall, 1999.
%
% Author: G. Levin, Apr. 26, 2004.

function [x,c]=CXCORR(a,b)
na=norm(a);
nb=norm(b);
a=a/na; %normalization
b=b/nb;
for k=1:length(b)
    c(k)=a*b';
    b=[b(end),b(1:end-1)]; %circular shift
end
x=[0:length(b)-1]; %lags


text:
function [cc] = l2cc(l)
% Convert volume from liters to cubic centimeters.
% Chad Greene 2012
cc = l*1000;


text:
function doTheCalcThingSLI(obj)
		%everything for the calculation itself goes in here
		import mapseis.projector.*;
		import mapseis.declustering.*;

		%get parameter into variables
		ShortCat = getShortCatalog(obj.Datastore,obj.SendedEvents);
		Sli_Type=obj.CalcParameter.Sli_Type;
		MainMag=obj.CalcParameter.MainMag;
		BackTime=obj.CalcParameter.BackTime;
		ForeTime=obj.CalcParameter.ForeTime;
		MinRad=obj.CalcParameter.MinRad;
		MagRadScale=obj.CalcParameter.MagRadScale;
		MagRadConst=obj.CalcParameter.MagRadConst;

		if Sli_Type==1
			%run declustering with magnitude
			disp('please wait....');
			tic;
			[clusterID,EventType,AlgoInfo] = clusterSLIDERmag(ShortCat,MainMag,BackTime,ForeTime,MinRad,MagRadScale,MagRadConst);
			obj.CalcTime=toc;
			disp('finished!');

		elseif Sli_Type==2
			%run declustering with time
			disp('please wait....');
			tic;
			[clusterID,EventType,AlgoInfo] = clusterSLIDERtime(ShortCat,MainMag,BackTime,ForeTime,MinRad,MagRadScale,MagRadConst);
			obj.CalcTime=toc;
			disp('finished!');

		end

		%store result
		obj.CalcRes.clusterID=clusterID;
		obj.CalcRes.EventType=EventType;
		obj.CalcRes.AlgoInfo=AlgoInfo;

		%write to datastore
		obj.Datastore.setDeclusterData(EventType,clusterID,obj.SendedEvents,[]);
		obj.Datastore.DeclusterMetaData=AlgoInfo;

		%correct for old datastore version, it makes sense to do this here, because it is needed only by the declustering
		NumberedUserData=getNumberedFields(obj.Datastore);
		setNumberedFields(obj.Datastore,union(NumberedUserData,{'Month','Day','Hour','Minute','Second','MilliSecond','DecYear'}));

end

author: CelsoReyes


text:
function avLL = GammaavLL(hmm,Gamma,Xi,T)
% average loglikelihood for state time course

% if isfield(hmm.train,'grouping') % DEPRECATED
%     Q = length(unique(hmm.train.grouping));
% else
%     Q = 1;
% end
Q = 1;
N = length(T); K = hmm.K;
order = (sum(T) - size(Gamma,1))/N;
avLL = 0;
do_clustering = isfield(hmm.train,'cluster') && hmm.train.cluster;

% avLL initial state % DEPRECATED
% if Q>1
%     for i = 1:Q
%         PsiDir_alphasum = psi(sum(hmm.Dir_alpha(:,i)));
%         ii = hmm.train.grouping==i;
%         for l = 1:K
%             if ~hmm.train.Pistructure(l), continue; end
%             avLL = avLL + sum(Gamma(jj(ii),l)) * (psi(hmm.Dir_alpha(l,i)) - PsiDir_alphasum);
%         end
%     end
% else
%     PsiDir_alphasum = psi(sum(hmm.Dir_alpha));
%     for l = 1:K
%         if ~hmm.train.Pistructure(l), continue; end
%         avLL = avLL + sum(Gamma(jj,l)) * (psi(hmm.Dir_alpha(l)) - PsiDir_alphasum);
%     end
% end

if ~isempty(Xi) && ~do_clustering % a proper HMM
    jj = zeros(N,1); % reference to first time point of the segments
    for in = 1:N
        jj(in) = sum(T(1:in-1)) - order*(in-1) + 1;
    end
    PsiDir_alphasum = psi(sum(hmm.Dir_alpha));
    % first time point
    for l = 1:K
        if ~hmm.train.Pistructure(l), continue; end
        avLL = avLL + sum(Gamma(jj,l)) * (psi(hmm.Dir_alpha(l)) - PsiDir_alphasum);
    end
    % avLL remaining time points
    for i = 1:Q
        if Q > 1
            ii = find(hmm.train.grouping==i)';
        else
            ii = 1:length(T);
        end
        PsiDir2d_alphasum = zeros(K,1);
        for l = 1:K, PsiDir2d_alphasum(l) = psi(sum(hmm.Dir2d_alpha(l,:,i))); end
        %PsiDir2d_alphasum = psi(sum(sum(hmm.Dir2d_alpha(:,:,i))));
        for k = 1:K
            for l = 1:K
                if ~hmm.train.Pstructure(l,k), continue; end
                if Q==1
                    avLL = avLL + sum(Xi(:,l,k)) * (psi(hmm.Dir2d_alpha(l,k))-PsiDir2d_alphasum(l));
                    %avLL = avLL + sum(Xi(:,l,k)) * (psi(hmm.Dir2d_alpha(l,k))-PsiDir2d_alphasum);
                    if isnan(avLL)
                        error(['Error computing log likelihood of the state time courses  - ' ...
                            'Out of precision?'])
                    end
                else
                    for n = ii
                        t = (1:T(n)-1-order) + sum(T(1:n-1)) - (order+1)*(n-1) ;
                        avLL = avLL + sum(Xi(t,l,k)) * (psi(hmm.Dir2d_alpha(l,k,i))-PsiDir2d_alphasum(l));
                    end
                end
            end
        end
    end
else % Simple mixture of distributions
    PsiDir_alphasum = psi(sum(hmm.Dir_alpha));
    for k = 1:K
        avLL = avLL + sum(Gamma(:,k)) * (psi(hmm.Dir_alpha(k)) - PsiDir_alphasum);
    end
end

end

author: OHBA-analysis

repo: HMM-MAR


text:
% VARGPLVMFEATURECLASSIFICATION Use vargplvm features to perform
% discriminative classification with log. regression.
%
%
% ARG: model   : The vargplvm model
% ARG: data    : A struct with:
%                 data.lbls      : Training labels in 1-of-K encoding
%                 data.lblsTs    : Test labels
%                 data.Xpred     : A test latent point (optional if data.Yts given)
%                 data.Xpred_var : A test latent variance (optional if data.Yts given)
%                 data.Yts       : Test output (optional if the 2 above given)
% ARG: options : A struct with options for the function: (Optional)
%                 options.samplesPerObserved: Populate tr. set with this amount of samples per q(x_i)
%                 options.samplesPerOutput  : When predicting FROM q(x*),
%                                             instead predict from many x*_j's sampled from q(x*) and
%                                             average predictions. samplesPerOutput says how many
%                                             x*_j's to sample.
%----
% RET: LogRegError      : Error when not populating training/predictions
% RET: LogRegErrorExt   : Error when populating training set
% RET: LogRegErrorExtOut: Error when populating trainign and predictions
%
% COPYRIGHT: Andreas C. Damianou, 2015
%
% SEEALSO: vargplvmPredictLatent.m
function [LogRegError, LogRegErrorExt, LogRegErrorExtOut, X_pred, varX_pred] = vargplvmFeatureClassification(model, data, options)

% Train classifiers from X to labels
labels = transformLabels(data.lbls)';
labelsTs = transformLabels(data.lblsTs)';

if ~isfield(data, 'X_pred')
    [X_pred, varX_pred] = ...
        vargplvmPredictLatent(model, data.Yts, [], false, model.globalOpt.reconstrIters,0,[],[],1);
else
    X_pred = data.X_pred;
    varX_pred = data.varX_pred;
end

if nargin < 3 || ~isfield(options, 'samplesPerObserved')
    samplesPerObserved = 5;
else
    samplesPerObserved = options.samplesPerObserved;
end
if nargin < 3 || ~isfield(options, 'samplesPerOutput')
    samplesPerOutput = 10;
else
    samplesPerOutput = options.samplesPerOutput;
end
% Use only SOME dimensions as features (e.g. the ones with best ARD weight)
if nargin < 3 || ~isfield(options, 'dims') || isempty(options.dims)
    dims = 1:model.q;
else
    dims = options.dims;
end

% Populate training set
%--- Take variance into account by sampling new data from the distribution
% Init sizes
Xnew = nan(size(model.vardist.means,1)*samplesPerObserved, size(model.vardist.means,2));
labelsNew = nan(size(Xnew,1), size(labels,2));
k=1;
% Take samples
for n=1:size(model.vardist.means,1)
    for kk = 1:samplesPerObserved
        Xnew(k,:) = model.vardist.means(n,:) + randn(size(model.vardist.means(n,:))).*sqrt(model.vardist.covars(n,:));
        labelsNew(k,:) = labels(n,:);
        k = k + 1;
    end
end
% Augment set with samples
Xext = [model.vardist.means; Xnew];
labelsExt = [labels; labelsNew];
clear 'Xnew' 'labelsNew';

% Training of logistic regression classifier. One for each label
% separately.
nClasses = length(unique(labels));
clear 'B' 'BExt'
for i=1:nClasses
    fprintf('\n # LogReg training for class # %d\n', i)
    lb = zeros(size(model.vardist.means(:,dims),1),1);
    lbExt = zeros(size(Xext(:,dims),1),1);
    lb(labels == i) = 1;
    lbExt(labelsExt == i) = 1;
    B{i} = glmfitWrapper(model.vardist.means(:,dims), lb,'binomial','logit',[],[],[],[],1000);
    BExt{i} = glmfitWrapper(Xext(:,dims), lbExt,'binomial','logit',[],[],[],[],1000);
end

svmmodel    = svmtrain(labels,                 model.vardist.means(:,dims),'-q');
%svmmodelExt = svmtrain(transformLabels(lbExt), model.layer{options.lOut}.vardist.means(:,dims),'-q');

[~, acc,~]    = svmpredict(labelsTs, X_pred(:,dims), svmmodel,'-q');
%[~, accExt] = svmpredict(labelsTs',X_pred, svmmodelExt);

% Prediction of each binary classifier
Ypred_logReg = zeros(size(data.lblsTs));
Ypred_logRegExtOut = zeros(size(data.lblsTs));
Ypred_logRegExt = zeros(size(data.lblsTs));
for i=1:nClasses
    for k=1:samplesPerOutput
        % Sample from the OUTPUT distribution, and then average predictions
        Xsamp = X_pred + randn(size(X_pred)).*sqrt(varX_pred);
        Ypred_logRegExtOut(:,i) = Ypred_logRegExtOut(:,i)+1/samplesPerOutput*glmval(BExt{i}, Xsamp(:,dims), 'logit');
    end
    Ypred_logReg(:,i) = glmval(B{i}, X_pred(:,dims), 'logit')';
    Ypred_logRegExt(:,i) = glmval(BExt{i}, X_pred(:,dims), 'logit')';
end
% Replace predictions with maximum probability (ie, make a decision)
[~,ind]=max(Ypred_logReg');
[~,indExt]=max(Ypred_logRegExt');
[~,indExtOut]=max(Ypred_logRegExtOut');
LogRegError = 0;
LogRegErrorExt = 0;
LogRegErrorExtOut = 0;
for i=1:size(X_pred,1)
    LogRegError = LogRegError + (ind(i) ~= labelsTs(i));
    LogRegErrorExt = LogRegErrorExt + (indExt(i) ~= labelsTs(i));
    LogRegErrorExtOut = LogRegErrorExtOut + (indExtOut(i) ~= labelsTs(i));
end

% fprintf('\n========================== REPORT ===========================\n')
% fprintf('# Error                                             : %d\n', LogRegError);
% fprintf('# Error with %d samp. PerObserved                    : %d\n', samplesPerObserved,LogRegErrorExt);
% fprintf('# Error with %d samp. PerObserved, %d samp.PerOutput : %d\n', samplesPerObserved, samplesPerOutput,LogRegErrorExtOut);
% fprintf('----------------------------------------------------------------\n\n');
%

N = size(X_pred,1);
fprintf('\n========================== REPORT ===========================\n')
fprintf('# Acc Reg:                                             : %.2f%%\n', (N-LogRegError)/N * 100);
fprintf('# Acc Reg: with %d samp. PerObserved                    : %.2f%%\n', samplesPerObserved,(N-LogRegErrorExt)/N*100);
fprintf('# Acc Reg: with %d samp. PerObserved, %d samp.PerOutput : %.2f%%\n', samplesPerObserved, samplesPerOutput,(N-LogRegErrorExtOut)/N*100);
fprintf('# Acc SVM:                                             : %.2f%%\n', acc(1));
%fprintf('# Acc SVM: with % samp. PerObserved                    : %.2f%%\n', samplesPerObserved,accExt(1));
fprintf('----------------------------------------------------------------\n\n');

repo: vargplvm


text:
% EX_MAXWELL_SRC_PERIODIC_SQUARE: solve Maxwell source problem in the unit square.

% 1) PHYSICAL DATA OF THE PROBLEM
clear problem_data
% Physical domain, defined as NURBS map given in a text file
problem_data.geo_name = 'geo_square.txt';
% The domain must also have the right continuity conditions on the periodic
%  sides, otherwise the convergence rate may deteriorate.
%  See the example in ex_laplace_square_periodic.

% Type of boundary conditions
problem_data.nmnn_sides     = [];
problem_data.drchlt_sides   = [3 4];
problem_data.periodic_directions = [1];

% Physical parameters
problem_data.c_stiff = @(x, y) ones(size(x));
problem_data.c_mass  = @(x, y) pi^2 * ones(size(x));

% Source and boundary terms
problem_data.f = @(x, y) cat(1, zeros (1, size (x, 1), size (x, 2)), ...
                             reshape (17*pi^2*sin(4*pi*x), [1, size(x)]));
problem_data.g = @(x, y, ind) zeros (2, size (x, 1), size (x, 2));
problem_data.h = @(x, y, ind) zeros (2, size (x, 1), size (x, 2));

% Exact solution (optional)
problem_data.uex     = @(x, y) cat(1, zeros (1, size (x, 1), size (x, 2)), ...
                                      reshape ( sin(4*pi*x), [1, size(x)]));
problem_data.curluex = @(x, y) 4*pi*cos(4*pi*x);

% 2) CHOICE OF THE DISCRETIZATION PARAMETERS
clear method_data
method_data.degree     = [3 3];     % Degree of the bsplines
method_data.regularity = [2 2];     % Regularity of the splines
method_data.nsub       = [9 9];   % Number of subdivisions
method_data.nquad      = [4 4];     % Points for the Gaussian quadrature rule

% 3) CALL TO THE SOLVER
[geometry, msh, space, u] = solve_maxwell_src (problem_data, method_data);

% 4) POST-PROCESSING
vtk_pts = {linspace(0, 1, 30), linspace(0, 1, 30)};
% 4.1) EXPORT TO PARAVIEW
output_file = 'maxwell_square_periodic_Deg3_Reg2_Sub8';
fprintf ('The result is saved in the file %s \n \n', output_file);
sp_to_vtk (u, space, geometry, vtk_pts, output_file, 'u')

% 4.2) Plot in Matlab. Comparison with the exact solution
[eu, F] = sp_eval (u, space, geometry, vtk_pts);
[X, Y]  = deal (squeeze(F(1,:,:)), squeeze(F(2,:,:)));
eu2     = problem_data.uex (X, Y);

subplot(1,2,1)
quiver (X, Y, squeeze(eu(1,:,:)), squeeze(eu(2,:,:)))
axis equal tight
title('Computed solution')
ylim([0,1]); xlim([0,1]);
subplot(1,2,2)
quiver (X, Y, squeeze(eu2(1,:,:)), squeeze(eu2(2,:,:)))
axis equal tight
title('Exact solution')
ylim([0,1]); xlim([0,1]);

[error_hcurl, error_l2] = ...
    sp_hcurl_error (space, msh, u, problem_data.uex, problem_data.curluex)

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function [res] = tt_mvdot(a_tt,x_tt,y_tt)
%Fast computation of the bilinear form (Ax,y) in the TT-format
%   [RES]=TT_MVDOT(A_TT,X_TT,Y_TT) Computes scalar product(AX,Y), where
%   A_TT is a matrix, X_TT is a vector, Y_TT  is a vector
%
%
% TT-Toolbox 2.2, 2009-2012
%
%This is TT Toolbox, written by Ivan Oseledets et al.
%Institute of Numerical Mathematics, Moscow, Russia
%webpage: http://spring.inm.ras.ru/osel
%
%For all questions, bugs and suggestions please mail
%ivan.oseledets@gmail.com
%---------------------------
d=size(a_tt,1);
%The c
core=a_tt{1}; matx=x_tt{1}; maty=y_tt{1};
n1=size(core,1); n2=size(core,2);
rc=size(core,3);
rx=size(matx,2); ry=size(maty,2);
core=maty'*reshape(core,[n1,n2*rc]); %core is ryxn2xrc
core=reshape(permute(reshape(core,[ry,n2,rc]),[3,1,2]),[rc*ry,n2]);
phi=permute(reshape(core*matx,[rc,ry,rx]),[1,3,2]); %core is rcxryxrx-> rcxrx*ry

for i=2:d-1
  core=a_tt{i}; core_x = x_tt{i}; core_y=y_tt{i};
  %Decipher sizes
 n1=size(core,1); n2=size(core,2); rc1=size(core,3); rc2=size(core,4);
 rx1=size(core_x,2); rx2=size(core_x,3); ry1=size(core_y,2); ry2=size(core_y,3);
 %Convolve phi & core over the first index
 %core is n1xn2xrc1xrc2 phi is rc1x rx1 x ry1
  core=reshape(permute(core,[3,1,2,4,5]),[rc1,n1*n2*rc2]);
  phi=reshape(phi,[rc1,rx1*ry1])'*core; %Is rx1 x ry1 x n1 x n2 x rc2
  %Convolve over n2 & rx1
  phi=permute(reshape(phi,[rx1,ry1,n1,n2,rc2]),[4,1,2,3,5]);
  %n2 x rx1 x ry1 x n1 x rc2;
  phi=reshape(core_x,[n2*rx1,rx2])'*reshape(phi,[n2*rx1,ry1*n1*rc2]);
  %phi is now rx2 x ry1 x n1 x rc2
  %Convolve over n1 & ry1
  phi=reshape(phi,[rx2,ry1*n1,rc2]); phi=permute(phi,[3,1,2]);
  %keyboard;
  phi=reshape(phi,[rc2*rx2,ry1*n1])*reshape(permute(core_y,[2,1,3]),[ry1*n1,ry2]);
  phi=reshape(phi,[rc2,rx2,ry2]);
end
  core=a_tt{d}; mat_x=x_tt{d}; mat_y=y_tt{d};
  n1=size(core,1); n2=size(core,2); rc=size(core,3);
  rx=size(mat_x,2); ry=size(mat_y,2);
  %Convolve phi & core over aux index
  %core is n1 x n2 x rc
  phi=reshape(core,[n1*n2,rc])*reshape(phi,[rc,rx*ry]);
  %phi is n1 x n2 x rx x ry
  %Convolve over n2 & rx
  phi=permute(reshape(phi,[n1,n2*rx,ry]),[1,3,2]);
  phi=reshape(phi,[n1*ry,n2*rx])*reshape(mat_x,[n2*rx,1]);
  %phi is a column of length
  res= phi'*reshape(mat_y,[n1*ry,1]);

return

end

author: oseledets

repo: TT-Toolbox


text:
function [event] = read_trigger(filename, varargin)

% READ_TRIGGER extracts the events from a continuous trigger channel
% This function is a helper function to read_event and can be used for all
% dataformats that have one or multiple continuously sampled TTL channels
% in the data.
%
% The optional trigshift (default is 0) causes the value of the
% trigger to be obtained from a sample that is shifted N samples away
% from the actual flank.
%
% This is a helper function for READ_EVENT
%
% TODO
%  - merge read_ctf_trigger into this function (requires trigshift and bitmasking option)
%  - merge biosemi code into this function (requires bitmasking option)

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

event = [];

% get the optional input arguments
hdr         = keyval('header',      varargin);
dataformat  = keyval('dataformat',  varargin);
begsample   = keyval('begsample',   varargin);
endsample   = keyval('endsample',   varargin);
chanindx    = keyval('chanindx',    varargin);
detectflank = keyval('detectflank', varargin);
denoise     = keyval('denoise',     varargin); if isempty(denoise),     denoise = 1;      end
trigshift   = keyval('trigshift',   varargin); if isempty(trigshift),   trigshift = 0;    end
trigpadding = keyval('trigpadding', varargin); if isempty(trigpadding), trigpadding = 1;  end
fixctf      = keyval('fixctf',      varargin); if isempty(fixctf),      fixctf = 0;       end
fixneuromag = keyval('fixneuromag', varargin); if isempty(fixneuromag), fixneuromag = 0;  end
fix4dglasgow= keyval('fix4dglasgow', varargin); if isempty(fix4dglasgow), fix4dglasgow = 0; end

if isempty(hdr)
  hdr = ft_read_header(filename);
end

if isempty(begsample)
  begsample = 1;
end

if isempty(endsample)
  endsample = hdr.nSamples*hdr.nTrials;
end

% read the trigger channel as raw data, can safely assume that it is continuous
dat = ft_read_data(filename, 'header', hdr, 'dataformat', dataformat, 'begsample', begsample, 'endsample', endsample, 'chanindx', chanindx, 'checkboundary', 0);

if isempty(dat)
  % there are no triggers to detect
  return
end

% Detect situations where the channel value changes almost at every time
% step which are likely to be noise
if denoise
  for i=1:length(chanindx)
    if (sum(diff(find(diff(dat(i,:))~=0)) == 1)/length(dat(i,:))) > 0.8
      warning(['trigger channel ' hdr.label{chanindx(i)} ' looks like noise and will be ignored']);
      dat(i,:) = 0;
    end
  end
end

if fixctf
  % correct for reading the data as signed 32-bit integer, whereas it should be interpreted as an unsigned int
  dat(dat<0) = dat(dat<0) + 2^32;
end

if fixneuromag
  % according to Joachim Gross, real events always have triggers > 5
  % this is probably to avoid the noisefloor
  dat(dat<5) = 0;
end

if fix4dglasgow
  % synchronization pulses have a value of 8192 and are set to 0
  dat = dat - bitand(dat, 8192);
  %% triggers containing the first bit assume a value of 4096 when sent by presentation
  %% this does not seem to hold for matlab; check this
  %dat = dat - bitand(dat, 4096)*4095/4096;
end

for i=1:length(chanindx)
  % process each trigger channel independently
  channel = hdr.label{chanindx(i)};
  trig    = dat(i,:);

  if trigpadding
    pad = trig(1);
  else
    pad = 0;
  end

  switch detectflank
    case 'up'
      % convert the trigger into an event with a value at a specific sample
      for j=find(diff([pad trig(:)'])>0)
        event(end+1).type   = channel;
        event(end  ).sample = j + begsample - 1;      % assign the sample at which the trigger has gone down
        event(end  ).value  = trig(j+trigshift);      % assign the trigger value just _after_ going up
      end
    case 'down'
      % convert the trigger into an event with a value at a specific sample
      for j=find(diff([pad trig(:)'])<0)
        event(end+1).type   = channel;
        event(end  ).sample = j + begsample - 1;      % assign the sample at which the trigger has gone down
        event(end  ).value  = trig(j-1-trigshift);    % assign the trigger value just _before_ going down
      end
    case 'both'
      % convert the trigger into an event with a value at a specific sample
      for j=find(diff([pad trig(:)'])>0)
        event(end+1).type   = [channel '_up'];        % distinguish between up and down flank
        event(end  ).sample = j + begsample - 1;      % assign the sample at which the trigger has gone down
        event(end  ).value  = trig(j+trigshift);      % assign the trigger value just _after_ going up
      end
      % convert the trigger into an event with a value at a specific sample
      for j=find(diff([pad trig(:)'])<0)
        event(end+1).type   = [channel '_down'];      % distinguish between up and down flank
        event(end  ).sample = j + begsample - 1;      % assign the sample at which the trigger has gone down
        event(end  ).value  = trig(j-1-trigshift);    % assign the trigger value just _before_ going down
      end
    otherwise
      error('incorrect specification of ''detectflank''');
  end
end

author: PatternRecognition

repo: OpenBMI


text:
function dirS = GetDirEdit(dirS)

% dirS = GetDirEdit(dirS);
%
% Allow user to edit the fields of the session directory structure
% [dirS].

% The following two cell arrays determine the field names to
% display, their descriptive label, and whether or not they are
% to be editable. Each row of the array has the form:
% {field name, label, edit flag}.
dirFields = { ...
'home', 'Home (rw)', 1; ...
'raw', 'Raw (r)', 1; ...
};

% Create the top-level fields
% if isunix
%
%   fontSize = 10;
%
% else
%   fontSize = 9;
%
% end
fontSize = 12;
titleFontSize = 12;
butWidth=10;
marginFields = 3.5;
topMargin = 1.5;
% scale factors for x and y axis coordinates
xs = 2.45;
ys = 1.6;
height = 1 * ys;
vSkip = 0.10 * ys;
dy = (height + vSkip);
x = 1 * xs;

dirData = CreateEditData(dirFields, dirS);
numFields = length(dirData);
maxWidth = max([dirData.width]);
pos = [35, 3, maxWidth, (numFields+marginFields)] .* [2, ys, 2.1, dy];

% Create the figure:
topH = figure(...
    'MenuBar', 'none', ...
    'Units', 'char', ...
    'Position', pos, ...
    'Resize','off', ...
    'Name', 'FMRI Directories', ...
    'NumberTitle','off' ...
    );
bkColor = get(topH, 'color');
topHS = num2str(topH);

% Label the session parameters:
y = (numFields + marginFields - topMargin) * dy;
titleStr = 'Choose/modify these directories';
uicontrol( ...
    'Style', 'text', ...
    'Units', 'char', ...
    'String', titleStr, ...
    'BackgroundColor', bkColor, ...
    'Position', [0, y, length(titleStr)*xs*1.5, height*1.2], ...
    'HorizontalAlignment', 'left', ...
    'FontWeight', 'bold', ...
    'FontSize', titleFontSize ...
    );
maxLabelWidth = 0.9 * max(size(char(dirData.label), 2)) * xs;
y = y - dy;
% Create the edit fields:
for iField=1:numFields
  labelPos = [x, y, length(dirData(iField).label)*xs, height];
  h = CreateEditRow(dirData(iField), labelPos, maxLabelWidth, topH, ...
    fontSize, 'UpdateDirEdit');
  dirData(iField).handle = h;
  y = y - dy;
end

% Build the UI data structure and attach it to the figure:
uiData.dirS = dirS;
uiData.origDirS = dirS;
uiData.dirData = dirData;
set(topH, 'UserData', uiData);

% Finally, install the file-control buttons
y = y - dy/2;
butWidth = 11;
bpos = [x, y, butWidth, height*1.2];
uicontrol( ...
    'Style', 'pushbutton', ...
    'String', 'Accept', ...
    'Units', 'char', ...
    'HorizontalAlignment', 'center', ...
    'Position', bpos, ...
    'Callback', 'uiresume', ...
    'FontSize', titleFontSize ...
    );
x = x + butWidth + 1;
bpos = [x, y, butWidth, height*1.2];
uicontrol( ...
    'Style', 'pushbutton', ...
    'String', 'Revert', ...
    'Units', 'char', ...
    'HorizontalAlignment', 'center', ...
    'Position', bpos, ...
    'FontSize', titleFontSize, ...
    'Callback', ['RevertDirEdit(', topHS, ');'] ...
    );
x = x + butWidth + 1;
bpos = [x, y, butWidth, height*1.2];
uicontrol( ...
    'Style', 'pushbutton', ...
    'String', 'Cancel', ...
    'Units', 'char', ...
    'HorizontalAlignment', 'center', ...
    'Position', bpos, ...
    'FontSize', titleFontSize, ...
    'Callback', ['CancelEdit(', topHS, ');'] ...
    );

% Wait until we get a uiresume, then perform an update. Repeat
% this cycle until the update reports no errors.
ok = 0;
while ~ok
  uiwait(topH);
  uiData = get(topH, 'UserData');
  if isfield(uiData, 'cancel')
    dirS = [];
    close(topH);
    return
  end
  ok = UpdateDirEdit(topH);
end

% After the update is successful, unpack the data and clean up.
uiData = get(topH, 'UserData');
dirS = uiData.dirS;
close(topH);

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function EP = updateEP(EP,Offsprings,nEP)
% Update the external population

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

% This function is written by Lucas Farias

    %% Select the non-dominated solutions
    EP = [EP,Offsprings];
    EP = EP(NDSort(EP.objs,1)==1);
    [N,M] = size(EP.objs);

    %% Delete the overcrowded solutions
    Dis = pdist2(EP.objs,EP.objs);
    Dis(logical(eye(length(Dis)))) = inf;
    Del = false(1,N);
    while sum(Del) < N-nEP
        Remain = find(~Del);
        subDis = sort(Dis(Remain,Remain),2);
        [~,worst] = min(prod(subDis(:,1:min(M,length(Remain))),2));
        Del(Remain(worst)) = true;
    end

    EP = EP(~Del);

end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [s_m_path, x_m_path, y_m_path, psi_rad_path,...
  d_m_veh, psi_rad_veh_p] = findPathPos(s_m_vec, x_m_vec, y_m_vec, psi_rad_vec,...
  x_m_veh, y_m_veh, psi_rad_veh_g)
%_________________________________________________________________
%% Documentation
%
% Authors:      Alexander Wischnewski (alexander.wischnewski@tum.de)
%
% Start Date:   01.02.2018
%
% Description:  Outputs the closest match for the actual position in the
%               given path vector. It is based on a brute force search.
%
% Inputs:
%   s_m_vec:        Vector with arc length parameter in meter
%   x_m_vec:        Vector with path x position in meter
%   y_m_vec:        Vector with path y position in meter
%   psi_rad_vec:    Vector with path orientation in radians
%   x_m_veh:        Actual vehicle position in x coordinates
%   y_m_veh:        Actual vehicle position in y coordinates
% Outputs:
%   s_m_path:       Arc length in meter of closest point
%   x_m_path:       x position in meter of closest point
%   y_m_path:       y position in meter of closest point
%   psi_rad_path:   Orientation in radians of closest point
%   d_m_veh:        Vehicle lateral deviation in meters from closest point
%   psi_rad_veh:    Vehicle heading deviation from closest point

% calculate difference vectors
diff_x_m = x_m_veh - x_m_vec;
diff_y_m = y_m_veh - y_m_vec;
dist_squared_trialpoints = diff_x_m.^2 + diff_y_m.^2;
% find minimum distance element
[mindist_squared_m, idx_mindist] = min(dist_squared_trialpoints);
% retrieve path points
s_m_path = s_m_vec(idx_mindist);
x_m_path = x_m_vec(idx_mindist);
y_m_path = y_m_vec(idx_mindist);
psi_rad_path = psi_rad_vec(idx_mindist);
% calculate sign of control error
diffvec_rot = [cos(-psi_rad_vec(idx_mindist)), -sin(-psi_rad_vec(idx_mindist))]*...
  [diff_x_m(idx_mindist);diff_y_m(idx_mindist)];
% calculate path orientation
psi_rad_veh_p = normalizeAngle(psi_rad_veh_g - psi_rad_path);
% project difference vector onto path
d_m_veh = -sign(diffvec_rot)*sqrt(mindist_squared_m)*cos(psi_rad_veh_p);

author: TUMFTM


text:
function [str] = struct2assoc(s, varargin)
    % Converts a Matlab structure into a string that describes an
    % expression for Mathematica association
    %
    % @author ayonga @date 2016-11-23
    %
    % Copyright (c) 2016, AMBER Lab
    % All right reserved.
    %
    % Redistribution and use in source and binary forms, with or without
    % modification, are permitted only in compliance with the BSD 3-Clause
    % license, see
    % http://www.opensource.org/licenses/bsd-license.php

    if length(s) == 1 % a scalar structure
        fields = fieldnames(s);
        nfield = length(fields);

        field_strs = cell(nfield, 1);
        for i = 1:nfield
            field = fields{i};
            value = s.(field);
            value_str = general2math(value, varargin{:});
            %         if iscell(value)
            %             % Use double-braces to prevent a struct array from being created
            %             value_str = ['{', value_str, '}'];
            %         end
            field_strs{i} = [str2mathstr(field), '->', value_str];
        end

        str = ['<| ', ...
            implode(field_strs,', '), ...
            ' |>'];
    elseif length(s) > 1 % structure array
        expr_s = arrayfun(@(y)struct2assoc(y, varargin{:}),s,'UniformOutput',false);
        str = cell2tensor(expr_s,'ConvertString',false);
    end

end

author: ayonga

repo: frost-dev

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


text:
function dUdabcd = msfmRotationDerivative(Q)
% wrapper around a mex to find the derivative of a rotation matrix
%
% Find the derivatives of a rotation matrix with respect to its quaternions
% For each quaternion in Q, it returns a [ 2 x 3 x 4 ] matrix: [ 2 x 3 ]
% is the derivative of the first two rows of the rotation, with respect to
% on of the quaternions
%
% USAGE
%  dUdabcd = msfmRotationDerivativeMex(Q);
%
% INPUTS
%  Q          - [ 4 x nFrame ] initial guess of the quaternions
%
% OUTPUTS
%  dUdabcd    - [ 2 x 3 x 4 x nFrame ] optimized quaternions
%
% EXAMPLE
%
% See also
%
% Vincent's Structure From Motion Toolbox      Version 3.0
% Copyright (C) 2008-2011 Vincent Rabaud.  [vrabaud-at-cs.ucsd.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the GPL [see external/gpl.txt]

dUdabcd = msfmRotationDerivativeMex(Q);

author: vrabaud


text:
function [fb_sim] = acc_gen (ref, imu)
% acc_gen: generates simulated accelerometers measurements from reference
%  data and IMU error profile.
%
% INPUT
%	ref: data structure with true trajectory.
%	imu: data structure with IMU error profile.
%
% OUTPUT
%	fb_sim: Nx3 matrix with simulated accelerations in the
%		body frame [X Y Z] (m/s^2, m/s^2, m/s^2).
%
%   Copyright (C) 2014, Rodrigo González, all rights reserved.
%
%   This file is part of NaveGo, an open-source MATLAB toolbox for
%   simulation of integrated navigation systems.
%
%   NaveGo is free software: you can redistribute it and/or modify
%   it under the terms of the GNU Lesser General Public License (LGPL)
%   version 3 as published by the Free Software Foundation.
%
%   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 Lesser General Public License for more details.
%
%   You should have received a copy of the GNU Lesser General Public
%   License along with this program. If not, see
%   <http://www.gnu.org/licenses/>.
%
% References:
%
%	R. Gonzalez, J. Giribet, and H. Patiño. NaveGo: a
% simulation framework for low-cost integrated navigation systems,
% Journal of Control Engineering and Applied Informatics, vol. 17,
% issue 2, pp. 110-120, 2015. Sec. 2.2.
%
%   Aggarwal, P. et al. MEMS-Based Integrated Navigation. Artech
% House. 2010.
%
%   Thinking about accelerometers and gravity by Dave Redell
% http://www.lunar.org/docs/LUNARclips/v5/v5n1/Accelerometers.html
%
% Version: 012
% Date:    2022/08/22
% Author:  Rodrigo Gonzalez <rodralez@frm.utn.edu.ar>
% URL:     https://github.com/rodralez/navego

N = max(size(ref.t));
M = [N, 3];

%% SIMULATION OF ACC

% If true, accelerations are provided...
if (isfield(ref, 'fb'))

    acc_b = ref.fb;

% If not, accelerations are obtained from velocity
elseif (isfield(ref, 'vel'))

    acc_raw = (diff(ref.vel)) ./ [diff(ref.t) diff(ref.t) diff(ref.t)];
    acc_raw = [ 0 0 0; acc_raw; ];

    % Noise introduced by derivatives should be smoothed
    acc_ned = my_sgolayfilt(acc_raw);
    acc_b = acc_nav2body(acc_ned, ref.DCMnb_m);

% If not, accelerations are obtained from position
else

    % Method: LLH > ECEF > NED
    [~, acc_ned] = pllh2vned (ref);
    acc_b = acc_nav2body(acc_ned, ref.DCMnb_m);
end

%% SIMULATION OF GRAVITY AND CORIOLIS

% Gravity and Coriolis in nav-ref
g_n = -gravity(ref.lat, ref.h);              % Accelerometer in Z axis senses an
                                                % acceleration of 1.0 G straight up
cor_n = coriolis(ref.lat, ref.vel, ref.h);

% Gravity and Coriolis from nav-ref to body-ref
g_b = zeros(M);
cor_b = zeros(M);
for i = 1:N
    dcmnb = reshape(ref.DCMnb_m(i,:), 3, 3);
    gb = dcmnb * g_n(i,:)';
    corb =  dcmnb * cor_n(i,:)';
    g_b(i,:) = gb';
    cor_b(i,:) = corb';
end

%% SIMULATION OF NOISES

% -------------------------------------------------------------------------
% Simulation of static bias as a constant random variable

ab_sta = noise_b_sta (imu.ab_sta, N);

% -------------------------------------------------------------------------
% Simulation of white noise

wn = randn(M);
acc_wn = zeros(M);

for i=1:3

    acc_wn(:, i) = imu.a_std(i).* wn(:,i);
end

% -------------------------------------------------------------------------
% Simulation of dynamic bias (bias instability) as a first-order Gauss-Markov model

dt = 1/ref.freq;
ab_dyn = noise_b_dyn (imu.ab_corr, imu.ab_dyn, dt, M);

% -------------------------------------------------------------------------
% Simulation of rate random walk

acc_rrw = noise_rrw (imu.vrrw, dt, M);

% -------------------------------------------------------------------------

fb_sim = acc_b + cor_b + g_b + acc_wn + ab_sta + ab_dyn + acc_rrw;

end

author: rodralez

repo: NaveGo

save_path: github-repos/MATLAB/rodralez-NaveGo


text:
function [Iex,Iexmask,Icb,Icbmask] = esvm_get_exemplar_icon(models, ...
                                                  index,flip,subind, ...
                                                  loadseg, VOCopts)
% Extract an exemplar visualization image (one from gt box, one from
% cb box) [and does flip if specified]
% Allows allows for the loading of a segmentation too
% NOTE(TJM): this function works, but needs cleanup
%
% Copyright (C) 2011-12 by Tomasz Malisiewicz
% All rights reserved.
%
% This file is part of the Exemplar-SVM library and is made
% available under the terms of the MIT license (see COPYING file).
% Project homepage: https://github.com/quantombone/exemplarsvm

if ~isfield(models{index}.model,'bb') || numel(models{index}.model.bb)==0
  Iex = ones(10,10,3);
  Iexmask = ones(10,10,1);
  Icb = ones(10,10,3);
  Icbmask = ones(10,10,1);
  return;
end

%Subind indicates which window to show (defaults to the base)
if ~exist('subind','var')
  subind = 1;
%else
%  flip = models{index}.model.bb(subind,7);
end

if ~exist('flip','var')
  flip = 0;
end

if ~exist('loadseg','var')
  loadseg = 1;
end

%cb = models{index}.model.bb(subind,1:4);
%TJM: changed
cb = models{index}.gt_box;
d1 = max(0,1 - cb(1));
d2 = max(0,1 - cb(2));
d3 = max(0,cb(3) - models{index}.sizeI(2));
d4 = max(0,cb(4) - models{index}.sizeI(1));
mypad = max([d1,d2,d3,d4]);
PADDER = round(mypad)+2;

if isfield(models{index},'I')
  I = convert_to_I(models{index}.I);
else
  I = ...
      convert_to_I(models{index}.train_set{models{index}.model.bb(subind,11)});
end
%pre-pad mask because GT region can be outside image
mask = pad_image(zeros(size(I,1),size(I,2)),PADDER);
g = models{index}.gt_box;
g2 = g + PADDER;
mask(g2(2):g2(4),g2(1):g2(3)) = 1;

if loadseg == 1 && exist('VOCopts','var')
  try
    %NOTE(TJM): this is voc-only and should not be here in the most
    %general scenario
    [I2, mask2] = load_seg(VOCopts,models{index});
    if numel(I2) > 0
      I = I2;
      mask = mask2;
    end
  catch
  end
end

cb = models{index}.gt_box;
Iex = pad_image(I, PADDER);
Iexmask = pad_image(mask,PADDER);
cb = round(cb + PADDER);

Iex = Iex(cb(2):cb(4),cb(1):cb(3),:);
Iexmask = Iexmask(cb(2):cb(4),cb(1):cb(3));

cb = models{index}.model.bb(subind,1:4);
Icb = pad_image(I, PADDER);
cb = round(cb + PADDER);
Icb = Icb(cb(2):cb(4),cb(1):cb(3),:);
Icbmask = mask(cb(2):cb(4),cb(1):cb(3));

if flip == 1
  Iex = flip_image(Iex);
  Icb = flip_image(Icb);
  Iexmask = flip_image(Iexmask);
  Icbmask = flip_image(Icbmask);
end

function [I, mask] = load_seg(VOCopts, model)

filer = sprintf('%s/%s/SegmentationObject/%s.png',VOCopts.datadir, ...
                VOCopts.dataset,model.curid);

filer_class = sprintf('%s/%s/SegmentationClass/%s.png',VOCopts.datadir, ...
                      VOCopts.dataset,model.curid);

cmap=VOClabelcolormap;

if ~fileexists(filer)
  I = [];
  mask = [];
  return;
end

classes = {'aeroplane','bicycle','bird','boat','bottle','bus', ...
           'car','cat','chair','cow','diningtable','dog','horse', ...
           'motorbike','person','pottedplant','sheep','sofa', ...
           'train','tvmonitor'};

clsid = find(ismember(classes,model.cls));

has_seg = 1;
res = imread(filer);
res_class = imread(filer_class);
res_class = reshape(cmap(res_class(:)+1,:),size(res_class,1), ...
                    size(res_class,2),3);

mask = double((res==model.objectid));
I= res_class;

author: quantombone

repo: exemplarsvm