text:
function h = plot_correlation(X,Y,varargin)
% :Usage:
% ::
%
%    handles = plot_correlation(X,Y,varargin)
%
% plots robust or OLS simple or partial correlations
% replaces prplot and plot_correlation_samefig. *Remove NaNs from X and Y
% before calling this function.*
%
% :Inputs:
%
%   **X:**
%        is matrix of columns of interest plus nuisance
%        default is to plot partial effect of 1st column
%
%   **Y:**
%        is one or more columns of data
%
% :Optional Inputs:
%
%   **'robust':**
%        robust IRLS plot
%
%   **'noprint':**
%        suppress text output
%
%   **'doquad':**
%        quadratic term; not tested, may not work
%
%   **'col':**
%        followed by column of interest
%
%   **'labels':**
%        followed by cell array of text labels for each obs.
%
%   **'colors':**
%        followed by cell array of colors for each column of Y
%
%   **'ylabel':**
%        followed by y-axis label string
%
%   **'xlabel':**
%        followed by x-axis label string
%
%   **'weights':**
%        followed by weights that override any computed ones
%
% :Examples: Plot robust partial corr. 2 of X against col. 17 of Y,
% controlling for other X
% ::
%
%    figure;
%    h = plot_correlation(X,Y(:,17),'col',2,'robust','ylabel','Brain
%                         data','xlabel','Order effect');
%
% Plot Col. 1 of X vs. Y in red squares
% ::
%
%    figure;
%    h = plot_correlation(X,Y(:,17),'robust','colors',{'rs'});
%
%    tor_fig;
%
% ..
%    tor wager, august 2006
% ..

% ---------------------------------------------
% default behaviors
% ---------------------------------------------
dorobust = 0;
doprint = 1;
wh_interest = 1;
doquad = 0;
mylabels = [];
mycol = {'ko' 'rv' 'bs' 'gd' 'y^' 'cv' 'mx'};
ylabelstr = 'Contrast beta value';
xlabelstr = 'Behavioral score';

% ---------------------------------------------
% Optional inputs
% ---------------------------------------------

for i = 1:length(varargin)
    if isstr(varargin{i})
        switch varargin{i}
            % reserved keywords
            case 'robust', dorobust = 1;
            case 'noprint', doprint = 0;
            case 'doquad', doquad = 1;

            % functional commands
            case 'labels', labels = varargin{i+1};
            case 'colors', mycol = varargin{i+1};
            case 'ylabel', ylabelstr = varargin{i+1}; varargin{i+1} = [];
            case 'xlabel', xlabelstr = varargin{i+1}; varargin{i+1} = [];
            case 'col', wh_interest = varargin{i+1};
            case 'weights', myweights = varargin{i+1};
            otherwise, warning(['Unknown input string option:' varargin{i}]);
        end
    end
end

% ---------------------------------------------
% drop NANs
% ---------------------------------------------
todrop = isnan(X) | isnan(Y);
if any(todrop)
    warning('Dropping %d observations due to NaNs', sum(todrop))
    X(todrop) = [];
    Y(todrop) = [];
end

% ---------------------------------------------
% intercept
% ---------------------------------------------
wh_intercept = find(all(diff(X) < eps));
if isempty(wh_intercept)
    wh_intercept = size(X,2) + 1;
    X(:,wh_intercept) = 1;
end

ny = size(Y,2);
while length(mycol) < ny, mycol = [mycol mycol]; end

for i = 1:ny
    [x,y,r,p,se,meany,stats] = partialcor(X,Y(:,i),wh_interest,doprint,dorobust);

    % b(1) is param, b(2) is intercept for partial plot
    b = stats.b([wh_interest wh_intercept]);

    % weights
    w = stats.w;
    if exist('myweights','var'), w = myweights; end

    h{i} = makefigure(x,y,mycol{i},mylabels,doquad,b,w,xlabelstr,ylabelstr);

    text(min(X(:,1)),max(y),sprintf('r = %3.2f',r),'FontSize',16,'Color',mycol{i}(1));

end

return

function h = makefigure(xvec,yvec,mycol,mylabels,doquad,b,varargin)
% h = makefigure(xvec,yvec,mycol,mylabels,doquad,b,[weights],[xlabel],[ylabel])

ylabelstr = 'Contrast beta value';
xlabelstr = 'Behavioral score';

hold on; grid on; set(gca,'FontSize',18)

if length(varargin) > 0
    % ROBUST (or we just have weights)
    w = varargin{1};
    for i = 1:length(xvec)

        if isnan(xvec(i)) | isnan(yvec(i)), i, continue, end

        h(i) = plot(xvec(i),yvec(i),mycol,'LineWidth',.5,'MarkerSize',8, ...
            'MarkerFaceColor',mycol(1));
        set(h(i),'MarkerFaceColor',[repmat(1-w(i),1,3)] )
    end

    % set axis
    xlims = [min(xvec) max(xvec)];
    ylims = [min(yvec) max(yvec)];
    xrange = (xlims(2) - xlims(1)) *.1; % % margin
    yrange = (ylims(2) - ylims(1)) *.1;
    xlims = xlims + [-xrange xrange];
    ylims = ylims + [-yrange yrange];

    set(gca,'Xlim',xlims,'YLim',ylims);

    % plot regression line by hand
    xl = get(gca,'Xlim'); yl = get(gca,'Ylim');
    x = xl(1)-5:xl(2)+5;
    h2 = plot(x,b(1) * x + b(2),'-','Color',mycol(1),'LineWidth',2);
    h = [h h2];
    set(gca,'XLim',xl,'YLim',yl)

else
    % not robust

    h = plot(xvec,yvec,mycol,'LineWidth',3,'MarkerSize',6,'MarkerFaceColor',mycol(1));

    if doquad
        %refcurve(doquad)
    else
        try
            refline
        catch
            tmp = get(gca,'XLim');
            x = min(xvec)-std(xvec):min(std(xvec),.01):max(xvec)+std(xvec);
            plot(x,b(1) * x + b(2),[mycol(1) '-'],'LineWidth',.5)
            set(gca,'XLim',tmp);

        end
    end

end

drawnow

if length(varargin) > 1
    xlabelstr = varargin{2};
end

if length(varargin) > 2
    ylabelstr = varargin{3};
end

if ~isempty(mylabels)
    for j = 1:length(xvec)
        text(xvec(j),yvec(j),mylabels{j},'FontWeight','b')
    end
end

ylabel(ylabelstr)
xlabel(xlabelstr)

return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function g = gaussianKernGradient(kern, x, varargin)

% GAUSSIANKERNGRADIENT Gradient of gaussian kernel's parameters.
% FORMAT
% DESC computes the gradient of
%	functions with respect to the gaussian kernel's
%	parameters. As well as the kernel structure and the input positions,
%	the user provides a matrix PARTIAL which gives the partial
%	derivatives of the function with respect to the relevant elements of
%	the kernel matrix.
% RETURN g:  gradients of the function of interest with respect to the
%	   kernel parameters. The ordering of the vector should match that
%	   provided by the function kernExtractParam.
% ARG kern : the kernel structure for which the gradients are being
%	   computed.
% ARG x : the input locations for which the gradients are being
%	   computed.
% ARG partial : matrix of partial derivatives of the function of
%	   interest with respect to the kernel matrix. The argument takes the
%	   form of a square matrix of dimension  numData, where numData is
%	   the number of rows in X.
%
% FORMAT
% DESC  computes the derivatives
%	as above, but input locations are now provided in two matrices
%	associated with rows and columns of the kernel matrix.
% RETURN g : gradients of the function of interest with respect to the
%	   kernel parameters.
% ARG kern : the kernel structure for which the gradients are being
%	   computed.
% ARG x1 : the input locations associated with the rows of the kernel
%	   matrix.
% ARG x2 : the input locations associated with the columns of the kernel
%	   matrix.
% ARG partial : matrix of partial derivatives of the function of
%	   interest with respect to the kernel matrix. The matrix should have
%	   the same number of rows as X1 and the same number of columns as X2
%	   has rows.
%
% SEEALSO : gaussianKernParamInit, kernGradient,
% gaussianKernDiagGradient, kernGradX
%
% COPYRIGHT : Mauricio A. Alvarez and Neil D. Lawrence, 2008
%
% MODIFICATIONS : Mauricio A. Alvarez, 2009.

% KERN

if nargin < 4
    x2 = x;
    covPar = varargin{1};
else
    x2 = varargin{1};
    covPar = varargin{2};
end

if kern.isArd
    [K, Kbase] = gaussianKernCompute(kern, x, x2);
    matGrad = zeros(kern.inputDimension,1);
    temp = 0.5*covPar.*K;
    for i = 1:kern.inputDimension,
        pX = x(:,i);
        X = pX(:, ones(1, size(x2,1)));
        pX2 = x2(:,i)';
        X2 = pX2(ones(size(x,1),1), :);
%         X = repmat(x(:,i),1, size(x2,1));
%         X2 = repmat(x2(:,i)',size(x,1),1);
        X_X2 = X - X2;
        matGrad(i) = -sum(sum(temp.*X_X2.*X_X2));
    end
else
    [K, Kbase, n2] = gaussianKernCompute(kern, x, x2);
    matGrad = - 0.5*sum(sum(covPar.*K.*n2,2));

end
g = [matGrad(:)' sum(sum(covPar.*Kbase))];

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [white, anatpath] = BuildWhiteVolume(hemisphere, fName)

% [white, path] = BuildWhiteVolume([hemisphere, fName]);
%
% Create a binary map of the white matter of the designated hemispheres.
% Optional input hemisphere is 0 for both [default], 1 for left, or 2 for
% right. Optional input fName allows specification of classFile name;
% default guesses file name from anatomy path.
%
% Ress, 07/05

mrGlobals

if ~exist('hemisphere', 'var') | isempty(hemisphere), hemisphere = 0; end

switch hemisphere
  case 0
    hemispheres = {'Left' 'Right'};
  case 1
    hemispheres = {'Left'};
  case 2
    hemispheres = {'Right'};
  otherwise, hemispheres = {'Left' 'Right'};
end

nHemi = length(hemispheres);

if ~exist('fName', 'var') | isempty(fName)
  if exist('vANATOMYPATH', 'var')
    [anatpath, name] = fileparts(vANATOMYPATH);
    for iH=1:nHemi
      dName = fullfile(anatpath, hemispheres{iH});
      if exist(dName, 'dir')
        fList = dir(fullfile(dName, '*.?lass'));
        if ~isempty(fList), fName{iH} = fullfile(dName, fList(1).name); end
      else
        disp(['Could not find ' hemispheres{iH} ' hemisphere class file!'])
        nHemi = nHemi - 1;
      end
    end
  end
end

if ~exist('fName', 'var') | isempty(fName)
  disp('Must provide class-file name');
  dist = [];
  return
else
  if ~iscell(fName), fName = {fName}; end
end

% Get classification volume and create binary map of white matter
for iH=1:nHemi
  class = readClassFile(fName{iH}, 0, 0);
  if exist('white', 'var')
    white = white | (class.data == class.type.white);
  else
    white = (class.data == class.type.white);
  end
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
classdef SMOP2 < PROBLEM
% <multi/many> <real> <large/none> <expensive/none> <sparse/none>
% Benchmark MOP with sparse Pareto optimal solutions
% theta --- 0.1 --- Sparsity of the Pareto set

%------------------------------- Reference --------------------------------
% Y. Tian, X. Zhang, C. Wang, and Y. Jin, An evolutionary algorithm for
% large-scale sparse multi-objective optimization problems, IEEE
% Transactions on Evolutionary Computation, 2020, 24(2): 380-393.
%------------------------------- 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".
%--------------------------------------------------------------------------

    properties(Access = private)
        theta = 0.1;    % Sparsity of the Pareto set
    end
    methods
        %% Default settings of the problem
        function Setting(obj)
            obj.theta = obj.ParameterSet(0.1);
            if isempty(obj.M); obj.M = 2; end
            if isempty(obj.D); obj.D = 100; end
            obj.lower    = [zeros(1,obj.M-1)+0,zeros(1,obj.D-obj.M+1)-1];
            obj.upper    = [zeros(1,obj.M-1)+1,zeros(1,obj.D-obj.M+1)+2];
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,X)
            K = ceil(obj.theta*(obj.D-obj.M+1));
            g = sum(g2(X(:,obj.M:obj.M+K-1),pi/3),2) + sum(g3(X(:,obj.M+K:end),0),2);
            PopObj = repmat(1+g/(obj.D-obj.M+1),1,obj.M).*fliplr(cumprod([ones(size(X,1),1),X(:,1:obj.M-1)],2)).*[ones(size(X,1),1),1-X(:,obj.M-1:-1:1)];
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(obj,N)
            R = UniformPoint(N,obj.M);
        end
        %% Generate the image of Pareto front
        function R = GetPF(obj)
            if obj.M == 2
                R = obj.GetOptimum(100);
            elseif obj.M == 3
                a = linspace(0,1,10)';
                R = {a*a',a*(1-a'),(1-a)*ones(size(a'))};
            else
                R = [];
            end
        end
    end
end

function g = g2(x,t)
    g = 2*(x-t).^2 + sin(2*pi*(x-t)).^2;
end

function g = g3(x,t)
    g = 4-(x-t)-4./exp(100*(x-t).^2);
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [sl, sh] = lowpass(s, lambda, npad)

% lowpass -- Lowpass filter image and return low and high frequency
%            components, consisting of the lowpass filtered image and
%            its difference with the input image. The lowpass filter
%            is equivalent to Tikhonov regularization with lambda as
%            the regularization parameter and a discrete gradient as
%            the operator in the regularization term.
%
% Usage:
%       [sl, sh] = lowpass(s, lambda, npad)
%
% Input:
%       s         Input image or 3d array of images
%       lambda    Regularization parameter controlling lowpass filtering
%       npad      Number of samples to pad at image boundaries
%
% Output:
%       sl        Lowpass component
%       sh        Highpass component
%
%
% Author: Brendt Wohlberg <brendt@lanl.gov>  Modified: 2015-04-09
%
% This file is part of the SPORCO library. Details of the copyright
% and user license can be found in the 'License' file distributed with
% the library.

if nargin < 3,
  npad = 16;
end

grv = [-1 1];
gcv = [-1 1]';
Gr = fft2(grv, size(s,1)+2*npad, size(s,2)+2*npad);
Gc = fft2(gcv, size(s,1)+2*npad, size(s,2)+2*npad);
A = 1 + lambda*conj(Gr).*Gr + lambda*conj(Gc).*Gc;
sp = padarray(s, [npad npad], 'symmetric', 'both');
slp = ifft2(bsxfun(@rdivide, fft2(sp), A), 'symmetric');
sl = slp((npad+1):(size(slp,1)-npad), (npad+1):(size(slp,2)-npad), :);
sh = s - sl;

return

author: xingchenzhang

repo: VIFB

save_path: github-repos/MATLAB/xingchenzhang-VIFB


text:
function n = p14_n ( )

%*****************************************************************************80
%
%% P14_N returns the number of variables for problem 14.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 March 2000
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer N.  If N is positive, N represents the
%    only legal value for N for this problem.  If N is
%    negative, the absolute value of N is the least legal
%    value of N, but other values are allowable.
%
  n = - 1;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [rxnList, rxnFormulaList] = findRxnsFromMets(model, metList, varargin)
% Returns a list of reactions in which at least one
% metabolite listed in metList participates.
%
% USAGE:
%
%    [rxnList, rxnFormulaList] = findRxnsFromMets(model, metList, printFlag)
%
% INPUTS:
%    model:             COBRA model structure
%    metList:           Metabolite list
%
% OPTIONAL INPUTS:
%    printFlag:          Print reaction formulas to screen (Default = false)
%    Property/Value:    Allowed Properties are:
%                        * `containsAll` - If true only reactions containing all metabolites in metList are returned (Default: false)
%                        * `printFlag` - as above, will overwrite a `printFlag` set individually (Default: false)
%                        * `producersOnly` - Only return reactions which produce any of the given metabolites. (Default: false)
%                        * `consumersOnly` - Only return reactions which consume any of the given metabolites. (Default: false)
%                        * `exclusive` - Only return those reactions, which do not contain any metabolites but those given (Default: false)
%
% OUTPUTS:
%    rxnList:           List of reactions
%    rxnFormulaList:    Reaction formulas coresponding to `rxnList`
%
% .. Authors:
%       - Richard Que (08/12/2010)
%       - Almut Heinken (09/25/2015)- made change so formulas are not printed if reaction list is empty.
%       - Thomas Pfau (21/1/2016) - Additional Options, and minimal speedup of the indexing, also updated behaviour of verbFlag to accurately reflect the description.

verbFlag = false;
containsAll = false;
if isa(metList,'char')
    metList = {metList};
end

if mod(numel(varargin), 2) == 1 % the first argument has to be verbFlag, the remaining are property/value pairs
    verbFlag = varargin{1};
    if numel(varargin) > 1
        varargin = varargin(2:end);
    end
end

parser = inputParser();
parser.addParameter('containsAll',false,@(x) (isnumeric(x) && (x==1 || x ==0)) || islogical(x));
parser.addParameter('producersOnly',false,@(x) (isnumeric(x) && (x==1 || x ==0)) || islogical(x));
parser.addParameter('consumersOnly',false,@(x) (isnumeric(x) && (x==1 || x ==0)) || islogical(x));
parser.addParameter('printFlag',verbFlag,@(x) isnumeric(x) );
parser.addParameter('verbFlag',verbFlag,@(x) isnumeric || islogical(x) ); % backward compatability
parser.addParameter('exclusive',false,@(x)  (isnumeric(x) && (x==1 || x ==0)) || islogical(x));

parser.parse(varargin{:});
%Verbosity flag backward compatability
if ismember('verbFlag',parser.UsingDefaults)
    verbFlag = parser.Results.printFlag;
elseif ismember('printFlag',parser.UsingDefaults)
    verbFlag = parser.Results.verbFag;
else
    verbFlag = max([parser.Results.printFlag,parser.Results.verbFlag]);
end

containsAll = parser.Results.containsAll;
producersOnly = parser.Results.producersOnly;
consumersOnly = parser.Results.consumersOnly;
exclusive = parser.Results.exclusive;

%Find met indicies
index = ismember(model.mets,metList);
if producersOnly && consumersOnly
    producers = sum(model.S(index,:) > 0,1)';
    consumers = sum(model.S(index,:) < 0,1)';
    rels = producers > 0 & model.ub > 0 | consumers > 0 & model.lb < 0 | producers > 0 & model.lb < 0 | consumers > 0 & model.ub > 0;
    totals = producers + consumers;
elseif producersOnly
    producers = sum(model.S(index,:) > 0,1)';
    consumers = sum(model.S(index,:) < 0,1)';
    rels = producers > 0 & model.ub > 0 | consumers > 0 & model.lb < 0;
    totals = producers + consumers;
elseif consumersOnly
    producers = sum(model.S(index,:) > 0,1)';
    consumers = sum(model.S(index,:) < 0,1)';
    rels = producers > 0 & model.lb < 0 | consumers > 0 & model.ub > 0;
    totals = producers + consumers;
else
    totals = sum(model.S(index,:) ~= 0,1);
    rels = totals > 0;
end

if exclusive
    others = sum(model.S(~index,:) ~= 0);
    %exclude anything that has other reactants.
    rels = rels & others == 0;
end

if containsAll
   rxnList = model.rxns(totals == numel(metList) & rels);
else
    %rxns = repmat(model.rxns,1,length(index));
    %find reactions i.e. all columns with at least one non zero value
    rxnList = model.rxns(totals > 0 & rels);
end

if (nargout > 1) || verbFlag
    if ~isempty(rxnList)
        rxnFormulaList = printRxnFormula(model,rxnList,verbFlag);
    else
        rxnFormulaList={};
    end
end

author: opencobra

repo: cobratoolbox


text:
function b = legcoeffs(f, varargin)
%LEGCOEFFS   Compute Legendre series coefficients of a BNDFUN object.
%   B = LEGCOEFFS(F) returns the Legendre series coefficients of BNDFUN F, so
%   that F = B(N+1)*P_N + ... + B(1)*P_0, where P_k is the kth Legendre
%   polynomial (scaled to the domain of F).
%
%   If F is an array-valued BNDFUN, then a matrix of coefficients is returned so
%   that F(:,k) = B(N+1,k)*P_N + ... + B(1,k)*P_0.
%
% See also CHEBCOEFFS.

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

b = legcoeffs(f.onefun, varargin{:});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [weights, th1, th2] = gncWeightsUpdate(weights, mu, residuals, barc2)
%GNC - GNC weights update for the TLS cost function
% Reference:
%    - "Graduated Non-Convexity for Robust Spatial Perception: From Non-Minimal Solvers to Global Outlier Rejection"
%      Yang, Antonante, Tzoumas, Carlone (2020). IEEE Robotics and Automation Letters (RA-L), 5(2), 1127–1134.
%      https://arxiv.org/pdf/1909.08605.pdf
%    - "Outlier-Robust Estimation: Hardness, Minimally-Tuned Algorithms, and Applications"
%      Antonante, Tzoumas, Yang, Carlone (2020).
%      https://arxiv.org/pdf/2007.15109.pdf
%

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

th1 = (mu+1)/mu * barc2;
th2 = (mu)/(mu+1) * barc2; % th1 > th2
for k = 1:length(residuals)
    if residuals(k) - th1 >= 0
        weights(k) = 0;
    elseif residuals(k) - th2 <= 0
        weights(k) = 1;
    else
        weights(k) = sqrt( barc2*mu*(mu+1)/residuals(k) ) - mu;
        assert(weights(k)>= 0 && weights(k) <=1, 'weights %g calculation wrong!', weights(k));
    end
end
end

author: MIT-SPARK

repo: GNC-and-ADAPT


text:
function TFM = createRotationVectorPoint3d(A,B,P)
%CREATEROTATIONVECTORPOINT3D Calculates the rotation between two vectors.
%   around a point
%
%   TFM = createRotationVectorPoint3d(A,B,P) returns the transformation
%   to rotate the vector A in the direction of vector B around point P
%
%   Example
%     A=-5+10.*rand(1,3);
%     B=-10+20.*rand(1,3);
%     P=-50+100.*rand(1,3);
%     ROT = createRotationVectorPoint3d(A,B,P);
%     C = transformVector3d(A,ROT);
%     figure('color','w'); hold on; view(3)
%     drawPoint3d(P,'k')
%     drawVector3d(P, A,'r')
%     drawVector3d(P, B,'g')
%     drawVector3d(P, C,'r')
%
%   See also
%   transformPoint3d, createRotationVector3d
%
% ---------
% Author: oqilipo
% Created: 2017-08-07
% Copyright 2017

P = reshape(P,3,1);

% Translation from P to origin
invtrans = [eye(3),-P; [0 0 0 1]];

% Rotation from A to B
rot = createRotationVector3d(A, B);

% Translation from origin to P
trans = [eye(3),P; [0 0 0 1]];

% Combine
TFM = trans*rot*invtrans;

end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [segments, fs] = detectVoiced(wavFileName,t)

%
% function [segments, fs] = detectVoiced(wavFileName)
%
% Theodoros Giannakopoulos
% http://www.di.uoa.gr/~tyiannak
%
% (c) 2010
%
% This function implements a simple voice detector. The algorithm is
% described in more detail, in the readme.pdf file
%
% ARGUMENTS:
%  - wavFileName: the path of the wav file to be analyzed
%  - t: if provided, the detected voiced segments are played and some
%  intermediate results are also ploted
%
% RETURNS:
%  - segments: a cell array of M elements. M is the total number of
%  detected segments. Each element of the cell array is a vector of audio
%  samples of the respective segment.
%  - fs: the sampling frequency of the audio signal
%
% EXECUTION EXAMPLE:
%
% [segments, fs] = detectVoiced('example.wav',1);
%

% Check if the given wav file exists:
fp = fopen(wavFileName, 'rb');
if (fp<0)
	fprintf('The file %s has not been found!\n', wavFileName);
	return;
end
fclose(fp);

% Check if .wav extension exists:
if  (strcmpi(wavFileName(end-3:end),'.wav'))
    % read the wav file name:
    [x,fs] = wavread(wavFileName);
else
    fprintf('Unknown file type!\n');
    return;
end

% Convert mono to stereo
if (size(x, 2)==2)
	x = mean(x')';
end

% Window length and step (in seconds):
win = 0.050;
step = 0.050;

%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  THRESHOLD ESTIMATION
%%%%%%%%%%%%%%%%%%%%%%%%%%%

Weight = 5; % used in the threshold estimation method

% Compute short-time energy and spectral centroid of the signal:
Eor = ShortTimeEnergy(x, win*fs, step*fs);
Cor = SpectralCentroid(x, win*fs, step*fs, fs);

% Apply median filtering in the feature sequences (twice), using 5 windows:
% (i.e., 250 mseconds)
E = medfilt1(Eor, 5); E = medfilt1(E, 5);
C = medfilt1(Cor, 5); C = medfilt1(C, 5);

% Get the average values of the smoothed feature sequences:
E_mean = mean(E);
Z_mean = mean(C);

% Find energy threshold:
[HistE, X_E] = hist(E, round(length(E) / 10));  % histogram computation
[MaximaE, countMaximaE] = findMaxima(HistE, 3); % find the local maxima of the histogram
if (size(MaximaE,2)>=2) % if at least two local maxima have been found in the histogram:
    T_E = (Weight*X_E(MaximaE(1,1))+X_E(MaximaE(1,2))) / (Weight+1); % ... then compute the threshold as the weighted average between the two first histogram's local maxima.
else
    T_E = E_mean / 2;
end

% Find spectral centroid threshold:
[HistC, X_C] = hist(C, round(length(C) / 10));
[MaximaC, countMaximaC] = findMaxima(HistC, 3);
if (size(MaximaC,2)>=2)
    T_C = (Weight*X_C(MaximaC(1,1))+X_C(MaximaC(1,2))) / (Weight+1);
else
    T_C = Z_mean / 2;
end

% Thresholding:
Flags1 = (E>=T_E);
Flags2 = (C>=T_C);
flags = Flags1 & Flags2;

if (nargin==2) % plot results:
	clf;
	subplot(3,1,1); plot(Eor, 'g'); hold on; plot(E, 'c'); legend({'Short time energy (original)', 'Short time energy (filtered)'});
    L = line([0 length(E)],[T_E T_E]); set(L,'Color',[0 0 0]); set(L, 'LineWidth', 2);
    axis([0 length(Eor) min(Eor) max(Eor)]);

    subplot(3,1,2); plot(Cor, 'g'); hold on; plot(C, 'c'); legend({'Spectral Centroid (original)', 'Spectral Centroid (filtered)'});
	L = line([0 length(C)],[T_C T_C]); set(L,'Color',[0 0 0]); set(L, 'LineWidth', 2);
    axis([0 length(Cor) min(Cor) max(Cor)]);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  SPEECH SEGMENTS DETECTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
count = 1;
WIN = 5;
Limits = [];
while (count < length(flags)) % while there are windows to be processed:
	% initilize:
	curX = [];
	countTemp = 1;
	% while flags=1:
	while ((flags(count)==1) && (count < length(flags)))
		if (countTemp==1) % if this is the first of the current speech segment:
			Limit1 = round((count-WIN)*step*fs)+1; % set start limit:
			if (Limit1<1)	Limit1 = 1; end
		end
		count = count + 1; 		% increase overall counter
		countTemp = countTemp + 1;	% increase counter of the CURRENT speech segment
	end

	if (countTemp>1) % if at least one segment has been found in the current loop:
		Limit2 = round((count+WIN)*step*fs);			% set end counter
		if (Limit2>length(x))
            Limit2 = length(x);
        end

        Limits(end+1, 1) = Limit1;
        Limits(end,   2) = Limit2;
    end
	count = count + 1; % increase overall counter
end

%%%%%%%%%%%%%%%%%%%%%%%
% POST - PROCESS      %
%%%%%%%%%%%%%%%%%%%%%%%

% A. MERGE OVERLAPPING SEGMENTS:
RUN = 1;
while (RUN==1)
    RUN = 0;
    for (i=1:size(Limits,1)-1) % for each segment
        if (Limits(i,2)>=Limits(i+1,1))
            RUN = 1;
            Limits(i,2) = Limits(i+1,2);
            Limits(i+1,:) = [];
            break;
        end
    end
end

% B. Get final segments:
segments = {};
for (i=1:size(Limits,1))
    segments{end+1} = x(Limits(i,1):Limits(i,2));
end

if (nargin==2)
    subplot(3,1,3);
    % Plot results and play segments:
    time = 0:1/fs:(length(x)-1) / fs;
    for (i=1:length(segments))
        hold off;
        P1 = plot(time, x); set(P1, 'Color', [0.7 0.7 0.7]);
        hold on;
        for (j=1:length(segments))
            if (i~=j)
                timeTemp = Limits(j,1)/fs:1/fs:Limits(j,2)/fs;
                P = plot(timeTemp, segments{j});
                set(P, 'Color', [0.4 0.1 0.1]);
            end
        end
        timeTemp = Limits(i,1)/fs:1/fs:Limits(i,2)/fs;
        P = plot(timeTemp, segments{i});
        set(P, 'Color', [0.9 0.0 0.0]);
        axis([0 time(end) min(x) max(x)]);
        sound(segments{i}, fs);
        clc;
        fprintf('Playing segment %d of %d. Press any key to continue...', i, length(segments));
        pause
    end
    clc
    hold off;
    P1 = plot(time, x); set(P1, 'Color', [0.7 0.7 0.7]);
    hold on;
    for (i=1:length(segments))
        for (j=1:length(segments))
            if (i~=j)
                timeTemp = Limits(j,1)/fs:1/fs:Limits(j,2)/fs;
                P = plot(timeTemp, segments{j});
                set(P, 'Color', [0.4 0.1 0.1]);
            end
        end
        axis([0 time(end) min(x) max(x)]);
    end
end

author: Sable

repo: mcbench-benchmarks


text:
% PURPOSE: Gets the baseline (mean) value
%
%
% FORMAT:
%
% blv = blvalue2(datax, timex, blcorr)
%
% INPUT:
%
% datax      - input data
% timex      - time vector
% blcorr     - time window for getting the mean value
%
% OUTPUT:
%
% blv        - mean value from window "blcorr"
%
%
% Example
% Get the baseline value for a window of -200 to 0 ms, at bin 4, channel 23
%
% blv = blvalue(ERP.bindata(23,:,4), ERP.times, [-200 0])
%
%
% See also blvalue2.m geterpvalues.m
%
%
% *** This function is part of ERPLAB Toolbox ***
% Author: Javier Lopez-Calderon
% Center for Mind and Brain
% University of California, Davis,
% Davis, CA
% 2013

function blv = blvalue2(datax, timex, blcorr)

%
% Baseline assessment
%
if ischar(blcorr)
      if ~strcmpi(blcorr,'no') && ~strcmpi(blcorr,'none')

            if strcmpi(blcorr,'pre')
                  [bbxx bb] = closest(timex, 0);    % zero-time locked
                  aa = 1;
            elseif strcmpi(blcorr,'post')
                  bb = length(timex);
                  [aax aa] = closest(timex, 0);
            elseif strcmpi(blcorr,'all') || strcmpi(blcorr,'whole')
                  bb = length(timex);
                  aa = 1;
            else
                  blcnum = str2num(blcorr); % in ms

                  %
                  % Check & fix baseline range
                  %
                  if blcnum(1)<min(timex)
                        blcnum(1) = min(timex); %ms
                  end
                  if blcnum(2)>max(timex)
                        blcnum(2) = max(timex); %ms
                  end

                  [xxx, cindex] = closest(timex, blcnum); % 04/21/2011
                  aa = cindex(1); % ms to sample pos
                  bb = cindex(2); % ms to sample pos
            end
            blv = mean(datax(aa:bb));
      else
            blv = 0;
      end
else
      %
      % Check & fix baseline range
      %
      if blcorr(1)<min(timex)
            blcorr(1) = min(timex); %ms
      end
      if blcorr(2)>max(timex)
            blcorr(2) = max(timex); %ms
      end
      [xxx, cindex] = closest(timex, blcorr); % 04/21/2011
      aa = cindex(1); % ms to sample pos
      bb = cindex(2); % ms to sample pos
      blv = mean(datax(aa:bb));
end

author: ucdavis

repo: erplab

save_path: github-repos/MATLAB/ucdavis-erplab


text:
function Q=QCoordTurn(T,x,sigmaV2,sigmaTurn2,sigmaLin2)
%%QCOORDTURN Get the discrete-time process noise covariance matrix for a
%            2D or 3D coordinated turn model with a Cartesian state. The
%            turn rate can be specified in terms of a turn rate in radians
%            per second, or in terms of a transversal acceleration.
%            Additionally, a linear acceleration can be given.
%
%INPUTS: T The time-duration of the propagation interval in seconds.
%        x The target state for 2D or 3D motion. If there is no linear
%          acceleration (acceleration along the direction of motion), then
%          x in 2D can either be x=[x;y;xdot;ydot;omega], where omega is
%          the turn rate estimate in radians per second counterclockwise
%          from the x-axis or x can be x=[x;y;xdot;ydot;at] where at is
%          the transversal acceleration, which is orthogonal to the
%          velocity and is defined such that positive values of at map to
%          positive values of omega. In 3D, the state has z and zdot
%          components after y and ydot. If there is a linear acceleration,
%          then in 2D the target state is either x=[x;y;xdot;ydot;omega;al]
%          where omega is the turn rate and al is the linear acceleration
%          or the target state is x=[x;y;xdot;ydot;at;al] if the turn is
%          expressed in terms of a transversal acceleration. Again, in 3D,
%          z and zdot components are respectively added after y and ydot. The
%          dimensionality of the state is used to determine the
%          dimensionality and whether a linear acceleration component is
%          present. The linear acceleration component changes the speed.
%          That means that it acts in the direction of the velocity vector.
%  sigmaV2 The variance driving the process noise affecting the velocity
%          components. This has units of m^2/s^4 and is assumed to be the
%          same in both the x and y dimensions.
% sigmaTurn2 If the turn is specified in terms of a turn rate in radians
%          per second, then this is the variance driving the process noise
%          of the turn rate having units of radians squared per second
%          squared. If the turn is expressed in terms of a transverse
%          acceleration, then this is the variance of the transverse
%          acceleration noise, having units of m^2/s^4.
% sigmaLin2 This parameter is only needed if a linear acceleration is
%          present. It is the variance of the linear acceleration noise
%          having units of m^2/s^4.
%
%OUTPUTS: Q The process noise covariance matrix under a direct discrete-
%           time coordinated turn model where the velocity is specified in
%           Cartesian coordinates. Note that Q is singular, due to the
%           direct-discrete modeling of the process noise.
%
%The basic 2D coordinated turn model is described in Section VA of [1]. It
%is assumed that the continuous-time turn rate model is
%omegaDot=-(1/tau)*Omega+noise, which discretizes to
%omega[k+1]=exp(-T/tau)*omega[k]+noise.
%The turn rate for an unknown process noise is taken from Equation 73 of
%the paper and is a discrete-time approximation, not an exact
%discretization. The ordering of the elements in the state has been changed
%from the paper.
%
%The velocity and position components of the process noise model are
%similar to a discrete white noise acceleration model. Thus, a starting
%point for the process noise parameter sigmaV2 can be found using the
%processNoiseSuggest function with 'PolyKalDirectDisc-ROT' as the model and
%order=1 to cover velocity changes that are not solely due to turns or
%linear accelerations. The changes in the turn rate/transverse acceleration
%that are not due to exponential decay can similarly be covered using
%process noise  starting with sigmaTurn2 being set based on
%processNoiseSuggest with 'PolyKalDirectDisc-ROT' and order=1/2. Similarly,
%a starting point for sigmaLin2 is from processNoiseSuggest with
%'PolyKalDirectDisc-ROT' and order=2.
%
%Note that the process noise added to the turn component and the linear
%acceleration term (if present) is not integrated into the velocity/
%position components. If the only changes allowed/ expected are due to
%modifications in the linear and transverse accelerations (or turn rate),
%then sigmaV2=0. However, regardless of the size of T, the uncertainty in
%the other components would not directly integrate into the position and
%velocity components over any one step. Thus, one would expect the model to
%be particularly poor when the step size is large.
%
%The corresponding transition matrix in 2D is given by the function
%FCoordTurn2D. The corresponding continuous-time functions are
%aCoordTurn2DOmega and aCoordTurn2DTrans with diffusion matrix
%DCoordTurn2D. However, note that the discrete-time functions with
%unknown noise is a direct-discrete model and not a discretization of the
%continuous-time model.
%
%REFERENCES:
%[1] X. R. Li and V. P. Jilkov, "Survey of maneuvering target tracking.
%   Part I: Dynamic models," IEEE Transactions on Aerospace and Electronic
%   Systems, vol. 39, no. 4, pp. 1333-1364, Oct. 2003.
%
%July 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

sigmaV=sqrt(sigmaV2);
sigmaTurn=sqrt(sigmaTurn2);

switch(length(x))
    case 5%2D, there is no linear acceleration.
        G=[T^2/2;%Position x
           T^2/2;%Position y
           T;%Velocity x
           T;%Velocity y
           T];%Turn rate/ linear acceleration
        G(1:4)=G(1:4)*sigmaV;
        G(5)=G(5)*sigmaTurn;
        Q=G*G';
    case 6%2D, there is linear acceleration.
        sigmaLin=sqrt(sigmaLin2);
        G=[T^2/2;%Position x
           T^2/2;%Position y
           T;%Velocity x
           T;%Velocity y
           T;%Turn rate/ linear acceleration
           T];%Transverse acceleration
        G(1:4)=G(1:4)*sigmaV;
        G(5)=G(5)*sigmaTurn;
        G(6)=G(6)*sigmaLin;
        Q=G*G';
    case 7%3D, there is no linear acceleration.
        G=[T^2/2;%Position x
           T^2/2;%Position y
           T^2/2;%Position z
           T;%Velocity x
           T;%Velocity y
           T;%Velocity z
           T];%Turn rate/ linear acceleration
        G(1:6)=G(1:6)*sigmaV;
        G(7)=G(7)*sigmaTurn;
        Q=G*G';
    case 8%3D, there is linear acceleration.
        sigmaLin=sqrt(sigmaLin2);
        G=[T^2/2;%Position x
           T^2/2;%Position y
           T^2/2;%Position z
           T;%Velocity x
           T;%Velocity y
           T;%Velocity z
           T;%Turn rate/ linear acceleration
           T];%Transverse acceleration
        G(1:6)=G(1:6)*sigmaV;
        G(7)=G(7)*sigmaTurn;
        G(8)=G(8)*sigmaLin;
        Q=G*G';
    otherwise
        error('The length of x is neither 5 nor 6 (for 2D) and neither 7 or 8 (for 3D).');
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 myrng()
% Author		: Tiep Vu (http://www.personal.psu.edu/thv102/)
% Time created	: Wed Jan 27 00:23:58 2016
% Last modified	: Wed Jan 27 00:23:59 2016
% Description	: make sure that random function does not return
%                 the same output

    c = clock; %get current time

    t = mod(floor(c(6))*13, 100); % c(6) is 'second'
    for i = 1: t
        randi(t);
    end
end

author: tiepvupsu

repo: DICTOL

save_path: github-repos/MATLAB/tiepvupsu-DICTOL


text:
function [GP dataOut] = iGPRFun(data, options)
GP = [];
dataOut = [];
PtsNumThr = options.PtsNumThr;
epsilon = options.epsilon;
MinPts = options.MinPts;
tData = options.tData;
tDist = options.tDist;
if size(data, 2) <= PtsNumThr
    return;
end
[~, ValidIdx] = Cluster(data, MinPts, epsilon, PtsNumThr);
if length(ValidIdx) <= options.PtsNumThr % point number is smaller than PtsNumThr
    return;
end
GP = GPR_CurbNew_IV(data, ValidIdx, options, []);
dataOut =  GP.data(GP.ValidIdx, : )';   % data0 is valid road edge points.
end

author: DrGabor

repo: LiDAR

save_path: github-repos/MATLAB/DrGabor-LiDAR


text:
%PRWAITBARINIT Low level routine to simplify PRWAITBAR init
%
%	[N,S,COUNT] = PRWAITBARINIT(STRING,N)
%
% INPUT
%   STRING  - String with text to be written in every waitbar,
%             e.g. 'Processing %i items: '. This will be parsed
%             by S = SPRINTF(STRING,N);
%   N       - Total number of items to be processed
%
% OUTPUT
%   N       - Resulting N (Input N may be expression)
%   S       - Resulting STRING (see above)
%   COUNT   - Counter initialisation, COUNT = 0
%
% This routine has to be used in combination with PRWAITBARNEXT, e.g.:
%
%  [n,s,count] = prwaitbarinit('Processing %i items:',size(a,1)*size(a,2));
%  for i=1:size(a,1)
%    for j=1:size(a,2)
%      < process a(i,j) >
%      count = prwaitbarnext(n,s,count);
%    end
%  end
%
% PRWAITBARINIT and PRWAITBARNEXT are written on top of PRWAITBAR. If it
% is possible that loops like the above are not fully processed, it is
% necessary to include as a final call PRWAITBAR(0)

% Copyright: R.P.W. Duin, r.p.w.duin@37steps.com
% Faculty EWI, Delft University of Technology
% P.O. Box 5031, 2600 GA Delft, The Netherlands

function [n,s,count] = prwaitbarinit(ss,n)

s = sprintf(ss,n);
if n > 1
	prwaitbar(n,s);
	prwaitbar(n,1,[s int2str(1)]);
end
count = 0;

return

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
%% flandmark_simple_example.m
% Test flandmark detector
%
% 27-03-12 Michal Uricar

clc;
clearvars; close all;

%% Add path

addpath('./mex/');

DIR = '../data/Images/';
IMGS = dir([DIR '*.jpg']);

% Load flandmark_model into MATLAB memory
model = flandmark_load_model('../data/flandmark_model.dat');

%% Run detector

% run detector 100 times for randomly selected image
for example = 1 : 100
    % pick random image
    idx = randi(numel(IMGS));
    filename = IMGS(idx).name;
    fname = filename(1:end-4);
    % load image and detected face bbox
    I = imread([DIR filename]);
    Ibw = rgb2gray(I);
    bbox = dlmread([DIR fname '.det']);

    % image output
    figure(1);
    imshow(I, [], 'Border', 'tight'); hold on;
    % plotbox(bbox);

    for i = 1 : size(bbox, 1)
        tic
        P = flandmark_detector(Ibw, int32(bbox(i, :)),  model);
        t1 = toc;
        fprintf('MEX:    Elapsed time %f ms\n', t1*1000);

        % show landmarks
        comps = ['S0'; 'S1'; 'S2'; 'S3'; 'S4'; 'S5'; 'S6'; 'S7'];
        plot(P(1, 1), P(2, 1), 'bs', 'LineWidth', 1, 'MarkerSize', 5, 'MarkerFaceColor', 'b');
        text(P(1, 1)+1, P(2, 1)+1, comps(1,:), 'color', 'b', 'FontSize', 12);
        plot(P(1, 2:end), P(2, 2:end), 'rs', 'LineWidth', 1, 'MarkerSize', 5, 'MarkerFaceColor', 'r');
        text(P(1, 2:end)+1, P(2, 2:end)+1, comps(2:end,:), 'color', 'r', 'FontSize', 12);
    end;
    pause(0.2);
end;

%%

clear mex;

author: uricamic

repo: flandmark

save_path: github-repos/MATLAB/uricamic-flandmark


text:
%data comes from the EUROC MAV dataset V1_01_medium
clear;
close all;
addpath myToolbox
addpath filters
addpath data

tMin = 1081; % starting IMU time
tImagesMin = 109; % starting camera time

%% Initialization
freqIMU = 200; %Hz
freqCam = 20; %Hz

ParamGlobal.dirImage = 'V1_02_medium/mav0/cam0/data/';
dirImage = 'V1_02_medium/mav0/cam0/data/';
fileData = 'DATA_MEDIUM.mat';
fileImages = 'fileImages_MEDIUM.mat';
offset = 881;

% IMU and frame date
load(fileData);
t = tIMU;
t = t/10^9; % ns -> s
g = 9.81*[0;0;-1]; %gravity field
omega = omega(:,tMin:end);
acc = acc(:,tMin:end);
t = t(tMin:end);
tIMU = tIMU(tMin:end);
ParamGlobal.tIMU = tIMU;
NbSteps = length(tIMU);
NbStepsMax = 8000;

load(fileImages);
tImages = tImages(tImagesMin:end);
fileImages = fileImages(tImagesMin:end);
ParamGlobal.fileImages = fileImages;
%corect offset
tReal(1:offset) = [];
trajReal.x(:,1:offset) = [];
trajReal.v(:,1:offset) = [];
trajReal.quat(1:offset) = [];
trajReal.psi(1:offset) = [];
trajReal.theta(1:offset) = [];
trajReal.phi(1:offset) = [];
trajReal.omega_b(:,1:offset) = [];
trajReal.a_b(:,1:offset) = [];
obsTimes = zeros(length(t),1);
obsTimes(1:10:end) = 1; % IMU is 10 times faster than camera

%% Init Filter
P0amers = diag([1;1;1]*1.e-3); %initial landmark covariance
R = 1.0^2*eye(2); %measurement noise for one landmark
ParamFilter.NbAmers = 30; %nominal number of landmarks in the state
ParamFilter.NbAmersMin = ParamFilter.NbAmers;
ParamFilter.EcartPixelMax = 20;

% init covariance
p0Rot = (0.01*pi/180)^2;
p0v =  1.e-4;
p0x =  1.e-8;
p0omegab = 1.e-6;
p0ab = 1.e-6;
P0 = diag([p0Rot*ones(3,1);p0v*ones(3,1);p0x*ones(3,1);...
    p0omegab*ones(3,1);p0ab*ones(3,1)]);

% process noises
q_omega = (1.6968e-4)^2*200;
q_a = (2e-3)^2*200;
q_omegab = (1.9393e-5)^2*200;
q_ab = (3e-3)^2*200;
Q = diag([q_omega*ones(3,1);q_a*ones(3,1);q_omegab* ...
    ones(3,1);q_ab*ones(3,1)]);
Qc = chol(Q);
P0 = blkdiag(P0,kron(eye(ParamFilter.NbAmers),P0amers));

%depending on the chosen camera %cam0
ParamFilter.chiC = [0.0148655429818, -0.999880929698, 0.00414029679422, -0.0216401454975;
    0.999557249008, 0.0149672133247, 0.025715529948, -0.064676986768;
    -0.0257744366974, 0.00375618835797, 0.999660727178, 0.00981073058949;
    0.0, 0.0, 0.0, 1.0]; % camera pose
ParamFilter.Pi = [458.654 0 0; 0 457.296 0; 367.215, 248.375 1]';%camera calibration matrix
ParamFilter.cameraParams = cameraParameters('IntrinsicMatrix', ParamFilter.Pi',...
    'RadialDistortion',[-0.28340811, 0.07395907],...
    'TangentialDistortion',[0.00019359, 1.76187114e-05]);

% Initialisation of the state is obtained following [Mur-Artal,2017],
load('data/ORB_SLAM_init.mat');
Rot0 = eul2rotm([trajReal.psi(1),trajReal.theta(1),trajReal.phi(1)]);
x0 = trajReal.x(:,1);
v0 = trajReal.v(:,1);
omega_b0 = orb_slam.omega_b;
a_b0 = orb_slam.a_b;
PosAmers0 = orb_slam.PosAmers;
trackerMain = orb_slam.trackerMain;
trackerBis = orb_slam.trackerBis;
myTracks = orb_slam.myTracks;

IdxImage = 2; % image index

trajL = initTraj(NbSteps);
RotL = Rot0;
vL = v0;
xL = x0;
omega_bL = omega_b0;
a_bL = a_b0;
PosAmersL = PosAmers0;
P_L = P0;
S_L = chol(P_L);
chiL = state2chi(RotL,vL,xL,PosAmersL);

trajR = initTraj(NbSteps);
RotR = Rot0;
vR = v0;
xR = x0;
omega_bR = omega_b0;
a_bR = a_b0;
PosAmersR = PosAmers0;
P_R = P0;
S_R = chol(P_R);
chiR = state2chi(RotR,vR,xR,PosAmersR);

trajRef = initTraj(NbSteps);
RotRef = Rot0;
vRef = v0;
xRef = x0;
omega_bRef = omega_b0;
a_bRef = a_b0;
PosAmersRef = PosAmers0;
P_Ref = blkdiag(P0);
S_Ref = chol(P_Ref);
chiRef = [RotRef xRef;0 0 0 1];
xidotRef = zeros(6,1);

trajU = initTraj(NbSteps);
RotU = Rot0;
vU = v0;
xU = x0;
omega_bU = omega_b0;
a_bU = a_b0;
PosAmersU = PosAmers0;
P_U = blkdiag(P0);
S_U = chol(P_U);

trajI = initTraj(NbSteps);
RotI = Rot0;
vI = v0;
xI = x0;
omega_bI = omega_b0;
a_bI = a_b0;
PosAmersI = PosAmers0;
P_I = blkdiag(P0);

%% Filtering
for i = 2:NbStepsMax
    % propagation
    dt = t(i)-t(i-1);
    omega_i = omega(:,i);
    acc_i = acc(:,i);

    chiAntR = chiR;
    RotR = RotR*expSO3((omega_i-omega_bR)*dt);
    vR = vR+(RotR*(acc_i-a_bR)+g)*dt;
    xR = xR+vR*dt;
    chiR = state2chi(RotR,vR,xR,PosAmersR);
    S_R = rukfPropagation(dt,chiR,chiAntR,omega_bR,a_bR,S_R,omega_i,...
        acc_i,Qc,g);

    % propagation for others filters
    chiAntL = chiL;
    RotL = RotL*expSO3((omega_i-omega_bL)*dt);
    vL = vL+(RotL*(acc_i-a_bL)+g)*dt;
    xL = xL+vL*dt;
    chiL = state2chi(RotL,vL,xL,PosAmersL);
    S_L = lukfPropagation(dt,chiL,chiAntL,omega_bL,a_bL,S_L,omega_i,...
        acc_i,Qc,g);

    chiAntRef = [RotRef,xRef;zeros(1,3) 1];
    vAnt = vRef;
    RotRef = RotRef*expSO3((omega_i-omega_bRef)*dt);
    vRef = vRef+(RotRef*(acc_i-a_bRef)+g)*dt;
    xRef = xRef+vRef*dt;
    xidotRef = xidotRef + [omega_i-omega_bRef;vRef]*dt;
    chiRef = [RotRef,xRef;zeros(1,3) 1];
    S_Ref = ukfRefPropagation(dt,chiRef,chiAntRef,vAnt,omega_bRef,...
        a_bRef,S_Ref,omega_i,acc_i,Qc,g);

    RotAntU = RotU;
    RotU = RotU*expSO3((omega_i-omega_bU)*dt);
    vU = vU+(RotU*(acc_i-a_bU)+g)*dt;
    xU = xU+vU*dt;
    chiU = [RotU,xU;zeros(1,3) 1];
    S_U = ukfPropagation(dt,RotU,RotAntU,vAnt,omega_bU,a_bU,S_U,omega_i,...
        acc_i,Qc,g);

    [RotI,vI,xI,PosAmersI,P_I] = iekfPropagation(dt,RotI,vI,xI,...
        omega_bI,a_bI,PosAmersI,P_I,omega_i,acc_i,Q,g);
    chiI = state2chi(RotI,vI,xI,PosAmersI);

    % if measurement
    if obsTimes(i) == 1
        % track points in image
        [y,yAmers,trackerMain,trackerBis,pointsMain,validityMain,...
            myTracks,pointsBis] = ...
            ObserveLandmarks(trackerMain,trackerBis,dirImage,IdxImage,...
            fileImages,ParamFilter,RotR,xR,PosAmersR,i,S_R,myTracks);

        % update state
        param.yAmers = yAmers;
        [chiR,omega_bR,a_bR,S_R] = rukfUpdate(chiR,omega_bR,...
            a_bR,S_R,y,param,R,ParamFilter);
        [RotR,vR,xR,PosAmersR] = chi2state(chiR);

        % update other filters
        chiL = state2chi(RotL,vL,xL,PosAmersL);
        [chiL,omega_bL,a_bL,S_L] = lukfUpdate(chiL,omega_bL,...
            a_bL,S_L,y,param,R,ParamFilter);
        [RotL,vL,xL,PosAmersL] = chi2state(chiL);

        param.PosAmers = PosAmersRef;
        chiRef = state2chi(RotRef,vRef,xRef,PosAmersRef);
        [chiRef,vRef,PosAmersRef,omega_bRef,a_bRef,S_Ref,xidotRef] = ukfRefUpdate(chiRef,vRef,omega_bRef,...
            a_bRef,S_Ref,y,param,R,ParamFilter,PosAmersRef,xidotRef);
        RotRef = chiRef(1:3,1:3);
        xRef = chiRef(1:3,4);

        param.PosAmers = PosAmersU;
        chiU = state2chi(RotU,vU,xU,PosAmersU);
        [RotU,vU,xU,PosAmersU,omega_bU,a_bU,S_U] = ukfUpdate(RotU,vU,xU,omega_bU,...
            a_bU,S_U,y,param,R,ParamFilter,PosAmersU);

        [chiI,omega_bI,a_bI,P_I] = iekfUpdate(chiI,omega_bI,a_bI,...
            P_I,y,R,ParamFilter,yAmers);
        [RotI,vI,xI,PosAmersI] = chi2state(chiI);

        % save trajectory
        trajR = updateTraj(trajR,RotR,vR,xR,omega_bR,a_bR,i);
        trajL = updateTraj(trajL,RotL,vL,xL,omega_bL,a_bL,i);
        trajRef = updateTraj(trajRef,RotRef,vRef,xRef,omega_bRef,a_bRef,i);
        trajU = updateTraj(trajU,RotU,vU,xU,omega_bU,a_bU,i);
        trajI = updateTraj(trajI,RotI,vI,xI,omega_bI,a_bI,i);

        % remplace non visible landmarks
        [S_R,PosAmersR,ParamFilter,trackerBis,myTracks,PosAmersNew,...
            IdxAmersNew,trackCov,pointsMain,validityMain] = manageAmers(S_R,...
            PosAmersR,ParamFilter,ParamGlobal,trackerBis,...
            trajR,i,pointsMain,validityMain,IdxImage,myTracks,pointsBis);
        chiR = state2chi(RotR,vR,xR,PosAmersR);

        setPoints(trackerMain,pointsMain,validityMain);

        % remplace non visible landmarks for others filters with same
        % landmarks
        if isempty(IdxAmersNew) == 0
            P_L = S_L'*S_L;
            P_Ref = S_Ref'*S_Ref;
            P_U = S_U'*S_U;
            for jj = 1:length(IdxAmersNew)
                idx = IdxAmersNew(jj);
                idxP = 15+(3*idx-2:3*idx);
                P_L(:,idxP) = 0;
                P_L(idxP,:) = 0;
                P_L(idxP,idxP) = trackCov{jj};
                PosAmersL(:,idx) = PosAmersNew(jj,:)';
                P_Ref(:,idxP) = 0;
                P_Ref(idxP,:) = 0;
                P_Ref(idxP,idxP) = trackCov{jj};
                PosAmersRef(:,idx) = PosAmersNew(jj,:)';
                P_U(:,idxP) = 0;
                P_U(idxP,:) = 0;
                P_U(idxP,idxP) = trackCov{jj};
                PosAmersU(:,idx) = PosAmersNew(jj,:)';
                P_I(:,idxP) = 0;
                P_I(idxP,:) = 0;
                P_I(idxP,idxP) = trackCov{jj};
                PosAmersI(:,idx) = PosAmersNew(jj,:)';
            end
            S_L = chol(P_L);
            chiL = state2chi(RotL,vL,xL,PosAmersL);
            S_Ref = chol(P_Ref);
            chiRef = [RotRef,xRef;zeros(1,3) 1];
            S_U = chol(P_U);
            chiI = state2chi(RotI,vI,xI,PosAmersI);
        end
        disp(i/200);
        IdxImage = IdxImage+1;
    else
        trajR = updateTraj(trajR,RotR,vR,xR,omega_bR,a_bR,i);
        trajL = updateTraj(trajL,RotL,vL,xL,omega_bL,a_bL,i);
        trajRef = updateTraj(trajRef,RotRef,vRef,xRef,omega_bRef,a_bRef,i);
        trajU = updateTraj(trajU,RotU,vU,xU,omega_bU,a_bU,i);
        trajI = updateTraj(trajI,RotI,vI,xI,omega_bI,a_bI,i);
    end
end

%% Plots
errorR = computeError(trajR,trajReal,i);
errorL = computeError(trajL,trajReal,i);
errorU = computeError(trajU,trajReal,i);
errorRef = computeError(trajRef,trajReal,i);
errorI = computeError(trajI,trajReal,i);

figure;hold on;
plot(t(2:i-1)-t(2),errorR.errorR);
plot(t(2:i-1)-t(2),errorL.errorR);
plot(t(2:i-1)-t(2),errorRef.errorR);
plot(t(2:i-1)-t(2),errorU.errorR);
plot(t(2:i-1)-t(2),errorI.errorR);
disp(sqrt(mean([errorR.errorR errorL.errorR  errorRef.errorR errorU.errorR errorI.errorR].^2)));
legend('R-UKF-LG','L-UKF-LG','SE(3)-UKF','UKF','IEKF')
xlabel('t (s)')
ylabel('RMSE attitdude (°)')
title('RMSE on attitude as function of time')
figure
hold on;
plot(t(2:i-1)-t(2),errorR.errorX);
plot(t(2:i-1)-t(2),errorL.errorX);
plot(t(2:i-1)-t(2),errorRef.errorX);
plot(t(2:i-1)-t(2),errorU.errorX);
plot(t(2:i-1)-t(2),errorI.errorX);
disp(sqrt(mean([errorR.errorX  errorL.errorX errorRef.errorX errorU.errorX errorI.errorX].^2)));
legend('R-UKF-LG','L-UKF-LG','SE(3)-UKF','UKF','IEKF')
xlabel('t (s)')
ylabel('RMSE position (m)')
title('RMSE on position as function of time')

author: mbrossar


text:
function [soln,eqn,info] = Poisson3WG(node,elem,bdFlag,pde,option,varargin)
%% POISSON3WG Poisson equation: P1 linear element in 3-D.
%
%   u = POISSON3WG(node,elem,bdFlag,pde) produces the P0-P0 Weak Galerkin
%   finite element approximation of the Poisson equation
%
%       -div(d*grad(u))=f  in \Omega, with
%       Dirichlet boundary condition u=g_D on \Gamma_D,
%       Neumann boundary condition   d*grad(u)*n=g_N on \Gamma_N,
%       Robin boundary condition     g_R*u + d*grad(u)*n=g_N on \Gamma _R
%
%   The mesh is given by node and elem and the boundary face is given by
%   bdFlag. See meshdoc, bddoc for details. The data is given by the
%   structure pde which contains function handles f, g_D, g_N, g_R, or d.
%   For general elliptic equations with convection and reaction
%   coefficients, see ellipticpde.
%
%   The function Poisson assembes the matrix equation AD*u = b and solves
%   it by the direct solver (small size <= 2e3) or the multigrid solver
%   (large size > 2e3). The Dirichlet boundary condition is built into the
%   matrix AD and the Neumann boundary condition is build into b.
%
%   The diffusion coefficient d is a scalar function or a column array with
%   the same length as the elem array.
%
%   u = Poisson3WG(node,elem,bdFlag,pde,option) specifies the options.
%     - option.solver == 'direct': the built in direct solver \ (mldivide)
%     - option.solver == 'mg':     multigrid-type solvers mg is used.
%     - option.solver == 'notsolve': the solution u = u_D.
%   The default setting is to use the direct solver for small size problems
%   and multigrid solvers for large size problems. For more options on the
%   multigrid solver mg, type help mg.
%
%   When only one type of boundary condition is imposed, the input argument
%   bdFlag can be skipped. The boundary condition is implicitly given in
%   the pde structure by specifying g_D or g_N only. See examples below.
%
%   [u,A] = Poisson3WG(node,elem,bdFlag,pde,HB) returns also the
%   non-modified stiffness matrix A, which is semi-definite. The kernel of
%   A consists of constant vectors. The matrix A can be used to evulate the
%   bilinear form A(u,v) = u'*A*v, especially the enery norm of a finite
%   element function u by sqrt(u'*A*u).
%
%   [u,A,eqn] = Poisson3WG(node,elem,bdFlag,pde,HB) returns also the equation
%   structure eqn, which includes:
%     - eqn.AD:  modified stiffness matrix AD;
%     - eqn.b:   the right hand side.
%   The solution u = AD\b. The output eqn can be used to test other solvers.
%
%   [u,A,eqn,info] = Poisson3WG(node,elem,bdFlag,pde) returns also the
%   information on the assembeling and solver, which includes:
%     - info.assembleTime: time to assemble the matrix equation
%     - info.solverTime:   time to solve the matrix equation
%     - info.itStep:       number of iteration steps for the mg solver
%     - info.error:        l2 norm of the residual b - A*u
%     - info.flag:         flag for the mg solver.
%       flag = 0: converge within max iteration
%       flag = 1: iterated maxIt times but did not converge
%       flag = 2: direct solver
%       flag = 3: no solve
%
%   Example
%     clear all
%     node = [0,0; 1,0; 1,1; 0,1];
%     elem = [2,3,1; 4,1,3];
%     for k = 1:4
%       [node,elem] = uniformbisect(node,elem);
%     end
%     % Homogenous Dirichlet boundary condition
%     pde.f = inline('ones(size(p,1),1)','p');
%     pde.g_D = 0;
%     u = Poisson3WG(node,elem,pde);
%     figure(1);
%     showresult(node,elem,u);
%     % Non-homogenous Dirichlet boundary condition
%     pde.f = inline('-4*ones(size(p,1),1)','p');
%     pde.g_D = inline('p(:,1).^2 + p(:,2).^2','p');
%     u = Poisson3WG(node,elem,pde);
%     figure(2);
%     showresult(node,elem,u);
%     % Homogenous Neumann boundary condition
%     clear pde
%     pde.f = inline('pi^2*cos(pi*p(:,1)).*cos(pi*p(:,2))','p');
%     u = Poisson3WG(node,elem,pde);
%     figure(3);
%     showresult(node,elem,u);
%
%   Example
%     cubePoisson;

%   Reference: Programming of Weak Galerkin Method. Long Chen
%   https://www.math.uci.edu/~chenlong/ifemdoc/fem/WGprogramming.pdf
%
%   See also Poisson3, Poisson3WGfemrate, Lshape, crack, mg
%
%   Copyright (C) Long Chen. See COPYRIGHT.txt for details.

if ~exist('bdFlag','var'), bdFlag = []; end
if ~exist('option','var'), option = []; end
NT = size(elem,1);

%% Diffusion coefficient
t = cputime;  % record assembling time
if ~isfield(pde,'d'), pde.d = []; end
if ~isfield(option,'dquadorder'), option.dquadorder = 1; end
if ~isempty(pde.d) && isnumeric(pde.d)
   K = pde.d;                   % d is an array
end
if ~isempty(pde.d) && ~isnumeric(pde.d)
    [lambda,weight] = quadpts3(option.dquadorder);
    nQuad = size(lambda,1);
    K = zeros(NT,1);
    for p = 1:nQuad
		pxy = lambda(p,1)*node(elem(:,1),:) ...
			+ lambda(p,2)*node(elem(:,2),:) ...
			+ lambda(p,3)*node(elem(:,3),:) ...
            + lambda(p,4)*node(elem(:,4),:);
        K = K + weight(p)*pde.d(pxy);           % d is a function
   end
end

%% Compute geometric quantities and gradient of local basis
[elem2face,face] = dof3face(elem);
NF = size(face,1);
Ndof = NT + NF;
elem2dof = NT + elem2face;

%% Assemble stiffness matrix
A = sparse(Ndof,Ndof);
% compute ct2 = 1/mean(||x-xc||^2)
center = (node(elem(:,1),:) + node(elem(:,2),:) + ...
          node(elem(:,3),:) + node(elem(:,4),:))/4;
[lambda,weight] = quadpts3(2);
nQuad = size(lambda,1);
ct2 = zeros(NT,1);
for p = 1:nQuad
    pxyz = lambda(p,1)*node(elem(:,1),:) ...
         + lambda(p,2)*node(elem(:,2),:) ...
         + lambda(p,3)*node(elem(:,3),:) ...
         + lambda(p,4)*node(elem(:,4),:);
    ct2 = ct2 + weight(p)*sum((pxyz-center).^2,2);
end
ct2 = 1./ct2;
[Dphi,volume] = gradbasis3(node,elem);
clear center pxyz

% Mbb: face - face
for i = 1:4
    for j = i:4
        % local to global index map
		ii = double(elem2dof(:,i));
		jj = double(elem2dof(:,j));
        % compuation of local stiffness matrix.
        Aij = 9*dot(Dphi(:,:,i),Dphi(:,:,j),2).*volume + 9/16*ct2.*volume;
        if ~isempty(pde.d)
            Aij = K.*Aij;
        end
        if (j==i)
            A = A + sparse(ii,jj,Aij,Ndof,Ndof);
        else
            A = A + sparse([ii,jj],[jj,ii],[Aij; Aij],Ndof,Ndof);
        end
    end
end

% Mob: interior - face
Aij = -9/4*ct2.*volume;
if ~isempty(pde.d)
    Aij = K.*Aij;
end
Mob = sparse([(1:NT)', (1:NT)', (1:NT)', (1:NT)'], ...
             double(elem2dof(:)), [Aij, Aij, Aij, Aij], Ndof, Ndof);
A = A + Mob + Mob';

% Moo: diagonal of interor
Aij = 9*ct2.*volume;
if ~isempty(pde.d)
    Aij = K.*Aij;
end
A =  A + sparse(1:NT, 1:NT, Aij, Ndof, Ndof);

clear K Aij

%% Assemble right hand side
b = zeros(Ndof,1);
if ~isfield(option,'fquadorder')
    option.fquadorder = 2;   % default order
end
if ~isfield(pde,'f') || (isreal(pde.f) && (pde.f==0))
    pde.f = [];
end
if ~isempty(pde.f)
    [lambda,weight] = quadpts3(option.fquadorder);
	nQuad = size(lambda,1);
    bt = zeros(NT,1);
    for p = 1:nQuad
		% quadrature points in the x-y-z coordinate
		pxyz = lambda(p,1)*node(elem(:,1),:) ...
			 + lambda(p,2)*node(elem(:,2),:) ...
			 + lambda(p,3)*node(elem(:,3),:) ...
             + lambda(p,4)*node(elem(:,4),:);
		fp = pde.f(pxyz);
        bt = bt + weight(p)*fp;
    end
    bt = bt.*volume;
    b(1:NT) = bt;
end
clear pxyz bt

%% Set up boundary conditions
[AD,b,u,freeDof,isPureNeumann] = getbdWG3(A,b);

%% Record assembling time
assembleTime = cputime - t;
if ~isfield(option,'printlevel'), option.printlevel = 1; end
if option.printlevel >= 2
    fprintf('Time to assemble matrix equation %4.2g s\n',assembleTime);
end

%% Solve the system of linear equations
if isempty(freeDof), return; end
% Set up solver
if isempty(option) || ~isfield(option,'solver')    % no option.solver
    if Ndof <= 2e3  % Direct solver for small size systems
        option.solver = 'direct';
    else            % MGCG  solver for large size systems
        option.solver = 'mg';
    end
end
solver = option.solver;
% solve
switch solver
    case 'direct'
        tic;
        u(freeDof) = AD(freeDof,freeDof)\b(freeDof);
        residual = norm(b - AD*u);
        info = struct('solverTime',toc,'itStep',0,'err',residual,'flag',2,'stopErr',residual);
    case 'none'
        info = struct('solverTime',[],'itStep',0,'error',[],'flag',3);
    case 'mg'
        option.solver = 'CG';
        if nargin>=6
            HB = varargin{1};
        else
            HB = [];
        end
        if isfield(option,'reducesystem') && (option.reducesystem == 0)
            % original system
            option.x0 = u;
            option.freeDof = freeDof;
            [u,info] = mg(AD,b,elem,option,face,HB);
        else  % reduced system
            option.x0 = u(NT+1:end);
            option.freeDof = freeDof(NT+1:end)-NT;
            % eleminate elementwise dof
            Aoinv = spdiags(1./diag(AD(1:NT,1:NT)),0,NT,NT);
            Aob = AD(1:NT,NT+1:end);
            Abo = Aob';
            Abb = AD(NT+1:end,NT+1:end);
            Abbm = Abb - Abo*Aoinv*Aob;
            bm = -Abo*Aoinv*b(1:NT) + b(NT+1:end);
            [ub,info] = mg(Abbm,bm,elem,option,face,HB);
            u(1:NT) = Aoinv*(b(1:NT) - Aob*ub);
            u(NT+1:end) = ub;
        end
    case 'amg'
        option.solver = 'CG';
        [u(freeDof),info] = amg(AD(freeDof,freeDof),b(freeDof),option);
end
% post-process for pure Neumann problem
if isPureNeumann
    uc = sum(u(1:NT).*volume);
    u = u - uc;   % normalization for pure Neumann problem
end

%% Compute Du
dudx = u(elem2dof(:,1)).*Dphi(:,1,1)+u(elem2dof(:,2)).*Dphi(:,1,2) ...
     + u(elem2dof(:,3)).*Dphi(:,1,3)+u(elem2dof(:,4)).*Dphi(:,1,4);
dudy = u(elem2dof(:,1)).*Dphi(:,2,1)+u(elem2dof(:,2)).*Dphi(:,2,2) ...
     + u(elem2dof(:,3)).*Dphi(:,2,3)+u(elem2dof(:,4)).*Dphi(:,2,4);
dudz = u(elem2dof(:,1)).*Dphi(:,3,1)+u(elem2dof(:,2)).*Dphi(:,3,2) ...
     + u(elem2dof(:,3)).*Dphi(:,3,3)+u(elem2dof(:,4)).*Dphi(:,3,4);
Du = -3*[dudx, dudy, dudz];

%% Output information
if nargout == 1
    soln = u;
else
    soln = struct('u',u,'Du',Du);
    eqn = struct('A',AD,'b',b,'face',face,'freeDof',freeDof,'Lap',A);
    info.assembleTime = assembleTime;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% subfunctions getbd3
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [AD,b,u,freeDof,isPureNeumann]= getbdWG3(A,b)
    %% GETBDCR Boundary conditions for Poisson equation: Crouzeix-Raviart element.

    u = zeros(Ndof,1);
    %% Initial check
    if ~isfield(pde,'g_D'), pde.g_D = []; end
    if ~isfield(pde,'g_N'), pde.g_N = []; end
    if ~isfield(pde,'g_R'), pde.g_R = []; end

    %% Part 1: Modify the matrix for Dirichlet and Robin condition
    % Robin boundary condition
    Robin = [];
    idxR = (bdFlag(:) == 3);      % index of Robin faces in bdFlag
    if any(idxR)
        isRobin = false(Ndof,1);
        isRobin(elem2face(idxR)) = true;
        Robin = face(isRobin,:);  % Robin faces
    end
    if ~isempty(Robin) && ~isempty(pde.g_R) && ~(isnumeric(pde.g_R) && (pde.g_R == 0))
        v12 = node(Robin(:,2),:)-node(Robin(:,1),:);
        v13 = node(Robin(:,3),:)-node(Robin(:,1),:);
        area = 0.5*sqrt(abs(sum(mycross(v12,v13,2).^2,2)));
        center = (node(Robin(:,1),:) + node(Robin(:,2),:)+node(Robin(:,3),:))/3;
        ii = NT + find(isRobin);
        ss = pde.g_R(center).*area; % exact for linear g_R
        A = A + sparse(ii,ii,ss,Ndof,Ndof);
    end

    % Dirichlet boundary nodes: fixedFace
    fixedFace = []; freeFace= [];
    if ~isempty(bdFlag) % find boundary faces
        idxD = (bdFlag(:)==1); % all Dirichlet faces in bdFlag
        isFixedFace = false(NF,1);
        isFixedFace(elem2face(idxD)) = true;
        fixedFace = find(isFixedFace);
        freeFace = find(~isFixedFace);
    end
    if isempty(bdFlag) && ~isempty(pde.g_D) && isempty(pde.g_N) && isempty(pde.g_R)
        % no bdFlag, only pde.g_D is given in the input
        s = accumarray(elem2face(:), 1, [NF 1]);
        fixedFace = find(s == 1);
        freeFace = find(s == 2);
    end
    isPureNeumann = false;
    if isempty(fixedFace) && isempty(Robin)  % pure Neumann boundary condition
        % pde.g_N could be empty which is homogenous Neumann boundary condition
        isPureNeumann = true;
        fixedFace = 1;
        freeFace = (2:NF)';    % eliminate the kernel by enforcing u(1) = 0;
    end
    % Modify the matrix
    % Build Dirichlet boundary condition into the matrix AD by enforcing
    % AD(fixedFace,fixedFace)=I, AD(fixedFace,freeFace)=0, AD(freeFace,fixedFace)=0.
    if ~isempty(fixedFace)
        bdidx = zeros(Ndof,1);
        bdidx(NT + fixedFace) = 1;
        Tbd = spdiags(bdidx,0,Ndof,Ndof);
        T = spdiags(1-bdidx,0,Ndof,Ndof);
        AD = T*A*T + Tbd;
    else
        AD = A;
    end

    %% Part 2: Find boundary faces and modify the right hand side b
    % Find boundary faces: Neumann
    Neumann = [];
    if ~isempty(bdFlag)  % bdFlag specifies different bd conditions
        idxN = (bdFlag(:) == 2);      % all Neumann faces in bdFlag
        isNeumann = elem2face(idxN | idxR); % index of Neumann faces
        % since boundary integral is also needed for Robin edges
        Neumann = face(isNeumann,:);      % Neumann faces
    end
    if isempty(bdFlag) && (~isempty(pde.g_N) || ~isempty(pde.g_R))
        % no bdFlag, only pde.g_N or pde.g_R is given in the input
        s = accumarray(elem2face(:), 1, [NF 1]);
        Neumann = face(s == 1,:);
    end

    % Neumann boundary condition
    if ~isempty(Neumann) && ~isempty(pde.g_N) && ~(isnumeric(pde.g_N) && (pde.g_N == 0))
        if ~isfield(option,'gNquadorder')
            option.gNquadorder = 3;   % default order exact for linear gN
        end
        [lambdagN,weightgN] = quadpts(option.gNquadorder);
        nQuadgN = size(lambdagN,1);
        gf = zeros(size(Neumann,1),1);
        v12 = node(Neumann(:,2),:)-node(Neumann(:,1),:);
        v13 = node(Neumann(:,3),:)-node(Neumann(:,1),:);
        area = 0.5*sqrt(abs(sum(mycross(v12,v13,2).^2,2)));
        for pp = 1:nQuadgN
            % quadrature points in the x-y coordinate
            ppxy = lambdagN(pp,1)*node(Neumann(:,1),:) ...
                 + lambdagN(pp,2)*node(Neumann(:,2),:)...
                 + lambdagN(pp,3)*node(Neumann(:,3),:);
            gNp = pde.g_N(ppxy);
            gf = gf+ weightgN(pp)*gNp;
        end
        gf = gf.*area;
        b(NT+isNeumann) = b(NT+isNeumann) + gf;
    end
    % The case with non-empty Neumann faces but g_N=0 or g_N=[] corresponds to
    % the zero flux boundary condition on Neumann faces and no modification of
    % A,u,b is needed.

    % Dirichlet boundary condition
    if ~isPureNeumann && ~isempty(fixedFace) && ...
       ~isempty(pde.g_D) && ~(isnumeric(pde.g_D) && all(pde.g_D == 0))    % nonzero g_D
        if isnumeric(pde.g_D)
            u(fixedFace) = pde.g_D(fixedFace);
        else % pde.g_D is a function handle
            center = (node(face(fixedFace,1),:)+node(face(fixedFace,2),:)+node(face(fixedFace,3),:))/3;
            u(NT+fixedFace) = pde.g_D(center);
        end
        b = b - A*u;
        b(NT+fixedFace) = u(NT+fixedFace);
    end
    % The case with non-empty Dirichlet nodes but g_D=0 or g_D=[] corresponds
    % to the zero Dirichlet boundary condition and no modification of u,b is
    % needed.

    % Pure Neumann boundary condition
    if isPureNeumann
        b = b - mean(b); % compatilbe condition: sum(b) = 0
        b(1) = 0;
    end

    freeDof = [(1:NT)'; NT+freeFace];
    end % end of getbd3WG
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end

repo: ifem


text:
function [w]=iwtt(v,wtt_transform)
%Inverse WTT transform with previously computed filters
%   W=IWTT(V,WTT_TRANSFORM) computes the inverse WTT tranformation of a
%   given tensor with filters, stored in the WTT_TRANSFORM structure.
%
%
% TT-Toolbox 2.2, 2009-2012
%
%This is TT Toolbox, written by Ivan Oseledets et al.
%Institute of Numerical Mathematics, Moscow, Russia
%webpage: http://spring.inm.ras.ru/osel
%
%For all questions, bugs and suggestions please mail
%ivan.oseledets@gmail.com
%---------------------------

u=wtt_transform.u;
r=wtt_transform.rks;
sz=wtt_transform.sz;
w=iwtt_loc(v,u,r,sz);
return
end
function [w]=iwtt_loc(v,u,r,sz)
%[W]=IWTT(V,U,R,SZ)
%Applies Inverse WTT transform to vector V with linear filters U
%and ranks R
%SZ is optional
if ( nargin == 1 )
  sz=size(v);% d=numel(sz);
else
 % d=numel(sz);
end
%Apply one transform back
  N=numel(v);
  v=reshape(v,[r(1)*sz(1),N/(r(1)*sz(1))]);
if ( numel(u) == 1 )
   w=u{1}*v;
    w=reshape(w,[r(1),N/(r(1))]);
else
  %Simplest one is recursion
  w=v;
  w0=v(1:r(2),:);
  unew=u(2:numel(u)); rnew=r(2:numel(r)); sznew=[sz(2),sz(3:numel(sz))];
  w0=iwtt_loc(w0,unew,rnew,sznew);
  w(1:r(2),:)=w0;
  m=size(u{1},1);
  w=reshape(w,[m,N/m]);
  w=u{1}*w;
  w=reshape(w,[r(1),N/r(1)]);
end
  return
end

author: oseledets

repo: TT-Toolbox


text:
%SUBSREF Subscripted reference. Dataset overload
%
%
% $Id: subsref.m,v 1.20 2010/06/01 08:47:05 duin Exp $


text:
%ISDATAIM Returns true if a dataset contains image objects or image features
%
% 	N = ISDATAIM(A)
%
% INPUT
%		A 	Dataset
%
% OUTPUT
%		N	  Scalar: 1 if A contains images as objects or features, otherwise 0
%
% DESCRIPTION
% If no output argument is given, the function will produce an error if A does
% not contain image objects or features (i.e. it will act as an assertion).
%
% SEE ALSO (PRTools Guide)
% ISOBJIM, ISFEATIM

% $Id: isdataim.m,v 1.2 2006/03/08 22:06:58 duin Exp $

function n = isdataim (a)

		if (nargin < 1)
		error('Insufficient number of arguments.');
	end

	% The FEATSIZE and OBJSIZE fields of a dataset indicate whether it contains
	% images.

	n = (isa(a,'prdataset')) & ((length(a.featsize) > 1) | (length(a.objsize) > 1)) ;

	% Generate and error if the input is not a dataset with image data and
	% no output is requested (assertion).

	if (nargout == 0) & (n == 0)
		error('Dataset with images expected.')
	end

return


text:
function bnet = fgraph_to_bnet(fg)
% FGRAPH_TO_BNET Convert a factor graph to a Bayes net
% bnet = fgraph_to_bnet(fg)
%
% We assume all factors are tabular_CPD.
% We create 1 dummy observed node for every factor.

N = fg.nvars + fg.nfactors;
vnodes = 1:fg.nvars;
fnodes = fg.nvars+1:N;
dag = zeros(N);
for x=1:fg.nvars
  dag(x, fnodes(fg.dep{x})) = 1;
end
ns = [fg.node_sizes ones(1, fg.nfactors)];
discrete = [fg.dnodes fnodes];
bnet = mk_bnet(dag, ns, 'discrete', discrete);
for x=1:fg.nvars
  bnet.CPD{x} = tabular_CPD(bnet, x, 'CPT', 'unif');
end
ev = cell(1, fg.nvars); % no evidence
for i=1:fg.nfactors
  f = fnodes(i);
  e = fg.equiv_class(i);
  pot = convert_to_pot(fg.factors{e}, 'd', fg.dom{i}, ev);
  m = pot_to_marginal(pot);
  bnet.CPD{f} = tabular_CPD(bnet, f, 'CPT', m.T);
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
%Program for Construction of a two-out-of-two Visual Cryptography Scheme

%Author : Athi Narayanan S
%M.E, Embedded Systems,
%K.S.R College of Engineering
%Erode, Tamil Nadu, India.
%http://sites.google.com/site/athisnarayanan/
%s_athi1983@yahoo.co.in

%Program Description
%This program generates a two-out-of-two Visual Cryptography Scheme shares.
%The input image for this program should be a binary image.
%The shares and the overlapping result of the shares are written as output.
%
%Usage:
%Input: inImg - A binary image
%Output: share1  - Generated share 1
%        share2  - Generated share 2
%        share12 - Overlapped result of shares 1 & 2

function [share1, share2, share12] = VisCrypt(inImg)

s = size(inImg);
share1 = zeros(s(1), (2 * s(2)));
share2 = zeros(s(1), (2 * s(2)));

%%White Pixel Processing
%White Pixel share combinations
disp('White Pixel Processing...');
s1a=[1 0];
s1b=[1 0];
[x y] = find(inImg == 1);
len = length(x);

for i=1:len
    a=x(i);b=y(i);
    pixShare=generateShare(s1a,s1b);
    share1((a),(2*b-1):(2*b))=pixShare(1,1:2);
    share2((a),(2*b-1):(2*b))=pixShare(2,1:2);
end

%Black Pixel Processing
%Black Pixel share combinations
disp('Black Pixel Processing...');
s0a=[1 0];
s0b=[0 1];
[x y] = find(inImg == 0);
len = length(x);

for i=1:len
    a=x(i);b=y(i);
    pixShare=generateShare(s0a,s0b);
    share1((a),(2*b-1):(2*b))=pixShare(1,1:2);
    share2((a),(2*b-1):(2*b))=pixShare(2,1:2);
end

share12=bitor(share1, share2);
share12 = ~share12;
disp('Share Generation Completed.');


text:
function model = ivmOptimiseIvm(model, display)

% IVMOPTIMISEIVM Selects the points for an IVM model.
% FORMAT
% DESC selects points in an IVM model for inclusion in the active
% set.
% ARG model : the model for which the points are to be selected.
% ARG display : the display level, should be 0, 1 or 2. The higher
% the display, the more verbose.
% RETURN model : the model with the active set selected.
%
% SEEALSO : ivmInit, ivmSelectVisualise, ivmAddPoint,
% ivmSelectPoint, ivmEpUpdatePoint
%
% COPYRIGHT : Neil D. Lawrence, 2003, 2004, 2005

% IVM

if nargin < 2
  display = 1;
end
tol = 1e-4;
model = ivmInit(model);
numData = size(model.X, 1);

if model.d > numData
  warning('Active set size is larger than data-set, resetting to data-set size');
  model.d = numData;
end

% Start with all data-points inactive.
model.J = 1:numData;

% Set first infoChange to NaN
infoChange(1) = NaN;
dVal = model.d;
for k = 1:dVal

  [indexSelect, infoChange(k)] = ivmSelectPoint(model);
  dataIndexSelect = model.J(indexSelect);
  model = ivmAddPoint(model, dataIndexSelect);
  if display
    ivmSelectVisualise(model, display, k, dataIndexSelect);
  end
end
if model.epUpdate
  lengthIndex = length(model.I);
  betaChange = 1;
  oldBeta = model.beta;
  counter = 0;
  while betaChange > tol
    I = model.I';
    counter = counter + 1;
    for i = I;
      model = ivmEpUpdatePoint(model, i);
    end
    betaChange = full(max(abs(model.beta - oldBeta)));
    fprintf('EP Update %d, beta change %2.4f\n', counter, betaChange)
    oldBeta =  model.beta;
  end
end

author: SheffieldML

repo: GPmat


text:
function k = exor ( i, j )

%% EXOR calculates the exclusive OR of two integers.
%
%  Modified:
%
%    31 March 2003
%
%  Author:
%
%   John Burkardt
%
%  Reference:
%
%    Bennett Fox,
%    Algorithm 647:
%    Implementation and Relative Efficiency of Quasirandom
%    Sequence Generators,
%    ACM Transactions on Mathematical Software,
%    Volume 12, Number 4, pages 362-376, 1986.
%
%  Parameters:
%
%    Input, integer I, J, two values whose exclusive OR is needed.
%
%    Output, integer K, the exclusive OR of I and J.
%
  k = 0;
  l = 1;
%
  i = floor ( i );
  j = floor ( j );

  while ( i ~= 0 | j ~= 0 )
%
%  Check the current right-hand bits of I and J.
%  If they differ, set the appropriate bit of K.
%
    i2 = floor ( i / 2 );
    j2 = floor ( j / 2 );

    if ( ...
      ( ( i == 2 * i2 ) & ( j ~= 2 * j2 ) ) | ...
      ( ( i ~= 2 * i2 ) & ( j == 2 * j2 ) ) )
      k = k + l;
    end

    i = i2;
    j = j2;
    l = 2 * l;

  end

author: canlab

repo: CanlabCore


text:
function rsamobj = loadbobfile(infile, snum, enum)
    % LOADBOBFILE Load an RSAM binary file
    % Usage:
    %   self = loadbobfile(infile, snum, enum)
    % Example:
    rsamobj = rsam();
    if nargin<3
        help rsam/loadbobfile
        warning('Not enough input parameters')
        return
    end
    f = struct;
    f.file = infile;
    f.found = false;
    f.snum = snum;
    f.enum = enum;

    if exist(infile,'file')
        f.found = true;
    end

    [dirname,basename] = fileparts(infile);
    if length(basename)==8
        rsamobj.sta = basename(1:4);
    end

    rsamobj.files = f;
    rsamobj = load(rsamobj);
    rsamobj.snum = snum;
    rsamobj.enum = enum;
end

author: geoscience-community-codes

repo: GISMO


text:
function vis_im = visualise_bboxes(id, objectNames, cmap, rgb_dir, detection_dir)

font_size = 25;

city = strtok(id, '_');

image_filename = fullfile(rgb_dir, city, [id '_leftImg8bit.png']);
detection_filename = fullfile(detection_dir, [id '_leftImg8bit.mat']);

rgb_im = imread(image_filename);
temp = load(detection_filename);
detections = temp.dets.annotation;
vis_im = rgb_im;

if ~isfield(detections, 'object') || isempty(detections.object)
    return;
end

rectangle_position = zeros(length(detections), 4);
colour = zeros(length(detections), 3);
text_strings = cell(length(detections), 1);
% we eventually want to display the bboxes differently depending on whether
% they are a group annotation or not. unfortunately it seems insertShape as
% it is now doesn't support changes to line style.
is_grp = zeros(length(detections), 1);

for k = length(detections.object):-1:1
    name = detections.object(k).name;
    train_id = find(strcmp(objectNames, name)) - 1;
    assert(~isempty(train_id));
    xmin = str2double(detections.object(k).bndbox.xmin) + 1;
    ymin = str2double(detections.object(k).bndbox.ymin) + 1;
    xmax = str2double(detections.object(k).bndbox.xmax) + 1;
    ymax = str2double(detections.object(k).bndbox.ymax) + 1;
    width = xmax - xmin + 1;
    height = ymax - ymin + 1;
    rectangle_position(k, :) = [xmin, ymin, width, height];

    colour(k, :) = cmap(train_id + 1,:);
    text_strings{k} = name;
    is_grp(k) = logical(str2double(detections.object(k).is_grp));
    % this is a temporary workaround
    if is_grp(k)
        text_strings{k} = [text_strings{k}, 'group'];
    end

end

text_position = rectangle_position(:,1:2);
vis_im = insertShape(vis_im, 'Rectangle', rectangle_position,...
    'LineWidth', 4, 'Color', colour*255);

% scale font size according to width of bbox for clearer visual
max_font_shrink_ratio = 2;
pixel_to_font_size_ratio = 2.6; % for three-letter word
step_size = 0.25;
width_threshes = [0, ...
    (pixel_to_font_size_ratio/max_font_shrink_ratio + step_size) : step_size : pixel_to_font_size_ratio,...
    Inf] * font_size;
for k = (numel(width_threshes) - 1):(-1):1
    dets_in_range = width_threshes(k) <= rectangle_position(:,3) & ...
        rectangle_position(:,3) < width_threshes(k+1);
    adjusted_font_size = floor(max(width_threshes(k)/pixel_to_font_size_ratio, ...
        font_size/max_font_shrink_ratio));
    if any(dets_in_range)
        vis_im = insertText(vis_im, text_position(dets_in_range, :), text_strings(dets_in_range),...
            'Font', 'UbuntuMono-R', 'FontSize', adjusted_font_size,...
            'BoxColor', colour(dets_in_range, :)*255, 'BoxOpacity', 0.5, 'TextColor', 'white', ...
            'AnchorPoint', 'LeftBottom');
    end
end
end

author: qizhuli

repo: Weakly-Supervised-Panoptic-Segmentation


text:
function stop = stopOutFun(x,optimValues,state)
stop = false;
if optimValues.constrviolation < 1e-9
    stop = true;
end
end

author: ayonga

repo: frost-dev

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


text:
function [areaErrorMap,meanX,meanY] = mfmAreaErrorMap(unfoldMesh,nFaces,unfolded2D,errorList)
%
%  [areaErrorMap,meanX,meanY] = mfmAreaErrorMap(unfoldMesh,nFaces,errorList)
%
%Author: Wandell
%Purpose:
%   Compute the area distortions.  Extracted from ARW code.
%

Xcogs=unfolded2D(unfoldMesh.uniqueFaceIndexList(:),1);
Xcogs=reshape(Xcogs,nFaces,3);
meanX=mean(Xcogs,2);

Ycogs=unfolded2D(unfoldMesh.uniqueFaceIndexList(:),2);
Ycogs=reshape(Ycogs,nFaces,3);
meanY=mean(Ycogs,2);
areaErrorMap=flipud(makeMeshImage([meanY(:),meanX(:)],errorList,128));

return;

author: vistalab

repo: vistasoft


text:
function [D, alpha, E, A] = inexact_alm_rosl_subsampling(X, K, L, lambda, tol, maxIter)

% This matlab code implements the inexact augmented Lagrange multiplier
% method for Robust Orthogonal Subspace Learning with random sampling.
%
% X - m x n matrix of observations/data (required input)
%
% lambda - weight on sparse error term in the cost function
%
% tol - tolerance for stopping criterion.
%
% maxIter - maximum number of iterations
%
%
% min \sum_{1<i<T}\|\alpha_i\|_2 + \lambda\|E\|_1, s.t. \|D_i\|_2= 1 ;
%
% The laplacian function:
% L(A,E,Y,u) =  \sum_1^T |alpha_i|_2 + \lambda*|E|_1 + \mu/2* |X- \sum_1^T D_i*alpha_i-E+Y/mu|_F^2;
%
% Xianbiao Shu (xshu2@illinois.edu)
% Copyright: Mitsubishi Electric Research Lab
% Reference: X. Shu, F. Porikli, N. Ahujia "Robust Orthonormal Subspace Learning: Efficient Recovery of Corrupted Low-rank Matrices". CVPR 2014;

[m, n] = size(X);

% Random subsampling (using randsample or randperm)

row_sampled = randsample(m,L);
row_notsampled = [1:m]';
row_notsampled(row_sampled) =[];
row_ind = [row_sampled; row_notsampled];
column_sampled = randsample(n,L);
column_notsampled = [1:n]';
column_notsampled(column_sampled) =[];
column_ind = [column_sampled; column_notsampled];

X_perm = X(row_ind, :);
X_perm = X_perm(:, column_ind);

% % solve the simplified RSL
% [Dict_hat1, alpha] = inexact_alm_rsl(X_perm(1:L,:), K, lambda, tol, maxIter);
%
% % robust linear regression
% [Dict_hat2] = inexact_alm_rlr(X_perm(L+1:m,1:L), K, alpha(:, 1:L), tol, maxIter) ;

% solve the simplified RSL
[D, alpha_1] = inexact_alm_rosl(X_perm(:,1:L), K, lambda, tol, maxIter);

% robust linear regression
[alpha] = inexact_alm_rlr(X_perm(1:L, :), K, D(1:L, :), tol, maxIter) ;

% compute the low-rank matrix
%  D = [Dict_hat1; Dict_hat2];
 A = D*alpha;

 %permute back
 A(:, column_ind) = A;
 A(row_ind, :) = A;

 E = X -A;

end

author: andrewssobral

repo: lrslibrary


text:
function error=Tsai2SwaminathanError(X,camera,uv_d,xy_u)

% JMM Montiel. Sept 2005
%
% Error between the Tsia undistortion model and the
%   Swaminathan model
% Input
%   k1      -Swaminathan distortion model
%   camera  -Camera parameters. camera.k1 is irrelevant
%   uv_d    -Distorted point coordinates
%   uv_u    -Undistorted coordinates according to Tsai model

k1=X(1);
f=X(2);
camera.k1=k1;
camera.f=f;
camera.K(1,3)=X(3);
camera.K(2,3)=X(4);
uv_u_SN=undistort(camera,uv_d);
xy_u_SN=[uv_u_SN(:,1)-camera.K(1,3), uv_u_SN(:,2)-camera.K(2,3)]*camera.dx/f;
if(1)
    cla;
    plot(xy_u_SN(:,1),xy_u_SN(:,2),'+r');
    hold on
    plot(xy_u(:,1),xy_u(:,2),'+g');
    xlabel('red, Swaminathan model.    green Tsai with 1 radial distorion parameter')
end
error=[xy_u(:,1)-xy_u_SN(:,1);xy_u(:,2)-xy_u_SN(:,2)];

repo: ekfmonoslam


text:
function yval = r8poly_der_val ( n, poly_cof, xval )

%*****************************************************************************80
%
%% R8POLY_DER_VAL evaluates the derivative of a polynomial in standard form.
%
%  Discussion:
%
%    A polynomial in standard form, with coefficients POLY_COF(*),
%    may be written:
%
%      P(X) = POLY_COF(1)
%           + POLY_COF(2) * X
%           ...
%           + POLY_COF(N) * X**(N-1)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the polynomial.
%
%    Input, real POLY_COF(1:N), the polynomial coefficients.
%    POLY_COF(1) is the constant term, and POLY_COF(N) is the coefficient of
%    X**(N-1).
%
%    Input, real XVAL, a value where the derivative of the
%    polynomial is to be evaluated.
%
%    Output, real YVAL, the value of the derivative of the
%    polynomial at XVAL.
%
  yval = ( n - 1 ) * poly_cof(n);

  for i = n-1 : -1 : 2
    yval = yval * xval + ( i - 1 ) * poly_cof(i);
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function F_cell = flip_vec(F_cell)

chkarg(istypesizeof(F_cell, 'complexcell', [1 Axis.count]) || istypesizeof(F_cell, 'complexcell', [Axis.count 1]), ...
	'"F_cell" should be length-%d row or column cell array whose each element is array with complex elements.', Axis.count);

for w = Axis.elems
	F_cell{w} = flip_array(F_cell{w});
end

author: wsshin

repo: maxwellfdfd

save_path: github-repos/MATLAB/wsshin-maxwellfdfd


text:
function outputImage = plus(this, otherImage)
% Adds data of other image to this one and writes output in new image
%
% NOTE: If voxel dimensions of 2nd image do not match - or a scalar is
% given as 2nd argument - data in the 2nd argument is automatically
% replicated to match this image geometry.
%
%
%   Y = MrImage()
%   outputImage = plus(Y, otherImage, ...
%   functionHandle)
%
% This is a method of class MrImage.
%
%
% IN
%   otherImage              image that will be added to this one
%
% OUT
%   outputImage             new MrImage, sum of this and otherImage
%
% EXAMPLE
%
%   % Compute sum of 2 images
%		Y = MrImage();
%		Z = MrImage();
%		X = Y.plus(Z);
%
%   % OR (cool overload!):
%       X = Y+Z
%
%   See also MrImage perform_binary_operation

% Author:   Saskia Bollmann & Lars Kasper
% Created:  2014-11-13
% Copyright (C) 2014 Institute for Biomedical Engineering
%                    University of Zurich and ETH Zurich
%
% This file is part of the TAPAS UniQC Toolbox, which is released
% under the terms of the GNU General Public Licence (GPL), version 3.
% You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version).
% For further details, see the file COPYING or
%  <http://www.gnu.org/licenses/>.

outputImage = this.perform_binary_operation(otherImage, @plus);

author: translationalneuromodeling

repo: tapas


text:
%% EXAMPLE AGENT CLASS (agent_example.m) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This agent is intended to describe how to create a new agent class for
% simulation under the Open Multi-Agent Simulation (OpenMAS) framework. It
% contains examples of the necessary properties and methods the class must
% have in order to interact with the environment.

% GENERAL NOTICE FOR AGENT DEFINITIONS
% - All objects(including agents) must be in '/objects' for dependancies.
% - The 'example_setup' demonstrates how the dependancies are setup.
% - All objects must have:
%   + A 'constructor' method (to create it).
%   + A 'main' method (to evaluate updates).
% - During the 'main' method, one cycle of the agents loop is computed,
%   updating the objects local state, and global properties.

% To call the agent class, your work directory must be within '/objects',
% otherwise it must be added using 'addpath('objects')' to include the
% complete inheritance tree.

% The class is called using the line:
% testAgent = agent_example();

% If have any questions or help getting started, please contact me directly
% at: jadouthwaite1@sheffield.ac.uk

% Author: James A. Douthwaite

classdef agent_example < agent
    % SPECIFY YOUR AGENT-SPECIFIC PROPERTIES HERE
    % These are in addition to those inherited from 'agent'.
    properties
        % EXAMPLE PROPERTIES (Unique to this agent)
        % noOfWings = 2;
    end
    %% //////////////////// MAIN (required) METHODS ///////////////////////
    methods
        % CONSTRUCTOR
        function obj = agent_example(varargin)
            % This function is called to create the 'agent_example' class,
            % which then takes on the new parameters specified.

            % INPUT HANDLING
            if length(varargin) == 1 && iscell(varargin)                   % Catch nested cell array inputs
                varargin = varargin{:};
            end

            % CALL THE SUPERCLASS CONSTRUCTOR
            obj@agent(varargin);                                           % Create the super class 'agent'
            % Set some of the parameters
            obj = obj.SetRadius(0.5);
            obj = obj.SetDetectionRadius(50);                              % Update the range attribute to the SIM VIRTUAL property
            % Parse any overrides
            [obj] = obj.configurationParser(obj,varargin);
        end
        % AGENT MAIN CYCLE
        function [obj] = main(obj,ENV,varargin)
            % This function is designed to contain everything your agent does
            % in a given simulation timestep. As an 'agent', a list of
            % detected entities is given if detected.
            % cycle that results in an acceleration vector in the global axis.
            % INPUTS:
            % ENV      - The simulation environmental update (currentTime,dt, etc..)
            % varargin - Cell array of optional inputs

            % OUTPUTS:
            % obj      - The updated project

            % GET THE TIMESTEP
            if isstruct(ENV)
                dt = ENV.dt;
            else
                error('OpenMAS environment structure is invalid.');
            end

            % PARSE THE ENTITY PROPERTIES FROM THE ENVIRONMENT
            [obj,obstacleSet,agentSet] = obj.GetAgentUpdate(dt,varargin{1});       % IDEAL INFORMATION UPDATE

            % \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
            %        INSERT ALGORITHM/DECISION MAKING PROCESS HERE
            % \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
            % DO NOTHING
            linearRates  = [0.5;0;0]; % (m/s)
            angularRates = [0;0;0.1]; % (rad/s)

            % \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
            %    DECIDE HOW THE AGENTS/ENTITIES LOCAL STATE IS UPDATED
            % \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
            % How each agent updates its state can be defined locally, as
            % a dedicated function.
            % CALL THE SEPERATE UPDATE FUNCTION
            [state_k_plus] = obj.updateLocalState(ENV,...
                                       obj.localState,...
                                          linearRates,...
                                         angularRates);

            % \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
            %  USE THE UPDATED LOCAL STATE TO UPDATE THE GLOBAL PROPERTIES
            % \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
            % How the objects global state is updated is dependant on how
            % the 'state_k_plus' is defined.
            obj = obj.updateGlobalProperties(dt,state_k_plus);
        end
    end
    %% /////////////////////// AUXILLARY METHODS //////////////////////////
    methods
        % ASSIGNED STATE UPDATE FUNCTION (USING ODE45)
        function [X] = updateLocalState(obj,TIME,X,velocity,omega)
            % This function computes the state update for the current agent
            % using the ode45 function.

            % DETERMINE THE INTEGRATION PERIOD
            if TIME.currentTime == TIME.timeVector(end)
                return
            else
                tspan = [TIME.currentTime TIME.timeVector(TIME.currentStep + 1)];
                opts = odeset('RelTol',1e-2,'AbsTol',1e-4);
                [~,Xset] = ode45(@(t,X) obj.dynamics_simple(X,velocity,omega), tspan, X, opts);
                X = Xset(end,:)';
            end
        end
    end
end

repo: OpenMAS


text:
function [anat, inplanes] = GetAnatomyFromAnalyze(fileName)
% [anat, inplanes] = GetAnatomyFromAnalyze(fileName)
%
% Similar to getAnatomy but you pass in a 3d Analyze / nifti file that you
% already made. Requires SPM5
% ARW 061407

anat = [];

if ~exist(fileName, 'file');
    disp(fileName);
    disp('File does not exist')

    return;
end

    V=spm_vol(fileName);

    [anat,XYZ] = spm_read_vols(V);
      % Do a rot90 on each slice of the anat. That's just the way it is...
    for thisSlice=1:size(anat,3)
        anat2(:,:,thisSlice)=rot90(anat(:,:,thisSlice),1);
    end
    anat=anat2;

    % Use imatrix to recover slice spacing etc.
    [inTrans]=spm_imatrix(V.mat);

    sz = size(anat);
    mm=abs(inTrans(7:9));

    inplanes.FOV = mm(1)*sz(1);
    inplanes.fullSize = sz(1:2);
    inplanes.voxelSize = mm;
    inplanes.spacing = 0;
    % We already checked hdr.image.slquant==nList, not again -- Junjie
    inplanes.nSlices = sz(3);

    inplanes.examNum = V.descrip;
    inplanes.crop = [];
    inplanes.cropSize = [];


text:
function [x, zone, variable, index_c1, index_c2, index_bc] = ...
         twod_readcadyf(file,num_fig)
%-----------------------------------------------------------------------
%  twod_readcadyf.m - Reads in mesh, connectivity and boundary node
%                     numbers from a *.cadyf file generated with
%                     amiral (preprocessor for cadyf).
%
%                     This code is quite problem specific, since the
%                     *.cadyf file is setup to solve the Navier-Stokes
%                     equations.
%
%                     ! NOTE: Matlab has difficulty with D+00 double
%                     precision notation, globally change to E+00
%
%  Copyright (c) 2001, Jeff Borggaard, Virginia Tech
%  Version: 1.0a
%
%  Usage:     [x, zone, variable, index_c1, index_c2, index_bc] = ...
%             twod_readcadyf(file,num_fig)
%
%  Variables:     file
%                        Name of the *.cadyf file.
%                 num_fig
%                        If included, plot the mesh
%
%                 x
%                        Nodal coordinates
%                 zone(#).e_conn
%                        Element connectivity
%                 variable
%                        Contains finite element information
%                 index_c1
%                        Node numbers corresponding to curve c1
%                 index_c2
%                        Node numbers corresponding to curve c2
%                 index_bc
%                        Node numbers corresponding to curve bc
%-----------------------------------------------------------------------
fid = fopen(file);

% match up with indata.f and contrl.f in cadyf

%-----------------------------------------------------------------------
%  From control blocks
%-----------------------------------------------------------------------
% card # 1
temp  = fgetl(fid);           % title line

% card # 2
A     = fscanf(fid,'%i',9);   % itim, lin, etc.    (cart43.f)
itemp = A(5);
iturb = A(8);

A     = fscanf(fid,'%g',9);   % icomp, itype, etc.

if (iturb == 3)
  temp = fgetl(fid);          % turbulent source terms
  temp = fgetl(fid);          % model coefficients
end

% card # 3
A = fscanf(fid,'%i %i %i %i %i %i %i %i %i',9); %  (cart44.f)

% get the number of nodes and elements
n_node  = A(1);
n_elem  = A(2);
n_zones = A(3);

% fprintf('twod_readcadyf: reading in mesh with %g nodes and %g elements\n',...
%          n_node,n_elem);

% skip over cards # 4 and # 5
temp = fscanf(fid,'%s',4);  % execution type, etc.  (cart45.f)
temp = fscanf(fid,'%s',4);  % body forces           (cart46.f)

% card # 6
A = fscanf(fid,'%i %i %i %i %i %i',6);    %         (cart47.f)

% read in number of specified boundary conditions
n_dbc = sum(A(1:6)); % number of bcs for u,v,w,T,k,e
% fprintf('                reading %g boundary conditions\n',n_dbc);

% card # 7 (assuming no neuman bc's)
temp = fscanf(fid,'%g',6);                %         (cart48.f)
% card # 8 (assuming no periodic bc's)
temp = fscanf(fid,'%g',1);                %         (cart49.f)

temp = fscanf(fid,'%g',12); % card #9     %         (ct4101.f)
temp = fscanf(fid,'%g',7);  % card #10    %         (ct4102.f)

% card #11
A = fscanf(fid,'%g',8);                   %         (ct411.f)
ninitu = A(5); % number of initial guess to read/skip
ninitt = A(6);
ninitk = A(7);
ninite = A(8);
ntoti  = ninitu + ninitt + ninitk + ninite;
% fprintf('                reading %g initial guess lines\n',ntoti);

temp = fscanf(fid,'%g',15); % card #12
temp = fscanf(fid,'%g',8);  % card #13
temp = fscanf(fid,'%g',2);  % card #14
temp = fscanf(fid,'%g',5);  % card #15
temp = fscanf(fid,'%g',6);  % card #16              (ct416.f)

%-----------------------------------------------------------------------
%  Now read in nodes and connectivity information
%-----------------------------------------------------------------------
% read in inblok data (skipped for now)
% read in time functions (intfns) (skipped for now)
% read in periodic nodes (skipped for now)

% read in nodal coordinates

for i=1:n_node
  A = fscanf(fid,'%g',12);    %      (cart7.f)
  if (length(A)<12)
    error('You need to change Ds to Es, Matlab doesnt recognize doubles')
  end
  x(i,1) = A(10);
  x(i,2) = A(11);
end

% read in blade normal information (inblad) (skipped for now)

% read in boundary conditions (incons)
%  this can be used to "tag" certain nodes...
index_c1 = [];
index_c2 = [];
index_bc = [];
for i=1:n_dbc
  A = fscanf(fid,'%i %g %i %g',4);      %   (cart9.f)

  if (A(2)==1)
    index_c1 = [index_c1 A(1)];
  elseif (A(2)==2)
    index_c2 = [index_c2 A(1)];
  else
    index_bc = [index_bc A(1)];
  end

  variable.dirichlet_node(i) = A(1);
  variable.dirichlet_value(i) = max(0,A(2));
end

% read in neuman data (inload, skipped for now)
% read in volumetric body forces (invol, skipped for now)

% read in, or skip, reading in the initial guess... (inital)
for i=1:ninitu
  temp = fscanf(fid,'%g',4);              %  (cart10.f)
end
for i=1:ninitt
  temp = fscanf(fid,'%g',2);
end
for i=1:ninitk
  temp = fscanf(fid,'%g',1);
end
for i=1:ninite
  temp = fscanf(fid,'%g',1);
end

for n=1:n_zones
  % skip over three lines
  A = fscanf(fid,'%i',8);
  n_elzone = A(2);

  A = fscanf(fid,'%i',8);
  A = fscanf(fid,'%i',8);

  % skip over material properties
  temp   = fscanf(fid,'%g',1); % rho
  temp   = fscanf(fid,'%g',1); % viscosity
  if ( itemp==1 | itemp==2 )
    temp = fscanf(fid,'%g',1); % cp
    temp = fscanf(fid,'%g',1); % cond
  end

  % read in element connectivity (assuming 6 noded elements)
  zone(n).e_conn = [];
  for i=1:n_elzone
    A = fscanf(fid,'%i',7);
    zone(n).e_conn(i,:) = A(2:7)';
  end
end

if ( nargin==2 )
  for n=1:n_zones
    twod_plotm2(num_fig,x,zone(n).e_conn)
  end
  figure(num_fig)

  hold on
  for i=1:length(index_c1)
    n = index_c1(i);
    plot(x(n,1),x(n,2),'g*')
  end
  for i=1:length(index_c2)
    n = index_c2(i);
    plot(x(n,1),x(n,2),'r*')
  end
  for i=1:length(index_bc)
    n = index_bc(i);
    plot(x(n,1),x(n,2),'b*')
  end
  pause(.001)
end


text:
N = 400;
A = gallery('poisson', N);
[L2,p,Ac] = achol(A);
% [L2,p,Ac] = acholold(A);
% save olddata L2 p Ac
disp(size(Ac,1));
    b = ones(size(A,1),1);
    tol = 1e-6; maxit = 100;
    tic;
%     Ap = A(p,p);
    [x2,fl2,rr2,it2,rv2] = pcg(A,b,tol,maxit,@(r)acholpre(r,A,L2,L2',p,Ac));
    toc;
    fprintf('#dof: %8.0u,  iter: %2.0u\n',size(A,1), it2)
    semilogy(0:it2,rv2./norm(b),'b.');
    hold on