text:
function printConstraints(model, minInf, maxInf, rxnSelection, modelAfter, printLevel)
% Print all network constraints that are between `-Inf (minInf)` or `+Inf
% (maxInf) inclusive`
%
% USAGE:
%
%    printConstraints(model, minInf, maxInf)
%
% INPUTS:
%    model:     COBRA model structure
%
% OPTIONAL INPUTS:
%   minInf:       value that is considered as -Inf (or desired minimum cutoff value)
%   maxInf:       value that is considered as +Inf (or desired maximum cutoff value)
%   rxnSelection: boolean vector or cell array of reaction abbreviations for the reactions to be printed
%   modelAfter:   model after some perturbation to the bounds
%   printLevel:
%

% .. Authors:
%       - Ines Thiele 02/09, Ronan Fleming 2020

if ~exist('minInf','var')|| isempty(minInf)
    minInf=-Inf;
end
if ~exist('maxInf','var') || isempty(maxInf)
    maxInf=Inf;
end
if ~isfield(model,'S')
    if isfield(model,'A')
        model.S = model.A(ismember(model.csense,'E'),:);
        model.dsense = model.csense(~ismember(model.csense,'E'));
        model.csense = model.csense(ismember(model.csense,'E'));
        model.C = model.A(~ismember(model.csense,'E'),:);
    else
        error('model.S missing')
    end
end
if ~exist('rxnSelection','var')
    rxnSelection=true(size(model.S,2),1);
end
if ~exist('printLevel','var')
    printLevel=0;
end
if exist('modelAfter','var')
    if isempty(modelAfter)
        clear modelAfter
    end
end

if ischar(rxnSelection) || iscell(rxnSelection)
    rxnSelection = ismember(model.rxns,rxnSelection);
end

if ~any(rxnSelection)
    return
end

closedRxnBool = model.lb == model.ub & model.lb == 0 & rxnSelection;
reversibleRxnBool = model.lb >= minInf & model.lb < 0 & model.ub <= maxInf & model.ub > 0 & rxnSelection;
%Forward and reverse reactions with NON-ZERO bounds
fwdRxnBoolNon0b = model.lb >= minInf & model.lb > 0 & ~reversibleRxnBool & ~closedRxnBool & model.ub <= maxInf & rxnSelection;
revRxnBoolNon0b = model.lb >= minInf & model.ub <= maxInf & model.ub < 0 & ~reversibleRxnBool & ~closedRxnBool & rxnSelection;
%Forward and reverse reactions with ZERO bounds (standard)
fwdRxnBool0b = model.lb >= minInf & model.lb == 0 & ~reversibleRxnBool & ~closedRxnBool & model.ub <= maxInf & rxnSelection;
revRxnBool0b = model.lb >= minInf & model.ub <= maxInf & model.ub == 0 & ~reversibleRxnBool & ~closedRxnBool & rxnSelection;

if ~any(closedRxnBool | reversibleRxnBool | fwdRxnBool0b | revRxnBool0b | fwdRxnBoolNon0b | revRxnBoolNon0b)
    boolRemainder = rxnSelection & ~(closedRxnBool | reversibleRxnBool | fwdRxnBool | revRxnBool);
    warning('no subset with bounds between [minInf maxInf]')
else
    boolRemainder=0;
end

if any(closedRxnBool & reversibleRxnBool) || any(closedRxnBool & fwdRxnBool0b) || any(closedRxnBool & revRxnBool0b) || any(fwdRxnBool0b & revRxnBool0b)
    warning('inconsistent boolean variables')
end

if isfield(model,'rxnNames')
rxnNames=model.rxnNames;
else
    rxnNames=cell(size(model.S,2),1);
end
for j=1:size(model.S,2)
    rxnNames{j}=rxnNames{j}(1:min(60,length(rxnNames{j})));
end

if ~any(closedRxnBool)
    if printLevel>0
        fprintf('%s\n','No closed reactions with non-default constraints.');
    end
else
    if printLevel>0
        fprintf('%s\n', ['...closed reaction constraints:']);
    end
    if exist('modelAfter','var')
        T = table(model.rxns(closedRxnBool),rxnNames(closedRxnBool),model.lb(closedRxnBool),modelAfter.lb(closedRxnBool),model.ub(closedRxnBool),modelAfter.ub(closedRxnBool),printRxnFormula(model, 'rxnAbbrList',model.rxns(closedRxnBool),'printFlag',0),'VariableNames',{'Closed_Reaction','Name','lb_before','lb_after','ub_before','ub_after','equation'});
    else
        T = table(model.rxns(closedRxnBool),rxnNames(closedRxnBool),model.lb(closedRxnBool),model.ub(closedRxnBool),printRxnFormula(model, 'rxnAbbrList', model.rxns(closedRxnBool),'printFlag',0),'VariableNames',{'Closed_reaction','Name','lb','ub','equation'});
    end
    disp(T);
end

if ~any(fwdRxnBool0b)
    if printLevel>0
        fprintf('%s\n','No forward reactions with non-default constraints.');
    end
else
    if printLevel>0
        fprintf('%s\n', ['...forward reactions with non-[0, ' num2str(maxInf) '] constraints:']);
    end
    if exist('modelAfter','var')
        T = table(model.rxns(fwdRxnBool0b),rxnNames(fwdRxnBool0b),model.lb(fwdRxnBool0b),modelAfter.lb(fwdRxnBool0b),model.ub(fwdRxnBool0b),modelAfter.ub(fwdRxnBool0b),printRxnFormula(model, 'rxnAbbrList',model.rxns(fwdRxnBool0b),'printFlag',0),'VariableNames',{'Forward_Reaction, 0 bound','Name','lb_before','lb_after','ub_before','ub_after','equation'});
    else
        T = table(model.rxns(fwdRxnBool0b),rxnNames(fwdRxnBool0b),model.lb(fwdRxnBool0b),model.ub(fwdRxnBool0b),printRxnFormula(model, 'rxnAbbrList', model.rxns(fwdRxnBool0b),'printFlag',0),'VariableNames',{'Forward_Reaction, 0 bound','Name','lb','ub','equation'});
    end
    disp(T);
end

if ~any(revRxnBool0b)
    if printLevel>0
        fprintf('%s\n','No reverse reactions with non-default constraints.');
    end
else
    if printLevel>0
        fprintf('%s\n',['...reverse reactions with non-[' num2str(minInf)  ', 0]  constraints:']);
    end
    if exist('modelAfter','var')
        T = table(model.rxns(revRxnBool0b),rxnNames(revRxnBool0b),model.lb(revRxnBool0b),modelAfter.lb(revRxnBool0b),model.ub(revRxnBool0b),modelAfter.ub(revRxnBool0b),printRxnFormula(model, 'rxnAbbrList', model.rxns(revRxnBool0b),'printFlag',0),'VariableNames',{'Reverse_Reaction, 0 bound','Name','lb_before','lb_after','ub_before','ub_after','equation'});
    else
        T = table(model.rxns(revRxnBool0b),rxnNames(revRxnBool0b),model.lb(revRxnBool0b),model.ub(revRxnBool0b),printRxnFormula(model, 'rxnAbbrList', model.rxns(revRxnBool0b),'printFlag',0),'VariableNames',{'Reverse_Reaction, 0 bound','Name','lb','ub','equation'});
    end
    disp(T);
end

if ~any(fwdRxnBoolNon0b)
    if printLevel>0
        fprintf('%s\n','No forward reactions with non-default constraints.');
    end
else
    if printLevel>0
        fprintf('%s\n', ['...forward reactions with non-[0, ' num2str(maxInf) '] constraints:']);
    end
    if exist('modelAfter','var')
        T = table(model.rxns(fwdRxnBoolNon0b),rxnNames(fwdRxnBoolNon0b),model.lb(fwdRxnBoolNon0b),modelAfter.lb(fwdRxnBoolNon0b),model.ub(fwdRxnBoolNon0b),modelAfter.ub(fwdRxnBoolNon0b),printRxnFormula(model, 'rxnAbbrList',model.rxns(fwdRxnBoolNon0b),'printFlag',0),'VariableNames',{'Forward_Reaction, non-0 bound','Name','lb_before','lb_after','ub_before','ub_after','equation'});
    else
        T = table(model.rxns(fwdRxnBoolNon0b),rxnNames(fwdRxnBoolNon0b),model.lb(fwdRxnBoolNon0b),model.ub(fwdRxnBoolNon0b),printRxnFormula(model, 'rxnAbbrList', model.rxns(fwdRxnBoolNon0b),'printFlag',0),'VariableNames',{'Forward_Reaction, non-0 bound','Name','lb','ub','equation'});
    end
    disp(T);
end

if ~any(revRxnBoolNon0b)
    if printLevel>0
        fprintf('%s\n','No reverse reactions with non-default constraints.');
    end
else
    if printLevel>0
        fprintf('%s\n',['...reverse reactions with non-[' num2str(minInf)  ', 0]  constraints:']);
    end
    if exist('modelAfter','var')
        T = table(model.rxns(revRxnBoolNon0b),rxnNames(revRxnBoolNon0b),model.lb(revRxnBoolNon0b),modelAfter.lb(revRxnBoolNon0b),model.ub(revRxnBoolNon0b),modelAfter.ub(revRxnBoolNon0b),printRxnFormula(model, 'rxnAbbrList', model.rxns(revRxnBoolNon0b),'printFlag',0),'VariableNames',{'Reverse_Reaction, non-0 bound','Name','lb_before','lb_after','ub_before','ub_after','equation'});
    else
        T = table(model.rxns(revRxnBoolNon0b),rxnNames(revRxnBoolNon0b),model.lb(revRxnBoolNon0b),model.ub(revRxnBoolNon0b),printRxnFormula(model, 'rxnAbbrList', model.rxns(revRxnBoolNon0b),'printFlag',0),'VariableNames',{'Reverse_Reaction, non-0 bound','Name','lb','ub','equation'});
    end
    disp(T);
end

if ~any(reversibleRxnBool)
    if printLevel>0
        fprintf('%s\n','No reversible reactions with non-default constraints.');
    end
else
    if printLevel>0
        fprintf('%s\n',['...reversible reactions with non-[' num2str(minInf)  ', ' num2str(maxInf) ']  constraints:']);
    end
    if exist('modelAfter','var')
        T = table(model.rxns(reversibleRxnBool),rxnNames(reversibleRxnBool),model.lb(reversibleRxnBool),modelAfter.lb(reversibleRxnBool),model.ub(reversibleRxnBool),modelAfter.ub(reversibleRxnBool),printRxnFormula(model, 'rxnAbbrList', model.rxns(reversibleRxnBool),'printFlag',0),'VariableNames',{'Reversible_Reaction','Name','lb_before','lb_after','ub_before','ub_after','equation'});
    else
        T = table(model.rxns(reversibleRxnBool),rxnNames(reversibleRxnBool),model.lb(reversibleRxnBool),model.ub(reversibleRxnBool),printRxnFormula(model, 'rxnAbbrList', model.rxns(reversibleRxnBool),'printFlag',0),'VariableNames',{'Reversible_Reaction','Name','lb','ub','equation'});
    end
    disp(T);
end

% if ~any(reversibleRxnBool)
%     if printLevel>0
%         fprintf('%s\n','No reversible reactions with non-default constraints.');
%     end
% else
%     if printLevel>0
%         fprintf('%s\n',['...reversible reactions with non-[' num2str(minInf)  ', ' num2str(maxInf) ']  constraints:']);
%     end
%     if exist('modelAfter','var')
%         T = table(model.rxns(reversibleRxnBool),rxnNames(reversibleRxnBool),model.lb(reversibleRxnBool),modelAfter.lb(reversibleRxnBool),model.ub(reversibleRxnBool),modelAfter.ub(reversibleRxnBool),printRxnFormula(model, 'rxnAbbrList', model.rxns(reversibleRxnBool),'printFlag',0),'VariableNames',{'Reversible_Reaction','Name','lb_before','lb_after','ub_before','ub_after','equation'});
%     else
%         T = table(model.rxns(reversibleRxnBool),rxnNames(reversibleRxnBool),model.lb(reversibleRxnBool),model.ub(reversibleRxnBool),printRxnFormula(model, 'rxnAbbrList', model.rxns(reversibleRxnBool),'printFlag',0),'VariableNames',{'Reversible_Reaction','Name','lb','ub','equation'});
%     end
%     disp(T);
% end

