text:
function h = draw_circle(x, r, outline_color, fill_color)
% draw filled circles at centers x with radii r.
% x is a matrix of columns.  r is a row vector.

n = 40;					% resolution
radians = [0:(2*pi)/(n-1):2*pi];
unitC = [sin(radians); cos(radians)];

% extend r if necessary
if length(r) < cols(x)
  r = [r repmat(r(length(r)), 1, cols(x)-length(r))];
end

h = [];
% hold is needed for fill()
held = ishold;
hold on
for i=1:cols(x)
  y = unitC*r(i) + repmat(x(:, i), 1, n);
  if nargin < 4
    h = [h line(y(1,:), y(2,:), 'Color', outline_color)];
  else
    h = [h fill(y(1,:), y(2,:), fill_color, 'EdgeColor', outline_color)];
  end
end
if ~held
  hold off
end

author: andrewssobral

repo: lrslibrary


text:
function [qp,QP_ep] = epose2qpose(ep)

% EPOSE2QPOSE  Euler-specified to quaternion-specified pose conversion.
%
%   QP = EPOSE2QPOSE(EP) returns a full 7-pose QP=[X;Q] from a full 6-pose
%   QE=[X;E], where X is 3D opsition and Q and E are 3D orientations.
%
%   [QP,Jep] = EPOSE2QPOSE(...) returns also the Jacobian matrix
%
%   See also QPOSE2EPOSE, EULERANGLES, QUATERNION, FRAME.

%   Copyright 2008-2009 Joan Sola @ LAAS-CNRS.

if any(size(ep) ~= [6,1])
    warning('Input Euler-pose should be a column 6-vector')
end

qp            = zeros(7,1);  % empty quaternion pose
qp(1:3)       = ep(1:3);     % position copy
P             = eye(3);      % Jacobian of copy function
[qp(4:7),Q_e] = e2q(ep(4:6));% orientation and Jacobian

QP_ep         = [P zeros(3);zeros(4,3) Q_e]; % Full Jacobian

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB 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.
%
%   SLAMTB 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 SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function SubmitHitsPlus( input_path, startID, endID, sandbox )
%SUBMITHITS Summary of this function goes here
%   Detailed explanation goes here
load('./data/HITData.mat');
Title = hit.Title.Text;
Description = hit.Description.Text;
Description = [Description ' If some of the images are not loaded, please do NOT continue working and email us to fix them immediately.'];
Reward = hit.Reward;
AssignmentDurationInSeconds = hit.AssignmentDurationInSeconds.Text;
Keywords = hit.Keywords.Text;
QualificationRequirement = hit.QualificationRequirement;

LifetimeInSeconds = '604800';
MaxAssignments = '2';
AutoApprovalDelayInSeconds = '172800';

name = regexp(input_path,'/','split');
name2 =regexp(name{end},'\.','split');
RequesterAnnotation = ['Batch:' name2{1}];
if strcmp(input_path(1:6),'/n/fs/')
    parameter = gendataLsunNew(name2{1});
else
    parameter = gendataLsunNew(input_path);
end

fprintf('%d hits pending submission: $%3.2f\n', length(parameter), length(parameter)*2*0.2);
pause(60);

aws_access_key = 'AKIAIAGDVEVNP3ML45RQ';
aws_secret_key = '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';
% sandbox = false;
% sandbox = true;
turk = InitializeTurk(aws_access_key, aws_secret_key, sandbox);