if any(boolRemainder)
    if printLevel>0
        fprintf('%s\n',['...reversible reactions with non-[' num2str(minInf)  ', ' num2str(maxInf) ']  constraints:']);
    end
    if exist('modelAfter','var')
        T = table(model.rxns(boolRemainder),rxnNames(boolRemainder),model.lb(boolRemainder),modelAfter.lb(boolRemainder),model.ub(boolRemainder),modelAfter.ub(boolRemainder),printRxnFormula(model, 'rxnAbbrList', model.rxns(boolRemainder),'printFlag',0),'VariableNames',{'Reversible_Reaction','Name','lb_before','lb_after','ub_before','ub_after','equation'});
    else
        T = table(model.rxns(boolRemainder),rxnNames(boolRemainder),model.lb(boolRemainder),model.ub(boolRemainder),printRxnFormula(model, 'rxnAbbrList', model.rxns(boolRemainder),'printFlag',0),'VariableNames',{'Reversible_Reaction','Name','lb','ub','equation'});
    end
    disp(T);
end

fprintf('\n');

author: opencobra

repo: cobratoolbox


text:
function [m_0, m_1, m_2] = quad_moments(fun, a, b, rtol, atol, minsubs)
% QUAD_MOMENTS Calculate the 0th, 1st and 2nd moment of a given
%              (unnormalized) probability distribution
%
%   [m_0, m_1, m_2] = quad_moments(fun, a, b, varargin)
%   Inputs:
%      fun  = Function handle to the unnormalized probability distribution
%      a,b  = integration limits [a,b]
%      rtol = relative tolerance for the integration (optional, default 1e-6)
%      atol = absolute tolerance for the integration (optional, default 1e-10)
%
%   Returns the first three moments:
%      m0  = int_a^b fun(x) dx
%      m1  = int_a^b x*fun(x) dx / m0
%      m2  = int_a^b x^2*fun(x) dx / m0
%
%   The function uses an adaptive Gauss-Kronrod quadrature. The same set of
%   integration points and intervals are used for each moment. This speeds up
%   the evaluations by factor 3, since the function evaluations are done only
%   once.
%
%   The quadrature method is described by:
%   L.F. Shampine, "Vectorized Adaptive Quadrature in Matlab",
%   Journal of Computational and Applied Mathematics, 211, 2008,
%   pp. 131-140.

%   Copyright (c) 2010 Jarno Vanhatalo, Jouni Hartikainen

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

    maxsubs = 650;

    if nargin < 4
        rtol = 1.e-6;
    end
    if nargin < 5
        atol = 1.e-10;
    end
    if nargin < 6
        minsubs = 10;
    end

    rtol = max(rtol,100*eps);
    atol = max(atol,0);
    minsubs = max(minsubs,2); % At least two subintervals are needed

    % points and weights
    points15 = [0.2077849550078985; 0.4058451513773972; 0.5860872354676911; ...
        0.7415311855993944; 0.8648644233597691; 0.9491079123427585; ...
        0.9914553711208126];
    points = [-points15(end:-1:1); 0; points15];

    w15 = [0.2044329400752989, 0.1903505780647854, 0.1690047266392679, ...
        0.1406532597155259, 0.1047900103222502, 0.06309209262997855, ...
        0.02293532201052922];
    w = [w15(end:-1:1), 0.2094821410847278, w15];

    w7 = [0,0.3818300505051189,0,0.2797053914892767,0,0.1294849661688697,0];
    ew = w - [w7(end:-1:1), 0.4179591836734694, w7];

    samples = numel(w);

    % split the interval.
    if b-a <= 0
        c = a; a = b; b=c;
        warning('The start of the integration interval was less than the end of it.')
    end
    apu = a + (1:(minsubs-1))./minsubs*(b-a);
    apu = [a,apu,b];
    subs = [apu(1:end-1);apu(2:end)];

    % Initialize partial sums.
    Ifx_ok = 0;
    Ifx1_ok = 0;
    Ifx2_ok = 0;
    % The main loop
    while true
        % subintervals and their midpoints
        midpoints = sum(subs)/2;
        halfh = diff(subs)/2;
        x = bsxfun(@plus,points*halfh,midpoints);
        x = reshape(x,1,[]);

        fx = fun(x);
        fx1 = fx.*x;
        fx2 = fx.*x.^2;

        fx = reshape(fx,samples,[]);
        fx1 = reshape(fx1,samples,[]);
        fx2 = reshape(fx2,samples,[]);

        % Subintegrals.
        Ifxsubs = (w*fx) .* halfh;
        errsubs = (ew*fx) .* halfh;
        Ifxsubs1 = (w*fx1) .* halfh;
        Ifxsubs2 = (w*fx2) .* halfh;

        % Ifx and tol.
        Ifx = sum(Ifxsubs) + Ifx_ok;
        Ifx1 = sum(Ifxsubs1) + Ifx1_ok;
        Ifx2 = sum(Ifxsubs2) + Ifx2_ok;
        tol = max(atol,rtol*abs(Ifx));

        % determine the indices ndx of Ifxsubs for which the
        % errors are acceptable and remove those from subs
        ndx = find(abs(errsubs) <= (2/(b-a)*halfh*tol));
        subs(:,ndx) = [];
        if isempty(subs)
            break
        end

        % Update the integral.
        Ifx_ok = Ifx_ok + sum(Ifxsubs(ndx));
        Ifx1_ok = Ifx1_ok + sum(Ifxsubs1(ndx));
        Ifx2_ok = Ifx2_ok + sum(Ifxsubs2(ndx));

        % Quit if too many subintervals.
        nsubs = 2*size(subs,2);
        if nsubs > maxsubs
            warning('quad_moments: Reached the limit on the maximum number of intervals in use.');
            break
        end
        midpoints(ndx) = [];
        subs = reshape([subs(1,:); midpoints; midpoints; subs(2,:)],2,[]); % Divide the remaining subintervals in half
    end

    % Scale moments
    m_0 = Ifx;
    m_1 = Ifx1./Ifx;
    m_2 = Ifx2./Ifx;
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function asa147_test ( )

%*****************************************************************************80
%
%% ASA147_TEST tests ASA147.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'ASA147_TEST\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the ASA147 library.\n' );

  asa147_test01 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'ASA147_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function CPDpot = convert_dbn_CPDs_to_tables(bnet, evidence)
% CONVERT_DBN_CPDS_TO_TABLES Convert CPDs of (possibly instantiated) DBN nodes to tables
% CPDpot = convert_dbn_CPDs_to_tables(bnet, evidence)
%
% CPDpot{n,t} is a table containing P(n,t|pa(n,t), ev)
% All hidden nodes are assumed to be discrete.
% We assume the observed nodes are the same in every slice.
%
% Evaluating the conditional likelihood of long evidence sequences can be very slow,
% so we take pains to vectorize where possible.

[ss T] = size(evidence);
%obs_bitv = ~isemptycell(evidence(:));
obs_bitv = zeros(1, 2*ss);
obs_bitv(bnet.observed) = 1;
obs_bitv(bnet.observed+ss) = 1;

ns = bnet.node_sizes(:);
CPDpot = cell(ss,T);

for n=1:ss
  % slice 1
  t = 1;
  ps = parents(bnet.dag, n);
  e = bnet.equiv_class(n, 1);
  if ~any(obs_bitv(ps))
    CPDpot{n,t} = convert_CPD_to_table_hidden_ps(bnet.CPD{e}, evidence{n,t});
  else
    CPDpot{n,t} = convert_to_table(bnet.CPD{e}, [ps n], evidence(:,1));
  end

% special cases: c=child, p=parents, d=discrete, h=hidden, 1sl=1slice
% if c=h=1 then c=d=1, since hidden nodes must be discrete
% c=h c=d p=h p=d 1sl method
% ---------------------------
% 1   1   1   1   -   replicate CPT
% -   1   -   1   -   evaluate CPT on evidence *
% 0   1   1   1   1   dhmm
% 0   0   1   1   1   ghmm
% other               loop
%
% * = any subset of the domain may be observed

% Example where all of the special cases occur - a hierarchical HMM
% where the top layer (G) and leaves (Y) are observed and
% all nodes are discrete except Y.
% (O turns on if Y is an outlier)

% G ---------> G
% |            |
% v            v
% S  --------> S
% |            |
% v            v
% Y            Y
% ^            ^
% |            |
% O            O

% Evaluating P(yt|St,Ot) is the ghmm case
% Evaluating P(St|S(t-1),gt) is the eval CPT case
% Evaluating P(gt|g(t-1) is the eval CPT case (hdom = [])
% Evaluating P(Ot) is the replicated CPT case

% Cts parents (e.g., inputs) would require an additional special case for speed

  % slices 2..T
  [ss T] = size(evidence);
  self = n+ss;
  ps = parents(bnet.dag, self);
  e = bnet.equiv_class(n, 2);

  if 1
  debug = 0;
  hidden_child = ~obs_bitv(n);
  discrete_child = myismember(n, bnet.dnodes);
  hidden_ps = all(~obs_bitv(ps));
  discrete_ps = mysubset(ps, bnet.dnodes);
  parents_in_same_slice = all(ps > ss);

  if hidden_child & discrete_child & hidden_ps & discrete_ps
    CPDpot = helper_repl(bnet, evidence, n, CPDpot, obs_bitv, debug);
  elseif discrete_child & discrete_ps
    CPDpot = helper_eval(bnet, evidence, n, CPDpot, obs_bitv, debug);
  elseif discrete_child & hidden_ps & discrete_ps & parents_in_same_slice
    CPDpot = helper_dhmm(bnet, evidence, n, CPDpot, obs_bitv, debug);
  elseif ~discrete_child & hidden_ps & discrete_ps & parents_in_same_slice
    CPDpot = helper_ghmm(bnet, evidence, n, CPDpot, obs_bitv, debug);
  else
    if debug, fprintf('node %d, slow\n', n); end
    for t=2:T
      CPDpot{n,t} = convert_to_table(bnet.CPD{e}, [ps self], evidence(:,t-1:t));
    end
  end
  end

  if 0
  for t=2:T
    CPDpot2{n,t} = convert_to_table(bnet.CPD{e}, [ps self], evidence(:,t-1:t));
    if ~approxeq(CPDpot{n,t}, CPDpot2{n,t})
      fprintf('CPDpot n=%d, t=%d\n',n,t);
      keyboard
    end
  end
  end

end

%%%%%%%
function CPDpot = helper_repl(bnet, evidence, n, CPDpot, obs_bitv, debug)

[ss T] = size(evidence);
if debug, fprintf('node %d, repl\n', n); end
e = bnet.equiv_class(n, 2);
CPT = convert_CPD_to_table_hidden_ps(bnet.CPD{e}, []);
CPDpot(n,2:T) = num2cell(repmat(CPT, [1 1 T-1]), [1 2]);

%%%%%%%
function CPDpot = helper_eval(bnet, evidence, n, CPDpot, obs_bitv, debug)

[ss T] = size(evidence);
self = n+ss;
ps = parents(bnet.dag, self);
e = bnet.equiv_class(n, 2);
ns = bnet.node_sizes(:);
% Example: given CPT(p1, p2, p3, p4, c), where p1,p3 are observed
% we create CPT([p2 p4 c], [p1 p3]).
% We then convert all observed p1,p3 into indices ndx
% and return CPT(:, ndx)
CPT = CPD_to_CPT(bnet.CPD{e});
domain = [ps self];
% if dom is [3 7 8] and 3,8 are observed, odom_rel = [1 3], hdom_rel = 2,
% odom = [3 8], hdom = 7
odom_rel = find(obs_bitv(domain));
hdom_rel = find(~obs_bitv(domain));
odom = domain(odom_rel);
hdom = domain(hdom_rel);
if isempty(hdom)
  CPT = CPT(:);
else
  CPT = permute(CPT, [hdom_rel odom_rel]);
  CPT = reshape(CPT, prod(ns(hdom)), prod(ns(odom)));
end
parents_in_same_slice = all(ps > ss);
if parents_in_same_slice
  if debug, fprintf('node %d eval 1 slice\n', n); end
  data = cell2num(evidence(odom-ss,2:T)); %data(i,t) = val of i'th obs parent at t+1
else
  if debug, fprintf('node %d eval 2 slice\n', n); end
  % there's probably a way of vectorizing this...
  data = zeros(length(odom), T-1);
  for t=2:T
    ev = evidence(:,t-1:t);
    ev = ev(:);
    ev2 = ev(odom);
    data(:,t-1) = cat(1, ev2{:});
    %data(:,t-1) = cell2num(ev2);
  end
end
ndx = subv2ind(ns(odom), data'); % ndx(t) encodes data(:,t)
if isempty(hdom)
  CPDpot(n,2:T) = num2cell(CPT(ndx)); % a cell array of floats
else
  CPDpot(n,2:T) = num2cell(CPT(:, ndx), 1); % a cell array of column vectors
end

%%%%%%%
function CPDpot = helper_dhmm(bnet, evidence, n, CPDpot, obs_bitv, debug)

if debug, fprintf('node %d, dhmm\n', n); end
[ss T] = size(evidence);
self = n+ss;
ps = parents(bnet.dag, self);
e = bnet.equiv_class(n, 2);
ns = bnet.node_sizes(:);
CPT = CPD_to_CPT(bnet.CPD{e});
CPT = reshape(CPT, [prod(ns(ps)) ns(self)]); % what if no parents?
%obslik = mk_dhmm_obs_lik(cell2num(evidence(n,2:T)), CPT);
obslik = eval_pdf_cond_multinomial(cell2num(evidence(n,2:T)), CPT);
CPDpot(n,2:T) = num2cell(obslik, 1);

%%%%%%%
function CPDpot = helper_ghmm(bnet, evidence, n, CPDpot, obs_bitv, debug)

if debug, fprintf('node %d, ghmm\n', n); end
[ss T] = size(evidence);
e = bnet.equiv_class(n, 2);
S = struct(bnet.CPD{e});
ev2 = cell2num(evidence(n,2:T));
%obslik = mk_ghmm_obs_lik(ev2, S.mean, S.cov);
obslik = eval_pdf_cond_gauss(ev2, S.mean, S.cov);
CPDpot(n,2:T) = num2cell(obslik, 1);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function blend_test06 ( )

%*****************************************************************************80
%
%% BLEND_TEST06 checks out BLEND_IJK_1D1
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 October 2008
%
%  Author:
%
%    John Burkardt
%
  m1 = 4;
  m2 = 3;
  m3 = 3;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'BLEND_TEST06\n' );
  fprintf ( 1, '  BLEND_IJK_1D1 interpolates data in a table,\n' );
  fprintf ( 1, '  from edge data.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The table is %d rows by %d columns by %d layers.\n', ...
    m1, m2, m3 );
%
%  Load data on the faces.
%
  for i = 1 : m1
    r = ( i - 1 ) / ( m1 - 1 );
    for j = 1 : m2
      s = ( j - 1 ) / ( m2 - 1 );
      for k = 1 : m3
        t = ( k - 1 ) / ( m3 - 1 );

        num_extreme = 0;
        if ( i == 1 | i == m1 )
          num_extreme = num_extreme + 1;
        end
        if ( j == 1 | j == m2 )
          num_extreme = num_extreme + 1;
        end
        if ( k == 1 | k == m3 )
          num_extreme = num_extreme + 1;
        end

        if ( 2 <= num_extreme )
          x(i,j,k) = quad_rst ( r, s, t, 1 );
        else
          x(i,j,k) = 0.0;
        end

      end
    end
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Data given to BLEND_IJK_1D1:\n' );
  fprintf ( 1, '\n' );

  for k = 1 : m3
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Layer K =  %d\n', k );
    fprintf ( 1, '\n' );
    for i = 1 : m1
      fprintf ( 1, '  %10f  %10f  %10f\n', x(i,1:m2,k) );
    end
  end

  x = blend_ijk_1d1 ( x, m1, m2, m3 );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Values interpolated by BLEND_IJK_1D1:\n' );
  fprintf ( 1, '\n' );

  for k = 1 : m3
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Layer K =  %d\n', k );
    fprintf ( 1, '\n' );
    for i = 1 : m1
      fprintf ( 1, '  %10f  %10f  %10f\n', x(i,1:m2,k) );
    end
  end
%
%  Load all data.
%
  for i = 1 : m1
    r = ( i - 1 ) / ( m1 - 1 );
    for j = 1 : m2
      s = ( j - 1 ) / ( m2 - 1 );
      for k = 1 : m3
        t = ( k - 1 ) / ( m3 - 1 );
        x(i,j,k) = quad_rst ( r, s, t, 1 );
      end
    end
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Exact data:\n' );
  fprintf ( 1, '\n' );

  for k = 1 : m3
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Layer K =  %d\n', k );
    fprintf ( 1, '\n' );
    for i = 1 : m1
      fprintf ( 1, '  %10f  %10f  %10f\n', x(i,1:m2,k) );
    end
  end

  return
end


text:
function stats = grayrlprops(varargin)

%GRAYCOPROPS Properties of gray-level run-length matrix.
%  -------------------------------------------
%  STATS = GRAYCOPROPS(GLRLM,PROPERTIES) Each element in  GLRLM, (r,c),
%   is the probability occurrence of pixel having gray level values r, run-length c in the image.
%   GRAYCOPROPS is to calculate PROPERTIES.
%  -------------------------------------------
%  Requirements:
%  -------------------------------------------
%   GLRLM mustbe an cell array of valid gray-level run-length
%   matrices.Recall that a valid glrlm must be logical or numerical.
%  -------------------------------------------
%  Current supported statistics include:
%  -------------------------------------------
%   Short Run Emphasis (SRE)
%   Long Run Emphasis (LRE)
%   Gray-Level Nonuniformity (GLN)
%   Run Length Nonuniformity (RLN)
%   Run Percentage (RP)
%   Low Gray-Level Run Emphasis (LGRE)
%   High Gray-Level Run Emphasis (HGRE)
%   Short Run Low Gray-Level Emphasis (SRLGE)
%   Short Run High Gray-Level Emphasis (SRHGE)
%   Long Run Low Gray-Level Emphasis (LRLGE)
%   Long Run High Gray-Level Emphasis (LRHGE)
%  --------------------------------------------
%  Reference:
%  --------------------------------------------
%   Xiaoou Tang,Texture Information in Run-Length Matrices
%   IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL.7, NO.11,NOVEMBER 1998
% ---------------------------------------------
%  See also GRAYRLMATRIX.
% ---------------------------------------------
% Author:
% ---------------------------------------------
%    (C)Xunkai Wei <xunkai.wei@gmail.com>
%    Beijing Aeronautical Technology Research Center
%    Beijing %9203-12,10076
% ---------------------------------------------
% History:
% ---------------------------------------------
% Creation: beta         Date: 01/10/2007
% Revision: 1.0          Date: 12/11/2007
% 1.Accept cell input now
% 2.Using MATLAB file style
% 3.Fully vectorized programming
% 4.Fully support the IEEE reference
% 5. ...

% Check GLRLM
[GLRLM numGLRLM] = ParseInputs(varargin{:});

% Initialize output stats structure.
% 11 statistics for each GLRLM
numStats = 11;

% % count number of GLRLM
% numGLRLM = length(GLRLM);

% Initialization default 4*11 matrix
stats = zeros(numGLRLM,numStats);

for p = 1 : numGLRLM
    %N-D indexing not allowed for sparse.

    if numGLRLM ~= 1
        % transfer to double matrix
        tGLRLM = GLRLM{p};
    else
        tGLRLM = GLRLM;
    end
    %     if numGLRLM ~= 1
    %         % transfer to double matrix
    %         tGLRLM = normalizeGLRL(GLRLM{p});
    %     else
    %         tGLRLM = normalizeGLRL(GLRLM);
    %     end
    % Get row and column subscripts of GLRLM.  These subscripts correspond to the
    % pixel values in the GLRLM.
    s = size(tGLRLM);
    % colum indicator
    c_vector =1:s(1);
    % row indicator
    r_vector =1:s(2);
    % matrix element indicator
    % Matrix form col and row: using meshgrid, you should transpose before using
    % i.e. if tGLRLM is m*n, then this function return c_matrix n*m,
    % r_matrix n*m.
    [c_matrix,r_matrix] = meshgrid(c_vector,r_vector);

    % Total number of runs
    N_runs = sum(sum(tGLRLM));

    % total number of elements
    N_tGLRLM = s(1)*s(2);

    %--------------------Prepare four matrix for speedup--------------
    % 1.Gray Level Run-Length Pixel Number Matrix
    %     p_p = calculate_p_p(tGLRLM,c_matrix');

    % 2.Gray-Level Run-Number Vector
    %   This vector represents the sum distribution of the number of runs
    %   with gray level i.
    p_g = sum(tGLRLM);

    % 3.Run-Length Run-Number Vector
    %   This vector represents the sum distribution of the number of runs
    %   with run length j.
    p_r = sum(tGLRLM,2)';

    % 4.Gray-Level Run-Length-One Vector
    %
    % p_o = tGLRLM(:,1); % Not used yet
    % ----------------------End four matrix---------------------------
    %
    %------------------------Statistics-------------------------------
    % 1. Short Run Emphasis (SRE)
    SRE = sum(p_r./(c_vector.^2))/N_runs;
    % 2. Long Run Emphasis (LRE)
    LRE = sum(p_r.*(c_vector.^2))/N_runs;
    % 3. Gray-Level Nonuniformity (GLN)
    GLN = sum(p_g.^2)/N_runs;
    % 4. Run Length Nonuniformity (RLN)
    RLN = sum(p_r.^2)/N_runs;
    % 5. Run Percentage (RP)
    RP = N_runs/N_tGLRLM;
    % 6. Low Gray-Level Run Emphasis (LGRE)
    LGRE = sum(p_g./(r_vector.^2))/N_runs;
    % 7. High Gray-Level Run Emphasis (HGRE)
    HGRE = sum(p_g.*r_vector.^2)/N_runs;
    % 8. Short Run Low Gray-Level Emphasis (SRLGE)
    SGLGE =calculate_SGLGE(tGLRLM,r_matrix',c_matrix',N_runs);
    % 9. Short Run High Gray-Level Emphasis (SRHGE)
    SRHGE =calculate_SRHGE(tGLRLM,r_matrix',c_matrix',N_runs);
    % 10. Long Run Low Gray-Level Emphasis (LRLGE)
    LRLGE =calculate_LRLGE(tGLRLM,r_matrix',c_matrix',N_runs);
    % 11.Long Run High Gray-Level Emphasis (LRHGE
    LRHGE =calculate_LRHGE(tGLRLM,r_matrix',c_matrix',N_runs);
    %----------------insert statistics----------------------------
    stats(p,:)=[SRE LRE GLN RLN  RP LGRE HGRE SGLGE SRHGE LRLGE  LRHGE ];
end % end all run length matrixs

%   ----------------------Utility functions--------------------
%-----------------------------------------------------------------------------
% function glrl = normalizeGLRL(glrl)
%
% % Normalize glcm so that sum(glcm(:)) is one.
% if any(glrl(:))
%   glrl = glrl ./ sum(glrl(:));
% end
% function p_p = calculate_p_p(GLRLM,c) % Note: currently not used
%
% % p_p(i; j) = GLRLM(i,j)*j
% % Each element of the matrix represents the number of pixels of run length
% % j and gray-level i. Compared to the original matrix, the new matrix gives
% % equal emphasis to all length of runs in an image.
%
% term1 =  c; % j index in matrix size
% term2 = GLRLM;
% p_p = term1 .* term2;
%---------------------------------
function SGLGE =calculate_SGLGE(tGLRLM,r_matrix,c_matrix,N_runs)
% Short Run Low Gray-Level Emphasis (SRLGE):

term = tGLRLM./((r_matrix.*c_matrix).^2);
SGLGE= sum(sum(term))./N_runs;

%------------------------------------
function  SRHGE =calculate_SRHGE(tGLRLM,r_matrix,c_matrix,N_runs)
% Short Run High Gray-Level Emphasis (SRHGE):
%
term  = tGLRLM.*(r_matrix.^2)./(c_matrix.^2);
SRHGE = sum(sum(term))/N_runs;
%------------------------------------
function   LRLGE =calculate_LRLGE(tGLRLM,r_matrix,c_matrix,N_runs)
% Long Run Low Gray-Level Emphasis (LRLGE):
%
term  = tGLRLM.*(c_matrix.^2)./(r_matrix.^2);
LRLGE = sum(sum(term))/N_runs;
%---------------------------------------
function  LRHGE =calculate_LRHGE(tGLRLM,r_matrix,c_matrix,N_runs)
% Long Run High Gray-Level Emphasis (LRHGE):
%
term  = tGLRLM.*(c_matrix.^2).*(r_matrix.^2);
LRHGE = sum(sum(term))/N_runs;
%----------------------------------------

%-----------------------------------------------------------------------------
function [glrlm num_glrlm] = ParseInputs(varargin)
% check stability of inputs
%
% first receive all inputs
glrlm = varargin{:};
% get numbers total
num_glrlm=length(glrlm);
% then for each element, check its stability
for i=1:num_glrlm
    % The 'nonnan' and 'finite' attributes are not added to iptcheckinput because the
    % 'integer' attribute takes care of these requirements.
    % iptcheckinput(glrlm,{'cell'},{'real','nonnegative','integer'}, ...
    % mfilename,'GLRLM',1);
    iptcheckinput(glrlm{i},{'logical','numeric'},{'real','nonnegative','integer'},...
        mfilename,'GLRLM',1);
    % Cast GLRLM to double to avoid truncation by data type. Note that GLRLM is not an
    % image.
    if ~isa(glrlm,'double')
        glrlm{i}= double(glrlm{i});
    end
end

author: mvallieres

repo: radiomics


text:
function rv = run_userfcn(userfcn, stage, varargin)
%RUN_USERFCN   Runs the userfcn callbacks for a given stage.
%   RV = RUN_USERFCN(USERFCN, STAGE, VARARGIN)
%
%   USERFCN : the 'userfcn' field of mpc, populated by ADD_USERFCN
%   STAGE   : the name of the callback stage being executed
%   (additional arguments) some stages require additional arguments.
%
%   Example:
%       mpc = om.get_mpc();
%       om = run_userfcn(mpc.userfcn, 'formulation', om);
%
%   See also ADD_USERFCN, REMOVE_USERFCN, TOGGLE_RESERVES, TOGGLE_IFLIMS,
%   RUNOPF_W_RES.

%   MATPOWER
%   Copyright (c) 2009-2016, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

rv = varargin{1};
if ~isempty(userfcn) && isfield(userfcn, stage)
    for k = 1:length(userfcn.(stage))
        if isfield(userfcn.(stage)(k), 'args')
            args = userfcn.(stage)(k).args;
        else
            args = [];
        end
        rv = feval(userfcn.(stage)(k).fcn, rv, varargin{2:end}, args);
        % mpc     = userfcn_*_ext2int(mpc, mpopt, args);
        % om      = userfcn_*_formulation(om, mpopt, args);
        % results = userfcn_*_int2ext(results, mpopt, args);
        % results = userfcn_*_printpf(results, fd, mpopt, args);
        % mpc     = userfcn_*_savecase(mpc, fd, prefix, args);
    end
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
% MAIN_passiveSimulate.m
%
% This script runs a passive simulation of the acrobot, to do a sanity
% check on the dynamics and plotting routines
%
% Things to try:
%
% 1) Adjust the parameters and states and see if you can guess what will
% happen when you run the simulation. What if one mass or length is much
% bigger than the other? Can you make the system behave like a single
% pendulum?
%
% 2) Type:      >> help animate
% to see the keyboard commands for controlling the animation. You can use
% then to see the system in slow motion, pause, or go back in time, to name
% a few possiblities. This is useful for understanding what the system is
% doing.
%

clc; clear;

% Physical parameters
p.m1 = 10;
p.m2 = 20;
p.g = 9;
p.l1 = 3;
p.l2 = 5;

% Initial state:
q1 = (pi/180)*120;
q2 = (pi/180)*50;
dq1 = 0;
dq2 = 0;
z0 = [q1;q2;dq1;dq2];  %Pack up initial state

tSpan = [0,8];  %time span for the simulation
dynFun = @(t,z)( acrobotDynamics(z,0,p) );  %passive dynamics function

% Run simulation:
sol = ode45(dynFun,tSpan,z0);
t = linspace(tSpan(1),tSpan(2),500);
z = deval(sol,t);
u = zeros(size(t));

% Animate the results:
A.plotFunc = @(t,z)( drawAcrobot(t,z,p) );
A.speed = 1.0;
A.figNum = 101;
A.verbose = true;
animate(t,z,A)

% Plot the results:
figure(1337); clf; plotAcrobot(t,z,u,p)

author: MatthewPeterKelly

repo: dscTutorials


text:
function [trg] = read_eep_trg(fn);

% READ_EEP_TRG reads triggers from an EEProbe *.trg file
%
% This function returns an Nx1 array with the N triggers
%
% [trg] = read_eep_trg(filename)
%
% trg(i).time   ... trigger latency in ms
% trg(i).offset ... byte offset
% trg(i).code   ... trigger code (string)
% trg(i).type   ... numeric value of trg.code
%
% where i as number between 1 and N (the number of triggers found)
%
% An EEProbe trigger file is formatted like
%   0.00195312 256
%   0.000     10350  __
%   17.033   2242926   1
%   20.535   2701934   5
%   21.096   2775406  13
%   21.098   2775662   8
%   ...
% where the first column is the trigger latency in seconds, the second
% column is the byte offset in the file and the third column is the triggercode.
% The triggercode can be numeric or a string. The first line of the file contains the
% sample duration.
%
% Author: Robert Oostenveld, Aalborg University, Denmark, 11 March 2003
%
% See also READ_EEP_CNT, READ_EEP_REJ, READ_EEP_AVR
%

% Copyright (C) 2002, Robert Oostenveld
%                     Aalborg University, Denmark
%                     http://www.smi.auc.dk/~roberto/
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

% $Log: read_eep_trg.m,v $
% Revision 1.2  2005/06/08 08:16:37  mvelde
% converted files to unix format
%
% Revision 1.1  2004/11/26 13:17:02  jwiskerke
% Added m-files without binary code in maple distribution.
%
% Revision 1.2  2003/10/24 13:34:41  Maarten-Jan Hoeve
% Added GNU Licence and updated revision history
%
% Revision 1.1.1.2  2003/10/17 09:55:20  mvelde
% updated: consistent copyrights, arguments/data labels, fixed some typos
%
% Revision 1.1.1.1  2003/03/11 15:24:51  roberto
% updated help and copyrights
% ANT Software BV, The Netherlands, www.ant-software.nl / info@ant-software.nl
%

trg = [];

fid = fopen(fn, 'rb');
if fid<0
   return
end

header = fgetl(fid);
while ~feof(fid)
  tmp = fscanf(fid, '%f %d %s', 3);
  if ~isempty(tmp)
    new.time   = tmp(1)*1000;			% in ms
    new.offset = tmp(2)+1;			% offset 1
    new.code   = char(tmp(3:end));		% string
    new.type   = str2double(new.code);		% numeric
    trg = [trg; new];
  end
end

fclose(fid);

author: PatternRecognition

repo: OpenBMI


text:
function pval_ttest = kbttest(X,y)

y = y(:)';

n0 = sum(~isnan(X(y==0,:)),1);
n1 = sum(~isnan(X(y==1,:)),1);
goodidx = find(n0 > 1 & n1 > 1); % AL: n0>=1 & n1>=1?

m0 = nanmean(X(y==0,:),1);
v0 = nanvar(X(y==0,:),1,1); % AL: nanvar(~,0,1)?
m1 = nanmean(X(y==1,:),1);
v1 = nanvar(X(y==1,:),1,1); % AL: nanvar(~,0,1)?
df = (v0./n0 + v1./n1).^2 ./ ...
  max(eps,( (v0./n0).^2./(n0-1) + (v1./n1).^2./(n1-1) ));

tstat = abs((m1-m0)./ max(eps,sqrt(v1./n1 + v0./n0)));
pval_ttest = 2 * tcdf(-abs(tstat),df);
pval_ttest(tstat==0) = 1;

author: kristinbranson

repo: JAABA


text:
% ------------------------------------------------------------------------
%  Copyright (C)
%  Universitat Politecnica de Catalunya BarcelonaTech (UPC) - Spain
%  University of California Berkeley (UCB) - USA
%
%  Jordi Pont-Tuset <jordi.pont@upc.edu>
%  Pablo Arbelaez <arbelaez@berkeley.edu>
%  June 2014
% ------------------------------------------------------------------------
% This file is part of the MCG package presented in:
%    Arbelaez P, Pont-Tuset J, Barron J, Marques F, Malik J,
%    "Multiscale Combinatorial Grouping,"
%    Computer Vision and Pattern Recognition (CVPR) 2014.
% Please consider citing the paper if you use this code.
% ------------------------------------------------------------------------
function params = nyud_params(varargin)
  params = [];
  % Clobber with overrides passed in as arguments
  for i = 1:2:length(varargin)
    key = varargin{i};
    val = varargin{i+1};
    eval(['params.' key ' = val;']);
  end

  params = cv(params, 'gt_set_test', 'test');
  params = cv(params, 'database', 'nyud40Obj');

  % GT sets for training
  params = cv(params, 'gt_set_pareto', 'trainval');
  params = cv(params, 'gt_set_ranking', 'trainval');

  % String identifying the raw hierarchy files
  params = cv(params, 'root_cache_dir', '.');
  params = cv(params, 'raw_hiers', 'hiers');

  % String identifying the current setup (hiers and n_r_cand)
  % to use at different file names
  params = cv(params, 'hiers_id', 'multi_3sc_u_3r');  % multi and 3 scales 'up'

  % String identifying the current Pareto chosen point
  params = cv(params, 'pareto_id', '20k');

  % Maximum number of regions per candidate
  params = cv(params, 'n_r_cand', 3);

  % To identify whether we want to use depth based features or not
  params = cv(params, 'feature_id', 'no_depth');
  params = cv(params, 'depth_features', false);
  params = cv(params, 'camera_matrix', []);

  % Hierarchies to combine
  params.hiers    = {'multi','scale_2.0', 'scale_1.0', 'scale_0.5'};

  % Level of overlap to erase duplicates
  params = cv(params, 'J_th', 0.95);

  % Number of samples per image that the ranking training will get
  params = cv(params, 'n_samples', 30);

  % Max margin parameter
  params = cv(params, 'theta', 0.7);

  %% Do not edit the lines below
  % Where to find the hierarchies
  for ii=1:length(params.hiers)
      params.hier_dirs{ii} = fullfile(params.root_cache_dir, 'output', 'ucm', params.hiers{ii});
  end

  % Cache directory
  params.cacheDir = fullfile(params.root_cache_dir, 'regions', 'cache');
  exists_or_mkdir(params.cacheDir);

  % String identifying the current final global setup at MCG
  params.mcg_id = [params.raw_hiers '_' params.hiers_id '_' params.pareto_id '_' params.feature_id];

  % ----- Results files -----
  % Global Pareto file
  params.files.pareto_full = fullfile(params.root_cache_dir, 'regions', 'pareto', [params.raw_hiers '_' params.hiers_id '_nr_' num2str(params.n_r_cand) '_' params.gt_set_pareto '_pareto_full.mat']);

  % Single-scale results
  for ii=1:length(params.hiers)
      params.files.pareto_singles{ii} = fullfile(params.root_cache_dir, 'regions', 'pareto', [params.raw_hiers '_' params.hiers{ii} '_nr_' num2str(params.n_r_cand) '_' params.gt_set_pareto '_pareto_single.mat']);
  end

  % Parameters of chosen Pareto point
  params.files.output_dir = fullfile(params.root_cache_dir, 'output', 'regions', params.mcg_id);
  exists_or_mkdir(params.files.output_dir);

  params.files.pareto_point = fullfile(params.root_cache_dir, 'regions', 'pareto', [params.mcg_id '_' params.gt_set_pareto '_pareto_point.mat']);
  exists_or_mkdir(fullfile(params.root_cache_dir, 'regions', 'pareto'));

  % File where full features will be stored
  params.files.features_file = fullfile(params.root_cache_dir, 'regions', 'features', [params.mcg_id '_features_' params.gt_set_ranking  '.mat']);
  exists_or_mkdir(fullfile(params.root_cache_dir, 'regions', 'features'));

  % File where the trained classifier will be stored
  params.files.trained_classifier = fullfile(params.root_cache_dir, 'regions', 'classifiers',[params.mcg_id '_rand_forest_' params.gt_set_ranking  '.mat']);
  exists_or_mkdir(fullfile(params.root_cache_dir, 'regions', 'classifiers'));
end

% -------------------------------------------------------------------
% Does nothing if conf.key exists, otherwise sets conf.key to val
function conf = cv(conf, key, val)
  try
    eval(['conf.' key ';']);
  catch
    eval(['conf.' key ' = val;']);
  end
end

author: s-gupta

repo: rcnn-depth

save_path: github-repos/MATLAB/s-gupta-rcnn-depth


text:
function w= rsam2waveform(s)
    w = [];
    for c=1:numel(s)
        if numel(s(c).data) > 0
            wc = waveform(s(c).ChannelTag, 1.0/s(c).sampling_interval, s(c).snum, s(c).data, s(c).units);
        else
            wc = waveform(s(c).ChannelTag, 1.0/s(c).sampling_interval, 0, [], s(c).units);
        end
        w = [w wc];
    end
end

author: geoscience-community-codes

repo: GISMO


text:
function r = isintval(a)
%ISINTVAL     Returns 1 if  a  is intval Taylor
%
%   r = isintval(a)

% written  05/21/09     S.M. Rump
%

  r = isa(a.t,'intval');

repo: OpenMAS


text:
% navego_example_canada_data: post-processing integration of
% IMU, GNSS and visual data for the Canada dataset.
%
% The main goal is to integrate IMU and GNSS measurements from Ekinox-D
% sensor which includes both IMU and GNSS sensors.
%
% Sensors dataset was generated driving a car through the streets of
% Turin city (Italy).
%
%   Copyright (C) 2014, Rodrigo Gonzalez, all rights reserved.
%
%   This file is part of NaveGo, an open-source MATLAB toolbox for
%   simulation of integrated navigation systems.
%
%   NaveGo is free software: you can redistribute it and/or modify
%   it under the terms of the GNU Lesser General Public License (LGPL)
%   version 3 as published by the Free Software Foundation.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU Lesser General Public License for more details.
%
%   You should have received a copy of the GNU Lesser General Public
%   License along with this program. If not, see
%   <http://www.gnu.org/licenses/>.
%
% References:
%
%  The Canadian Planetary Emulation Terrain Energy-Aware Rover Navigation
% Dataset, Run 3. https://starslab.ca/enav-planetary-dataset/
%
%  Johann Diep et al. (2022). Investigating the Performance of LCNS with
% Visual-Inertial Odometry for Lunar Rover Navigation. NAVITEC 2022, April
% 2022.
%
%  Johann Diep (2022). Investigating the Performance of LCNS with
% Visual-Inertial Odometry for Lunar Rover Navigation. Youtube video.
% https://youtu.be/sPnPMBs1gSU
%
% Version: 002
% Date:    2022/07/19
% Author:  Johann Diep <johanndiep@gmail.com>
% URL:     https://github.com/rodralez/navego
%
% NOTE: NaveGo assumes that IMU is aligned with respect to body-frame as
% X-forward, Y-right and Z-down.
%
% NOTE: NaveGo assumes that yaw angle (heading) is positive clockwise.

clc
close all
clear
matlabrc

addpath ../../ins/
addpath ../../ins-gnss/
addpath ../../ins-visual/
addpath ../../conversions/
addpath ../../performance-analysis/
addpath ../../plot/
addpath ../../misc/
addpath ../../simulation/
addpath ./canada-planetary-data/data/
addpath ./canada-planetary-data/scripts/

navego_print_version;

fprintf('\nNaveGo: starting Canada planetary data integration... \n')

%% PARAMETERS

% Comment any of the following parameters in order to NOT execute a
% particular portion of code

GEN_DATA     = 'ON';

if (~exist('GEN_DATA','var')), GEN_DATA = 'OFF'; end

% FusionCase = "inertial_gnss";
% FusionCase = "inertial_visual";
FusionCase = "inertial_visual_gnss";

% Method to remove GNSS measurements for a certain time window. This is
% done to showcase that visual measurements provide a failsafe option in
% case satellite data is not available.
Sparse = "true";
% Sparse = "false";

fprintf('NaveGo: parameter FusionCase = %s \n', FusionCase)
fprintf('NaveGo: parameter Sparse = %s \n', Sparse)

%% Generating Data

if strcmp(GEN_DATA, 'ON')

    if FusionCase == "inertial_gnss"
        fprintf('NaveGo: generating IMU data... \n')
        imu_structure;

        fprintf('NaveGo: generating GNSS data... \n')
        if Sparse == "true"
            gnss_sparse_structure;
        else
            gnss_structure;
        end
    else
        imu_structure;
        fprintf('NaveGo: generating IMU data... \n')
        visual_structure;
        fprintf('NaveGo: generating VISUAL data... \n')
        if Sparse == "true"
            fprintf('NaveGo: generating GNSS data... \n')
            gnss_sparse_structure;
        else
            gnss_structure;
        end
    end
else
    fprintf('NaveGo: loading data... \n')

    if FusionCase == "inertial_gnss"
        load imu_planetary
        if Sparse == "true"
            load gnss_planetary_r
            load gnss_planetary_sparse_r
            load gnss_planetary
        else
            load gnss_planetary_r
            load gnss_planetary
        end
    else
        load imu_planetary
        load visual_planetary
        if Sparse == "true"
            load gnss_planetary_r
            load gnss_planetary_sparse_r
            load gnss_planetary
        else
            load gnss_planetary_r
            load gnss_planetary
        end
    end
end
%% Estimation

switch FusionCase
    case "inertial_gnss"
        fprintf('NaveGo: processing INS/GNSS integration... \n')
        nav_e = ins_gnss(imu_planetary,gnss_planetary,'dcm');
    case "inertial_visual"
        fprintf('NaveGo: processing INS/VISUAL integration... \n')
        nav_e = ins_visual(imu_planetary,gnss_planetary_r,visual_planetary,'dcm'); % note: figure out why we input gnss_planetary_r
    case "inertial_visual_gnss"
        fprintf('NaveGo: processing INS/GNSS/VISUAL integration... \n')
        nav_e = ins_visual_gnss(imu_planetary,gnss_planetary,visual_planetary,'dcm');
end

[nav_i, nav_ref] = navego_interpolation (nav_e, gnss_planetary_r);

if Sparse == "true"
    [gnss_i,gnss_ref] = navego_interpolation(gnss_planetary,gnss_planetary_sparse_r);
else
    [gnss_i,gnss_ref] = navego_interpolation(gnss_planetary,gnss_planetary_r);
end

%% Plotting

fprintf('NaveGo: plotting... \n')

switch FusionCase
    %% Plotting: IMU + GNSS
    case "inertial_gnss"

        % Position
        figure();
        hold on;
        plot(rad2deg(gnss_planetary.lon),rad2deg(gnss_planetary.lat),'.', 'Color', ones(1,3) * 0.75, 'LineWidth', 1.5);
        scatter(rad2deg(gnss_planetary_r.lon),rad2deg(gnss_planetary_r.lat),'y.','MarkerEdgeAlpha',.6);
        plot(rad2deg(nav_e.lon),rad2deg(nav_e.lat),'Color', [0, 0, 0], 'LineWidth', 1.5);
        grid on;
        xlabel('Longitude [deg]');
        ylabel('Latitude [deg]');
        legend('degraded GNSS','RTK','IMU + degraded GNSS','Location','Southeast');
        axis equal;

        % Position Errors
        [RN,RE]  = radius(nav_i.lat);
        LAT2M = RN + nav_i.h;
        LON2M = (RE + nav_i.h).*cos(nav_i.lat);
        [RN,RE]  = radius(gnss_i.lat);
        LAT2M_GR = RN + gnss_i.h;
        LON2M_GR = (RE + gnss_i.h).*cos(gnss_i.lat);

        figure();
        subplot(2,1,1);
        hold on;
        plot(gnss_i.t,LAT2M_GR.*(gnss_i.lat - gnss_ref.lat), '.', 'Color', ones(1,3) * 0.75, 'LineWidth', 1.5)
        plot(nav_i.t, LAT2M.*(nav_i.lat - nav_ref.lat),'Color', [0, 0, 0], 'LineWidth', 1.5)
        grid on;
        xlabel('Time [s]')
        ylabel('[m]')
        legend('GNSS', 'IMU + degraded GNSS', 'Location', 'northoutside');
        title('Latitude Error');
        xlim([0,max(gnss_planetary.t)]);
        subplot(2,1,2);
        hold on;
        plot(gnss_i.t, LON2M_GR.*(gnss_i.lon - gnss_ref.lon), '.', 'Color', ones(1,3) * 0.75, 'LineWidth', 1.5)
        plot(nav_i.t, LON2M.*(nav_i.lon - nav_ref.lon),'Color', [0, 0, 0], 'LineWidth', 1.5)
        grid on;
        xlabel('Time [s]')
        ylabel('[m]')
        legend('GNSS', 'IMU + degraded GNSS', 'Location', 'northoutside');
        title('Longitude Error');
        xlim([0,max(gnss_planetary.t)]);

        %% Plotting: Vision + IMU
    case "inertial_visual"

        % Position
        figure();
        hold on;
        plot(rad2deg(visual_planetary.lon),rad2deg(visual_planetary.lat),'.', 'Color', ones(1,3) * 0.75, 'LineWidth', 1.5);
        scatter(rad2deg(gnss_planetary_r.lon),rad2deg(gnss_planetary_r.lat),'y.','MarkerEdgeAlpha',.6);
        plot(rad2deg(nav_e.lon),rad2deg(nav_e.lat),'Color', [0, 0, 0], 'LineWidth', 1.5);
        grid on;
        xlabel('Longitude');
        ylabel('Latitude');
        legend('OpenVINS','RTK','IMU + OpenVINS','Location','Southeast');
        axis equal;
        % Position Errors
        [RN,RE]  = radius(nav_i.lat);
        LAT2M = RN + nav_i.h;
        LON2M = (RE + nav_i.h).*cos(nav_i.lat);

        [RN,RE]  = radius(gnss_i.lat);
        LAT2M_GR = RN + gnss_i.h;
        LON2M_GR = (RE + gnss_i.h).*cos(gnss_i.lat);

        figure();
        subplot(2,1,1);
        hold on;
        plot(gnss_i.t,  LAT2M_GR.*(gnss_i.lat - gnss_ref.lat), '.', 'Color', ones(1,3) * 0.75, 'LineWidth', 1.5)
        plot(nav_i.t, LAT2M.*(nav_i.lat - nav_ref.lat),'Color', [0, 0, 0], 'LineWidth', 1.5)
        grid on;
        xlabel('Time [s]')
        ylabel('[m]')
        legend('GNSS', 'IMU + OpenVINS', 'Location', 'northoutside');
        title('Latitude Error');
        xlim([0,max(gnss_planetary.t)]);
        subplot(2,1,2);
        hold on;
        plot(gnss_i.t, LON2M_GR.*(gnss_i.lon - gnss_ref.lon), '.', 'Color', ones(1,3) * 0.75, 'LineWidth', 1.5)
        plot(nav_i.t, LON2M.*(nav_i.lon - nav_ref.lon),'Color', [0, 0, 0], 'LineWidth', 1.5)
        grid on;
        xlabel('Time [s]')
        ylabel('[m]')
        legend('GNSS', 'IMU + OpenVINS', 'Location', 'northoutside');
        title('Longitude Error');
        xlim([0,max(gnss_planetary.t)]);

        %% Plotting: Vision + GNSS + INS
    case "inertial_visual_gnss"

        % Position
        figure();
        hold on;
        plot(rad2deg(gnss_planetary.lon),rad2deg(gnss_planetary.lat),'.', 'Color', ones(1,3) * 0.75, 'LineWidth', 1.5);
        scatter(rad2deg(gnss_planetary_r.lon),rad2deg(gnss_planetary_r.lat),'y.','MarkerEdgeAlpha',.6);
        plot(rad2deg(nav_e.lon),rad2deg(nav_e.lat),'Color', [0, 0, 0], 'LineWidth', 1.5);
        grid on;
        xlabel('Longitude [deg]');
        ylabel('Latitude [deg]');
        legend('degraded GNSS','RTK','IMU + degraded GNSS + OpenVINS','Location','Southeast');
        axis equal;

        % Position Errors
        [RN,RE]  = radius(nav_i.lat);
        LAT2M = RN + nav_i.h;
        LON2M = (RE + nav_i.h).*cos(nav_i.lat);

        [RN,RE]  = radius(gnss_i.lat);
        LAT2M_GR = RN + gnss_i.h;
        LON2M_GR = (RE + gnss_i.h).*cos(gnss_i.lat);

        figure();
        subplot(2,1,1);
        hold on;
        plot(gnss_i.t,  LAT2M_GR.*(gnss_i.lat - gnss_ref.lat), '.', 'Color', ones(1,3) * 0.75, 'LineWidth', 1.5)
        plot(nav_i.t, LAT2M.*(nav_i.lat - nav_ref.lat),'Color', [0, 0, 0], 'LineWidth', 1.5)
        grid on;
        xlabel('Time [s]')
        ylabel('[m]')
        legend('GNSS', 'IMU + degraded GNSS + OpenVINS', 'Location', 'northoutside');
        title('Latitude Error');
        xlim([0,max(gnss_planetary.t)]);
        subplot(2,1,2);
        hold on;
        plot(gnss_i.t, LON2M_GR.*(gnss_i.lon - gnss_ref.lon), '.', 'Color', ones(1,3) * 0.75, 'LineWidth', 1.5)
        plot(nav_i.t, LON2M.*(nav_i.lon - nav_ref.lon),'Color', [0, 0, 0], 'LineWidth', 1.5)
        grid on;
        xlabel('Time [s]')
        ylabel('[m]')
        legend('GNSS', 'IMU + degraded GNSS + OpenVINS', 'Location', 'northoutside');
        title('Longitude Error');
        xlim([0,max(gnss_planetary.t)]);
end

author: rodralez

repo: NaveGo

save_path: github-repos/MATLAB/rodralez-NaveGo


text:
function sparse_count_test025 ( dim_min, dim_max, level_max_min, ...
  level_max_max )

%*****************************************************************************80
%
%% SPARSE_COUNT_TEST025 tests CFN_E_SIZE_TOTAL.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 May 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_MIN, the minimum spatial dimension to consider.
%
%    Input, integer DIM_MAX, the maximum spatial dimension to consider.
%
%    Input, integer LEVEL_MAX_MIN, the minimum value of LEVEL_MAX to consider.
%
%    Input, integer LEVEL_MAX_MAX, the maximum value of LEVEL_MAX to consider.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST025\n' );
  fprintf ( 1, '  CFN_E_SIZE_TOTAL returns the number of\n' );
  fprintf ( 1, '  points in a CFN_E sparse grid made from \n' );
  fprintf ( 1, '  any closed fully nested family of 1D quadrature\n' );
  fprintf ( 1, '  rules with exponential growth, including:\n' );
  fprintf ( 1, '  * CC_E, the Clenshaw Curtis Exponential Growth family;\n' );
  fprintf ( 1, '  * NCC_E, the Newton Cotes Closed Exponential Growth family.\n' );
  fprintf ( 1, '  No reduction is made because of repeated points.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '   DIM: ' );

  for dim_num = dim_min : dim_max
    fprintf ( 1, '  %10d', dim_num );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '   LEVEL_MAX\n' );
  fprintf ( 1, '\n' );

  for level_max = level_max_min : level_max_max
    fprintf ( 1, '    %4d', level_max );
    for dim_num = dim_min : dim_max
      point_num = cfn_e_size_total ( dim_num, level_max );
      fprintf ( 1, '  %10d', point_num );
    end
    fprintf ( 1, '\n' );
  end

  return
end


text:
function circle_inout ( )

%*****************************************************************************80
%
%% CIRCLE_INOUT plots data in and out of a circle.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 April 2011
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CIRCLE_INOUT:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Make a scatterplot of two sets of data representing\n' );
  fprintf ( 1, '  points in the unit square that are also in or not in\n' );
  fprintf ( 1, '  the unit circle.\n' );
%
%  Read the data.
%
  xy_in = load ( 'circle_in.txt' );
  [ n_in, dim ] = size ( xy_in );

  xy_out = load ( 'circle_out.txt' );
  [ n_out, dim ] = size ( xy_out );

  n = n_in + n_out;
%
%  Report on the estimate for pi:
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of points inside the circle = %d\n', n_in );
  fprintf ( 1, '  Number outside                     = %d\n', n_out );
  fprintf ( 1, '  Total                              = %d\n', n );
  fprintf ( 1, '  Estimate for PI                    = %d\n', 4 * n_in / n );
%
%  Set up points on a circle
%
  n = 50;
  center = [ 0.0, 0.0 ];
  radius = 1.0;
  theta1 = 0.0;
  theta2 = 90.0;

  xy_circ = circle_arc ( n, center, radius, theta1, theta2 );
%
%  Plot the data.
%
  plot ( xy_in(:,1),  xy_in(:,2),  'b.', ...
         xy_out(:,1), xy_out(:,2), 'r.', ...
         xy_circ(:,1), xy_circ(:,2), 'r-', 'LineWidth', 2 );
  grid on
  axis ( [ 0.0, 1.0, 0.0, 1.0 ] )
  axis equal
  axis square
  xlabel ( '<--- X --->' );
  ylabel ( '<--- Y --->' );
  title ( 'Random points inside/outside the unit circle' );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'CIRCLE_INOUT:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  return
end
function xy_circ = circle_arc ( n, center, radius, theta1, theta2 )

%*****************************************************************************80
%
%% CIRCLE_ARC samples points on a circular arc.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 April 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of sample points.
%
%    Input, real CENTER(2), the center of the circle.
%
%    Input, real RADIUS, the radius of the circle.
%
%    Input, real THETA1, THETA2, the angular coordinates of the first and
%    last points on the arc, in degrees.
%
%    Output, real XY_CIRC(N,2), points along the arc.
%
  t = pi * linspace ( theta1, theta2, n ) / 180;
  t = t';

  xy_circ = zeros ( n, 2 );

  xy_circ(:,1) = center(1) + radius * cos ( t(:,1) );
  xy_circ(:,2) = center(2) + radius * sin ( t(:,1) );

  return
end


text:
%
% This is a calling routine to test & check out the power spectrum &
% spectrogram routines for unequal segment lengths. In addition, use it
% to compare with Chronux routines when segments are of equal length.
%
clear all;

if 0
    dir = 'G:\ravi\Chrowser\Pass~ Tioga_0e927741-9673-46e5-9050-ca1d7541bf22\';
    xfile = 'Pass~ Tioga_0e927741-9673-46e5-9050-ca1d7541bf22'
    %dir = 'G:\ravi\Chrowser\sample~ data_8ef647e3-e5ea-43a6-8c69-fb848b8db7c2\';
    %xfile = 'sample~ data_8ef647e3-e5ea-43a6-8c69-fb848b8db7c2'
else
    dir = 'Z:\xltekRawData\Wallis~ Terry_c3f44891-afa7-4fa7-a643-55c772a05241\'
    xfile = 'Wallis~ Terry_c3f44891-afa7-4fa7-a643-55c772a05241'
end

% Get header info
% Channels are labelled from C1 through C127 and ''
% total of 128 channels
hdr = eegMex( dir, xfile);

gram = 1 ; % 0=spectra, 1=coherence
chronux = 0 ; % 0=no comparison with Chronux; 1=compare with chronux

%nSamples = 4210;
nChannels = 2;
nSegments = 1;
movingwin = [25, 25];

%
% Spectral Parameters
%
params.fpass = [ 0 0.5 ];
params.pad = 2;
params.err = [2 0.05];  % err(1)=0 is no err estimates; err(1)=1 is asymptotic estimates; err(1)=2 is Jacknife
params.trialave = 1;
params.Fs = 1;

%
% Tapers issues
%
halfBandWidth = 2.5;
kCap = 2*halfBandWidth - 1;
%params.tapers = [ halfBandWidth, kCap ];
params.tapers = [ halfBandWidth, 2 ];

%
% Basic checks on inputs
%
if (gram==1) && (nChannels < 2), error( 'Coherence requires at least 2 channels...' ); end
%if (nSegments==1) && (params.err(1)==2), error( 'Jacknife requires more than 1 segment'); end

%
% Generate segments endpoints randomly
% myrandint is a 3rd party routine (from matlab site)
%
% Randomly generated segment end points
sMarkers = reshape( sort( myrandint( 2*nSegments, 1, [ ceil(hdr.nSamples/500) : ceil(hdr.nSamples/50) ], 'noreplace' ) )', 2, nSegments )';
%sMarkers = [ ceil(hdr.nSamples/80), ceil(hdr.nSamples/65) ];

%
% Randomly select a few channels
%
if ~chronux
    %chIndices = sort( myrandint( nChannels, 1, [ round(hdr.nChans/4) : round(3*hdr.nChans/4) ], 'noreplace' ) );
    chIndices = [ 3 : 3+nChannels-1 ];
else
    %chIndices = [ 10 : 10+nChannels-1 ];
    chIndices = [ 3, 7 ];
end

%
% Randomly generate the time series
%
fulldata = eegMex( dir, xfile, chIndices, [ 1 hdr.nSamples/50 1 ] );
mDiscardBits = 0;
conversionFactor = ( 8711 / (2^21 - 0.5) ) * 2^mDiscardBits;
fulldata{:} = fulldata{:} * conversionFactor;

%
% Create a data matrix with all the segments aligned one after another
%
totalSegmentLength = sum( sMarkers(:,2) - sMarkers(:,1) + 1 );
data = zeros( totalSegmentLength, length(chIndices) ); % preallocate to ensure contiguous memory
newMarkers(1,1) = 1;
newMarkers(1,2) = sMarkers(1,2) - sMarkers(1,1) + 1;
data( newMarkers(1,1):newMarkers(1,2), : ) = detrend( fulldata{1}( sMarkers(1,1):sMarkers(1,2), :) );
for sg = 2:size( sMarkers, 1 )
    newMarkers(sg,1) = newMarkers(sg-1,2) + 1;
    newMarkers(sg,2) = newMarkers(sg,1) + sMarkers(sg,2) - sMarkers(sg,1);
    data( newMarkers(sg,1):newMarkers(sg,2), : ) = detrend( fulldata{1}( sMarkers(sg,1):sMarkers(sg,2), :) );
end

% To ensure that we check results from array indices beyond 1
if nChannels > 1
    ix = sort( myrandint( 1, 2, [1:length(chIndices)], 'noreplace' ) ); % Arbitrarily pick two indices from selected channels for testing results
    i1=ix(1); i2=ix(2);
    % iC = m + (n-1)*(n-2)/2, for elements of the the coherence matrix, Cmn
    iC = ix(1) + (ix(2)-1)*(ix(2)-2)/2;
else
    ix = sort( myrandint( 1, 1, [1:length(chIndices)], 'noreplace' ) ); % Arbitrarily pick 1 indices from selected channels for testing results
    i1=ix(1);
end

%
% Power spectrum/spectrogram/coherence/coherogram
%
if gram==0
    [ S, f, Serr ] = avgSpectrum( data, movingwin, params, newMarkers );
    figure; plot( f, 10*log10( S(:,i1) ), 'k', f, 10*log10( Serr(2,:,i1) ), 'g--', f, 10*log10( Serr(1,:,i1)), 'g--' ); title('Avg. Routine:: Spectrum');
    %figure; plot( f, 10*log10( S(:,i1) )); title('Avg. Routine:: Spectrum');
elseif gram==1
    [Cmn,Phimn,Smn,Smm,f,ConfC,PhiStd,Cerr] = avgCoherence( data, movingwin, params, newMarkers );
    %  C(i,j) = Cmn(:,k) where k = j + (1/2)*(i-1)*(i-2)
    figure; plot( f, Cmn(:,iC), 'k', f, Cerr(2,:,iC), 'g--', f, Cerr(1,:,iC), 'g--' );
    title('Avg. Routine:: Coherence'); ylim([0 1])
    %figure; plot( f, 10*log10( Cmn(:,iC) ) ); title('Avg. Routine:: Coherence-Magnitude');
    %figure; plot( f, phimn(:,iC) ); title('Avg. Routine:: Coherence-Phase');
    disp( ['Confidence level for C (confC) at (1-p) level: ' num2str( ConfC(iC)) ] );
end

%
% Use to check against Chronux: only for equal length segments
%
if chronux

    win = floor( newMarkers(1,2) / movingwin(1) );
    newMarkers(1,2) = newMarkers(1,2) - mod( newMarkers(1,2), win );
    cdata = data( [1:newMarkers(1,2)], i1 );
    cdata = detrend( reshape( cdata, [ newMarkers(1,2)/win, win ] ) );
    cdata2 = data( [1:newMarkers(1,2)], i2 );
    cdata2 = detrend( reshape( cdata2, [ newMarkers(1,2)/win, win ] ) );
    params.trialave = 1;
    if gram==0
        [ cS, cf, cSerr ] = mtspectrumc( cdata, params );
        figure; plotvector( cS, cf, 'l', cSerr );
        %figure; plot( cf, 10*log10( cS )); title('Chronux:: Spectrum');
        figure; plot( cf, 10*log10(cSerr(1,:)), cf, 10*log10(cSerr(2,:)) ); title('Chronux Error-Bar Computations');
        figure; plot( cf, 10*log10( cS ) - 10*log10( S(:,i1) )); title('Error in Spectrum = |New Routines - Chronux|');
        figure; plot( cf, 10*log10(cSerr(1,:)) - 10*log10(Serr(1,:,i1)), cf, 10*log10(cSerr(2,:)) - 10*log10(Serr(2,:,i1)) );title('Error in Error-Bar Computations = |New Routines - Chronux| ');
    elseif gram==1

        [cC,cphi,cS12,cS1,cS2,cf,cconfC,cphistd,cCerr]=coherencyc( cdata, cdata2, params );

        %figure; plotvector( cC(:,1), cf, 'n', cCerr );
        figure; plot( cf, cC(:,iC), 'k', cf, cCerr(2,:,iC), 'g--', cf, cCerr(1,:,iC), 'g--' );
        title('Chronux:: Coherence'); ylim([0 1])
        %figure; plot( cf, 10*log10( cC(:,1) ) ); title('Chronux:: Coherence-Magnitude');
        figure; plot( cf, 10*log10( cC(:,1) ) - 10*log10( Cmn(:,iC) ) ); title('Error in Coherence = |New Routines - Chronux|');
        % Phase may give a problem of 2pi difference... look into it.
        figure; plot( cf, cphi(:,1) -  Phimn(:,iC) ); title('Error in Phase = |New Routines - Chronux|');
        %
        % Note the remaining quantities do not really need to checked since
        % coherence = cross-spectrum/power spectra* power spectra, ie C = S12/(S1*S2)
        % so unlikely that S12, S1, S2 are incorrect if C is correct.
        if 1
            figure; plot( cf, 10*log10( cS1(:,1) ) - 10*log10( Smm(:,ix(1)) ) ); title('Error in Power Spectrogram-1 = |New Routines - Chronux|');
            figure; plot( cf, 10*log10( cS2(:,1) ) - 10*log10( Smm(:,ix(2)) ) ); title('Error in Power Spectrogram-2 = |New Routines - Chronux|');
        end
        %
        % Error-Bars & Confidence Levels
        disp( ['Confidence levelfor C (confC) at (1-p) level: ' num2str( cconfC) ' (Chronux)' ] );
        disp( ['Error in confidence level, confC: ' num2str( ConfC(iC) - cconfC ) ] );
        %figure; plot( cf, cphistd(:,1), f, phistd(:,iC) ); title('Phase-Error-Bar Computations');
        figure; plot( cf, cphistd(:,1) -  PhiStd(:,iC) ); title('Error in PhiStd-1');
        figure; plot( cf, cphistd(:,1) -  PhiStd(:,iC) ); title('Error in PhiStd-2');
        figure; plot( cf, abs(cCerr(1,:,1) -  Cerr(1,:,iC)), cf, abs(cCerr(2,:,1) -  Cerr(2,:,iC))  ); title('Error in Abs(Coherence)-Error-Bar Computations = |New Routines - Chronux|');
    end
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [ x, w ] = rule05 ( n )

%*****************************************************************************80
%
%% RULE05 returns the rule of degree 5.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 July 2014
%
%  Author:
%
%    Original FORTRAN77 version by Hong Xiao, Zydrunas Gimbutas.
%    This MATLAB version by John Burkardt.
%
%  Reference:
%
%    Hong Xiao, Zydrunas Gimbutas,
%    A numerical algorithm for the construction of efficient quadrature
%    rules in two and higher dimensions,
%    Computers and Mathematics with Applications,
%    Volume 59, 2010, pages 663-676.
%
%  Parameters:
%
%    Input, integer N, the number of nodes.
%
%    Output, real X(2,N), the coordinates of the nodes.
%
%    Output, real W(N), the weights.
%
  xs = [ ...
    0.1775868202077551E-01,-.1775868202077539E-01, ...
    0.7788710544649639,-.7788710544649639, ...
    -.7703781288541645,0.7703781288541645, ...
    -.7490353914168658D-33 ];
  ys = [ ...
    -.9659285494001192,0.9659285494001192, ...
    -.5715708301251639,0.5715708301251639, ...
    -.5829672991828014,0.5829672991828014, ...
    0.1356144833394667D-33 ];
  ws = [ ...
    0.2246199725165690,0.2246199725165690, ...
    0.3901817339168917,0.3901817339168917, ...
    0.3953508381187504,0.3953508381187504, ...
    0.8081220356417684 ];

  x(1,1:n) = xs(1:n);
  x(2,1:n) = ys(1:n);
  w(1:n) = ws(1:n);

  return
end


text:
function r8mat_transpose_print_some ( m, n, a, ilo, jlo, ihi, jhi, title )

%*****************************************************************************80
%
%% R8MAT_TRANSPOSE_PRINT_SOME prints some of an R8MAT, transposed.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 May 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns.
%
%    Input, real A(M,N), an M by N matrix to be printed.
%
%    Input, integer ILO, JLO, the first row and column to print.
%
%    Input, integer IHI, JHI, the last row and column to print.
%
%    Input, string TITLE, an optional title.
%
  incx = 5;

  if ( 0 < s_len_trim ( title ) )
    fprintf ( 1, '\n' );
    fprintf ( 1, '%s\n', title );
  end

  for i2lo = max ( ilo, 1 ) : incx : min ( ihi, m )

    i2hi = i2lo + incx - 1;
    i2hi = min ( i2hi, m );
    i2hi = min ( i2hi, ihi );

    inc = i2hi + 1 - i2lo;

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Row: ' );
    for i = i2lo : i2hi
      fprintf ( 1, '%7d       ', i );
    end
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Col\n' );

    j2lo = max ( jlo, 1 );
    j2hi = min ( jhi, n );

    for j = j2lo : j2hi

      fprintf ( 1, '%5d ', j );
      for i2 = 1 : inc
        i = i2lo - 1 + i2;
        fprintf ( 1, '%12f', a(i,j) );
      end
      fprintf ( 1, '\n' );

    end

  end

  return
end


text:
classdef StudentTD
%%STUDENTD Functions to handle the multivariate Student-t distribution.
%    Note that the multivariate Student-t distribution with one degree of
%    freedom is the same as the multivariate Cauchy distribution.
%    Also the scalar student-t distribution is the distribution of
%    x/sqrt(y/n) where x is a normal random variable, y is a central
%    chi squared random variable and n is the number of degrees of freedom
%    of y. In terms of tracking, the distribution can be of interest,
%    because, as noted in [1], angular measurements in the presence of
%    glint have been fit to Cauchy distributions in some instances.
%Implemented methods are: mean, cov, PDF,
%   Methods only for scalar distributions: CDF, invCDF
%   Methods only for scalar, central, non-scaled distributions: rand,
%                                                               entropy
%
%REFERENCES:
%[1] U. Nickel, "Angular superresolution with phased array radar: A review
%    of algorithms and operational constraints," IEEE Proceedings, vol.
%    134, Part F, no. 1, pp. 53-59, Feb. 1987.
%
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

methods(Static)

    function val=mean(mu,nu)
    %%MEAN  Obtain the mean of the multivariate Student-t distribution for
    %       given location and scale parameters and degrees of freedom.
    %
    %INPUTS: mu The DX1 location vector of the student's-t distribution.
    %        nu The scalar number of degrees of freedom of the Student-t
    %           distribution. nu>=0.
    %
    %OUTPUTS: val The mean of the multivariate Student-t distribution.
    %
    %The mean is undefined if nu<=1.
    %
    %October 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.

        if(nu>1)
            val=mu;
        else
            val=NaN;
        end
    end

    function val=cov(Sigma,nu)
    %%MEAN  Obtain the covariance matrix of the multivariate Student-t
    %       distribution for given location and scale parameters and degrees
    %       of freedom.
    %
    %INPUTS: Sigma The DXD symmetric, positive definite scale matrix.
    %           nu The scalar number of degrees of freedom of the Student-t
    %              distribution. nu>=0.
    %
    %OUTPUTS: val The covariance matrix of the multivariate Student-t
    %             distribution.
    %
    %The covariance matrix is undefined if nu<=2.
    %
    %October 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.

        if(nu>2)
            val=nu/(nu-2)*Sigma;
        else
            val=NaN;
        end
    end

    function val=PDF(x,mu,Sigma,nu)
    %%PDF Evaluate the monovariate or multivariate Student-t distribution at
    %     the desired point.
    %
    %INPUTS: x The DX1 vector at which the (possibly multivariate) Student's-t
    %          distribution is to be evaluated.
    %       mu The DX1 location vector of the student's-t distribution.
    %    Sigma The DXD symmetric, positive definite scale matrix.
    %       nu The scalar number of degrees of freedom of the Student-t
    %          distribution. nu>=0.
    %
    %OUTPUTS: val The PDF of the Student's-t distribution at x with the given
    %             parameters.
    %
    %The vector version of the Student-t Distribution is given in Appendix B of
    %[1]. As nu->Inf, the distribution reduces to a multivariate Gaussian
    %distribution with mean mu and covariance matrix Sigma.
    %
    %Logarithms are used in the implementation to reduce the effect of
    %precision problems that can arise if nu is large. The problems arise due
    %to the ratio of gamma functions.
    %
    %REFERENCES:
    %[1] C. M. Bishop, Pattern Recognition and Machine Learning. Cambridge,
    %    United Kingdom: Springer, 2007.
    %
    %October 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.

        D=size(x,1);

        diff=x-mu;
        Delta2=diff'*pinv(Sigma)*diff;

        num=gammaln((nu+D)/2);
        denom=gammaln(nu/2)+log(sqrt(det(nu*pi*Sigma)))+((D+nu)/2)*log(1+Delta2/nu);

        val=num-denom;
        val=exp(val);
    end

    function val=CDF(x,nu,mu,sigma)
    %%CDF Evaluate cumulative distribution function (CDF) of a a scalar
    %     Student-t distribution at a specified points given the mean
    %     and the variance, or for a StudentT(nu,0,1) distribution if the
    %     mean and variance are omitted.
    %
    %INPUTS: x: a vector of values at which the CDF will be evaluated
    %       nu: scalar degrees of freedom
    %       mu: scalar mean parameter
    %    sigma: scalar scale parameter such that t = (x-mu)/sigma is a standard
    %           t-distributed variable. Note this is not the standard
    %           deviation.
    %
    %OUTPUTS: val: The scalar value(s) of the CDF with mean mu and
    %              evaluated at the point(s) x.
    %
    %Note: This function relies on the built-in betainc functions which
    %      computes the incomplete beta function.
    %
    %October 2020 Codie T. Lewis, Naval Research Laboratory, Washington D.C.

        if(nargin<3||isempty(mu))
           mu = 0;
        end

        if(nargin<4||isempty(sigma))
           sigma = 1;
        end

        %Standardize inputs
        t = (x-mu)/sigma;

        down = t<0;

        s = nu./(t.^2+nu);
        val = 1-0.5*betainc(s,nu/2,1/2);
        val(down) = 1-val(down);

    end

    function val=invCDF(prob,nu,mu,sigma)
    %%invCDF Evaluate the inverse cumulative distribution function (CDF) of
    %        a scalar Student-t distribution at a specified points given the
    %        mean and the variance, or for a StudentT(nu,0,1) distribution if
    %        the mean and variance are omitted.
    %
    %INPUTS: prob: a vector of values at which the invCDF will be evaluated
    %          nu: scalar degrees of freedom
    %          mu: scalar mean parameter
    %       sigma: scalar scale parameter such that t = (x-mu)/sigma is a
    %              standard t-distributed variable. Note this is not the
    %              standard deviation.
    %
    %OUTPUTS: val: The scalar value(s) of the invCDF with mean mu and
    %              evaluated at the point(s) x.
    %
    %Note: This function relies on the built-in betaincinv function which
    %      computes the incomplete beta function.
    %
    %October 2020 Codie T. Lewis, Naval Research Laboratory, Washington D.C.

        if(nargin<3||isempty(mu))
           mu = 0;
        end

        if(nargin<4||isempty(sigma))
           sigma = 1;
        end

        up = prob>=1/2;
        down = prob<1/2;
        t = zeros(size(prob));

        t(up) = sign(prob(up)-0.5).*sqrt((nu./betaincinv(2*(1-prob(up)),nu/2,1/2))-nu);
        t(down) = sign(prob(down)-0.5).*sqrt((nu./betaincinv(2*(prob(down)),nu/2,1/2))-nu);

        val = mu+t*sigma;
    end

    function vals=rand(N,nu)
    %%RAND Generate scalar Student-t random variables with unit scale factor
    %      and a given number of degrees of freedom.
    %
    %INPUTS: N If N is a scalar, then rand returns an NXN matrix of random
    %          variables. If N=[M,N1] is a two-element row vector, then rand
    %          returns an MXN1 matrix of random variables.
    %       nu The scalar number of degrees of freedom of the Student-t
    %          distribution. nu>=0.
    %
    %OUTPUTS: vals A matrix whose dimensions are determined by N of the
    %              generated scalar Student-t random variables.
    %
    %The algorithm implemented is the TIR algorithm in [1].
    %
    %REFERENCES:
    %[1] A. J. Kinderman, J. F. Monahan, and J. G. Ramage, "Computer methods
    %    for sampling from student's t distribution," Mathematics of
    %    Computation, vol. 31, no. 140, pp. 1009-1018, Oct. 1977.
    %
    %August 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.

    if(isscalar(N))
        dims=[N, N];
    else
        dims=N;
    end

    vals=zeros(dims);
    numVals=numel(vals);

    b=sqrt(2*exp(-1/2)-1);
    alpha=nu;
    for curVal=1:numVals
        while(1)
            %Step 1
            u=rand(1);

            if(u<b/2)
                x=4*u-b;
                %Step 2
                v=rand(1);
                if(v<=1-abs(x)/2)
                   vals(curVal)=x;
                   break;
                end

                uAlpha=(1+x^2/alpha)^(-(alpha+1)/2);

                if(v<=uAlpha)
                    vals(curVal)=x;
                    break;
                end
                continue;
            end

            if(u<0.5)
                %Step 3
                temp=4*u-1-b;
                x=(abs(temp)+b)*sign(temp);
                v=rand(1);

                %Step 4
                if(v<=1-abs(x)/2)
                    vals(curVal)=x;
                    break;
                end

                if(v>=(1+b^2)/(1+x^2))
                    continue;
                end

                uAlpha=(1+x^2/alpha)^(-(alpha+1)/2);
                if(v<=uAlpha)
                    vals(curVal)=x;
                    break;
                end
                continue;
            end

            if(u<0.75)
                %Step 5
                temp=8*u-5;
                x=2/((abs(temp)+1)*sign(temp));
                u1=rand(1);
                v=x^(-2)*u1;

                %Step 4 again
                if(v<=1-abs(x)/2)
                    vals(curVal)=x;
                    break;
                end

                if(v>=(1+b^2)/(1+x^2))
                    continue;
                end

                uAlpha=(1+x^2/alpha)^(-(alpha+1)/2);
                if(v<=uAlpha)
                    vals(curVal)=x;
                    break;
                end
                continue;
            end

            %Step 6
            x=2/(8*u-7);
            v=rand(1);

            uAlpha=(1+x^2/alpha)^(-(alpha+1)/2);
            if(v<x^2*uAlpha)
                vals(curVal)=x;
                break;
            end
        end
    end
    end

    function entropyVal=entropy(nu)
    %%ENTROPY Obtain the differential entropy in nats of the scalar Student-t
    %         distribution with unit scale factor and a given number of degrees
    %         of freedom.  The differential entropy of a continuous
    %         distribution is entropy=-int_x p(x)*log(p(x)) dx where the
    %         integral is over all values of x. Units of nats mean that the
    %         natural logarithm is used in the definition. Unlike the Shannon
    %         entropy for discrete variables, the differential entropy of
    %         continuous variables can be both positive and negative.
    %
    %INPUTS: nu The scalar number of degrees of freedom of the Student-t
    %           distribution. nu>=0.
    %
    %OUTPUTS: entropyVal The value of the differential entropy in nats.
    %
    %Differential entropy is defined in Chapter 8 of [1].
    %
    %REFERENCES:
    %[1] T. M. Cover and J. A. Thomas, Elements of Information Theory, 2nd ed.
    %    Hoboken, NJ: Wiley-Interscience, 2006.
    %
    %April 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.

        entropyVal=((nu+1)/2)*(psi((nu+1)/2)-psi(nu/2))+(1/2)*log(nu)+betaln(nu/2,1/2);
    end
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


function C=vec(M);

C=M(:);


text:
function acq_checkTriggers(bbci)
%ACQ_CHECKTRIGGERS - Check whether triggers are operational
%
%Synopsis:
% acq_checkTriggers(<BBCI>)
%
%If no BBCI structure is specified, the default bbci.source is assumed
%which is (currently) BrainVision.

misc_checkTypeIfExists('bbci', 'STRUCT');

bbci= bbci_apply_setDefaults(bbci);
bbci.source.record_signals= 0;
[data, bbci]= bbci_apply_initData(bbci);
pause(0.5);

inp= 2.^[0:7];
trigger_ok= 1;
for k= 1:length(inp),
  fprintf('sending trigger: %3d -> received: ', inp(k));
  outp{k}= acq_checkTrigger(bbci, data, inp(k), 500);
  if isempty(outp{k}),
    fprintf('nothing.\n');
  else
    fprintf('%3d.\n', outp{k});
  end
  if ~isequal(inp(k), outp{k}),
    trigger_ok= 0;
  end
end

bbci_apply_close(bbci);
if trigger_ok,
  fprintf('Trigger test successful.\n');
else
  error('!!! Trigger test unsuccessful !!!\n');
end
return

function outp= acq_checkTrigger(bbci, data, inp, timeout)

outp= [];
trigger_received= false;
t0= clock;
bbci_trigger(bbci, inp);
while isempty(outp) && etime(clock, t0)*1000<timeout,
  [source, marker]= ...
        bbci_apply_acquireData(data.source, bbci.source, data.marker);
  if ~source.state.running,
    break;
  end
  if length(marker.desc)>1,
    outp= marker.desc(end);
  end
end

author: bbci


text:
function [fr,th_val,phi_val] = tess_parametrize(vert, th_val, phi_val, p)
% TESS_PARAMETRIZE_NEW: Get a parametric representation of a closed and non-overlapping surface.
%
% USAGE:  tess_parametrize(vert, th_val, phi_val, p)
%         tess_parametrize(vert, th_val, phi_val)
%
% INPUTS:
%    - vert    : x,y,z coordinates of the surface to parametrize
%    - th_val  : Row array of theta values
%    - phi_val : Row array of phi values
%    - p       : Pad the function of p values in both direction with circular values

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2011

% Parse inputs
if (nargin < 4) || isempty(p)
    p = 0;
end

% Convert surface to spherical coordinates
[th,phi,r] = cart2sph(vert(:,1), vert(:,2), vert(:,3));
% Replicate the values to get a full coverage at the edges
th  = [th;  th;        th;        th+2*pi;  th+2*pi;  th+2*pi;   th-2*pi;  th-2*pi;  th-2*pi  ];
phi = [phi; phi+2*pi;  phi-2*pi;  phi;      phi+2*pi; phi-2*pi;  phi;      phi+2*pi; phi-2*pi ];
r   = [r;   r;         r;         r;        r;        r;         r;        r;        r        ];

% Build grid of (theta,phi) values
[fth, fphi] = meshgrid(th_val, phi_val);
% Estimate radius values for all those points
fr = griddata(th, phi, r, fth, fphi);

% Pad function with circular values, to avoid edge irregularities
if (p > 0)
    fr_pad  = bst_pad(fr, [0  p], 'circular');
    fr_pad  = bst_pad(fr_pad, [p  0], 'mirror');
    fr = fr_pad;
    th_val  = [th_val(end-p+1:end)-2*pi,  th_val,  th_val(1:p)+2*pi];
    phi_val = [phi_val(end-p+1:end)-2*pi, phi_val, phi_val(1:p)+2*pi];
end

author: brainstorm-tools


text:
function cvt_test14 ( )

%*****************************************************************************80
%
%% CVT_TEST14 generates a 10 point CVT on [0,1].
%
%  Discussion:
%
%    Generate 10 CVT points on the interval [0,1].
%    We expect them to be at 1/20, 3/20, 5/20, ..., 19/20.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 September 2010
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST14\n' );
  fprintf ( 1, '  Generate a CVT in the interval [0,1] using 10 points.\n' );
  fprintf ( 1, '  Exact answer: { 0.05, 0.15, 0.25, ..., 0.85, 0.95 }\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  It turns out that, for a fixed number of points N,\n' );
  fprintf ( 1, '  a 1D problem will converge much more slowly than for\n' );
  fprintf ( 1, '  cases where the dimension is higher.\n' );

  dim_num = 1;
  n = 10;
  batch = 10000;
  init = 1;
  init_string = 'uniform';
  it_max = 50;
  it_fixed = 1;
  sample = 1;
  sample_num = 100000;
  sample_string = 'uniform';
  seed = 123456789;
  r = [];

  seed_init = seed;

  [ r, seed, it_num, it_diff, energy ] = cvt ( dim_num, n, batch, init, ...
    sample, sample_num, it_max, it_fixed, seed, r );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Dimension DIM_NUM =        %12d\n', dim_num );
  fprintf ( 1, '  Number of points N =       %12d\n', n );
  fprintf ( 1, '  Initial SEED =             %12d\n', seed_init );
  fprintf ( 1, '  Current SEED =             %12d\n', seed );
  fprintf ( 1, '  INIT =                    "%s".\n', init_string );
  fprintf ( 1, '  Max iterations IT_MAX =    %12d\n', it_max );
  fprintf ( 1, '  IT_FIXED (fixed samples) = %12d\n', it_fixed );
  fprintf ( 1, '  Iterations IT_NUM =        %12d\n', it_num );
  fprintf ( 1, '  Difference IT_DIFF =       %14f\n', it_diff );
  fprintf ( 1, '  CVT ENERGY =               %14f\n', energy );
  fprintf ( 1, '  SAMPLE =                  "%s".\n', sample_string );
  fprintf ( 1, '  Samples SAMPLE_NUM    =    %12d\n', sample_num );
  fprintf ( 1, '  Sampling BATCH size =      %12d\n', batch );
  fprintf ( 1, '  EPSILON (unit roundoff) =  %12e\n', eps );

  r = sort ( r );

  r8mat_transpose_print ( dim_num, n, r, '  Generators (rows):' );

  return
end


text:
%%
 %  Copyright (c) 2014, Facebook, Inc.
 %  All rights reserved.
 %
 %  This source code is licensed under the BSD-style license found in the
 %  LICENSE file in the root directory of this source tree. An additional grant
 %  of patent rights can be found in the PATENTS file in the same directory.
 %
 %%
 function level2_features = generate_level2_features_cnn(level1_features, gfeatures, config)
 %combine deep learning features/scores for a second level feature
 %is_score = 1 means combining scores, otherwise combining the features.
 assert(~isempty(level1_features));

 basis = size(level1_features{1}.features,2);
 is_score = (basis==1);
 num_poselets = length(level1_features);

 num_examples = size(gfeatures,1);

 level2_features = zeros(num_examples, num_poselets*basis, 'single');

 for i = 1:num_poselets
     %if features is empty, meaning no cnn model of this poselet
     if isempty(level1_features{i})
         continue;
     end

     if is_score
         level2_features(level1_features{i}.used,i) = (level1_features{i}.features-0.5);
     else
         level2_features(level1_features{i}.used,(i-1)*basis+1:i*basis) = level1_features{i}.features;
     end
 end

 switch config.CNN_MODEL_TYPE
     case 'combined'
         level2_features = [level2_features gfeatures];
         fprintf('Combined model\n');
     case 'holistic'
         level2_features = gfeatures;
         fprintf('Global-only model\n');
     case 'parts'
         fprintf('Parts-only model\n');
 end

 %normalization
 switch config.CNN_NORMALIZATION_TYPE
     case 'max'
         fprintf('Normalization: max\n');
         level2_features = level2_features ./repmat(max(level2_features,[],2),1,size(level2_features,2));
     case 'L1'
         fprintf('Normalization: L1\n');
         level2_features = level2_features ./repmat(sum(abs(level2_features),2),1,size(level2_features,2));
     case 'power'
         ppp = 0.3;
         for i = 1:size(level2_features,1)
             level2_features(i,:) = sign(level2_features(i,:)).*abs(level2_features(i,:)).^ppp;
         end
         fprintf('Normalization: power %f\n',ppp);
 end
 assert(~any(isnan(level2_features(:))));
 end

author: facebookarchive

repo: pose-aligned-deep-networks