%%
for bid = startID:min(endID,length(parameter))
    fprintf('%d/%d\n', bid, length(parameter));
    success = false;
    Question2 = sprintf(QuestionPattern2, parameter(bid).DataStr(2:end-1), parameter(bid).exampleStr(2:end-1), ['''' parameter(bid).definitionStr '''']);
    Question = [QuestionPattern1 Question2 QuestionPattern3 parameter(bid).QuestionStr QuestionPattern4];

    while ~success
        try
            if ~sandbox
            status = RequestTurk(turk, 'CreateHIT', {'Title',Title,'Description',Description,...
                        'Reward.Amount',Reward.Amount.Text,...
                        'Reward.CurrencyCode',Reward.CurrencyCode.Text,...
                        'Reward.FormattedPrice',Reward.FormattedPrice.Text,...
                        'AssignmentDurationInSeconds',AssignmentDurationInSeconds,...
                        'LifetimeInSeconds',LifetimeInSeconds,'MaxAssignments',MaxAssignments,...
                        'Keywords',Keywords,'AutoApprovalDelayInSeconds',AutoApprovalDelayInSeconds,...
                        'RequesterAnnotation',RequesterAnnotation,...
                        'QualificationRequirement.1.QualificationTypeId',QualificationRequirement{1}.QualificationTypeId.Text, ...
                        'QualificationRequirement.1.Comparator',QualificationRequirement{1}.Comparator.Text, ...
                        'QualificationRequirement.1.IntegerValue',QualificationRequirement{1}.IntegerValue.Text, ...
                        'QualificationRequirement.1.RequiredToPreview',QualificationRequirement{1}.RequiredToPreview.Text, ...
                        'QualificationRequirement.2.QualificationTypeId',QualificationRequirement{2}.QualificationTypeId.Text, ...
                        'QualificationRequirement.2.Comparator',QualificationRequirement{2}.Comparator.Text, ...
                        'QualificationRequirement.2.IntegerValue',QualificationRequirement{2}.IntegerValue.Text, ...
                        'QualificationRequirement.2.RequiredToPreview',QualificationRequirement{2}.RequiredToPreview.Text, ...
                        'Question',Question});
            else
                status = RequestTurk(turk, 'CreateHIT', {'Title',Title,'Description',Description,...
                            'Reward.Amount',Reward.Amount.Text,...
                            'Reward.CurrencyCode',Reward.CurrencyCode.Text,...
                            'Reward.FormattedPrice',Reward.FormattedPrice.Text,...
                            'AssignmentDurationInSeconds',AssignmentDurationInSeconds,...
                            'LifetimeInSeconds',LifetimeInSeconds,'MaxAssignments',MaxAssignments,...
                            'Keywords',Keywords,'AutoApprovalDelayInSeconds',AutoApprovalDelayInSeconds,...
                            'RequesterAnnotation',RequesterAnnotation,...
                            'Question',Question});
            end

            if strcmp(status.CreateHITResponse.HIT.Request.IsValid.Text, 'True')
                success = true;
            else
                fprintf('%s\n', status.CreateHITResponse.HIT.Request.Errors.Error.Message.Text);
            end
        end
    end
end

%%
if ~sandbox
    if ~strcmp(input_path(end-3:end),'.txt')
        input_path = [input_path '.txt'];
    end
    new_batch_path = './data/New_Batches.txt';
    fp = fopen(new_batch_path,'a');
    fprintf(fp,'%s %s\n', RequesterAnnotation, input_path);
    fclose(fp);
end

end

author: jianxiongxiao

repo: ProfXkit


text:
function cstr= str_tail(cstr,varargin)
%tail= strtail(str,<delim>)
%
% returns the string after (not including) the first space
% works also for cells of strings
% 'delim' -  delimiter of the string head, e.g. '_' (default ' ')
%
% See also strhead

if nargin<2
    delim = ' ';
else
    delim = varargin{1};
end

if ~iscell(cstr),
  cstr= str_tail({cstr},delim);
  cstr= cstr{1};
  return
end

for il= 1:length(cstr),
  cstr{il}= cstr{il}(max([1 findstr(delim, cstr{il})+1]):end);
end

author: bbci


text:
function LU = getbounds(F,avoidequalitybounds,LU)

K.f = 0;
K.l = 0;
L = [];
U = [];
if nargin < 3
    LU =  yalmip('getbounds',1:yalmip('nvars'));
end
binary = yalmip('binvariables');
LU(binary,1) = 0;
LU(binary,2) = 1;
F = flatten(F);
simplex = [];
is_interval = is(F,'interval');
for i = 1:length(F.LMIid)
    if F.clauses{i}.type == 2
        X = F.clauses{i}.data;
        AB = getbase(X);
        K.l = prod(size(X));
        variables = getvariables(X);
        if is_interval(i)
            [lb,ub,cand_rows] = find_lp_bounds_interval(AB,K);
        else
            [lb,ub,cand_rows] = find_lp_bounds(AB,K);
        end
        LU(variables,1) = max([lb LU(variables,1)]')';
        LU(variables,2) = min([ub LU(variables,2)]')';
    elseif  F.clauses{i}.type == 55
        X = F.clauses{i}.data;X = X(:);
        AB = getbase(X);
        K.l = prod(size(X));
        variables = getvariables(X);
        if is_interval(i)
            [lb,ub,cand_rows] = find_lp_bounds_interval(AB,K);
        else
            [lb,ub,cand_rows] = find_lp_bounds(AB,K);
        end
        LU(variables,1) = max([lb LU(variables,1)]')';
        LU(variables,2) = min([ub LU(variables,2)]')';

    elseif F.clauses{i}.type == 3 && (nargin==1 || (nargin > 1 && isempty(avoidequalitybounds)))
        % FIX : Extract from equalities and binary constraints
        X = F.clauses{i}.data;
        AB = getbase(X);
        if size(AB,1)==1 && (AB(1) > 0) && all(AB(2:end)<0)
            % Save this simplex sum a_ix_i == b for later
            % and extract when as many variables as possible have been
            % fixed. Common in portfolio models
            simplex = [simplex i];
        end
        AB = [AB;-AB];
        K.l = 2*prod(size(X));
        variables = getvariables(X);
        if is_interval(i)
            [lb,ub,cand_rows] = find_lp_bounds_interval(AB,K);
        else
            [lb,ub,cand_rows] = find_lp_bounds(AB,K);
        end
        LU(variables,1) = max([lb LU(variables,1)]')';
        LU(variables,2) = min([ub LU(variables,2)]')';
    end
end

semicont = yalmip('semicontvariables');
if ~isempty(semicont)
    for i = 1:length(semicont)
        if LU(semicont(i),1) > 0
            LU(semicont(i),1) = 0;
        end
        if LU(semicont(i),2) < 0
            LU(semicont(i),2) = 0;
        end
    end
end

% Go through the simplex models with as much bounded as possible
for i = simplex
    X = F.clauses{i}.data;
    AB = getbase(X);
    variables = getvariables(X);
    if all(LU(variables,1)>=0)
        simplex_bound = AB(1)./(-AB(2:end));
        LU(variables,2) = min([simplex_bound(:) LU(variables,2)]')';
    end
end

% Try to bound some nonlinear terms
% FIX: complete code
[mt,variable_type] = yalmip('monomtable');
quadratic = find(variable_type == 2);
if ~isempty(quadratic)
    M = mt(quadratic,:);
    [ii,jj,kk] = find(M);
    LU(quadratic,1) = max([LU(quadratic,1) zeros(length(quadratic),1)],[],2);
    LU(quadratic,2) = max([LU(quadratic,1).^2 LU(quadratic,2).^2],[],2);
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [labels, distances] = elec_distance(elec_labels,X,Y,Z,xo,yo,zo)

% elec_distance - Calculates spherical interelectrode distances (arc length).
%
% Usage: [labels, distances] = elec_distance(elec_labels,X,Y,Z,xo,yo,zo)
%
% Notes:    Arc length method assumes a centroid at (xo,yo,zo) = (0,0,0)
%           and input values are in rectangular Cartesian coordinates.
%
%           Sphere radius is estimated by the average radius from
%           (xo,yo,zo) to any 2 pairs of electrodes.  It will vary
%           from one pair to another, but this method works well for
%           small theta (eg, nearest neighbours).
%
%           Returns 2 matrices, one for paired electrode labels and
%           another for spherical arc length estimates.
%

% $Revision: 1.1 $ $Date: 2009-04-28 22:13:54 $

% Author:   Darren.Weber_at_radiology.ucsf.edu
% Created:  18/05/00 - linear distance
% Modified: 24/06/01 - spherical arc length (should be OK for small theta
%                      but for large theta, elliptical arc length may be
%                      preferable).
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    % initialise centroid, unless input parameters defined
    if ~exist('xo','var')  xo = 0;   end
    if ~exist('yo','var')  yo = 0;   end
    if ~exist('zo','var')  zo = 0;   end

    fprintf('%s\n', 'Calculating inter-electrode spherical arc length.');

    labels = cell(1);
    distances = [];

    % convert numerical electrode labels to cell strings
    if isnumeric(elec_labels) > 0
        if size(elec_labels,1) > size(elec_labels,2)
            elec_labels = cellstr(strcat(num2str(elec_labels)));
        else
            elec_labels = cellstr(strcat(num2str(elec_labels')));
        end
    end

    A = [ (X-xo) (Y-yo) (Z-zo) ];  B = A;              % define electrode vectors A,B, given centroid (xo,yo,zo)
    rows = 50;                                         % progress indicator
    for a =   1:length(X)
        fprintf('.'); if( a == rows ) fprintf('\n'); rows = rows + 50; end  % progress indicator

        Aa = A(a,:);
        A_len = sqrt( sum( Aa.^2 ) );                        % length of electrode vector A

        for b = 1:length(X)

            Bb = B(b,:);
            B_len = sqrt ( sum(Bb.^2) );                     % length of electrode vector B

            if( Aa == Bb )
                arc_len = 0;                                    % no distance from electrode to itself
            else
                r = (A_len + B_len)/2;                          % estimate sphere radius from A_len and B_len
                theta = acos( dot(Aa,Bb) / (A_len * B_len) );   % Angle between A & B, in radians
                arc_len = r * theta;                            % arc length = radius * theta
            end

            distances(a, b) = arc_len;
            labels(a, b) = strcat(elec_labels(a), ', ', elec_labels(b));
        end
    end
    fprintf('\n');

author: PatternRecognition

repo: OpenBMI


text:
function eventLog = ma_executeCoast_goto_func_value(event, initialState, eventNum, considerSoITransitions, soiSkipIds, massLoss, orbitDecay, maData, celBodyData)
%ma_executeCoast_goto_func_value Summary of this function goes here
%   Detailed explanation goes here
    global number_state_log_entries_per_coast;

    old_number_state_log_entries_per_coast = number_state_log_entries_per_coast;
    number_state_log_entries_per_coast = 5; %#ok<NASGU>

    bnds = [0, event.maxPropTime];
    funcNoAbs = @(dt) getCoastFuncValue(dt, event.funcHandle, event.coastToValue, initialState, eventNum, considerSoITransitions, soiSkipIds, massLoss, orbitDecay, maData, celBodyData, false);
    funcAbs = @(dt) getCoastFuncValue(dt, event.funcHandle, event.coastToValue, initialState, eventNum, considerSoITransitions, soiSkipIds, massLoss, orbitDecay, maData, celBodyData, true);

    odeFunc = @(t,y) odefun(t,y);
    evtFnc = @(t,y) odeEventsFun(t,y,funcNoAbs);
    options = odeset('RelTol',1E-6, 'AbsTol',1E-6, 'Events',evtFnc, 'MaxStep',abs(diff(bnds))/50);
    [~,~,te,~,ie] = ode45(odeFunc, bnds, bnds(1), options);

    if(isempty(ie)) %fall back to old way of doing things
        options = optimset('TolX',1E-6);
        [x, fval, exitflag] = fminbnd(funcAbs, bnds(1), bnds(2), options);

        if(exitflag == 0 || abs(fval) > 1E-4)
            x = event.maxPropTime;
        end
    else
        x = te(1);
    end

    number_state_log_entries_per_coast = old_number_state_log_entries_per_coast;
    eventLog = ma_executeCoast_goto_dt(x, initialState, eventNum, considerSoITransitions, soiSkipIds, massLoss, orbitDecay, celBodyData);
end

function ydot = odefun(~,~)
    ydot = 1;
end

function [value,isterminal,direction] = odeEventsFun(t,~, func)
    value = func(t);
    isterminal = 1;
    direction = 0;
end

function value = getCoastFuncValue(dt, funcHandle, desiredValue, initialState, eventNum, considerSoITransitions, soiSkipIds, massLoss, orbitDecay, maData, celBodyData, useAbs)
    eventLog = ma_executeCoast_goto_dt(dt, initialState, eventNum, considerSoITransitions, soiSkipIds, massLoss, orbitDecay, celBodyData);
    finalState = eventLog(end,:);

    [datapt, ~, ~, taskStr, refBodyInfo, otherSC, station] = funcHandle(finalState, false, maData);

    if(useAbs)
        value = abs(datapt - desiredValue);
    else
        value = datapt - desiredValue;
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
% Like a durational HMM, except we use soft evidence on the observed nodes.
% Should give the same results as HSMM/test_mgram2.

past = 1;
% If past=1, P(Yt|Qt=j,Dt=d) = P(y_{t-d+1:t}|j)
% If past=0, P(Yt|Qt=j,Dt=d) = P(y_{t:t+d-1}|j) - future evidence

words = {'the', 't', 'h', 'e'};
data = 'the';
nwords = length(words);
word_len = zeros(1, nwords);
word_prob = normalise(ones(1,nwords));
word_logprob = log(word_prob);
for wi=1:nwords
  word_len(wi)=length(words{wi});
end
D = max(word_len);

alphasize = 26*2;
data = letter2num(data);
T = length(data);

% node numbers
W = 1; % top level state = word id
L = 2; % bottom level state = letter position within word
F = 3;
O = 4;

ss = 4;
intra = zeros(ss,ss);
intra(W,[F L O])=1;
intra(L,[O F])=1;

inter = zeros(ss,ss);
inter(W,W)=1;
inter(L,L)=1;
inter(F,[W L O])=1;

% node sizes
ns = zeros(1,ss);
ns(W) = nwords;
ns(L) = D;
ns(F) = 2;
ns(O) = alphasize;
ns2 = [ns ns];

% Make the DBN
bnet = mk_dbn(intra, inter, ns, 'observed', O);
eclass = bnet.equiv_class;

% uniform start distrib over words, uniform trans mat
Wstart = normalise(ones(1,nwords));
Wtrans = mk_stochastic(ones(nwords,nwords));
%Wtrans = ones(nwords,nwords);

% always start in state d = length(word) for each bottom level HMM
Lstart = zeros(nwords, D);
for i=1:nwords
  l = length(words{i});
  Lstart(i,l)=1;
end

% make downcounters
RLtrans = mk_rightleft_transmat(D, 0); % 0 self loop prob
Ltrans = repmat(RLtrans, [1 1 nwords]);

% Finish when downcoutner = 1
Fprob = zeros(nwords, D, 2);
Fprob(:,1,2)=1;
Fprob(:,2:end,1)=1;

% Define CPDs for slice 1
bnet.CPD{eclass(W,1)} = tabular_CPD(bnet, W, 'CPT', Wstart);
bnet.CPD{eclass(L,1)} = tabular_CPD(bnet, L, 'CPT', Lstart);
bnet.CPD{eclass(F,1)} = tabular_CPD(bnet, F, 'CPT', Fprob);

% Define CPDs for slice 2
bnet.CPD{eclass(W,2)} = hhmmQ_CPD(bnet, W+ss, 'Fbelow', F, 'startprob', Wstart,  'transprob', Wtrans);
bnet.CPD{eclass(L,2)} = hhmmQ_CPD(bnet, L+ss, 'Fself', F, 'Qps', W+ss, 'startprob', Lstart, 'transprob', Ltrans);

if 0
% To test it is generating correctly, we create an artificial
% observation process that capitalizes at the start of a new segment
% Oprob(Ft-1,Qt,Dt,Yt)
Oprob = zeros(2,nwords,D,alphasize);
Oprob(1,1,3,letter2num('t'),1)=1;
Oprob(1,1,2,letter2num('h'),1)=1;
Oprob(1,1,1,letter2num('e'),1)=1;
Oprob(2,1,3,letter2num('T'),1)=1;
Oprob(2,1,2,letter2num('H'),1)=1;
Oprob(2,1,1,letter2num('E'),1)=1;
Oprob(1,2,1,letter2num('a'),1)=1;
Oprob(2,2,1,letter2num('A'),1)=1;
Oprob(1,3,1,letter2num('b'),1)=1;
Oprob(2,3,1,letter2num('B'),1)=1;
Oprob(1,4,1,letter2num('c'),1)=1;
Oprob(2,4,1,letter2num('C'),1)=1;

% Oprob1(Qt,Dt,Yt)
Oprob1 = zeros(nwords,D,alphasize);
Oprob1(1,3,letter2num('t'),1)=1;
Oprob1(1,2,letter2num('h'),1)=1;
Oprob1(1,1,letter2num('e'),1)=1;
Oprob1(2,1,letter2num('a'),1)=1;
Oprob1(3,1,letter2num('b'),1)=1;
Oprob1(4,1,letter2num('c'),1)=1;

bnet.CPD{eclass(O,2)} = tabular_CPD(bnet, O+ss, 'CPT', Oprob);
bnet.CPD{eclass(O,1)} = tabular_CPD(bnet, O, 'CPT', Oprob1);

evidence = cell(ss,T);
%evidence{W,1}=1;
sample = cell2num(sample_dbn(bnet, 'length', T, 'evidence', evidence));
str = num2letter(sample(4,:))
end

if 1

[log_obslik, obslik, match] = mk_mgram_obslik(lower(data), words, word_len, word_prob);
% obslik(j,t,d)
softCPDpot = cell(ss,T);
ens = ns;
ens(O)=1;
ens2 = [ens ens];
for t=2:T
  dom = [F W+ss L+ss O+ss];
  % tab(Ft-1, Q2, Dt)
  tab = ones(2, nwords, D);
  if past
    tab(1,:,:)=1; % if haven't finished previous word, likelihood is 1
    %tab(2,:,:) = squeeze(obslik(:,t,:)); % otherwise likelihood of this segment
    for d=1:min(t,D)
      tab(2,:,d) = squeeze(obslik(:,t,d));
    end
  else
    for d=1:max(1,min(D,T+1-t))
      tab(2,:,d) = squeeze(obslik(:,t+d-1,d));
    end
  end
  softCPDpot{O,t} = dpot(dom, ens2(dom), tab);
end
t = 1;
dom = [W L O];
% tab(Q2, Dt)
tab = ones(nwords, D);
if past
  %tab = squeeze(obslik(:,t,:));
  tab(:,1) = squeeze(obslik(:,t,1));
else
  for d=1:min(D,T-t)
    tab(:,d) = squeeze(obslik(:,t+d-1,d));
  end
end
softCPDpot{O,t} = dpot(dom, ens(dom), tab);

%bnet.observed = [];
% uniformative observations
%bnet.CPD{eclass(O,2)} = tabular_CPD(bnet, O+ss, 'CPT', mk_stochastic(ones(2,nwords,D,alphasize)));
%bnet.CPD{eclass(O,1)} = tabular_CPD(bnet, O, 'CPT', mk_stochastic(ones(nwords,D,alphasize)));

engine = jtree_dbn_inf_engine(bnet);
evidence = cell(ss,T);
% we add dummy data to O to force its effective size to be 1.
% The actual values have already been incorporated into softCPDpot
evidence(O,:) = num2cell(ones(1,T));
[engine, ll_dbn] = enter_evidence(engine, evidence, 'softCPDpot', softCPDpot);

%evidence(F,:) = num2cell(2*ones(1,T));
%[engine, ll_dbn] = enter_evidence(engine, evidence);

gamma = zeros(nwords, T);
for t=1:T
  m = marginal_nodes(engine, [W F], t);
  gamma(:,t) = m.T(:,2);
end

gamma

xidbn = zeros(nwords, nwords);
for t=1:T-1
  m = marginal_nodes(engine, [W F W+ss], t);
  xidbn = xidbn + squeeze(m.T(:,2,:));
end

% thee
% xidbn(1,4)  = 0.9412  the->e
% (2,3)=0.0588 t->h
% (3,4)=0.0588 h-e
% (4,4)=0.0588 e-e

end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
% evaluate MPJPE up to aligning the reference body joint (pelvis)
function error = MPJPE_h36(joint,j_p)
        % joint: ground truth 3D pose
        % j_p: prediction 3D pose

        %%
        %h36m

        jroot=[(joint(9,1)+joint(12,1))/2, (joint(9,2)+joint(12,2))/2,(joint(9,3)+joint(12,3))/2];

        %jroot=mean(joint);

        lsum = joint-repmat(jroot,14,1);
        %Lsum = mean(sqrt(lsum(:,1).^2+lsum(:,2).^2+lsum(:,3).^2));

        jroot_p = [(j_p(9,1)+j_p(12,1))/2, (j_p(9,2)+j_p(12,2))/2,(j_p(9,3)+j_p(12,3))/2];
        %jroot_p = mean(j_p);
        lsu = j_p - repmat(jroot_p,14,1);
        %Lsu = mean(sqrt(lsu(:,1).^2+lsu(:,2).^2+lsu(:,3).^2));

        %scale = Lsum/Lsu;
        temp = lsum - lsu;%*scale;

        error = mean(sqrt(temp(:,1).^2+temp(:,2).^2+temp(:,3).^2));

end

author: flyawaychase


text:
% function Z = ZG_rdiv(X,Y)
%
% row division: Z = X / Y row-wise
% Y must have one column
%
% Machine Learning Toolbox
% Version 1.0  01-Apr-96
% Copyright (c) by Zoubin Ghahramani
% http://mlg.eng.cam.ac.uk/zoubin/software.html
%
% ------------------------------------------------------------------------------
% The MIT License (MIT)
%
% Copyright (c) 1996, Zoubin Ghahramani
%
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, including without limitation the rights
% to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
% copies of the Software, and to permit persons to whom the Software is
% furnished to do so, subject to the following conditions:
%
% The above copyright notice and this permission notice shall be included in
% all copies or substantial portions of the Software.
%
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
% IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
% FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
% AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
% LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
% OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
% THE SOFTWARE.
% ------------------------------------------------------------------------------

function Z = ZG_rdiv(X,Y)

if (length(X(:,1)) ~= length(Y(:,1))) || (length(Y(1,:)) ~= 1)
    error('Error in ZG_RDIV');
end

Z = zeros(size(X));

for i = 1:length(X(1,:))
    Z(:,i) = X(:,i)./Y;
end

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function [p1,p2,a,a1,a2,a3,h]=pigeon()
% PIGEON Draw a pigeon
% numandina@gmail.com
%

l=1; % pigeon size
hold on
axis([-5 5 -5 5 -5 5])
view(3)
grid on
rotate3d

% body
[a,b,c]=sphere(25);
aa1=surf(a.*l,b.*l,c.*l);
set(aa1,'edgecolor','none','facecolor',[.5 .5 .5])

% head
[a,b,c]=sphere(25);
aa2=surf(a.*.5.*l,b.*.5.*l+.5.*l,c.*.5.*l+1.*l);
set(aa2,'edgecolor','none','facecolor',[.5 .5 .5])

% eyes
[a,b,c]=sphere(25);
aa2=surf(a.*.1.*l,b.*.1.*l+.875.*l,c.*.1.*l+1.25.*l);
set(aa2,'edgecolor','none','facecolor',[0 0 0])
[a,b,c]=sphere(25);
aa2=surf(a.*.1.*l+.2,b.*.1.*l+.85.*l,c.*.1.*l+1.25.*l);
set(aa2,'edgecolor','none','facecolor',[0 0 0])

% beak
Cone2 = Cone([0 0.5 1].*l ,[0.5 1.5 0.5].*l,[.25 0].*l,25,[.9 .5 .5]);

% legs
l1=line([0 0]-.1,[-.2 -.2].*l,[-.5 -1.5].*l);
set(l1,'linewidth',2,'color',[.9 .5 .5])
l2=line([0 0]+.1,[-.2 -.2].*l,[-.5 -1.5].*l);
set(l2,'linewidth',2,'color',[.9 .5 .5])

a= [.5 0 0;
	.5 0 0;
	1.5 1 0;
	1.5 1 0;
	.5 0 1;
	.5 0 1;
	.5 1 1;
	.5 1 1].*l;

a(:,2)=a(:,2)-.5.*l;
a(:,3)=a(:,3)-.5.*l;
a(:,1)=a(:,1)+0.2.*l;

b= [1 2 6 5;
	2 3 7 6;
	3 4 8 7;
	4 1 5 8;
	1 2 3 4;
	5 6 7 8];

% wings, first position
p1=patch('vertices',a,'faces',b,'edgecolor',[.5 .5 .5],'facecolor',[.9 .9 .9]);
a1=a;
a1(:,1)=-a1(:,1);
p2=patch('vertices',a1,'faces',b,'edgecolor',[.5 .5 .5],'facecolor',[.9 .9 .9]);

% second position of wings
a2= [.5 0 2;
	.5 0 2;
	1.5 1 2;
	1.5 1 2;
	.5 0 1;
	.5 0 1;
	.5 1 1;
	.5 1 1].*l;
a2(:,2)=a2(:,2)-.5.*l;
a2(:,3)=a2(:,3)-.5.*l;
a2(:,1)=a2(:,1)+.2.*l;
a3=a2;
a3(:,1)=-a3(:,1);

% animate wings
while true
	set(p1,'vertices',a2)
	set(p2,'vertices',a3)
	pause(.25)
	set(p1,'vertices',a)
	set(p2,'vertices',a1)
	pause(.25)
end

% get everything's handle
h=[aa1,aa3,Cone2,EndPlate1,EndPlate2,l1,l2];

% edited from a file on the FEX, create a cone
	function Cone2 = Cone(X1,X2,R,n,cyl_color)
		length_cyl=norm(X2-X1);
		t=linspace(0,2*pi,n)';
		x1=[0 length_cyl];
		xx1=repmat(x1,n,1);
		xx2=[R(1)*cos(t) R(2)*cos(t)];
		xx3=[R(1)*sin(t) R(2)*sin(t)];
		Cone2=mesh(xx1,xx2,xx3);
		unit_Vx=[1 0 0];
		angle_X1X2=acos( dot( unit_Vx,(X2-X1) )/( norm(unit_Vx)*norm(X2-X1)) )*180/pi;
		axis_rot=cross([1 0 0],(X2-X1) );
		if angle_X1X2~=0
			rotate(Cone2,axis_rot,angle_X1X2,[0 0 0])
		end
		set(Cone2,'XData',get(Cone2,'XData')+X1(1),'facecolor',cyl_color)
		set(Cone2,'YData',get(Cone2,'YData')+X1(2))
		set(Cone2,'ZData',get(Cone2,'ZData')+X1(3))
		set(Cone2,'EdgeAlpha',0)
	end
end

author: Sable

repo: mcbench-benchmarks


text:
function GPIB
hMainFigure = figure(...       % The main GUI figure
                    'MenuBar','none', ...
                    'Toolbar','none', ...
                    'HandleVisibility','callback', ...
                    'Color', get(0,...
                             'defaultuicontrolbackgroundcolor'));

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

global eje;
global inicial;
global final;
global paso;
global hcontrol;
global ESP300;
global hmanual;
global hautomatico;
global hvar;

clc

h = uicontrol(hMainFigure,'Style', 'pushbutton', 'String', 'Mover',...
    'Position', [400 50 100 20],'Callback',@Mover);
%--------------------------------------------------------------------------
hcerrar = uicontrol(hMainFigure,'Style', 'pushbutton', 'String', 'Cerrar',...
    'Position', [50 50 100 20],'Callback',@Cerrar);
%--------------------------------------------------------------------------
uicontrol(hMainFigure,'Style', 'text','String', 'Eje',...
    'Position', [20 365 80 20])
heje = uicontrol(hMainFigure,'Style', 'edit','String', '1',...
    'Position', [20 350 80 20],'Callback',@heje_call);
eje = str2double(get(heje,'String'));
%--------------------------------------------------------------------------
uicontrol(hMainFigure,'Style', 'text','String', 'Inicial',...
    'Position', [20 315 80 20])
hinicial = uicontrol(hMainFigure,'Style', 'edit', 'String', '0',...
    'Position', [20 300 80 20],'Callback',@hinicial_call);
inicial = str2double(get(hinicial,'String'));
%--------------------------------------------------------------------------
uicontrol(hMainFigure,'Style', 'text','String', 'Final',...
    'Position', [20 265 80 20])
hfinal = uicontrol(hMainFigure,'Style', 'edit', 'String', '0',...
    'Position', [20 250 80 20],'Callback',@hfinal_call);
final = str2double(get(hfinal,'String'));
%--------------------------------------------------------------------------
uicontrol(hMainFigure,'Style', 'text','String', 'Paso',...
    'Position', [20 215 80 20])
hpaso = uicontrol(hMainFigure,'Style', 'edit', 'String', '0.05',...
    'Position', [20 200 80 20],'Callback',@hpaso_call);
paso = str2double(get(hpaso,'String'));
%--------------------------------------------------------------------------
uicontrol(hMainFigure,'Style', 'text','String', 'Control',...
    'Position', [20 165 80 20])
hcontrol = uicontrol(hMainFigure,'Style', 'edit',...
    'Position', [20 150 80 20]);
%--------------------------------------------------------------------------

%--------------------------------------------------------------------------
hmanual = uicontrol(hMainFigure,'Style', 'radiobutton',...
    'String', 'Manual',...
    'Position', [200 350 100 20],...
    'Callback',@Manual);
%--------------------------------------------------------------------------
hautomatico = uicontrol(hMainFigure,'Style', 'radiobutton',...
    'String', 'Automatico',...
    'Position', [200 320 100 20],...
    'Callback',@Automatico);
%--------------------------------------------------------------------------
heje1 = uicontrol(hMainFigure,'Style', 'radiobutton',...
    'String', 'Eje 1',...
    'Position', [300 350 100 20],...
    'Callback',@Eje1);
%--------------------------------------------------------------------------
heje2 = uicontrol(hMainFigure,'Style', 'radiobutton',...
    'String', 'Eje 2',...
    'Position', [300 320 100 20],...
    'Callback',@Eje2);
%--------------------------------------------------------------------------
heje3 = uicontrol(hMainFigure,'Style', 'radiobutton',...
    'String', 'Eje 3',...
    'Position', [300 290 100 20],...
    'Callback',@Eje3);
%--------------------------------------------------------------------------
ESP300 = gpib('ni',0,1);
ESP300.Timeout = 0.5;
fopen(ESP300);

function Mover(hObject, eventdata)

global eje
global inicial
global final
global paso
global hcontrol;
global ESP300;
global hmanual;
global hautomatico;

%--------------------------------------------
if (get(hmanual,'Value') == get(hmanual,'Max'))
    MovManual(eje,inicial,final,paso,ESP300,hcontrol);
end
if (get(hautomatico,'Value') == get(hautomatico,'Max'))
    MovAutomatico(eje,inicial,final,paso,ESP300,hcontrol);
end

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

function [PosFinal] = MovManual(eje,inicial,final,paso,ESP300,hcontrol)
ESP300.EOIMode = 'on';
ESP300.EOSMode = 'read';
ESP300.EOSCharCode = 'X';

Send = ([num2str(eje),'PA',num2str(final)]);
fprintf(ESP300,Send);

%--------------------------------------------
fprintf(ESP300,'TS')
Dato1a  = fscanf(ESP300,'%s');

while ~(strcmp(Dato1a,'P'))
    fprintf(ESP300,'TS')
    Dato1a  = fscanf(ESP300,'%s')
end

if Dato1a == 'P'

    Send1 = ([num2str(eje),'TP']);
    fprintf(ESP300,Send1);
    Dato1 = fscanf(ESP300,'%s')
   % [Dato1 Resi1] = query(ESP300,Send1,'%s','%s');
    set(hcontrol,'String',num2str(Dato1));
    pause(0.5)
end

function MovAutomatico(eje,inicial,final,paso,ESP300,hcontrol)

var = (inicial)

while final+paso >= var

    MovManual(eje,inicial,var,paso,ESP300,hcontrol)
    var = var+paso
end

function Eje1(hObject, eventdata)

global ESP300;

ESP300.EOIMode = 'on';
ESP300.EOSMode = 'read';
ESP300.EOSCharCode = 'X';

if (get(hObject,'Value') == get(hObject,'Max'))
	Send1 = (['1','MO']);
    fprintf(ESP300,Send1);
    set(hObject,'String','Eje 1 ON');
else
    Send1 = (['1','MF']);
    fprintf(ESP300,Send1);
	set(hObject,'String','Eje 1');
end

function Eje2(hObject, eventdata)

global ESP300;

ESP300.EOIMode = 'on';
ESP300.EOSMode = 'read';
ESP300.EOSCharCode = 'X';

if (get(hObject,'Value') == get(hObject,'Max'))
	Send1 = (['2','MO']);
    fprintf(ESP300,Send1);
    set(hObject,'String','Eje 2 ON');
else
    Send1 = (['2','MF']);
    fprintf(ESP300,Send1);
	set(hObject,'String','Eje 2');
end

function Eje3(hObject, eventdata)

global ESP300;

ESP300.EOIMode = 'on';
ESP300.EOSMode = 'read';
ESP300.EOSCharCode = 'X';

if (get(hObject,'Value') == get(hObject,'Max'))
	Send1 = (['3','MO']);
    fprintf(ESP300,Send1);
    set(hObject,'String','Eje 3 ON');
else
    Send1 = (['3','MF']);
    fprintf(ESP300,Send1);
	set(hObject,'String','Eje 3');
end

function Cerrar(hObject, eventdata)

global ESP300;

fclose(ESP300);
delete(ESP300);
close all

function Manual(hObject, eventdata)

global hautomatico;

if (get(hObject,'Value') == get(hObject,'Max'))
	set(hautomatico,'Value',0);
else
	% Radio button is not selected-take approriate action
end

function Automatico(hObject, eventdata)

global hmanual;

if (get(hObject,'Value') == get(hObject,'Max'))
	set(hmanual,'Value',0);
else
	% Radio button is not selected-take approriate action
end

function heje_call(hObject, eventdata)
global eje
eje = str2double(get(hObject,'String'));

function hinicial_call(hObject, eventdata)
global inicial
inicial = str2double(get(hObject,'String'));

function hfinal_call(hObject, eventdata)
global final
final = str2double(get(hObject,'String'));

function hpaso_call(hObject, eventdata)
global paso
paso = str2double(get(hObject,'String'));


text:
function L = dist2poly ( p, edgexy, lim )

%*****************************************************************************80
%
% Find the minimum distance from the points in P to the polygon defined by
% the edges in EDGEXY. LIM is an optional argument that defines an upper
% bound on the distance for each point.

% Uses (something like?) a double sweep-line approach to reduce the number
% of edges that are required to be tested in order to determine the closest
% edge for each point. On average only size(EDGEXY)/4 comparisons need to
% be made for each point.
%
%  Author:
%
%    Darren Engwirda
%
if nargin<3
   lim = [];
end
np = size(p,1);
ne = size(edgexy,1);
if isempty(lim)
   lim = inf*ones(np,1);
end

% Choose the direction with the biggest range as the "y-coordinate" for the
% test. This should ensure that the sorting is done along the best
% direction for long and skinny problems wrt either the x or y axes.
dxy = max(p)-min(p);
if dxy(1)>dxy(2)
    % Flip co-ords if x range is bigger
    p       = p(:,[2,1]);
    edgexy  = edgexy(:,[2,1,4,3]);
end

% Ensure edgexy(:,[1,2]) contains the lower y value
swap           = edgexy(:,4)<edgexy(:,2);
edgexy(swap,:) = edgexy(swap,[3,4,1,2]);

% Sort edges
[i,i]          = sort(edgexy(:,2));                                        % Sort edges by lower y value
edgexy_lower   = edgexy(i,:);
[i,i]          = sort(edgexy(:,4));                                        % Sort edges by upper y value
edgexy_upper   = edgexy(i,:);

% Mean edge y value
ymean = 0.5*( sum(sum(edgexy(:,[2,4]))) )/ne;

% Alloc output
L = zeros(np,1);

% Loop through points
tol = 1000.0*eps*max(dxy);
for k = 1:np

   x = p(k,1);
   y = p(k,2);
   d = lim(k);

   if y<ymean

      % Loop through edges bottom up
      for j = 1:ne
         y2 = edgexy_lower(j,4);
         if y2>=(y-d)
            y1 = edgexy_lower(j,2);
            if y1<=(y+d)

               x1 = edgexy_lower(j,1);
               x2 = edgexy_lower(j,3);

               if x1<x2
                  xmin = x1;
                  xmax = x2;
               else
                  xmin = x2;
                  xmax = x1;
               end

               if xmin<=(x+d) && xmax>=(x-d)
                  % Calculate the distance along the normal projection from [x,y] to the jth edge
                  x2mx1 = x2-x1;
                  y2my1 = y2-y1;

                  r = ((x-x1)*x2mx1+(y-y1)*y2my1)/(x2mx1^2+y2my1^2);
                  if r>1.0                                                 % Limit to wall endpoints
                     r = 1.0;
                  elseif r<0.0
                     r = 0.0;
                  end

                  dj = (x1+r*x2mx1-x)^2+(y1+r*y2my1-y)^2;
                  if (dj<d^2) && (dj>tol)
                     d = sqrt(dj);
                  end

               end

            else
               break
            end
         end
      end

   else

      % Loop through edges top down
      for j = ne:-1:1
         y1 = edgexy_upper(j,2);
         if y1<=(y+d)
            y2 = edgexy_upper(j,4);
            if y2>=(y-d)

               x1 = edgexy_upper(j,1);
               x2 = edgexy_upper(j,3);

               if x1<x2
                  xmin = x1;
                  xmax = x2;
               else
                  xmin = x2;
                  xmax = x1;
               end

               if xmin<=(x+d) && xmax>=(x-d)
                  % Calculate the distance along the normal projection from [x,y] to the jth edge
                  x2mx1 = x2-x1;
                  y2my1 = y2-y1;

                  r = ((x-x1)*x2mx1+(y-y1)*y2my1)/(x2mx1^2+y2my1^2);
                  if r>1.0                                                 % Limit to wall endpoints
                     r = 1.0;
                  elseif r<0.0
                     r = 0.0;
                  end

                  dj = (x1+r*x2mx1-x)^2+(y1+r*y2my1-y)^2;
                  if (dj<d^2) && (dj>tol)
                     d = sqrt(dj);
                  end

               end

            else
               break
            end
         end
      end

   end

   L(k) = d;

end

end      % dist2poly()

author: johannesgerer

repo: jburkardt-m


text:
function [ readerobj ] = open_ph_reader( filename, varargin )
%OPEN_PH_READER Open a generic reader for multiple formats of phase history data
%
% Returns a reader object structure with the following methods:
%
% [wbvectors, per_pulse_metadata] = readerobj.read_cphd(pulse_indices, sample_indices, channels)
%      If file format holds raw phase history, the pulse data in the file
%      is read and then the correct processing is done so that motion
%      compensated pulses and metadata are returned.  If the file format
%      holds motion compensated phase history, this data is returned
%      directly.  Either way, this function always returns data at a
%      consistent stage in the processing chain.
%      PULSE_INDICES: Vector with indices to all desired pulses.  Can
%         also use the string 'all' to specify all pulses.  Default is all
%         pulses.
%      SAMPLES_INDICES: Vector with indices to all desired pulses.  Can
%         also use the string 'all' to specify all samples.  Default is
%         all samples.
%      CHANNELS: Vector with indices to all desired channels.  Can
%         also use the string 'all' to specify all pulses.  Default is 1.
%      WBVECTORS: The returned phase history (in 'wbvectors') is a matrix
%         of requested pulses wherein each pulse is a column with the first
%         requested pulse in column 1, second requested pulse in column 2,
%         etc.
%                      pulses ->
%         samples |  | 1 2             |
%                 v  | 1 2 ...         |
%                    | 1 2             |
%      PER_PULSE_METADATA: MATLAB structure contain metadata for each
%         individual pulse.
%
%      Examples:
%         Read all data from all pulses, samples, and channels
%         [wbvectors, nbdata] = readerobj.read_cphd(filename,'all','all','all');
%
%         Read every other pulse and samples 1-1000 from channel number 4
%         [wbvectors, nbdata] = readerobj.read_cphd(filename,...
%                                                   first_pulse:2:last_pulse,
%                                                   1:1000,
%                                                   4);
%
%         Read only the narrowband data.  Will read preamble and all per-pulse
%         metadata, but WBVECTORS will be empty since no samples were requested.
%         [wbvectors, nbdata] = readerobj.read_cphd(filename,'all',[]);
%
% [wbvectors, per_pulse_metadata] = readerobj.read_raw(pulse_indices, sample_indices, channels)
%      Reads raw pulse data without any phase stabilization/motion
%      compensation.  Not available for file formats that area already
%      motion compensated.  Input and output arguments are the same as for
%      the READ_CPHD method.
% readerobj.get_meta() % Returns a MATLAB structure with the CPHD
%         "preamble" metadata
% readerobj.close() % Should always be run after all reading is finished
%
% Written by: Wade Schwartzkopf, NGA/IDT
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

%% Determine format and open
format_string=guess_ph_format(filename);
if isempty(format_string)
    error('OPEN_PH_READER:UnrecognizedFormat','Unable to determine format of phase history file.');
end
readerobj = eval(['open_' format_string '_reader(''' filename ''', varargin{:});']);

end

% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

author: ngageoint


text:
function vTrunc= util_trunc(v, digits, policy)
%vTrunc= util_trunc(v, <digits=4, policy>)
%
% Truncates the number v, digits gives the amount of digits after the
% comma.
%
% policy is 'floor', 'ceil', or 'round'

if ~exist('digits', 'var'), digits=4; end
if ~exist('policy','var'), policy='round'; end

a= 10^digits;
vTrunc= feval(policy, a*v)/a;


text:
classdef LHSintegrator_StokesD < handle

    properties (Access = private)
        pressureFun
        velocityFun
        mesh
        quadrature
    end

    methods (Access = public)

        function obj = LHSintegrator_StokesD(cParams)
            %             obj.init(cParams);
            %             obj.createQuadrature();
            %             obj.createInterpolation();
            %             obj.createGeometry();
            obj.initStokesD(cParams);
            obj.createQuadrature();
        end

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

    end

    methods (Access = protected)

        function lhs = computeElementalLHS(obj)
            dNdxV = obj.velocityFun.computeCartesianDerivatives(obj.quadrature);
            dvolV = obj.mesh.computeDvolume(obj.quadrature)';
            shpeP = obj.pressureFun.computeShapeFunctions(obj.quadrature);

            nElem = obj.mesh.nelem;
            nDimfV = size(dNdxV,1);
            nNodeV = size(dNdxV,2);
            nNodeP = size(shpeP,1);

            nGaus = size(dNdxV,4);

            D = zeros(nDimfV*nNodeV,nNodeP,nElem);
            for igaus = 1:nGaus
                for inode_var = 1:nNodeP
                    for inode_test = 1:nNodeV
                        for idime = 1:nDimfV
                            dof_test = inode_test*nDimfV - nDimfV + idime;
                            v = squeeze(dNdxV(idime,inode_test,:,igaus));
                            D(dof_test,inode_var,:)= squeeze(D(dof_test,inode_var,:)) - v(:).*shpeP(inode_var,igaus)...
                                .*dvolV(:,igaus);
                        end
                    end
                end
            end
            lhs = D;
        end

    end

    methods (Access = private)

        function initStokesD(obj, cParams)
            obj.mesh     = cParams.mesh;
            obj.pressureFun = cParams.pressureFun;
            obj.velocityFun = cParams.velocityFun;
%             obj.material = cParams.material;
        end

        function createQuadrature(obj)
            q = Quadrature.set(obj.mesh.type);
            q.computeQuadrature('QUADRATIC'); % ehhh
            obj.quadrature = q;
        end

        function LHS = assembleMatrix(obj, lhs)
            s.fun    = []; % !!!
            assembler = AssemblerFun(s);
            LHS = assembler.assembleFunctions(lhs, obj.velocityFun, obj.pressureFun);
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [wRed,xRed,PRed,exitCode]=kMeansGaussReduction(wOrig,xOrig,POrig,param4,xInit,PInit,maxIter)
%%KMEANSGAUSSREDUCTION Perform Gaussian mixture reduction by iterating a k-
%              means algorithm given an inital reduced estimator (or using
%              Runnals' algorithm for the initial estimate if none is
%              given.
%
%INPUTS: wOrig An NX1 or 1XN vector of weights of the components of the
%          original Gaussian mixture.
%   muOrig An xDimXN matrix of the means of the vector components of the
%          original Gaussian mixture.
%    POrig An xDimXxDim XN hypermatrix of the covariance matrices for the
%          components of the original Gaussian mixture.
%   param4 If not initial estimate of a reduced distribution is provided,
%          then this is k, the number of components to which the
%          distribution should be reduced (and xInit and PInit are either
%          not provided or are empty matrices) and RunnalsGaussMixRed will
%          be used to obtain an initial estimate. Otherwise, if an initial
%          reduced distribution is provided, this is wInit, the kX1 or 1Xk
%          set of weights for the initial estimate of the reduced
%          distribution.
%    xInit If an initial reduced distribution is provided, this is the
%          xDimXk set of Gaussian mean vectors. Otherwise, this parameter
%          can be omitted or an empty matrix passed.
%    PInit If an initial reduced distribution is provided, this is the
%          xDimXxDimXk set of covariance matrices. Otherwise, this
%          parameter can be omitted or an empty matrix passed.
%  maxIter The maximum number of iterations of the k means algorithm to
%          perform. The default if omitted or an empty matrix is passed is
%          50.
%
%OUTPUTS: wRed The KX1 weights of the mixture after reduction.
%        muRed The xDimXK means of the mixture after reduction.
%         PRed The xDimXxDimXK covariance matrices of the mixture after
%              reduction.
%     exitCode A parameter indicating how the algorithm terminated.
%              Possible values are:
%              0 The algorithm converged.
%              1 the maximum number of iterations elapsed.
%
%This implements the algorithm of [1] and [2], which is also described in
%[3]. TGlobal convergence is not guaranteed and the performance depends a
%lot on the quality of the initial estiamte. The Kullback-Leiblber
%divergence is used as a closeness criterion between components via the
%KLDivGauss function. Additional algorithmic optimizations precomputing
%various determinants are possible beyond what was done here.
%
%EXAMPLE:
%We reduce a 10-component mixture to five components. We first plot
%Runnals' algorithm. We then use a perturbed version of the Runnals
%solution as the initial estimate to the k-means estimator and we find that
%it converges back to the Runnals solution.
% w=[0.03,0.18,0.12,0.19,0.02,0.16,0.06,0.1,0.08,0.06];
% n=length(w);
% mu=[1.45,2.20,0.67,0.48,1.49,0.91,1.01,1.42,2.77,0.89];
% P=[0.0487,0.0305,0.1171,0.0174,0.0295,0.0102, 0.0323, 0.0380, 0.0115, 0.0679];
% P=reshape(P,[1,1,n]);
%
% k=6;%Number of reduced components.
%
% %Runnal's reduction
% [wRedR,muRedR,PRedR]=RunnalsGaussMixRed(w,mu,P,k);
%
% %Perturb Runnal's solution so that the kMeansReduction has a bad
% %initialization.
% wRedRP=wRedR;
% wRedRP(1)=wRedR(1)+0.1;
% wRedRP=wRedRP/sum(wRedRP);
% muRedRP=muRedR+0.1;
% [wRedkM,muRedkM,PRedkM]=kMeansGaussReduction(w,mu,P,wRedRP,muRedRP,PRedR);
%
% numPoints=500;
% xVals=linspace(0,3,numPoints);
% PDFVals1=GaussianMixtureD.PDF(xVals,w,mu,P);
% PDFVals2=GaussianMixtureD.PDF(xVals,wRedR,muRedR,PRedR);
% PDFVals3=GaussianMixtureD.PDF(xVals,wRedRP,muRedRP,PRedR);
% PDFVals4=GaussianMixtureD.PDF(xVals,wRedkM,muRedkM,PRedkM);
% figure(1)
% clf
% hold on
% plot(xVals,PDFVals1,'-k','linewidth',4)
% plot(xVals,PDFVals2,'-r','linewidth',4)
% plot(xVals,PDFVals3,'-.c','linewidth',2)
% plot(xVals,PDFVals4,'--g','linewidth',2)
% legend('Full Mixture','Algorithm of Runnals','Initialization of k-Means','k-Means Algorithm')
% set(gca,'FontSize',14,'FontWeight','bold','FontName','Times')
%
%REFERENCES:
%[1] D. Schieferdecker and M. F. Huber, "Gaussian mixture reduction via
%    clustering," in Proceedings of the 12th International Conference on
%    Information Fusion, Seattle, WA, 6-9 Jul. 2009, pp. 1536-1543.
%[2] A. Nikseresht and M. Gelgon, "Gossip-based computation of a Gaussian
%    mixture model for distributed multimedia indexing," IEEE Transactions
%    on Multimedia, vol. 10, no. 3, pp. 385-392, Apr. 2008.
%[3] D. F. Crouse, P. Willett, K. Pattipati, and L. Svensson, "A look at
%    Gaussian mixture reduction algorithms," in Proceedings of the 14th
%    International Conference on Information Fusion, Chicago, IL, 5-8 Jul.
%    2011.
%
%May 2019 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<7||isempty(maxIter))
   maxIter=20;
end

if(nargin<5||isempty(xInit))
%If no initial estimate is given but the desired number of components is
%given, then get an initial estimate using Runnall's algorithm.
    [wInit,xInit,PInit]=RunnalsGaussMixRed(wOrig,xOrig,POrig,param4);
else
    wInit=param4;
end

%The original number of mixture parameters.
N=length(wOrig);
%The reduced number of mixture parameters.
k=length(wInit);

cAssign=zeros(N,1);
cAssignNew=zeros(N,1);

%We assign the Gaussians in the full mixture to cluster centers in the
%current reduced mixture. Then, we recompute the cluster centers. These two
%steps are looped until the assignments no longer change or a maximum
%number of iterations has passed.

wRed=wInit;
xRed=xInit;
PRed=PInit;

exitCode=1;
for curIter=1:maxIter
    %For each point, find the closest center.
    for curP=1:N
        minCost=Inf;
        for curC=1:k
            cost=KLDistGauss(xOrig(:,curP),POrig(:,:,curP),xRed(:,curC),PRed(:,:,curC));
            if(cost<minCost)
                minCost=cost;
                cAssignNew(curP)=curC;
            end
        end
    end

    %Terminate when the assignments no longer change.
    if(sum(cAssignNew~=cAssign)==0)
        exitCode=0;%Convergence attained.
        break;
    end
    cAssign=cAssignNew;

    %Now, merge all points having common centers.
    [cSorted,idx]=sort(cAssign,'ascend');
    [~,numReps]=runLenEncode(cSorted);
    startIdx=1;
    for curC=1:k
        selIdx=idx(startIdx:(startIdx+numReps(curC)-1));

        wSelNorm=wOrig(selIdx);
        wSum=sum(wSelNorm);
        wSelNorm=wSelNorm/wSum;

        [xMerged,PMerged]=calcMixtureMoments(xOrig(:,selIdx),wSelNorm,POrig(:,:,selIdx));

        wRed(curC)=wSum;
        xRed(:,curC)=xMerged;
        PRed(:,:,curC)=PMerged;

        startIdx=startIdx+numReps(curC);
    end

    wRed=wRed/sum(wRed);
end
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [f] = spm_fp_cmc_tfm(x,u,P,M)
% parameter equations for a neural mass model (canonical microcircuit)
% FORMAT [f] = spm_fp_cmc_tfm(x,u,P,M)
%
% x      - state vector
%   x(:,1) - voltage     (spiny stellate cells)
%   x(:,2) - conductance (spiny stellate cells)
%   x(:,3) - voltage     (superficial pyramidal cells)
%   x(:,4) - conductance (superficial pyramidal cells)
%   x(:,5) - voltage     (inhibitory interneurons)
%   x(:,6) - conductance (inhibitory interneurons)
%   x(:,7) - voltage     (deep pyramidal cells)
%   x(:,8) - conductance (deep pyramidal cells)
%
% f        - dP = h(x(t),u(t),P,M)
%
% Prior fixed parameter scaling
%
% G  = intrinsic rates
% D  = propagation delays (intrinsic, extrinsic)
% T  = synaptic time constants
% R  = slope of sigmoid activation function
%
%__________________________________________________________________________
% David O, Friston KJ (2003) A neural mass model for MEG/EEG: coupling and
% neuronal dynamics. NeuroImage 20: 1743-1755
%__________________________________________________________________________
% Copyright (C) 2005 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_fp_cmc_tfm.m 6234 2014-10-12 09:59:10Z karl $

% Neuronal states (deviations from baseline firing)
%--------------------------------------------------------------------------
%   x(:,1) - voltage     (spiny stellate cells)
%   x(:,2) - conductance (spiny stellate cells)
%   x(:,3) - voltage     (superficial pyramidal cells)
%   x(:,4) - conductance (superficial pyramidal cells)
%   x(:,5) - voltage     (inhibitory interneurons)
%   x(:,6) - conductance (inhibitory interneurons)
%   x(:,7) - voltage     (deep pyramidal cells)
%   x(:,8) - conductance (deep pyramidal cells)
%--------------------------------------------------------------------------
persistent iG nP Ca G
if isempty(iG)
    Ca  = zeros(size(P.G));
    G   = zeros(size(P.G));                    % parameter (deviates)
    iG  = spm_fieldindices(P,'G');
    nP  = spm_length(P);
end

% get dimensions and configure state variables
%--------------------------------------------------------------------------
f  = zeros(nP,1);                              % flow
x  = spm_unvec(x,M.x);                         % neuronal states
x  = x(:,1:2:end);                             % depolarisation

% neuronal populations with Voltage-dependent connectivity
%==========================================================================
%                  ss sp ii dp                 % neuronal populations
%--------------------------------------------------------------------------
a     = [1 8 2 1]*64;                          % potentiation rate
b     = [4 2 2 1]*4;                           % decay rate

NMDA  = @(x)1./(1 + exp(-x)) - 1/2;            % depolarisation CDF

% NMDA-like Voltage-dependent changes in (recurrent) synaptic efficacy
%--------------------------------------------------------------------------
A     = exp(P.E)*diag(a);
B     = exp(P.F)*diag(b);
dC    = (A.*NMDA(8*x) - Ca).*B;
dG    = Ca.*(2 - G)/2  - G.*B;
Ca    = Ca + dC*M.dt;
G     = G  + dG*M.dt;
f(iG) = G(:);

author: spm


text:
function disp( x, prefix ) %#ok
if nargin < 2, prefix = ''; end
disp( [ prefix, 'cvx base object' ] );

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

repo: REST


text:
function rxcallback(src, msg)
    disp([char(msg.Data()),sprintf('\n  Message received: %s',datestr(now))]);

author: petercorke

repo: robotics-toolbox-matlab


text:
%% RATE OF CONVERGENCE OF BILINEAR FINITE ELEMENT METHOD
%
% This example is to show the rate of convergence of bilinear finite
% element approximation of the Poisson equation on the unit square:
%
% $$- \Delta u = f \; \hbox{in } (0,1)^2$$
%
% for the following boundary condition:
%
% # Non-empty Dirichlet boundary condition. $u=g_D \hbox{ on }\Gamma_D, \quad \nabla u\cdot n=g_N \hbox{ on }\Gamma_N. \Gamma _D = \{(x,y): x=0, y\in [0,1]\}, \; \Gamma _N = \partial \Omega \backslash \Gamma _D$.
% # Pure Neumann boundary condition. $\Gamma _N = \partial \Omega$.
% # Robin boundary condition. $g_R u + \nabla u\cdot n=g_N \hbox{ on }\partial \Omega$

%%
clear all; close all; clc;
[node,elem] = squarequadmesh([0,1,0,1],1/2^5);
option.L0 = 1;
option.maxIt = 4;
option.printlevel = 1;
option.plotflag = 0;
option.elemType = 'Q1';

%% Non-empty Dirichlet boundary condition.
pde = sincosdata;
bdFlag = setboundary(node,elem,'Dirichlet','~(x==0)','Neumann','x==0');
% bdFlag = setboundary(node,elem,'Dirichlet');
femPoisson(node,elem,pde,bdFlag,option);
return

%% Pure Neumann boundary condition.
pde = sincosNeumanndata;
bdFlag = setboundary(node,elem,'Neumann');
femPoisson(node,elem,pde,bdFlag,option);

%% Pure Robin boundary condition.
option.plotflag = 0;
pdeRobin = sincosRobindata;
bdFlag = setboundary(node,elem,'Robin');
femPoisson(node,elem,pdeRobin,bdFlag,option);

%% Conclusion
%
% The optimal rate of convergence of the H1-norm (2nd order) and L2-norm
% (3rd order) is observed. The order of ||DuI-Duh|| is almost 3rd order and
% thus superconvergence exists.
%
% MGCG converges uniformly in all cases.

repo: Aero-matlab


text:
function F_L = CIECAM02_F_L(L_A)
%
%
%       F_L = CIECAM02_F_L(L_A)
%
%       Input:
%           -L_A: is the luminance of the adapting field in cd/m^2.
%
%       Output:
%           -F_L: is a predictor a variety of luminance-dependent
%            appearance effect.
%
%     Copyright (C) 2015  Francesco Banterle
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
%     The paper describing this technique is:
%     "The CIECAM02 color appearance model"
% 	  by Nathan Moroney , Mark D. Fairchild , Robert W. G. Hunt ,
%     Changjun Li , M. Ronnier Luo , Todd Newman
%     in IS&T/SID 10 th Color Imaging Conference
%

k = 1.0 ./ (5 * L_A + 1.0); %Equation 1

F_L = 0.2 * k.^4 .* (5 * L_A) + ...
      0.1 * (1.0 - k.^4).^2 .* ((5 * L_A).^(1.0 /  3.0)); %Equation 2

end

author: banterle


text:
function out = InsertDate(in)

%InsertDate - Insert date into string.
%
% Replaces %y, %m, %d and %t with current year, month, day and time.
%
%  USAGE
%
%    out = InsertDate(in)
%
%    in             string to modify

% Copyright (C) 2013 by Michaël Zugaro
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 3 of the License, or
% (at your option) any later version.

n = clock;
year = int2str(n(1));
month = int2zstr(n(2),2);
day = int2zstr(n(3),2);
time = [int2zstr(n(4)) '' int2zstr(n(5))];

out = in;
out = strrep(out,'%y',year);
out = strrep(out,'%m',month);
out = strrep(out,'%d',day);
out = strrep(out,'%t',time);

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function StepThru(mfile,ws)
%STEPTHRU automatically steps through specified m-file
%
%   STEPTHRU(MFILE) sets breakpoints at all but first line comment and
%   executes MFILE in order to step through it using debugger.
%
%       Line comments begin with % character.
%
%       No breakpoint is generated if the text 'nobkpt' is
%       anywhere in the line comment.
%
%   STEPTHRU(MFILE,WS) executes MFILE in specified
%   workspace (ie, 'caller' or 'base' valid).

%This utility intended to simplify MATLAB script demos
%v2 Apr2004 rbemis@mathworks.com (added npbkpt option)
%v1 Oct2002 rbemis@mathworks.com

% Copyright 2004-2010 The MathWorks, Inc.

% find pure comment lines to set breakpoints
fid=fopen([mfile '.m']);
line_num=0;
while ~feof(fid) %read 1 line at-a-time
  dummy=NoLeadingBlanks(fgetl(fid));
  line_num=line_num+1;
  if length(dummy)>0 %nonempty line
    if dummy(1)=='%' & line_num>1 %comment line but not line1
      if strfind(lower(dummy),'nobkpt') %no breakpoint this line
        eval(['dbclear in ' mfile ' at ' num2str(line_num)])
      else
        eval(['dbstop in ' mfile ' at ' num2str(line_num)])
      end
    end
  end
end %file empty
fclose(fid);

%open file in Editor
edit(mfile)

% execute m-file (in specified workspace)
if nargin>1 %workspace specified
  evalin(ws,mfile)
else %default=local workspace (ie, this function)
  eval(mfile)
end
dbclear all

% helper function to remove any leading blanks
%--------------------------------------------------
function str1 = NoLeadingBlanks(str1)
Done=false;
while ~Done %keep looking
  if length(str1)==0 %empty string
    Done=true; %get out
  else %check first character
    if str1(1)==' ' %leading blank detected
      str1=str1(2:end); %remove it
    else %any other character
      Done=true; %leave it, done
    end
  end
end %stripped!


text:
function expoTimes = ReadExpoTimes(scenePath)

[~, expPaths, ~] = GetFolderContent(scenePath, '.txt');

fid = fopen(expPaths{1});
expoTimes = 2.^cell2mat(textscan(fid, '%f'));
fclose(fid);

author: qingsenyangit

repo: AHDRNet


text:
function value = lietzke_condition ( n )

%*****************************************************************************80
%
%% LIETZKE_CONDITION returns the L1 condition of the LIETZKE matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 April 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Output, real VALUE, the L1 condition.
%
  s = 0;
  k = n;
  for i = 1 : n
    s = s + k;
    if ( mod ( i, 2 ) == 1 )
      k = k - 1;
    end
  end
  a_norm = s;
  if ( n == 1 )
    b_norm = 0.25;
  elseif ( n == 2 )
    b_norm = 5.0 / 6.0;
  else
    b_norm = 2.0;
  end
  value = a_norm * b_norm;

  return
end


text:
function pass = test_equiOption( pref )
% Test funqi in 2D

if ( nargin < 1 )
    pref = chebfunpref;
end
tol = 100*pref.cheb2Prefs.chebfun2eps;

% Canonical domain:
dom = [-1 1 -1 1];
f = @(x,y) cos(x+y);
x = linspace(dom(1),dom(2),100);
y = linspace(dom(3),dom(4),100);
[xx, yy] = meshgrid( x, y );
A = f(xx, yy) ;
g = chebfun2( A , dom, 'equi' );
h = chebfun2( f, dom);
pass(1) = norm( h - g ) < tol ;

% Rectangular domain:
dom = [-1 2 -2 1];
f = @(x,y) cos(x+y);
x = linspace(dom(1),dom(2),100);
y = linspace(dom(3),dom(4),100);
[xx, yy] = meshgrid( x, y );
A = f(xx, yy) ;
g = chebfun2( A , dom, 'equi' );
h = chebfun2( f, dom);
pass(2) = norm( h - g ) < tol ;

% Nonsymmetric function:
dom = [-1 2 -2 1];
f = @(x,y) cos(x+2*y);
x = linspace(dom(1),dom(2),100);
y = linspace(dom(3),dom(4),100);
[xx, yy] = meshgrid( x, y );
A = f(xx, yy) ;
g = chebfun2( A , dom, 'equi' );
h = chebfun2( f, dom);
pass(3) = norm( h - g ) < tol ;

% Small domain;
h = 1e-3;
dom = [1-h 1+h 1-2*h 1+2*h];
f = @(x,y) cos(x+2*y);
x = linspace(dom(1),dom(2),100);
y = linspace(dom(3),dom(4),100);
[xx, yy] = meshgrid( x, y );
A = f(xx, yy) ;
g = chebfun2( A , dom, 'equi' );
h = chebfun2( f, dom);
pass(4) = norm( h - g ) < tol ;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function alpha = complexNormalAngle(varargin)
%COMPLEXNORMALANGLE compute normal angle of a vertex of a cellular complex
%
%   ALPHA = complexNormalAngle(NODES, EDGES, FACES, INDEX)
%   ALPHA = complexNormalAngle(NODES, EDGES, FACES, CELLS, INDEX)
%   Compute the nortmal angle of the polyhedral reconstruction defined be
%   nodes NODES, edges EDGES and faces FACES. For 3D reconstructions, it
%   can also contain cells CELLS. INDEX is the index of NODES for which the
%   normal angle ALPHA is computed.
%   Result is normalised between 0 and 2*PI.
%
%   ALPHA = complexNormalAngle(GRAPH, INDEX)
%   Internal data are stored in a structure GRAPH, with fields : 'nodes',
%   'edges', 'faces', and eventually 'cells'.
%
%
%   ALPHA = complexNormalAngle(..., INDICES)
%   If INDICES is an array of indices, the normal angle is computed for
%   each element of NODES(INDICES,:). The result ALPHA has the same size
%   than INDICES.
%
%
% ------
% Author: David Legland
% e-mail: david.legland@jouy.inra.fr
% Created: 2005-12-19
% Copyright 2005 INRA - CEPIA Nantes - MIAJ (Jouy-en-Josas).

%   2006-04-19 fix bug for small number of faces
%   2006-04-26 returns 0 and not [] for null complexes
%   2006-10-25 revert to return value=[] for null complex
%   2008-08-11 code clean up

cells = [];
if length(varargin)==4
    % no cells in cellular complex
    nodes = varargin{1};
    edges = varargin{2};
    faces = varargin{3};
    ind   = varargin{4};

elseif length(varargin)==5
    % cells are given
    nodes = varargin{1};
    edges = varargin{2};
    faces = varargin{3};
    cells = varargin{4};
    ind   = varargin{5};

elseif length(varargin)==2
    % data stored as structure
    graph = varargin{1};
    nodes = graph.nodes;
    edges = graph.edges;
    faces = graph.faces;
    if isfield(graph, 'cells')
        cells = graph.cells;
    end
    ind   = varargin{2};
else
    error('wrong number of arguments');
end

alpha0 = zeros([length(ind) 1]);
alpha1 = zeros([length(ind) 1]);
alpha2 = zeros([length(ind) 1]);
alpha3 = zeros([length(ind) 1]);
alpha = [];

if size(nodes, 2)==2
    % 2 dimensions

    if iscell(faces)

        % process faces as cell array
        for i=1:length(ind)
            % check that vertex is contained in the complex
            if ind(i)>size(nodes, 1)
                continue;
            end

            % normal angle of vertex
            alpha0(i) = 2*pi;

            % normal angle of edges
            alpha1(i) = length(find(sum(edges==ind(i), 2)))*pi;

            % normal angle of faces
            alpha2(i) = 0;
            for j=1:length(faces)
                face = faces{j};
                indf = find(face==ind(i));

                if ~isempty(indf)
                    alpha2(i) = alpha2(i) + polygonNormalAngle(nodes(face,:), indf);
                end
            end
        end

    else
        % process faces as arrays
        for i=1:length(ind)

            % check that vertex is contained in the complex
            if ind(i)>size(nodes, 1)
                continue;
            end

            % normal angle of vertex
            alpha0(i) = 2*pi;

            % normal angle of edges
            alpha1(i) = length(find(sum(edges==ind(i), 2)))*pi;

            % normal angle of faces
            alpha2(i) = 0;
            for j=1:size(faces, 1)
                face = faces(j,:);
                indf = find(face==ind(i));

                if ~isempty(indf)
                    alpha2(i) = alpha2(i) + polygonNormalAngle(nodes(face,:), indf);
                end
            end

        end
    end

    % compute total normal angle of reconstruction
    alpha = alpha0 - alpha1 + alpha2;

elseif size(nodes, 2)==3
    % 3 dimensions
    for i=1:length(ind)

        % check that vertex is contained in the complex
        if ind(i)>size(nodes, 1)
            continue;
        end

        % normal angle of vertex
        alpha0(i) = 4*pi;

        % normal angle of edges
        alpha1(i) = length(find(sum(edges==ind(i), 2)))*2*pi;

        % normal angle of faces
        alpha2(i) = 0;
        if iscell(faces)
            % process faces as cell array
            for j=1:length(faces)
                face = faces{j};
                indf = find(face==ind(i));

                if ~isempty(indf)
                    alpha2(i) = alpha2(i) + polygon3dNormalAngle(nodes(face,:), indf);
                end
            end

        else
            % process faces as array of double
            for j=1:size(faces, 1)
                face = faces(j,:);
                indf = find(face==ind(i));

                if ~isempty(indf)
                    alpha2(i) = alpha2(i) + polygon3dNormalAngle(nodes(face,:), indf);
                end
            end
        end

        % normal angle of cells
        alpha3(i) = 0;
        for j=1:length(cells)
            cell = cells{j};
            if iscell(faces)
                cellFaces = faces(cell);
            else
                cellFaces = faces(cell, :);
            end

            alpha3(i) = alpha3(i) + polyhedronNormalAngle(nodes, cellFaces, ind(i));
        end

        % compute total normal angle of reconstruction
        alpha = alpha0 - alpha1 + alpha2 - alpha3;
   end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% set of internal functions, for angle computations in 2D and 3D

function theta = polygonNormalAngle(points, ind)
%POLYGONNORMALANGLE compute normal angle at a vertex of the polygon

% number of points
np = size(points, 1);

% number of angles to compute
nv = length(ind);

theta = zeros(nv, 1);

for i=1:nv
    p0 = points(ind(i), :);

    % previous vertex
    if ind(i)==1
        p1 = points(np, :);
    else
        p1 = points(ind(i)-1, :);
    end

    % next vertex
    if ind(i)==np
        p2 = points(1, :);
    else
        p2 = points(ind(i)+1, :);
    end

    % compute angles
    theta1  = mod(atan2(p1(2)-p0(2), p1(1)-p0(1)) + 2*pi, 2*pi);
    theta2  = mod(atan2(p2(2)-p0(2), p2(1)-p0(1)) + 2*pi, 2*pi);
    dtheta  = mod(theta2-theta1+2*pi, 2*pi);

    % use simplification due to the fact that cells are convex
    dtheta  = min(dtheta, 2*pi-dtheta);
    theta(i)= pi - dtheta;
end

return;

function theta = polyhedronNormalAngle(nodes, faces, ind)
%POLYHEDRONNORMALANGLE compute normal angle at a vertex of a 3D polyhedron
%
%   THETA = polyhedronNormalAngle(NODES, FACES, IND);
%   where NODES is a set of 3D points, and FACES a set of faces, whose
%   elements are indices to NODES array, compute the normal angle at the
%   vertex whose index is given by IND.

% number of angles to compute
na = length(ind);

theta = zeros(na, 1);
for i=1:na

    % find faces containing given vertex,
    % and compute normal angle at each face containing vertex
    if iscell(faces)
        for j=1:length(faces)
            if ismember(ind(i), faces{j})
                % create 3D polygon
                face = nodes(faces{j}, :);

                % index of point in polygon
                indp = find(faces{j}==ind(i));

                % compute face angle
                thetaf = [thetaf polygon3dInnerAngle(face, indp)];
            end
        end
    else
        indf = find(sum(ismember(faces, ind(i)), 2));

        thetaf = zeros(length(indf), 1);
        for j=1:length(indf)
            ind2 = faces(indf(j), :);
            face = nodes(ind2, :);
            indp = find(ind2==ind(i));
            thetaf(j) = polygon3dInnerAngle(face, indp);
        end
    end

    % compute normal angle of polyhedron, by use of angle defect formula
    if ~isempty(thetaf)
        theta(i) = 2*pi - sum(thetaf);
    end
end

return;

function theta = polygon3dNormalAngle(points, ind)
%POLYGON3DNORMALANGLE compute normal angle at a vertex of the 3D polygon

theta = 2*pi - 2*polygon3dInnerAngle(points, ind);

return;

function theta = polygon3dInnerAngle(points, ind)
%POLYGON3DNORMALANGLE compute normal angle at a vertex of the 3D polygon

% number of points
np = size(points, 1);

% number of angles to compute
nv = length(ind);

theta = zeros(nv, 1);

for i=1:nv
    p0 = points(ind(i), :);

    % previous vertex
    if ind(i)==1
        p1 = points(np, :);
    else
        p1 = points(ind(i)-1, :);
    end

    % next vertex
    if ind(i)==np
        p2 = points(1, :);
    else
        p2 = points(ind(i)+1, :);
    end

    theta(i) = angle3d(p1, p0, p2);
    theta(i) = min(theta(i), 2*pi-theta(i));
    % todo: solve case for CW oriented polygons
end

return;

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function [x,info] = mgS(A,b,elem,option,varargin)
%% MGS multigrid-type solvers
%
%   x = MGS(A,b,elem) attempts to solve the system of linear equations A*x =
%   b for x using multigrid type solver. To acheive multigrid efficiency,
%   the mesh must be of bisection type. Inside mg, an coarsening algorithm
%   on bisection grids is applied. See doc coarsen for the coarsening algorithm.
%
%   It is a simplified version of |mg| such that it can serve as a template
%   for coding multigrid methods using iFEM.
%
% Example:
%
%   [node,elem] = squaremesh([0 1 0 1],0.5);
%   for k = 1:8
%     [node,elem] = uniformrefine(node,elem);
%   end
%   pde.f = inline('p(:,1).*p(:,2)','p');
%   pde.g_D = inline('zeros(size(p,1),1)','p');
%   option.solver = 'notsolve';
%   [u,Du,eqn] = Poisson(node,elem,pde,[],option);
%   fprintf('\n Number of unknowns: %8.0u\n',length(eqn.b))
%   tic; display('Direct solver'); u = eqn.A\eqn.b; toc;
%   tic; x = mgS(eqn.A,eqn.b,elem); toc;
%   format shorte
%   fprintf('Difference between direct and mg solvers %0.2g \n',norm(u-x));
%
% See also mg
%
% Documentation in Help browser ifem mgdoc
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

tic;

%% Size of systems
Ndof = length(b);                  % number of dof
N = max(elem(:));                  % number of nodes
% NT = size(elem,1);                 % number of elements
dim = size(elem,2)-1;

%% Options
% Assign default values to unspecified parameters
if ~exist('option','var')
    option = [];
end
option = mgoptions(option,Ndof);    % parameters
x0 = option.x0;
N0 = option.N0;
tol = option.tol;
maxIt = option.solvermaxit;
mu = option.smoothingstep;
solver = option.solver;
preconditioner = option.preconditioner;
coarsegridsolver = option.coarsegridsolver;
printlevel = option.printlevel;
% setupflag = option.setupflag;

%% Fixed dof
isFixDof = [];
if isfield(option,'freeDof') % freeDof is given
   isFreeDof = false(N,1);
   isFreeDof(option.freeDof) = true;
   NA = size(A,1);
   if NA > length(option.freeDof)
       isFixDof = true(NA,1);
       isFixDof(isFreeDof) = false;
   end
else % Find free dof and eliminate isolated dof
    deg = sum(spones(A));  % degree of the matrix graph
    isFreeDof = false(Ndof,1);
    isFreeDof(deg>1) = true;
    isFixDof = ~isFreeDof; % degree = 1 is isolated vertex
    isFixDof(deg == 0) = false; % degree = 0 the matrix is singular
end
if any(isFixDof) % a bigger matrix is given
    xD = zeros(Ndof,1);
    xD(isFixDof) = b(isFixDof)./diag(A(isFixDof,isFixDof));
    A = A(isFreeDof,isFreeDof);% shrink to freedof only
end
if length(b) > sum(isFreeDof)  % a bigger b is given
    b = b(isFreeDof);          % shrink to freedof only
end
if length(x0) > sum(isFreeDof) % a bigger x0 is given
    x0 = x0(isFreeDof);        % shrink to freedof only
    option.x0 = x0;
end
isFreeNode = isFreeDof(1:N); % free nodes (for P1)

%% Hierarchical Structure of Mesh
%  hierarchical structure for linear element.
if dim == 2  % 2-D
    [HB, NL, level] = HBstructure(elem,N0);
end
if dim == 3  % 3-D
    if nargin > 4
        HBmesh = varargin{end}; % need HB for bisection refinement
        [HB, NL, level] = HBstructure3(elem,N0,HBmesh);
    else % no HBmesh is given. only work for the red uniform refinement
        [HB, NL, level] = HBstructure3(elem,N0);
    end
end

%% Transfer operators between multilevel meshes for P1 element
% standard prolongation and restriction operator for P1 element
[Pro,Res] = transferoperator(HB,NL,isFreeNode);
if Ndof > N
    if ~exist('auxPro','var')
        disp('The current element is not supported by mg');
    else
        Pro{level} = auxPro;
        Res{level+1} = Pro{level}';
        level = level + 1; % add one more level from P1 to current element
    end
end
clear HB auxPro

%% Matrices in each level
Ai = cell(level,1);
Ai{level} = A;
for j = level:-1:2
    Ai{j-1} = Res{j}*Ai{j}*Pro{j-1};           % Ac = Res*Af*Pro
    switch option.smoother
        case 'GS'
            Bi{j} = tril(Ai{j});        % Forward Gauss-Seidel   B = D+L
            BBi{j} = triu(Ai{j});       % Backward Gauss-Seidel BB = D+U
        case 'JAC'
            Bi{j} = spdiags(diag(Ai{j}),0,size(Ai{j},1),size(Ai{j},1));        %
            BBi{j} = Bi{j};       % Jacobi iteration
    end
    if option.smoothingparameter~=1
        Bi{j} = Bi{j}/option.smoothingparameter;
        BBi{j} = BBi{j}/option.smoothingparameter;
    end
end

%% MG cycles
% initial set up
k = 1;
x = x0;
r = b - A*x;
nb = norm(b);
err = zeros(maxIt,2);
if nb > eps  % nb is non-zero
    err(1,:) = norm(r)/nb;
else
    err(1,:) = norm(r);
end
switch solver
    case 'VCYCLE'
        if printlevel >= 1
            fprintf('Multigrid Vcycle Iteration \n')
        end
        while (max(err(k,:)) > tol) && (k <= maxIt)
            k = k + 1;
            % Step 2: Compute Br by one Vcylce MG
            Br = vcycle(r);
            % Step 3: Correct the solution
            x = x + Br;
            % Step 1: Form residual r
            r = r - A*Br;
            err(k,1) = sqrt(abs(Br'*r/(x'*b))); % approximate relative error in energy norm
            err(k,2) = norm(r)/nb; % relative error of the residual in L2-norm
            if printlevel >= 2
                fprintf('#dof: %8.0u,  #nnz: %8.0u, MG Vcycle iter: %2.0u, err = %8.4e\n',...
                         Ndof, nnz(A), k-1, max(err(k,:)));
            end
        end
        err = err(1:k,:);
        itStep = k-1;
    case 'WCYCLE'
        if printlevel >= 1
            fprintf('Multigrid Wcycle Iteration \n')
        end
        while (max(err(k,:)) > tol) && (k <= maxIt)
            k = k + 1;
            % Step 2: Compute Br by one Vcylce MG
            Br = wcycle(r);
            % Step 3: Correct the solution
            x = x + Br;
            err(k,1) = sqrt(abs(Br'*r/(x'*b))); % approximate relative error in energy norm
            % Step 1: Form residual r
            r = r - A*Br;
            err(k,2) = norm(r)/nb; % relative error of the residual in L2-norm
            if printlevel >= 2
                fprintf('#dof: %8.0u,  #nnz: %8.0u, MG Wcycle iter: %2.0u, err = %8.4e\n',...
                         Ndof, nnz(A), k-1, max(err(k,:)));
            end
        end
        err = err(1:k,:);
        itStep = k-1;
end

%% Krylov space method
% set up preconditioner
switch preconditioner
    case 'V'
        prefunc = @vcycle;
        if printlevel >= 1
            fprintf('Multigrid V-cycle Preconditioner with ')
        end
    case 'W'
        prefunc = @wcycle;
        if printlevel >= 1
            fprintf('Multigrid W-cycle Preconditioner with ')
        end
end
switch solver
    case 'CG'
        if printlevel >= 1
            fprintf('Conjugate Gradient Method\n')
        end
        [x,flag,err,itStep] = pcg(A,b,tol,maxIt,prefunc,[],x0);
    case 'MINRES'
        if printlevel >= 1
            fprintf('Minimum Residual Method \n')
        end
        [x,flag,err,itStep] = minres(A,b,tol,maxIt,prefunc,[],x0);
    case 'GMRES'
        if printlevel >= 1
            fprintf('General Minimum Residual Method\n')
        end
        if isfield(option,'restart')
            restart = option.restart;
        else
            restart = min(N,10);
        end
        [x,flag,err,itStep] = gmres(A,b,restart,tol,maxIt,prefunc,[],x0);
        itStep = (itStep(1)-1)*restart + itStep(2);
end

%% Modify x to include fix dof
if any(isFixDof)
    xD(isFreeDof) = x;
    x = xD;
end

%% Output
if k > maxIt
    flag = 1;
else
    flag = 0;
end
time = toc;
if printlevel >= 2
    fprintf('#dof: %8.0u, level: %2.0u,   coarse grid %2.0u, #nnz: %8.0u\n',...
              Ndof, level, size(Ai{1},1), nnz(Ai{1}))
end
if printlevel >= 1
    fprintf('#dof: %8.0u,  #nnz: %8.0u, smoothing: %2.0u, iter: %2.0u,   err = %8.4e,   time = %4.2g s\n',...
                 Ndof, nnz(A), mu, itStep, max(err(end,:)), time)
end
if (flag == 1) && (printlevel>0)
   fprintf('NOTE: the iterative method does not converge! \n');
end
info = struct('solverTime',time,'itStep',itStep,'error',err,'flag',flag,'stopErr',max(err(end,:)));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% subfunctions vcycle, wcycle, fcycle, bpx
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Vcycle MG
    function Br = vcycle(r,J)        % solve equations Ae = r in each level
    if nargin<=1
        J = level;
    end
    ri = cell(J,1);            % record residual in each level
    ei = cell(J,1);            % record err in each level
    ri{J} = r;
    for i = J:-1:2
        ei{i} = Bi{i}\ri{i};   % pre-smoothing
        for s = 1:mu-1           % extra mu-1 steps smoothing
            ei{i} = ei{i} + Bi{i}\(ri{i}-Ai{i}*ei{i});
        end
        ri{i-1} = Res{i}*(ri{i} - Ai{i}*ei{i});
    end
    if strcmp(coarsegridsolver,'direct')
        ei{1} = Ai{1}\ri{1}; % direct solver in the coarest level
    else                         % iterative solver in the coarest level
        D = spdiags(diag(Ai{1}),0,size(Ai{1},1),size(Ai{1},1));
        [ei{1},flag] = pcg(Ai{1},ri{1},1/size(Ai{1},1),1000,D);
    end
    for i = 2:J
        ei{i} = ei{i} + Pro{i-1}*ei{i-1};
        ei{i} = ei{i} + BBi{i}\(ri{i}-Ai{i}*ei{i});
        for s = 1:mu-1
            ei{i} = ei{i} + BBi{i}\(ri{i}-Ai{i}*ei{i}); % post-smoothing
        end
    end
    Br = ei{J};
    end

%% Wcycle MG
    function e = wcycle(r,J)        % solve equations Ae = r in each level
    if nargin<=1
        J = level;
    end
    if J == 1
        e = Ai{J}\r;   % exact solver in the coaresest grid
        return
    end
    % fine grid pre-smoothing
    e = Bi{J}\r;        % pre-smoothing
    for s = 1:mu-1        % extra mu-1 steps smoothing
        e = e + Bi{J}\(r-Ai{J}*e);
    end
    % coarse grid correction twice
    rc = Res{J}*(r - Ai{J}*e);
    ec = wcycle(rc,J-1);
    ec = ec + wcycle(rc - Ai{J-1}*ec,J-1);
    e = e + Pro{J-1}*ec;
    % fine grid post-smoothing
    e = e + BBi{J}\(r-Ai{J}*e);
    for s = 1:mu-1
        e = e + BBi{J}\(r-Ai{J}*e); % post-smoothing
    end
    end

end

repo: ifem


text:
function [score] = FRIQUEE(I)

testFriqueeFeats = extractFRIQUEEFeatures(I);

% Load a learned model
load('support_functions/FRIQUEE/data/friqueeLearnedModel.mat');

% Scale the features of the test image accordingly.
% The minimum and the range are computed on features of all the images of
% LIVE Challenge Database
testFriqueeALL = testFeatNormalize(testFriqueeFeats.friqueeALL, friqueeLearnedModel.trainDataMinVals, friqueeLearnedModel.trainDataRange);

score = svmpredict (0, double(testFriqueeALL), friqueeLearnedModel.trainModel, '-b 1 -q');

end

author: HuiZeng


text:
function D = spm_diag_array(X)
% Extracts diagonal from 3-D arrays
% FORMAT D = spm_diag_array(X)
%
% X(:,i,i) -> D(:,i);
%__________________________________________________________________________
% Copyright (C) 2014 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_diag_array.m 5956 2014-04-16 14:34:25Z guillaume $

% extract diagnonals
%--------------------------------------------------------------------------
if iscell(X)
    D     = cell(size(X));
    for i = 1:numel(X)
        D{i} = spm_diag_array(X{i});
    end
    return
end
D     = zeros(size(X,1),size(X,2));
for i = 1:size(X,3)
    D(:,i) = X(:,i,i);
end


text:
function [x] = tt_meshgrid(varargin)
%Analogue of the meshgrid function for the TT-format
%   X = TT_MESHGRID(A,B,C,...) Computes the meshgrid based on "1d"
%       representations
%   X = TT_MESHGRID(A) Computes the meshgrid based on the cell array A of
%       the representations
%   X = TT_MESHGRID(T,D) Computes the d-dimensional meshgrid, using T as a
%       one-dimensional grid
if ( numel(varargin) == 2 )
    if ( ismatrix(varargin{2}) )
        d = varargin{2};
        t = varargin{1};
        z = cell(d,1);
        for i = 1:d
            z{i} = t;
        end
    else
        z = varargin{1};
    end
else
    z = varargin{1};
end
if ~iscell(z)
    z = varargin;
end
d = numel(z);
x = cell(d,1);
e = cell(d,1);
for i = 1:d
    e{i} = tt_ones(size(z{i}));
end
for i = 1:d
    v = z{i};
    for j = i-1:-1:1
        v = kron(e{j},v);
    end
    for j = i+1:d
        v = kron(v,e{j});
    end
    x{i} = v;
end
end

author: oseledets

repo: TT-Toolbox


text:
%SerialLink.cinertia Cartesian inertia matrix
%
% M = R.cinertia(Q) is the NxN Cartesian (operational space) inertia matrix which relates
% Cartesian force/torque to Cartesian acceleration at the joint configuration Q, and N
% is the number of robot joints.
%
% See also SerialLink.inertia, SerialLink.rne.

% MOD HISTORY
% 	4/99 add object support
% $Log: not supported by cvs2svn $
% $Revision: 1.2 $

% Copyright (C) 1993-2017, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB 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 Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

function Mx = cinertia(robot, q)
	J = jacob0(robot, q);
	Ji = inv(J);                %#ok<*MINV>
	M = inertia(robot, q);
	Mx = Ji' * M * Ji;


text:
function fv= proc_rSquare(fv, varargin)
%PROC_RSQUARE - computes r^2 values (measure for discriminance)
%
%Synopsis:
% 	fv = proc_rSquare(fv, <opt>)
%
%Returns:
% FV_RVAL - data structure of squared biserial correlation coefficients
%  .x     - squared biserial correlation between each feature and the class label
%  .se    - standard error of atanh(r), if opt.Stats==1
%  .p     - p value of null hypothesis that there is zero
%           correlation between feature and class-label, if opt.Stats==1
%  .sgnlogp - contains the signed log10 p-value, if opt.Stats==1
%             if opt.Bonferroni==1, the p-value is multiplied by
%             fv_rval.corrfac
%  .sgnlogp - contains the signed log10 p-value, if opt.Stats==1
%           if opt.Bonferroni==1, the p-value is multiplied by
%           fv_rval.corrfac, cropped, and then logarithmized
%  .sigmask - binary array indicating significance at alpha level
%             opt.Alphalevel, if opt.Stats==1 and opt.Alphalevel > 0
%  .corrfac - Bonferroni correction factor (number of simultaneous tests),
%             if opt.Bonferroni==1
%
%Properties:
% 'TolerateNans': observations with NaN value are skipped
%    (nanmean/nanstd are used instead of mean/std). Deafult: 0
% 'ValueForConst': constant feauture dimensions are assigned this
%    value. Default: NaN.
% 'MulticlassPolicy': possible options: 'pairwise' (default),
%    'all-against-last', 'each-against-rest', or provide specified
%    pairs as an [nPairs x 2] sized matrix. ('specified_pairs' is obsolete)
% 'Stats' - if true, additional statistics are calculated, including the
%           standard error of atanh(r), the p-value for the null
%           Hypothesis that the correlation is zero,
%           and the "signed log p-value"
% 'Bonferroni' - if true, Bonferroni corrected is used to adjust p-values
%                and their logarithms
% 'Alphalevel' - if provided, a binary indicator of the significance to the
%                alpha level is returned for each feature in fv_rval.sigmask
%
%Description
% Computes the r^2 value for each feature. The r^2 value is a measure
% of how much variance of the joint distribution can be explained by
% class membership.
%
% See also proc_classmeanDiff, proc_rValues, proc_rSquareSigned
%
% 03-03 Benjamin Blankertz
% 09-2012 stefan.haufe@tu-berlin.de

if nargin==0,
  fv=proc_rValues; return
end

fv= proc_rValues(fv, varargin{:});
fv.x= fv.x.^2;
for cc= 1:length(fv.className),
  fv.className{cc}= ['r^2' fv.className{cc}(2:end)];
end
fv.yUnit= 'r^2';

author: bbci


text:
function out = isempty( F )
%ISEMPTY empty boolean check for a CHEBFUN2V object.
%   ISEMPTY(F) returns 1 if every component of F is an empty CHEBFUN2, and
%   return 0 otherwise.

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

if ( isempty( F.components ) )
    out = 1;
    return
end

% Take isempty of each component:
out = cellfun( @isempty, F.components, 'UniformOutput', false );
out = all(cell2mat( out ) );

end


text:
% DEMOIL5 Oil data with partially independent training conditional.

% FGPLVM

% Fix seeds
randn('seed', 1e5);
rand('seed', 1e5);

dataSetName = 'oil';
experimentNo = 5;

% load data
[Y, lbls] = lvmLoadData(dataSetName);

% Set up model
options = fgplvmOptions('pitc');
options.optimiser = 'scg';
latentDim = 2;
d = size(Y, 2);

model = fgplvmCreate(latentDim, d, Y, options);

% Optimise the model.
iters = 1000;
display = 1;

model = fgplvmOptimise(model, display, iters);

% Save the results.
capName = dataSetName;;
capName(1) = upper(capName(1));
save(['dem' capName num2str(experimentNo) '.mat'], 'model');

if exist('printDiagram') & printDiagram
  fgplvmPrintPlot(model, lbls, capName, experimentNo);
end

% Load the results and display dynamically.
fgplvmResultsDynamic(dataSetName, experimentNo, 'vector')

% compute the nearest neighbours errors in latent space.
errors = fgplvmNearestNeighbour(model, lbls);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function test18
%TEST18 test cholmod2 on a few large matrices
% Example:
%   test18
% See also cholmod_test

% Copyright 2006-2007, Timothy A. Davis, University of Florida

fprintf ('=================================================================\n');
fprintf ('test18: test cholmod2 on a few large matrices\n') ;

rand ('state',1) ;
randn ('state',1) ;

Prob = UFget (887)							    %#ok
A = Prob.A ;
n = size (A,1) ;
b = rand (n,1) ;
x = cholmod2 (A,b) ;
norm (A*x-b,1)

clear all

Prob = UFget (936)							    %#ok
A = Prob.A ;
n = size (A,1) ;
b = rand (n,1) ;
x = cholmod2 (A,b) ;
norm (A*x-b,1)

clear all

Prob = UFget (887)							    %#ok
A = Prob.A ;
n = size (A,1) ;
b = rand (n,1) ;
x = cholmod2 (A,b) ;
norm (A*x-b,1)

fprintf ('test18 passed\n') ;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function ft = f03_title ( )

%*****************************************************************************80
%
%% F03_TITLE returns the title for function 3.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 January 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string FT, the function title.
%
  ft = 'Saddle';

  return
end


text:
function cg_rc_test02 ( )

%*****************************************************************************80
%
%% CG_RC_TEST02 tests CG_RC with the Wathen matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 January 2013
%
%  Author:
%
%    John Burkardt
%
  n = 79;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'CG_RC_TEST02\n' );
  fprintf ( 1, '  Use CG_RC to solve a linear system\n' );
  fprintf ( 1, '  involving the Wathen matrix.\n' );

  nx = 5;
  ny = 4;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  NX = %d\n', nx );
  fprintf ( 1, '  NY = %d\n', ny );
  fprintf ( 1, '  N  = %d\n', n );

  a = wathen ( nx, ny, n );

  seed = 123456789;
  [ x_exact, seed ] = r8vec_uniform_01 ( n, seed );

  size ( x_exact )
  size ( a )
  b(1:n) = a(1:n,1:n) * x_exact(1:n);

  x = zeros ( n, 1 );
%
%  Parameters we need for the stopping test.
%
  it = 0;
  it_max = 30;
  tol = 1.0E-05;
  bnrm2 = norm ( b(1:n) );
%
%  Set parameters for CG_RC.
%
  r = zeros ( n, 1 );
  z = zeros ( n, 1 );
  p = zeros ( n, 1 );
  q = zeros ( n, 1 );
  job = 1;
%
%  Repeatedly call CG_RC, and on return, do what JOB tells you.
%
  while ( 1 )

    [ x, r, z, p, q, job ] = cg_rc ( n, b, x, r, z, p, q, job );
%
%  Compute q = A * p.
%
    if ( job == 1 )

      q = a * p;
%
%  Solve M * z = r;
%
    elseif ( job == 2 )

      for i = 1 : n
        z(i) = r(i) / a(i,i);
      end
%
%  Compute r = r - A * x.
%
    elseif ( job == 3 )

      r = r - a * x;
%
%  Stopping test.
%
    elseif ( job == 4 )

      rnrm2 = norm ( r );

      if ( bnrm2 == 0.0 )
        if ( rnrm2 <= tol )
          break
        end
      else
        if ( rnrm2 <= tol * bnrm2 )
          break
        end
      end

      it = it + 1;

      if ( it_max <= it )
        fprintf ( 1, '\n' );
        fprintf ( 1, '  Iteration limit exceeded.\n' );
        fprintf ( 1, '  Terminating early.\n' );
        break
      end

    end

    job = 2;

  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of iterations was %d\n', it );
  fprintf ( 1, '  Estimated error is %g\n', rnrm2 );
  err = max ( abs ( x_exact(1:n) - x(1:n) ) );
  fprintf ( 1, '  Loo error is %g\n', err );

  fprintf ( 1, '\n' );
  fprintf ( 1, '     I      X(I)         X_EXACT(I)        B(I)\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '  %4d  %14f  %14f  %14f\n', ...
      i, x(i), x_exact(i), b(i) );
  end

  return
  end


text:
function Fitness = CalFitness(PopObj)
% Calculate the fitness of each solution

%------------------------------- 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".
%--------------------------------------------------------------------------

    N = size(PopObj,1);

    %% Detect the dominance relation between each two solutions
    Dominate = false(N);
    for i = 1 : N-1
        for j = i+1 : N
            k = any(PopObj(i,:)<PopObj(j,:)) - any(PopObj(i,:)>PopObj(j,:));
            if k == 1
                Dominate(i,j) = true;
            elseif k == -1
                Dominate(j,i) = true;
            end
        end
    end

    %% Calculate S(i)
    S = sum(Dominate,2);

    %% Calculate R(i)
    R = zeros(1,N);
    for i = 1 : N
        R(i) = sum(S(Dominate(:,i)));
    end

    %% Calculate D(i)
    Distance = pdist2(PopObj,PopObj);
    Distance(logical(eye(length(Distance)))) = inf;
    Distance = sort(Distance,2);
    D = 1./(Distance(:,floor(sqrt(N)))+2);

    %% Calculate the fitnesses
    Fitness = R + D';
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function determ = aegerter_determinant ( n )

%*****************************************************************************80
%
%% AEGERTER_DETERMINANT returns the determinant of the AEGERTER matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of rows and columns of A.
%
%    Output, real DETERM, the determinant.
%
  determ = ( n - ( ( n - 1 ) * n * ( 2 * n - 1 ) ) / 6 );

  return
end

author: johannesgerer


text:
function output = F_mean(input_layer, curr_layer)
input = input_layer.a;
if isfield(curr_layer, 'pool_idx')
    output = mean(input, curr_layer.pool_idx);
else
    output = mean(input,2);
end
end

author: singaxiong

repo: SignalGraph


text:
function out = shuffledMatrix(x)

    N = length(x(:,1));
    L = length(x(1,:));

    out = zeros(N,L);
    for i=1:L
        q = randperm(N);
        out(:,i) = x(q,i);
    end

author: gordonberman

repo: MotionMapper


text:
function [lon,lat]=nj_lonlat(ncRef,var);
% NJ_LONLAT  find longitude and latitude variables
% [lon,lat]=nj_lonlat(fname,var);
%    Inputs: fname = file name, url, or netcdf file object
%            var   = variable name (e.g. 'U10')
%         idouble  = [0] default, [1] convert output to double
%    Output: lon   = longitude array corresponding to "var"
%            lat   = latitude array corresponding to "var"

% import the NetCDF-Java methods we need for this example

% Rich Signell (rsignell@usgs.gov);

if nargin < 2, help(mfilename), return, end

try
    if (isa(ncRef, 'ncgeodataset')) %check for ncgeodataset Object
        nc = ncRef;
    else
        % open CF-compliant NetCDF File as a Common Data Model (CDM) "Grid Dataset"
        nc = ncgeodataset(ncRef);
    end
    gvar=nc.geovariable(var);
    lon=nc{gvar.getlonname}(:);
    lat=nc{gvar.getlatname}(:);
catch
    %gets the last error generated
    err = lasterror();
    disp(err.message);
end

author: nctoolbox

repo: nctoolbox

save_path: github-repos/MATLAB/nctoolbox-nctoolbox